Network Reconstruction with Ambient Noise
NNetwork Reconstruction with Ambient Noise
Melvyn Tyloo, Robin Delabays, and Philippe Jacquod School of Engineering, University of Applied Sciences of Western Switzerland HES-SO, CH-1951 Sion, Switzerland. Automatic Control Laboratory, ETH Zurich, CH-8092 Zurich, Switzerland Department of Quantum Matter Physics, University of Geneva, CH-1211 Geneva, Switzerland (Dated: August 3, 2020)The dynamics of systems of coupled agents is determined by the structure of their couplingnetwork. Often, the latter is not directly observable and a fundamental, open question is how toreconstruct it from system measurements. We develop a novel approach to identify the networkstructure underlying dynamical systems of coupled agents based on their response to homogeneous,ambient noise. We show that two-point frequency signal correlators contain all the informationon the network Laplacian matrix. Accordingly, when all agents are observable, the full Laplacianmatrix can be reconstructed. Furthermore, when only a fraction of the agents can be observed, pairsof observable agents can be ranked in order of their geodesic distance when the noise correlation timeis short enough. The method is computationally light and we show numerically that it is accurateand scalable to large networks.
Introduction.
Network science – the field that studiescomplex, networked systems [1] – has seen an enormousgrowth of activity in recent years. More and more di-verse systems of physical, life and human sciences areanalyzed through larger and larger models of agents con-nected to one another [2], thanks in large part to the ever-increasing capacity for data mining and processing [3].Network science draws both on analytical methods, forinstance from graph theory or statistical mechanics, andon purely data-based approaches. Analytical approachesare arguably much more powerful than data-based ap-proaches, however they rely on the precise knowledge ofthe underlying network and become harder to apply inlarge networks. Data-based approaches rely only on ob-servations, but the collected data are algorithmically pro-cessed with little gain in physical knowledge. Combiningthe two approaches may compensate for the weaknessesof one with the strengths of the other. In this manuscriptwe connect data-based to analytical approaches. Weshow that the topology of a coupling network can beprecisely reconstructed from sufficiently long sets of mea-surement data on systems subjected to uncontrolled, ho-mogeneous noise from their environment. Accordingly,one big advantage of our reconstruction method is thatit requires only the ability to measure the dynamics ofthe n agents, but does not require to precisely controlits input signal. The approach will be precious to inferthe structure of unknown, noisy networks such as, e.g.social networks which change over short time scales [4],interconnected power grids whose topology is determinedby line faults and disconnections that are not systemat-ically reported [5], or gene regulatory networks made ofsuch huge numbers of proteins and genes that the exactstructure of their interaction network cannot be knowna priori [6].A widely used method to reconstruct unknown inter-action networks is to inject a probe signal somewhere inthe network and to measure the response dynamics of the agents [7–13]. The successful reconstruction of thenetwork topology, through e.g., the Laplacian or adja-cency matrix, requires however that one is able not onlyto measure agents everywhere in the network, but addi-tionally that one can control and inject specially tailoredprobe signals. Another difficulty is that the probe signalsmay modify the network dynamics. Soft approaches havebeen used in power grids, with injected signals with smallamplitude and in frequency ranges outside the networkbandwidth [13–15]. Because of the restricted probe fre-quency range, the method is then limited and identifiesonly certain network modes [14, 15] or requires a largenumber of probings [13].A different approach relies on passive observation –i.e., without probe signal injection – of the agents dy-namics. A brute-force approach is then to minimize acost function over all the network parameters to identifythe network edges [16–18]. While the approach worksin principle, the required computation time scales atleast as O ( n ), with the number n of agents [17, 18],and it quickly becomes computationally prohibitive inlarge networks. Lighter approaches identify edges be-tween pairs of agents through trajectory correlators andGranger causality [19, 20], or leverage a Bayesian ap-proach to determine the most likely network structure,given a set of data [21, 22]. Ref. [23] extracts spectralmoments of the network from its dynamics, but cannotdirectly reconstruct its network matrix. Recent worksleverage the response of a system to some noisy signal toimprove the inference accuracy [24]. Rather interestingly,Refs. [25, 26] noticed that the accuracy of their machinelearning network inference was improved by dynamicalnoise. So far, this is however only an observation andno formal understanding has been provided. Reviews ofnetwork reconstruction in different scientific domains aregiven in Refs. [27] and [28]. It is commonly acceptedthat measurements on all nodes are prerequisite for fullnetwork reconstruction. a r X i v : . [ n li n . AO ] J u l All existing reconstruction methods we know of suf-fer from at least one of several weaknesses. They eitherrequire a large degree of control over the investigatedsystem, are approximate or computationally prohibitive,or they deliver only limited information on the networkstructure. In this paper, we propose a novel technique fornetwork reconstruction which does not suffer from any ofthese deficiencies. To fully reconstruct the network, itrequires the passive observation of the dynamics of the n agents subjected to unavoidable ambient noise. Becauseit relies on the computation of two-point correlators ofthe agent dynamics, the method requires sets of measure-ments on all n nodes and a computation time scaling as O ( n ) for high frequency ambient noise, or O ( n ) for lowfrequency ambient noise for which the method requires aLaplacian matrix inversion.Measuring all nodes may turn out to be impractical oreven unreachable in large networks, and one often has tosettle for partial measurements on only a fraction of theagents. In that case, our method is still able to identifythe geodesic distance between pairs of such observablenodes, when the noise has a sufficiently small correlationtime. Network-coupled dynamical systems.
We consider anensemble of i = 1 , . . . , n agents with coordinates x i ∈ R ,with a dynamics governed by a set of coupled ordinarydifferential equations,˙ x i = ω i − (cid:88) j a ij f ij ( x i − x j ) + δω i ( t ) , (1)where ω i ∈ R are constant natural frequencies with (cid:80) i ω i = 0, and δω i ( t ) is the unavoidable ambient noisedue to the environment. The interaction between agentsis a differentiable function f ij : R → R , that is even in itsindices i and j and odd in its argument, and a ij ≥ a ij are sufficientlylarge and numerous, Eq. (1) with δω i ( t ) = 0 has a stablefixed point which we denote x ∗ ∈ R n . We observe thedynamics generated by δω i ( t ) (cid:54) = 0 in the vicinity of thisfixed point. Accordingly, we consider small deviations δx = x − x ∗ and linearize Eq. (1) in δx to obtain δ ˙ x = − J ( x ∗ ) δx + δω , (2)where the Jacobian matrix reads J ij ( x ∗ ) = (cid:40) − a ij ∂ x f ij ( x ) (cid:12)(cid:12) x = x ∗ i − x ∗ j , i (cid:54) = j , (cid:80) k a ik ∂ x f ik ( x ) (cid:12)(cid:12) x = x ∗ i − x ∗ k , i = j . (3)The matrix J is a real, symmetric Laplacian matrix ofthe interaction network, that is weighted by the deriva-tive of f ij at the fixed point. It contains information onboth the interaction network and the fixed point. ThisLaplacian is the matrix we want to reconstruct. Un-der our assumption of a stable fixed point, it has real, nonnegative eigenvalues, 0 = λ < λ ≤ ... ≤ λ n andan orthonormal basis of eigenvectors { u α } nα =1 , with aconstant-component zero-mode u i = n − / , ∀ i .Equation (2) is solved by expanding δx over the eigen-basis { u α } , δx ( t ) = (cid:88) α c α ( t ) u α , (4)which yields a set of coupled Langevin equations for c α ( t )whose solution reads [29] c α ( t ) = e − λ α t (cid:90) t e λ α t (cid:48) u α · δω ( t (cid:48) )d t (cid:48) . (5) Network reconstruction from ambient noise.
It is rea-sonable to expect that the network’s randomly fluctuat-ing environment generates noise with spatial correlationsdecaying fast with distance. Accordingly we define δω i ( t )in Eq. (1) as a spatially uncorrelated noise with its firsttwo moments given by (cid:104) δω i ( t ) (cid:105) = 0 , (6a) (cid:104) δω i ( t ) δω j ( t ) (cid:105) = δω δ ij exp ( −| t − t | /τ ) , (6b)where δω is the noise standard deviation, τ its cor-relation time, and the brackets denote averaging overnoise realizations or a large enough observation time, (cid:104) . . . (cid:105) = T − (cid:82) T . . . d t , T (cid:29) τ .The equal time two-point frequency correlator isstraightforwardly calculated in the limit of long obser-vation times ( λ α T (cid:29) (cid:104) δ ˙ x i ( t ) δ ˙ x j ( t ) (cid:105) = δω (cid:16) δ ij − (cid:88) α ≥ u α,i u α,j λ α τ λ α τ (cid:17) . (7)Two regimes are of particular interest. First, in the limitof short correlation time, λ α τ <
1, Taylor-expandingEq. (7) and realizing that the matrix elements of the k th power of the Laplacian read ( J k ) ij = (cid:80) α λ kα u α,i u α,j gives (cid:104) δ ˙ x i δ ˙ x j (cid:105) = δω (cid:34) δ ij + ∞ (cid:88) k =1 ( − τ ) k ( J k ) ij (cid:35) . (8)Second, in the other limit of long correlation time, λ α τ >
1, a similar series expansion of Eq. (7) gives thistime (cid:104) δ ˙ x i δ ˙ x j (cid:105) = δω (cid:34) n − − ∞ (cid:88) k =1 ( − τ ) − k ( J − k ) ij (cid:35) , (9)where J − k stands for the k th power of the pseudo-inverse J † of the Laplacian.Equations (8) and (9) form the basis of our networkreconstruction approach. In both deep asymptotic limits Figure 1. (a): Random interaction network used for reconstruction, with color-coded geodesic distances from node a ij ∈ [0 . , λ n τ = 0 .
014 and (d) Eq. (10b) with λ τ = 4 . T with λ T = 500 and over 40 different noise realizations. of either very short noise correlation time τ or very long τ the network Laplacian is reconstructed asˆ J ij = ( δ ij − (cid:104) δ ˙ x i δ ˙ x j (cid:105) /δω ) τ − , λ α τ → , (10a)ˆ J † ij = ( (cid:104) δ ˙ x i δ ˙ x j (cid:105) /δω − n − ) τ , λ α τ → ∞ . (10b)The short correlation time asymptotics of Eq. (10a) al-lows for a direct reconstruction of the Laplacian matrix.Based on the calculation of two-point correlation func-tions it requires a computation time scaling as O ( n )and therefore allows to reconstruct large networks as wenumerically illustrate below. The long correlation timeasymptotic of Eq. (10b) requires on the other hand to firstcompletely reconstruct the pseudo-inverse of the Lapla-cian and second to invert that matrix to obtain the net-work Laplacian, thereby requiring a computation timescaling as O ( n ).More generally, ambient noise may be given by a su-perposition of different, uncorrelated noises δω i ( t ) = (cid:80) kα =1 δω ( α ) i ( t ), with noise sequences δω ( α ) i each with itsown standard deviation δω α and correlation time τ α . Ourapproach remains applicable in that case as long as ei-ther max α,j λ α τ j < α,j λ α τ j >
1. Equation (10)becomesˆ J ij = (cid:16) δ ij (cid:88) α δω α − (cid:104) δ ˙ x i δ ˙ x j (cid:105) (cid:17)(cid:46) (cid:88) α δω α τ α , (11a)ˆ J † ij = (cid:16) (cid:104) δ ˙ x i δ ˙ x j (cid:105) − n − (cid:88) α δω α (cid:17)(cid:46) (cid:88) α δω α τ − α , (11b)in the short and and long correlation time asymptotics,respectively. Up to noise-dependent but spatially homo-geneous factors (cid:80) α δω α , (cid:80) α δω α τ α , and (cid:80) α δω α τ − α , theLaplacian matrix can be reconstructed as before. Geodesic distances.
In the short correlation time limitof Eq. (8), it is furthermore possible to partially recon-struct the network even with a limited number of mea-surements on only a fraction of the nodes. This is so,because the geodesic distance between any two nodes, i and j , is given by the minimal exponent q for which( J q ) ij (cid:54) = 0. Therefore, (cid:104) δ ˙ x i δ ˙ x j (cid:105) = δω δ ij + ∞ (cid:88) k = q ( − τ ) k (cid:0) J k (cid:1) ij , (12)which makes it possible to determine the geodesic dis-tance q between any measurable pair of nodes ( i, j ) aslong asmin l,m ( J q − ) lm τ − (cid:29) ( J q ) ij (cid:29) max l,m ( J q +1 ) lm τ , (13)where the minimum (resp. maximum) is taken over pairs( l, m ) of nodes with geodesic distance ≤ q − ≥ q + 1). When Eq. (13) holds, pairs of nodes with geodesicdistance q have noise correlators sufficiently away fromthose with geodesic distances q − q + 1 that onecan easily identify them. Numerical illustrations.
We first validate our ap-proach on a random network with n = 20 agents. Thetime evolution is given by Eq. (1), with interaction net-work edges shown in Fig. 1(a), randomly distributed as a ij ∈ [0 . , f ij ( x ) = sin( x ),corresponding to Kuramoto oscillators [30].Figures 1(c) and (d) display the Laplacian matrix withelements ˆ J ij inferred through Eqs. (10a) and (10b) re-spectively. Figure 1(b) shows the original Laplacian ma-trix J . The relation between inferred and real elementsof the Laplacian is shown in Fig. 2(a). The agreementbetween real and reconstructed network Laplacian is al-most perfect, with both short and long correlation timemethods.We next demonstrate the scalability of our method andapply it to the PanTaGruEl network model of the Eu-ropean electric power grid [31, 32]. The network has n = 3809 agents and m = 4944 edges. Figure 2(b)shows the inferred vs. real Laplacian matrix elements forthis large network. The short correlation time method Figure 2. (a) Comparison between inferred, ˆ J , and true, J ,Laplacian of the weighted network shown in Fig. 1(a). Darkblue crosses are obtained using Eq. (10a) with λ n τ = 0 . λ τ = 4 .
97. (b)Comparison between inferred and true Laplacian for a net-work corresponding to the European high-voltage electricalgrid [31, 32]. Dark blue crosses are obtained using Eq. (10a)with λ n τ = 0 .
029 and yellow crosses using Eq. (10b) λ τ =94. of Eq. (10a) correctly reconstructs the network matrix,with a systematic rescaling factor of ∼ τ > τ →
0. Including the next order correction from Eq. (8)into Eq. (10) givesˆ J ij = J ij − ( J ) ij τ + O ( τ ) . (14)Because for most networks with J ij > i (cid:54) = j , and inparticular for the PanTaGruEl network J ij and ( J ) ij have the same sign, the error in reconstructing off-diagonal matrix elements goes systematically towardssmaller absolute values, in agreement with the datashown in Figs. 2(a) and 2(b). It is furthermore propor-tional to τ as we numerically checked, but do not show.The error between the reconstructed and the real ma-trix elements in the long correlation time also appearsto be systematic. A similar argument as just made in-dicates that it is proportional to τ − . We attribute itsclearly nonlinear behavior to the fact that reconstructionin the large correlation time case requires a matrix inver-sion. We have found that, in the limit τ → τ → ∞ , errors in de-tecting vanishing matrix elements occur due to the nec-essary matrix inversion. For the data shown in Fig. 2(b),we found only ten spurious edges, compared to 4944 realones and more than seven millions possible ones. Thisis particularly satisfying, given that coupling strengthsvary by an order of magnitude in PanTaGruEl [31, 32].One of our main results is that partial network recon-struction can still occur when one has access to partialmeasurements over a fraction of the network agents. Themethod is illustrated in Fig. 3 which plots noise correla-tors vs. τ for the network of Fig. 1(a) with homoge- Figure 3. Two point correlator of Eq. (7) between agent a ij = 1 if and only if i and j are connected.(b) Inhomogeneous edge weights with randomly distributed a ij ∈ [0 . , . τ = λ − n (left)and τ = λ − (right). neous (left panel) and inhomogeneous (right panel) edgeweights. For small enough τ , so that λ n τ <
1, indicatedby the left vertical dashed line in Fig. 3, noise correlatorscome in well defined bunches, indicating that Eq. (13) issatisfied. Geodesic distances are then clearly identified,as is seen in Fig. 3 where data colors correspond to thecolor-coded geodesic distances in Fig. 1(a).Figure 4 further illustrates partial network reconstruc-tion by showing a colorplot histogram of Eq. (8) for eachagent i = 1 , . . .
20. Bands corresponding to geodesicdistances are clearly identified, at least up to the thirdneighbor. While full reconstruction of the network re-quires observability of all n nodes, partial structures canbe inferred from geodesic distances between the observ-able nodes via the two-point noise correlator. The partialreconstruction method works well, however it is in prac-tice limited to identifying the first few neighbors. Thisis so because it requires short correlation times, and isbased on Eq. (8) which states that pairs of k th neighboragents ( i, j ) have a noise correlator given by (cid:104) δ ˙ x i δ ˙ x j (cid:105) = δω ( − τ ) k ( J k ) ij + O [ τ k +10 ( J k +1 ) ij ] . (15)For distant pairs of agents with large values of k , the two-point correlator therefore becomes smaller and smaller,until it becomes smaller than its statistical standard de-viation. For the networks we investigated we have foundthat geodesic distances at least up to k = 4 can be in-ferred in practice. Conclusion.
We have presented a network reconstruc-tion method based on the dynamics of its agents underunavoidable ambient noise from the network’s environ-ment. Our approach is scalable to large networks, it hasgood accuracy, and has the significant advantage of be-ing nonintrusive. In particular, compared to earlier ap-proaches [7–13], it does not require the ability to injectsignal at each network node. To fully reconstruct the
Figure 4. Reconstruction of the 1st to 4th neighbors fromthe two point correlator of Eq. (7) for the network shown inFig. 1(a). For each agent index i , and each correlator value,the color plot gives the number of agents j (cid:54) = i with thatcorrelator value at the smallest value of τ in Fig. 3(a). network, the approach requires that one is able to mea-sure the dynamics of all agents, however, when not allthose measurements are accessible, the method is stillable to determine the geodesic distance between pairs ofmeasurable agents. This still provides precious informa-tion on the network structure. We trust that our methodcould significantly improve existing approach of networkreconstruction with partial measurements based on mostlikely network estimates [21, 22]. Acknowledgments.
This work has been supportedby the Swiss National Science Foundation under grant2000020 − [1] M. E. J. Newman, SIAM Review , 167 (2003).[2] A.-L. Barab´asi, Network science (Cambridge UniversityPress, Cambridge, England, 2016).[3] M. Hilbert and P. Lopez, Science , 60 (2011).[4] V. Sekara, A. Stopczynski, and S. Lehmann, Proc. Natl.Acad. Sci. USA , 9977 (2016).[5] J. Machowski, J. W. Bialek, and J. R. Bumby,
PowerSystem Dynamics , 2nd ed. (Wiley, Chichester, U.K,2008).[6] D. Bray, Science , 1864 (2003).[7] D. Yu, M. Righero, and L. Kocarev, Phys. Rev. Lett. , 188701 (2006).[8] M. Timme, Phys. Rev. Lett. , 224101 (2007).[9] D. Yu and U. Parlitz, Europhys. Lett. , 48007 (2008).[10] Z. Dong, T. Song, and C. Yuan, PLoS ONE , e83263(2013).[11] F. Basiri, J. Casadiego, M. Timme, and D. Witthaut,Phys. Rev. E , 012305 (2018). [12] S. Furutani, C. Takano, and M. Aida, IEICE Trans.Commun. E102-B , 799 (2019).[13] R. Delabays and M. Tyloo, (2020), arXiv:2002.00490[math.DS].[14] J. W. Pierre, N. Zhou, F. K. Tuffner, J. F. Hauer, D. J.Trudnowski, and W. A. Mittelstadt, IEEE Trans. PowerSyst. , 835 (2010).[15] L. Dosiek, N. Zhou, J. W. Pierre, Z. Huang, and D. J.Trudnowski, IEEE Trans. Power Syst. , 779 (2013).[16] V. A. Makarov, F. Panetsos, and O. de Febo, J. Neurosci.Methods , 265 (2005).[17] S. G. Shandilya and M. Timme, New J. Phys. , 013004(2011).[18] M. J. Panaggio, M.-V. Ciocanel, L. Lazarus, C. M. Topaz,and B. Xu, Chaos , 103116 (2019).[19] R. Dahlhaus, M. Eichler, and J. Sandk¨uhler, J. Neurosci.Methods , 93 (1997).[20] K. Sameshima and L. A. Baccal´a, J. Neurosci. Methods , 93 (1999).[21] M. E. J. Newman, Nature Physics , 542 (2018).[22] T. P. Peixoto, Phys. Rev. Lett. , 128301 (2019).[23] A. Mauroy and J. Hendrickx, SIAM J. Dyn. Syst. ,479 (2017).[24] R. Shi, W. Jiang, and S. Wang, Chaos , 013138 (2020).[25] M. G. Leguia, C. G. B. Mart´ınez, I. Malvestio, A. T.Campo, R. Rocamora, Z. Levnaji´c, and R. G. Andrzejak,Phys. Rev. E , 012319 (2019).[26] A. Banerjee, J. Pathak, R. Roy, J. G. Restrepo, andE. Ott, Chaos , 121104 (2019).[27] W.-X. Wang, Y.-C. Lai, and C. Grebogi, Phys. Rep. , 1 (2016).[28] I. Brugere, B. Gallagher, and T. Y. Berger-Wolf, ACMComput. Surv. , 24 (2018).[29] M. Tyloo, T. Coletta, and P. Jacquod, Phys. Rev. Lett. , 084101 (2018).[30] Y. Kuramoto, Prog. Theor. Phys. Suppl. , 223 (1984).[31] L. Pagnier and P. Jacquod, “PanTaGruEl - a pan-European transmission grid and electricity genera-tion model (Zenodo Rep.),” https://doi.org/10.5281/zenodo.2642175 (2019).[32] M. Tyloo, L. Pagnier, and P. Jacquod, Sci. Adv.5