Global Network Prediction from Local Node Dynamics
GGlobal Network Prediction from Local NodeDynamics
Neave O’Clery a , Ye Yuan b , Guy-Bart Stan c , and Mauricio Barahona c a University of Oxford; b Huazhong University of Science and Technology; c Imperial College LondonThis manuscript was compiled on September 5, 2018
The study of dynamical systems on networks, de-scribing complex interactive processes, provides in-sight into how network structure affects global be-haviour. Yet many methods for network dynamicsfail to cope with large or partially-known networks, aubiquitous situation in real-world applications. Herewe propose a localised method, applicable to a broadclass of dynamical models on networks, whereby indi-vidual nodes monitor and store the evolution of theirown state and use these values to approximate, via asimple computation, their own steady state solution.Hence the nodes predict their own final state with-out actually reaching it. Furthermore, the localisedformulation enables nodes to compute global networkmetrics without knowledge of the full network struc-ture. The method can be used to compute globalrankings in the network from local information; todetect community detection from fast, local transientdynamics; and to identify key nodes that computeglobal network metrics ahead of others. We illus-trate some of the applications of the algorithm byefficiently performing web-page ranking for a largeinternet network and identifying the dynamic roles ofinter-neurons in the
C. Elegans neural network. Themathematical formulation is simple, widely applicableand easily scalable to real-world datasets suggestinghow local computation can provide an approach tothe study of large-scale network dynamics.
We live in an interconnected world. From social networksto environmental sensors and autonomous robots, networkedsystems have become ubiquitous in our everyday lives. Inorder to describe such complexity, we can employ mathemati-cal tools to model the behaviour of agents—whether they behuman, organism or machine—as they communicate, processand integrate information in a networked environment. For in-stance, opinion formation evolves over time based on the belieflevel, with each individual updating their own opinion basedthat of their contacts. Eventually, as information spreads, aconsensus emerges for the entire network as confidants embracethe emerging common opinion. Simple consensus models (16),whereby nodes converge to a common value or opinion, fallwithin the general class of dynamical linear systems (1) andare employed in the study of social, biological, robotic andmanufacturing systems (2–4) with applications in distributedsensors (26), social networks (27, 28) and synchronisation(5, 29, 30).A key focus of recent research has been understanding andexploiting the connection between network structure and dy-namics (5–7). A well-known application that exploits dynamicsto elucidate information about network structure is that ofnode ranking, which forms the basis for Google’s PageRank(8). Given a network of webpages connected by links, nodes iteratively update their own value by averaging the rankings ofneighbouring nodes under a dynamical model. The final noderanking is derived from computation of the long-run steadystate solution of this dynamical system. Other examples in theengineering and control literature focus on the ability of cer-tain nodes to influence or decide the outcome of a dynamicalsystem on a network (9). Such analysis is critical in the fieldof distributed and cooperative robotics, where teams of au-tonomous robots or vehicles must co-ordinate their movement,and respond to external stimuli (10). Related models are alsoused extensively in the analysis of social networks and opinionformation (11). These examples can be seen as applicationsof generalised linear dynamics on a network, where all nodesevolve towards a final steady state that is determined by globalinformation contained in the network (12).Despite much progress in this field, many existing tools andalgorithms at the network-dynamics interface are ill-equippedto cope with the current explosion in data availability and net-work size, as well as practical, design and security challenges.For example, researchers today are regularly faced with net-work data containing information on millions, if not billions,of nodes. Examples include social media networks (Facebookcurrently has 1.5 billion active users) (13), the internet (2), andbiological systems such as protein interaction networks (14)and neural networks (15). In order to compute the steady statevalue of the dynamics on these networks, the full dynamicalsystem is typically simulated until its equilibrium behaviouris attained. However, due to the fact that the convergenceprocess can be very slow (16), dynamical simulation is not a To whom correspondence should be addressed. E-mail: [email protected]
SIGNIFICANCE STATEMENT
Many technologies today rely on networks to describecomplex interactions between groups of agents. Dynamicalprocesses on networks, whereby each node evolves over timebased on the state of its neighbours, have a large number ofapplications, including opinion models, environmental sensingand internet search. However, the analysis of large networkeddynamical systems is computationally intensive, and requirespotentially unavailable knowledge of the global networkstructure. Here we propose an efficient localised approach,whereby each node computes its own final equilibriumvalue using only local information collected from its ownhistory. We show how this methodology can be used forweb-page ranking, community detection, and identification ofkey communicator nodes in neural networks in a local manner.
September 5, 2018 | a r X i v : . [ phy s i c s . s o c - ph ] S e p lways computationally feasible. Additionally, for many real-world practical applications, such as the internet, biological,finance and sensor networks, the full network structure is oftenunknown due to measurement difficulties, security concerns,and/or physical constraints. In these cases, individual agentsor nodes may only have access to a limited amount of locallycollected data. Hence it is often desirable, if not necessary,to employ limited local information (i.e., information on asingle node, or subset of nodes) to extract both local and global network characteristics.Here we propose a fully localised method, whereby indi-vidual nodes monitor the evolution of their own state as itchanges, and use these values of its dynamical history to per-form a simple computation that allows them to approximatetheir own long-term dynamics. This approach enables nodesto ‘predict’ their own final state, and in some cases that of thewhole network, well before the dynamics actually converge toan equilibrium state.Our method builds on work of Sundaram and Hadjicostis(17) who originally proposed a method to compute the finalstate using a sequence of initial state values of length equalto the rank of the observability matrix (a matrix constructedfrom powers of the adjacency matrix, see (18) for an overview).Yuan et al. (19) proposed an analogous but localised approachby replacing the observability matrix with a local Hankelmatrix containing only the history of state values of each node.While theoretically interesting, both approaches suffer fromthe fact that the observability (or, equivalently, Hankel) rankis typically large (usually equal to the network size) for mostreal-world networks, limiting their practical usefulness andapplicability.Here we show that a relaxation of the localised Hankelapproach generates a sequence of approximations to the steadystate value for each node. Specifically, for each node and ateach time step, we compute the singular value decompositionof a Hankel matrix composed of preceding state values, anduse it to approximate the Hankel nullspace vector to computean approximate steady state value for each node. Critically,this method, termed below as the ‘Hankel method’, does notrequire knowledge of the full network topology, and enablesnodes to not only ‘predict’ their own long-term dynamicalequilibrium or steady state, but also in many cases that ofthe full network. The localised and dynamical nature of thealgorithm has several desirable features. The Hankel method is fast compared to traditional methods,and scalable to large systems.
The sequence of Hankel approx-imations is guaranteed to converge in a number of steps lessthan or equal to the Hankel rank, but in practice typicallyaccurately approximates the steady state value in very fewsteps. For example, for a random (Erdős-Rényi) network ofone thousand nodes, only three to four steps are typicallyneeded to approximate the final value for each node. In con-trast, thousands of steps may be required for convergence ofthe iterative dynamics to the same accuracy. We show thatthis approach can be used to construct a localised algorithmto compute the Google PageRank vector (8), and we illustratethis for both the well-known Karate Club social network (20),and a larger network of over 10 ,
000 webpages within the 2002
Stanford.edu domain (21).
The Hankel method is node-specific.
The process of predict-ing global information from local information is not equally possible for all nodes, i.e., some nodes can compute their ownlong-run equilibrium value using fewer of their own initialvalues. Under a consensus model (16), whereby nodes con-verge to a common value or opinion, certain nodes are thefirst to be able to compute the full network steady state (andcan communicate this information to other nodes if needed).These nodes can been seen as knowledgeable ‘network insiders’.We illustrate this phenomenon by considering the well-known
C. Elegans neural network, which describes the chemical andelectrical wiring structure of sensory, motor and inter-neuronsfor this worm (15, 22). Assuming a simple dynamics underwhich neurons receive and assimilate information from theirneighbours, we observe that inter-neurons tend to requirefewer values of their history to compute the global outcomeof the system dynamics, consistent with their role as keycommunicators in the network.
The Hankel method is generalisable and adaptable.
Its math-ematical formulation is simple, and encompasses a large classof widely-used linear models for network dynamics. Hence,we can exploit its localised and computational efficiency fora range of applications. As illustrations, we show below thateach node can not only predict its own final value or rank-ing, but also detect clustering of intermediate values beforean equilibrium state is reached. This behaviour is relevantto the detection of modular structure in networks based onthe dynamics of random walker travelling from node to node.Transient clustering of the dynamics is indicative of modularcommunity structure, as random walkers become trapped indynamical ‘basins’ (6, 23, 24). Such analysis can illuminate, forexample, functional groups of proteins in proteome networks,which can be associated with different diseases or biologicalprocesses (25).
Localised, Finite-time Computation of Network SteadyStates
We will consider a generalised model of consensus dynamicswhereby nodes can both store their own ‘value’ or opinionat any point in time, and periodically learn their neighboursvalues in order to ‘update’ their own opinion. Mathematically,if x k ∈ R n is a vector containing the value or state of eachnode at time-step k , all future node values may be describedby the iterative process x k +1 = W x k [1]where x ∈ R n contains the initial condition for the dynamics,and the weight matrix W ∈ R n × n describes the inter-noderelationships giving rise to the dynamical system. All nodesasymptotically (i.e. eventually) converge to the same (con-sensus) steady state solution under two conditions: W = (where corresponds to a vector of ones) and all other eigen-values of W are within the unit disk.Depending on the situation, W can take several forms. Oneof the most common forms is that of the Laplacian consensusdynamics (16), which describes the process by which a node ofa network updates its values based on that of its neighbours.In this case, W = I − ωL , where L = D − A is the Laplacianmatrix with A being the adjacency matrix of the networkwith entries i, j containing the weight of the edge betweennode i and node j , and D being a matrix of zeros everywhere,except on the diagonal where it contains the correspondingnode out-degree. For a directed network, if (i) the graph is ig. 1. Predicting the final value of Laplacian consensus dynamics. (A)
Laplacianconsensus dynamics of the network shown in the inset for a random initial condition(time traces correspond to node colour). The iterative dynamics reaches the consen-sus value slowly beyond ∼
60 steps. (B)
The sequence of approximate consensusvalues computed via the Hankel method
Eq. (3) for each node converges to thetrue consensus value in less than 10 steps (detected with tolerance (cid:15) = 10 − ).The inset is a blow-up of the main figure showing that each node approximates itsfinal value in a different number of steps; in all cases less than r + 1) with ∆ r = { , , , , , , , } for nodes r = 1 , ..., , respectively Eq. (2) . balanced and strongly connected, and (ii) 0 < ω < / max( D ),then all variables in x k asymptotically reach a shared averageconsensus value. Another key example is that of random walkerprobability dynamics (24), where W = AD − and the finalstate vector represents the long-term transition probabilitiesof a random walker transitioning through each node. Randomwalker dynamics on a network forms the basis of a number ofalgorithms including Google’s PageRank (8), and communitydetection algorithms (6, 23).An example network is shown in the inset of Figure 1 A,where we simulate Laplacian consensus dynamics as describedabove, and observe that, despite the small size of the network( n = 8 nodes), the node dynamics do not converge to a con-sensus value (for a random initial condition) until step k ≈ W above, the rateof asymptotic convergence for such dynamics is upper bounded(16), yet there is no lower bound to this rate of convergence,i.e., there is no mathematical limit on the number of iterationsor steps this could take.Sundaram and Hadjicostis (17, 31) have shown that indi-vidual nodes, under certain conditions, can compute their ownfinal state value in a finite number of steps equal to the rankof a node-specific observability matrix (composed of matrixpowers of W , see SI and (18) for a review). However, in mostcases, the observability rank is equal to the size of the network,and nodes require knowledge of the full network structure—anassumption not likely to be met for many real-world applica-tions. Yuan et al. (19) proposed an analogous but localised methodology, employing a ‘local’ Hankel matrix in place ofthe ‘global’ observability matrix. Hankel matrices are squarematrices in which each ascending skew-diagonal from left toright is constant, and are frequently employed in control the-ory. Specifically, if x i ( r ) for i = 0 , , , . . . are the state valuesassociated with node r , then the Hankel matrix H ( r ) k ∈ R k × k of size k for node r is defined as H ( r ) k = x ( r ) − x ( r ) . . . x k ( r ) − x k − ( r ) x ( r ) − x ( r ) x k − ( r ) − x k − ( r )... ... x k ( r ) − x k − ( r ) . . . x k − ( r ) − x k − ( r ) . This matrix is increased in size (via addition of rows andcolumns) until its singular value decomposition (32) includes asingle zero eigenvalue and it has dropped rank at k = ∆ r + 1,i.e., rank (cid:16) H ( r ) k (cid:17) = ∆ r . The steady state value can be thencomputed by node r using the formula (see (19)): h ∗ ( r ) = [ x ( r ) . . . x ∆ r ( r )] v ( r )∆ r +1 T v ( r )∆ r +1 . [2]where v ( r )∆ r +1 ∈ R ∆ r +1 is the single nullspace vector of H ( r )∆ r +1 and x ( r ) , . . . , x ∆ r ( r ) are successive values of node r . A proofmay be found in (19), and is based on a Jordan decomposition(32) of W and a Vandermonde decomposition (32) of H ( r )∆ r +1 .When identical symmetry exists both in the graph and initialcondition, the rank of H ( r ) k does not accurately reflect thenumber of steps to compute the final value (19, 33). Yet, for arandom initial condition this is highly unlikely to be an issue.As mentioned above for the observability method, the num-ber of steps required for each node, given by the correspondingHankel rank, of most common classes of complex networksis usually trivial (i.e., it equals the graph size n ) (34). Thisimplies, firstly, that there is a prohibitive cost for very largenetworks with millions of nodes in computing the steady statevalue using the above approach, and, secondly, that all nodescompute the final value in the same number of state values—there is no node-specific predictive advantage. There is alsoa significant limitation to the accuracy of this approach asthe graph size grows due to the fact that the entries of theHankel matrix tend towards zero as the system converges tothe consensus value, i.e., x k ( r ) − x k − ( r ) → k → ∞ . De-termining the matrix rank accurately is therefore numericallyunstable for large matrices (7), and particularly so for thosecomposed of small numbers.We exploit the properties of the singular value decompo-sition to propose a relaxation of this approach. Specifically,we use the singular vector corresponding to the smallest sin-gular value σ to approximate the Hankel nullspace vectorfor increasing Hankel size k . More specifically, we compute asequence of approximations to the steady state solution fornode r for steps k ≤ ∆ r + 1 such that h ∗ k ( r ) = [ x ( r ) . . . x k − ( r )] v ( r ) k T v ( r ) k [3]where v ( r ) k is the singular vector corresponding to the smallestsingular value of the Hankel matrix H ( r ) k . As the distance tothe singular matrix decreases (35), this sequence of approxi-mates approaches the true solution. By construction we areuaranteed (in the worst case) to compute the steady statevalue in a maximum number of steps given by k = ∆ r + 1,but, in practice, we typically have k (cid:28) ∆ r + 1.Figure 1 highlights the contrast between the actual numberof steps needed for (cid:15) -convergence of the dynamical iteration(i.e., the dynamic simulation is within an (cid:15) -distance of thetrue value, or | x k ( r ) − x ∗ ( r ) | ≤ (cid:15) ), and the significantly fewernumber of steps needed for (cid:15) -convergence of the Hankel ap-proximation to the final value (i.e., | h k ( r ) − x ∗ ( r ) | ≤ (cid:15) ). Forpractical applications, these state values could be acquiredfrom sensor measurements. For most illustrations below, wecompute these initial state values via a short simulation of thedynamical system given by Eq. (1). We then show that thenumber of steps required to compute the final value of eachnode, in almost all cases, is significantly less than the numberof steps needed for convergence of dynamical system.We highlight several important features of this example:• Convergence of the Hankel approximation is guaranteedand upper-bounded by the Hankel rank of each node,yet in the consensus dynamics example of Figure 1, thenumber of steps given by the upper bound was neverrequired—and was significantly less than the number ofsteps needed for the corresponding dynamical iteration;• Each node exhibits a distinct number of steps requiredto compute the common consensus value—hence somenodes are in a sense more ‘knowledgeable’ than others,and can compute or predict the final value sooner;• In the case of consensus models such as the one consid-ered in Figure 1, a single node, without any informationother than a short sequence of its own state values, candetermine the final consensus value of the entire network.These powerful features, and their potential wide-rangingapplications, will be explored in the following section. Node-specific Predictive Capability in Complex Networks .
The example in the previous section highlights the fact thateach node of the network requires a unique number of initialstate values to compute its own steady state solution, andthat some nodes are better predictors in the sense that theyrequire fewer values. Here we investigate this result by ap-plying our approach to the well-known
C. Elegans neuronalnetwork which describes the wiring structure of sensory, motorand inter-neurons for the
Caenorhabditis Elegans worm.
C.Elegans has been used extensively as a model organism tostudy neuronal development (15) as it is one of the simplestorganisms with a nervous system, and is easy to grow in bulkpopulations.The
C. Elegans network is shown in Figure 2 A, with sen-sory (green), motor (grey) and inter- (orange) neurons locatedin distinct regions of the network. Key neurons, such as in-terneurons responsible for mediating signals between inputsensory neurons and motor neurons, are labelled explicitlyon the network. The edges represent a combination of gapjunction (electrical) connections and chemical synapse con-nections. Following (22) nodes have been positioned withprocessing depth (i.e., the number of synapses from sensoryto motor neurons) on the y-axis, and the node’s respectiveentry in the Fiedler eigenvector on the x-axis. The Fiedlervector is the eigenvector corresponding to the second smallest
Fig. 2.
Node-specific predictive capability. (A)
The well-known neuronal network ofthe
Caenorhabditis Elegans worm represents the synaptic connections and chemicaljunctions between sensory (green), motor (grey) and inter- (orange) neurons (22). Thelayout follows Varshney et al (22): the y -coordinate is given by the processing depth,and the x -coordinate is the node coordinate of the normalised Fiedler eigenvector. (B) Each node of the network is coloured according to the number of Hankel steps itneeds to approximate its own steady state value to a small tolerance ( (cid:15) = 10 − ). Thevalues are averaged over dynamics started from random initial conditions in theinterval [ − , . We observe that inter-neurons typically compute their steady-statevalue first. (C) There is a statistically significant difference in the number of Hankelsteps needed by the three neuron groups. eigenvalue of the Laplacian matrix of a graph. The entries ofthis vector have previously been used to partition the networkinto two densely connected groups of nodes (36) (positivevalues of the Fiedler vector in one group; negative values inthe other), and for spectral embedding, i.e., for finding anoptimum one-dimensional representation of a graph (37).We model the communication abilities of the neurons viaa simple consensus model, where neurons (nodes) communi-cate with their neighbours in the network via chemical andlectrical signalling (22, 38). We seek to identify nodes which,after receiving a limited number of signals from their neigh-bours, can predict the final consensus value of all nodes orneurons in the fewest steps. These nodes, we propose, arehighly knowledgeable—in the sense that they are influentialcommunicators—about the network dynamics given their po-sition in the network. Figure 2 B shows the network colouredby the number of Hankel steps needed by each node to approx-imate its own steady state value (the consensus value of thenetwork). These values are averaged over 5,000 random initialconditions in the interval [ − , (cid:15) = 10 − . We observe that inter-neurons, such as AVER/L,AVKL and RIGL (associated with locomotion in response tostimuli) and sensory neurons ADEL/R located in the head,typically compute their steady-state value first. Overall, inter-neurons exhibit a statistically significant decrease in numberof consecutive values needed to approximate the final steady-state value, shown via comparing the three distinct neurongroups in Figure 2 C.Hence, key inter-neurons, having collected a small numberof their own state values based on signals from their neighbours,can compute the common consensus value of the whole network,and broadcast or communicate it to other nodes who have notyet been able to compute this consensus value. In the light ofthis, it appears that inter-neurons have an increased predictiveability compared to more peripheral nodes, consistent withtheir functional role as key communicators between sensoryand motor neurons in the network. Global Node Ranking Based on Local Information.
Beyondthe consensus framework, the Hankel approximation methodmay be applied to a wide range of linear models. Here weconsider node ranking, which forms the basis of the Googlesearch engine technology (8). Analogous to the previous exam-ple, using our approach, individual nodes can compute theirown ranking while employing significantly fewer iteration stepsthan traditionally needed—and without the requirement ofknowing the full network structure.The classic PageRank algorithm (39) is based on a model ofa random walker moving from node to node along the edges ofa network formed by hyperlinked websites. A node’s rankingcan be seen as the long-run probability of the walker traversingthat node relative to other nodes. The PageRank vector (i.e.,the node ranking vector) is obtained as the solution of thelinear system given in Eq. (1) with W = αAD − + − αn E where E is a matrix of ones. The second term may be seen asa cost or damping term for which choices of α close to 1 yieldthe most accurate ranking—yet are more computationallyexpensive in the sense that convergence is slower. In essence,a value of α less than 1 adds constant edges to the networkenabling information to diffuse and spread more quickly. Forsmaller values of α , these edges are more heavily weighted andthe system reaches an equilibrium state faster.While details of the state-of-the-art Google ranking algo-rithms are unavailable, for large networks (the largest beingthe whole internet), webpage ranking employing this classicalgorithm is normally approximated via 50 −
100 iterationsof the full system, and has been reported to take a few daysto complete each month (40). Here we show that, for botha well-known social network and a large web network, wecan use our Hankel method to compute the ranking valuefor each node in a localised and efficient manner compared to the traditional approach of iteration of the full dynamics.Figure 3 A illustrates the well-known Karate Club network(20). This is a social network of friendship links between 34members of a karate club at a US university in the 1970s. Adisagreement between the club’s president and main instructorcreated a split between the members, and ultimately led tothe breakup of the club. The nodes are coloured according tothe true ranking (dark blue corresponds to high ranking). InFigure 3 B, for a random initial condition (corresponding to arandom initial ranking shown on the left y-axis), we show thenode ranking for each step of the Hankel ranking algorithm.We observe that the final ranking is achieved by most nodesby step k = 6, and all nodes by step k = 8. The presidentand instructor are the top-ranked nodes, using their densefriendship links to other club members to cement their positionin rival camps.In order to illustrate the efficiency of our method on alarger scale, we compare the number of steps needed to ob-tain the node ranking of a large web network, the undirectedlargest connected component of the Stanford web networkwhich describes hyperlinks between almost 10 ,
000 web pagesin the domain stanford.edu for the year 2002 with n = 8929nodes, and over 25 ,
000 edges (21). Using a damping factorvalue of α = 0 .
9, we compute the node ranking via both theHankel approximation method and, for comparison, iterationof the full dynamical system. In order to compare the speedof each of these approaches, we compute the Spearman rankcorrelation (41) between the approximate ranking (Hankel ordynamical iteration) and the true ranking (previously com-puted via a long iteration) at each step. Figure 3 C showsthat our Hankel approximation can obtain the node ranking(i.e., the correlation reaches 1) in significantly fewer steps thanthe linear dynamic iteration.For large networks, the quality of this approximation de-pends on the choice of damping factor α (42)—the closer α is to 1, the higher the accuracy of the approximation butalso the slower the convergence rate (43) ( α = 0 .
85 is themost commonly used value in the literature (39)). Figure 3 Dshows that, as we increase α and, thereby, the accuracy of theranking, both the iteration and the Hankel method requiremore steps to converge. The inset shows that the ratio ofHankel steps to dynamics steps decreases with increasing α ,implying that the Hankel method gains in relative efficiencyas the damping increases.These results are important as many applications todaycontain millions, or sometimes even billions, of nodes. Thepotential to both locally and efficiently determine metrics suchas node rankings, or equivalently a variety of centrality ordynamics-based measures, using our newly proposed method,could render previously intractable problems solvable. Community Detection from Local Node Transients.
Beyondconsensus models and node ranking for networks, linear sys-tems models form the basis for a large class of computationalalgorithms for the analysis of network structure. Due to thelarge-scale, complex nature of real-world networks, the de-tection of communities or groups of nodes, typically tightlyconnected, can yield powerful insights into network behaviourby revealing the underlying organisation of the network, orproviding insight into its function. Furthermore, network sizemay be reduced via aggregation of nodes in communities, oftenyielding more a tractable and informative topology (6, 7, 44). ig. 3.
Node ranking from local information. (A)
The network of the
Karate Club social network (20) represents a friendship network within a US karate club dominated by twocliques—one allied to the
President , and the other to the
Instructor . The nodes are coloured by their node ranking (darker nodes are ranked higher) confirming the influence ofthese two key actors. (B)
Starting from a random initial condition (the random initial ranking on the left y-axis), a new node ranking can be obtained for each Hankel step(x-axis). By step 8, all nodes converge to their final ranking (colours of the lines as in A). (C)
The node ranking of the largest connected component of the
Stanford webnetwork ( n = 8929 nodes, inset) (21) is computed using the Hankel approximation (with damping factor α = 0 . , blue circles), and with the full iteration method (grey circles)starting from a random initial condition. Spearman correlation of both rankings against the ground truth (i.e., the Pagerank achieved asymptotically). The Hankel approachproduces an accurate ranking (i.e., the correlation reaches ) in significantly fewer steps than the convergent dynamics. (D) The number of steps required by the Hankelmethod compared to the full iteration method when the damping factor α is varied between . and . (shown is an average over random initial conditions for eachvalue of α with tolerance (cid:15) = 10 − ). For higher values of α the approximation is closer to the ‘true’ ranking for an undamped system with α = 1 . Hence both the iteration andthe Hankel method require more steps to converge. (Inset) The ratio of Hankel steps to full dynamical iteration steps decreases with increasing α , i.e., the Hankel method isincreasingly more efficient as we approximate the undamped network ranking. Many algorithms have been proposed for the analysis ofcommunity structure in graphs (3, 45) including normalisedcut (46), modularity (47, 48) and stability (6, 24). Many ofthese algorithms are based on the idea that a random walkeron a network (i.e., a walker that travels from node to node)becomes ‘trapped in wells’, circulating within sub-regions ofthe network. This can be due, for example, to a region ofhigh connectivity leading the walker to repeatedly traverse thesame set of nodes for an extended period of time.The probabilistic dynamics of a random walker on a graphmay be modelled via the construction of a linear system simi-lar to that introduced above in the webpage ranking example(6, 24), i.e., W = AD − . Transient clustering of these dynam-ics provides evidence of community structure in the networkas the dynamics of nodes in the same community (cid:15) -convergetemporarily (i.e., they exhibit similar state values) beforereaching a (possibly but not necessarily common) steady state.We find that beyond approximating the steady state solution,our Hankel approach can also detect these transient states.This important feature of our method means that, even before computing the final state for each node, Hankel approxima-tions converge for nodes within the same community, therebyallowing communities to be very quickly identified during theHankel iterations. Hence, as we will see, applying the Hankelmethod enables us to also detect community structure in alocalised manner, and using significantly less successive statevalues than existing methods.To illustrate this result, we generate an ensemble of net-works with community structure defined by a matrix of prob-abilities P such that P l,m is the probability of connectionbetween any node in community l and any node in community m . Figure 4 A shows the entries in the adjacency matrix ofsuch a network, with n = 200 nodes split into four equal sizecommunities (each with 50 nodes). In this case, the proba-bilities of node connection between communities is given by P l,m = 0 . l = m (i.e., nodes in the same community)and P l,m = 0 .
01 for all l = m .For each network we compute the iterative dynamics (using W = AD − and a random initial condition in the interval[0 ,
1] for Eq. (1)), and compute a sequence of Hankel values as ig. 4.
Community Detection. (A)
Adjacency matrix of a network with n = 200 nodes and four equal sized communities (the probability of connection for node pairs within thesame community is P i,j = 0 . , tand across different communities is P i,j = 0 . ). (B) The Hankel approximation method detects transient clustering of node dynamics,corresponding to the detection of community structure. We visualise the distance matrix D k , where entries correspond to the difference between the Hankel approximations foreach node pair at step k . The blocks of small values (dark shading) emerging at steps , and correspond to the transient clustering of the Hankel approximation withincommunities. By step 5, the iteration has converged and all nodes have the same value ( D k has only small values). (C) The corresponding dynamical matrices S k , whereentries correspond to the distance between dynamical iteration values for node pairs at step k , also detects communities, but at much longer times than the Hankel approach(e.g., community structure only starts to appear at around step ). While the Hankel method computes the final consensus value in just steps, the full system dynamicstakes up to steps to converge. defined by Eq. (2). We seek to detect ‘distances’ between theHankel values, and the iterative dynamics, for pairs of nodesat each step. We compute distance matrices D k and S k ateach step k of the iteration, with entries D k ( i, j ) = | h ∗ k ( i ) − h ∗ k ( j ) | S k ( i, j ) = | x k ( i ) − x k ( j ) | for all pairs of nodes i and j , where h ∗ k ( i ) is defined in Eq. (3)and x k ( i ) is the ith component of the vector x k in Eq. (1).Groups of nodes exhibiting similar dynamics (as captured bysmall distance values) signal community structure.In Figure 4 B and C, we visualise the matrices D k and S k for increasing values of the iteration number k . In both cases,dark shading indicates small distance values—and hence thedetection of community structure. We observe that the Hankelapproximation transiently detects four communities within just k = 3 steps, before closely approximating the final steady statevalue at k = 5. In contrast, the full dynamics only detect thecommunity structure after around k = 20 steps, and does notconverge to the consensus value until k > Discussion
Using local information to obtain global properties from in-terconnected dynamical systems. our work exhibits a numberof distinguishing features. Its localised formulation, wherebynodes need only have access to a limited number of their ownsuccessive state values, enables nodes or sensors to indepen-dently compute variables of interest, without knowledge ofthe network structure or of their neighbouring nodes’ state values. Its formulation is simple, and is easily adaptable toa large class of network-based dynamical models. It can bescaled to cope with large real-world systems, with significantefficiency compared to commonly used iterative methods as itavoids the need for convergence. Finally, it enables us to iden-tify functional attributes of nodes in networks, i.e., ‘predictor’or ‘communicator’ nodes that can estimate the full networkdynamics ahead of other nodes.Our algorithm performs well for a range of networks andapplications; yet alternative but related approaches couldbe explored for future work. While we have employed thesingular vector corresponding to the smallest singular valueas an approximate nullspace vector, it may be possible to usea combination of singular vectors defining a closest subspacein order to obtain ‘smoother’ convergence to the consensusvalue. Furthermore, computing the SVD or rank updatefor the Hankel matrix as new rows and columns of data areadded (without re-computing the decomposition) in the lineof thought explored by (49, 50) would be beneficial. A closed-form solution for the updated singular vector dependent onlyon the previous step and the new data does not appear possibleusing current techniques, but could be feasible under a re-formulation or further relaxation of the problem.For any such algorithm, it is desirable to have a fully lo-calised convergence criterion that can be computed by eachnode locally to ascertain when the consensus value is wellapproximated. In the SI, we show that the difference betweensuccessive steps of the algorithm is highly correlated with the‘distance’ to the true consensus value; hence this criterion canbe used as an effective stopping criterion for large graphs. De-veloping a theoretical convergence criterion would be beneficial,although this is a non-trivial task due to the non-monotonicityf the sequence of steady-state value approximations.The number of steps needed for an individual agent toapproximate its own final value is node-specific, and nodesthat require the fewest steps can be seen as key communicatorsin a network. In the case of consensus dynamics, a subsetof influential nodes can ‘predict’ the final value earlier, andpotentially communicate that result to other agents in thenetwork. There is much scope to investigate further howthese results could be used dynamically by nodes to aid ordisrupt consensus, and to develop further applications withinthe context of network design and autonomous sensing.
Materials and Methods • The adjacency matrix, node type and node position (xand y axis co-ordinates) for the
C. Elegans ∼ mejn/netdata/, and visualisedwith Gephi (https://gephi.org/).• The adjacency matrix for the 2002 Stanford web network (21) inFig. 3C downloaded from https://snap.stanford.edu/data/web-Stanford.html, and visualised with Gephi (https://gephi.org/). ACKNOWLEDGMENTS.
N.O’C. was funded by a Wellcome TrustDoctoral Studentship. G-B.S. acknowledges the support of EP-SRC through the EPSRC Fellowship for Growth EP/M002187/1.M.B. acknowledges support from EPSRC Grants EP/I017267/1 andEP/N014529/1.
1. Newman M, Barabasi AL, Watts DJ (2011)
The structure and dynamics of networks . (Prince-ton University Press).2. Strogatz SH (2001) Exploring complex networks.
Nature
SIAM Review
IEEE ControlSystems
Physical ReviewLetters
Proceedings of the National Academy of Sciences
Physical Review E
Com-puter networks and ISDN systems
Nature
Graph theoretic methods in multiagent networks . (PrincetonUniversity Press).11. Schaub MT, Delvenne JC, Lambiotte R, Barahona M (2018) Structured networks and coarse-grained descriptions: a dynamical perspective.
ArXiv e-prints .12. Schaub MT, Delvenne JC, Lambiotte R, Barahona M (2018) Multiscale dynamical embed-dings of complex networks.
ArXiv e-prints .13. Ellison NB, Steinfield C, Lampe C (2007) The benefits of facebook “friends:” social capitaland college students’ use of online social network sites.
Journal of Computer-Mediated Com-munication
Physical Biology
Nature
Proceedings of the IEEE
IEEE American Control Conference pp. 711–716.18. Liu YY, Slotine JJ, Barabási AL (2013) Observability of complex systems.
Proceedings of theNational Academy of Sciences
Automatica
Journalof Anthropological Research
Internet Mathe-matics
PLoS Computational Biology
Physics reports
ArXiv: 0812.1770 .25. Jonsson PF, Cavanna T, Zicha D, Bates PA (2006) Cluster analysis of networks generatedthrough homology: automatic identification of important protein communities involved in can-cer metastasis.
BMC bioinformatics
IEEE Fourth International Symposium on Information Processing in SensorNetworks .27. Krause U (1997) Soziale dynamiken mit vielen interakteuren. eine problemskizze.
Model-lierung und Simulation von Dynamiken mit vielen interagierenden Akteuren, Bremen Univer-sity pp. 37–51.28. Sood V, Redner S (2005) Voter model on heterogeneous graphs.
Physical Review Letters
Proceedings of the 42nd IEEE Conference on Decision and Control pp. 2996–3000.31. Sundaram S, Hadjicostis CN (2007) Distributed consensus and linear functional calculationin networks: an observability perspective.
Proceedings of the 6th international conference oninformation processing in sensor networks pp. 99–108.32. Golub GH, Van Loan CF (2012)
Matrix computations . (JHU Press) Vol. 3.33. Hendrickx JM, Olshevsky A, Tsitsiklis JN (2011) Distributed anonymous discrete functioncomputation.
IEEE Transactions on Automatic Control
LINPACK users’ guide . (Siam) No. 8.36. Newman ME (2006) Finding community structure in networks using the eigenvectors of ma-trices.
Physical Review E
Discrete Ap-plied Mathematics
Nature
Princeton University Press .40. Bryan K, Leise T (2006) The $25,000,000,000 eigenvector: The linear algebra behind google.
Siam Review
Numerical recipes in C . (Cam-bridge Univ Press) Vol. 2.42. Boldi P, Santini M, Vigna S (2005) Pagerank as a function of the damping factor. .43. Kamvar SD, Haveliwala TH, Manning CD, Golub GH (2003) Extrapolation methods for accel-erating pagerank computations.
Proceedings of the Twelfth International World Wide WebConference .44. Yaliraki SN, Barahona M (2007) Chemistry across scales: from molecules to cells.
Philosoph-ical Transactions of the Royal Society of London A: Mathematical, Physical and EngineeringSciences
Physical ReviewE
IEEE Transactions on PatternAnalysis and Machine Intelligence
PhysRev E
Proceedings of theNational Academy of Sciences
SIAM Journal on Numerical Analysis
Computer Vision pp. 707–720.51. Ogata K (2009)
Modern Control Engineering . (Prentice Hall), 5th edition. upplementary Information
1. Theoretical Framework
We consider linear dynamics where x k ∈ R n is a vector containingthe state of each node at step k , x k +1 = W x k [4]for initial condition x . When W has a simple eigenvalue at 1and all the other eigenvalues within the unit disk this system isguaranteed to reach a steady state solution.Under certain conditions the final value or steady state may becomputed by each variable from a finite number of initial steps orstate values (17, 31). In this case, the number of steps needed isgiven by the polynomial degree—in particular the coefficients ofthe minimal polynomial for the graph. Below we briefly review themathematics behind this approach.The minimal polynomial of matrix W ∈ R n × n with respectto node r is the unique monic polynomial q r with minimal degree∆ r + 1 q r ( z ) = z ∆ r +1 + ∆ r X i =0 α ri z i = 0 [5]such that e r T q r ( W ) = , where e r ∈ R n is the vector of zeros witha single 1 in position r .By application of the Final Value Theorem: x ∗ ( r ) = lim k →∞ x k ( r ) = lim z → ( z − X r ( z ) = F r (1) p r (1)where X r ( z ), F r ( z ) and p r ( z ) are given by X r ( z ) = ∞ X k =0 x k ( r ) z − k F r ( z ) = ∆ r X j =0 x j ( r ) z ∆ r +1 − j + ∆ r X i =1 α ri i − X j =0 x j ( r ) z i − j , and p r ( z ) = 1 z − q r ( z ) = z ∆ r + ∆ r − X i =0 (1 + ∆ r X j = i +1 α rj ) z i . In this case x ∗ ( r ) is given by x ∗ ( r ) = [ x ( r ) . . . x ∆ r ( r )] w r T w r [6]where w r = P ∆ r i =1 α ri P ∆ r i =2 α ri ...1 + α r ∆ r ∈ R ∆ r +1 . All that is needed in this case to compute the final value x ∗ ( r )for node r is the state values for node r up to step ∆ r (∆ r + 1 stepsin total) and the coefficients α ri for i = 1 , .., ∆ r . The coefficients α ri can be computed by considering the system (cid:2) α r . . . α r ∆ r (cid:3) e r I e r W ... e r W ∆ r = 0 . This matrix, denoted Q (∆ r +1) r , is the discrete observability matrixwith ∆ r + 1 rows (51). By forming the discrete observability matrixwith respect to node r with k rows Q ( r ) k = e r I e r W ... e r W k − ∈ R k × n . Fig. 5. (A)
For an Erd˝os-Rényi random graph of size n = 100 and edge density f = 0 . , we observe the convergence of the Hankel approximation for each node r to the final value x ∗ ( r ) . (B) While the Hankel approximation converges to the finalvalue, so do the differences between successive steps. and increasing k until the matrix loses rank at k = ∆ r + 1, thecoefficients α ri can be obtained from the left kernel vector (17, 31). ∗ This approach however, is not decentralised as knowledge of thenetwork structure (e.g., W ) is required to compute the final value.Yuan et al. (19) developed a new approach that employs a Hankelmatrix—which does not rely on W —to compute the α ri terms andderive a fully decentralised analogous method for the approximationof x ∗ ( r ) (see (19) for more details).
2. Hankel Method Algorithm
In the main text we proposed an algorithm to compute a sequenceof approximations to the final value of a individual variable or nodeby considering the singular vector corresponding to the smallestsingular value of a Hankel matrix of increasing size.Specifically, if H ( r ) k = x ( r ) − x ( r ) . . . x k ( r ) − x k − ( r ) x ( r ) − x ( r ) x k − ( r ) − x k − ( r )... ... x k ( r ) − x k − ( r ) . . . x k − ( r ) − x k − ( r ) , [7]and v ( r ) k is the singular vector corresponding to the smallest singularvalue of H ( r ) k , then the Hankel approximation for the final value ofnode r at step k is given by h ∗ k ( r ) = [ x ( r ) . . . x k − ( r )] v ( r ) k T v ( r ) k . [8]The steps to compute a sequence of final value approximations fornode r include:1. Initialisation:
Set counter k = 1.2. Iteration:
For each k • Build the Hankel matrix H ( r ) k ∈ R k × k given by Eq. (7).• Compute the SVD and set the singular vector correspond-ing to the smallest singular value equal to v ( r ) k .• Calculate the approximation h ∗ k ( r ) given by Eq. (8).3. Termination:
If 1 < k ≤ n and | h ∗ k ( r ) − h ∗ k − ( r ) | < ν for sometolerance ν (see below for a discussion), then an approximatefinal value has been found, otherwise increment k = k + 1 andreturn to step 2 † .By construction we are guaranteed to compute the consensus valuein a maximum number of steps k = ∆ r + 1 ≤ n .
3. Rate of Convergence
For an arbitrary graph, we do not know a priori the final value of anynode, and hence we must define a convergence or stopping criterionfor the algorithm. In order to assess the convergence properties ofthis approach, we examine the relationship between the ’true’ error, ∗ Due to the Cayley-Hamilton Theorem, if Q ( r )∆ r = ∆ r then rank ( Q ( r )∆ r +1 ) = ∆ r . † For many of the examples in the main paper, where x ∗ was known, the criterion for countingthe number of Hankel steps required for the method to converge for the full network was || h ∗ k − x ∗ || < (cid:15) . ig. 6. We explore the the relationship between the approximated vs. true consensus value, | h ∗ k ( r ) − x ∗ ( r ) | , on the one hand, and the change in the approximatedconsensus value from one step to the next, | h ∗ k +1 ( r ) − h ∗ k ( r ) | , on the other. Here we observe good correlation between these quantities over variation in edge density forboth Erd˝os-Rényi and scale free graphs, illustrating the appropriateness of the latter as stopping criteria for our algorithm for various kinds of large graphs, and a range of keygraph parameters. Note, the paler colour dots correspond to the first few approximation steps (corresponding to higher values of the distances on both the x and y axis). A-C
Erd˝os-Rényi graphs of size n = 100 with edge densities f equal to . , . and . respectively. D-F
Scale-free graphs, which display variation in degree heterogeneity ascontrolled by parameter γ , of size n = 100 (with γ = 2 . ) with edge densities of . , . and . respectively. G-I
Scale-free graphs of constant size ( n = 100 ) and edgedensity ( f = 0 . ). We vary the degree heterogeneity γ = 2 . , . , . (most heterogenous to least heterogeneous). | h ∗ k ( r ) − x ∗ ( r ) | , and the convergence error, | h ∗ k +1 ( r ) − h ∗ k ( r ) | , for arange of large graphs with differing sizes, edge densities and degreedistributions as seen in Figures 5 and 6.In Figure 5, we consider an Erdős-Rényi random graph of size n= 100 and edge density f = 0 . γ (a parameter used to control thedistribution/heterogeniety of high and low degree nodes, middlerow), and scale free graphs with varying γ (bottom row). We observe a consistently strong correlation between the distance of the Hankelvalue to the true value, and its preceding value.We can conclude that a small convergence error is stronglyindicative that the Hankel approximation has reached a similarlyclose approximation of the true final value. Hence, we deem thealgorithm(s) to have converged if | h ∗ k ( r ) − h ∗ k − ( r ) | < ν for some tolerance ν . Figure 6 tells us that, in the case of anErdős-Rényi or scale graph of size n = 100, a tolerance of νν