Spectral dimension reduction of complex dynamical networks

Predicting the global dynamical states of complex networks is essential to anticipate critical breakdowns and to study the underlying resilience. Yet, we often rely solely on approximations of the average activity, which often misrepresent the particular structure of the network on which the dynamics takes place. We propose an approach to construct dynamical observables as weighted averages, based entirely on a linear combination of the dominant eigenvectors of the adjacency matrix. The observables form a low-dimensional system that can be solved to predict each global dynamical state of the network. Moreover, our approach is able to predict the multiple activation of modular networks as well as the tipping-points of random networks with arbitrary degree distributions, including scale-free networks.


I. INTRODUCTION
Critical breakdowns generally arise unexpectedly in dynamical systems [1].Particularly noteworthy examples are financial crisis [2,3], epileptic seizures [4], climate changes [5], and species extinction [6].While a lot of effort is done to anticipate the appearance of these breakdowns [7], there is still no simple and intuitive method for predicting breakdowns at the large-scale level.This is mostly due to the fact that dynamical systems are typically of high dimensions with complex architecture of interactions.
Counter-intuitively perhaps, it remains difficult to explain the origin of emergent mechanisms at the macro-scale of large systems, even though precise models have been developed at the microscopic levels [8].This poses a serious threat when searching for specific interventions aimed at preventing global system breakdowns.Thus, to design intervention strategies, it is essential to understand the global and local dynamical behaviors of complex systems exposed to structural stress, that is a mechanism that affects the network connectivity.A simple approach is to use dimension reduction methods to construct n-dimensional dynamical systems whose variables describe the global characteristics of the N -dimensional original systems under study [9], with n N .By having a thorough description of the reduced systems, we then gain a powerful tool for the characterization of large systems.
Constructing low-dimensional representations with good predictive power is, however, a challenging task.The main difficulty is that the dynamical properties strongly depend on the architecture of the interactions between the systems' components.For instance, food webs are known to have, by nature, a redundant and modular architecture of interactions, which is a less efficient structure but makes them more resilient to random fluctuations of food availability [10].Nonetheless, some approaches have been developed and have been proven to work well on many network structures.One well-known approach is network sparsification [11].This technique consists in removing nonessential edges so that certain core properties are approximately preserved [9,12].While this is useful to extract the backbone of a network and for the design of nearly linear-time algorithms for solving sys-tems of linear equations [13], it does not exactly provide a low-dimensional representation of the global states of a complex network that can be used to predict activity breakdowns.
Recently, Gao et al. have presented a dimension reduction formalism based on the degree-weighted average activity [14].It reduces any mutualistic system, i.e. one in which the activity of individuals benefits the neighbors, to a 1-dimensional equation that predicts the evolution of the degree-weighted average activity in terms of a single structural parameter, the degree-weighted mean connectivity.They have found that many dynamical systems exhibit the same universal bifurcations, independently of the local network structure.However, the formalism is based on a restrictive approximation: the network must not have degree correlations, i.e. it should belong to the particular family named uncorrelated random networks.Nevertheless, the results are surprisingly accurate to predict the tipping-points of real bipartite plant-pollinator networks.This spectacular outcome is not totally satisfactory since no fundamental reason is provided that would explain why the activity should be weighted by degrees rather than other structural parameters.
In an attempt to determine the tipping-points of 59 bipartite mutualistic plant-pollinator networks [15], Jiang et al. have proposed a 2-dimensional reduction.A reduced system is constructed of two variables that each describes the average activity of a distinct population: plants and pollinators.They have numerically explored three different averages: eigenvector-based, unweighted and degree-weighted.Although the eigenvector-based and degree-weighted averages are more accurately described by a 2-dimensional system than the uniform approximation, no theoretical insight is provided to guide the use of one average scheme over the others.Their results nonetheless suggest that, on certain types of networks, increasing the number of effective variables can lead to better predictions than the 1-dimensional formalism of Gao et al..
For arbitrary network structures, it remains unclear if the dimension reduction procedure of Gao et al. allows accurate predictions and why some dynamical systems require more effective dimensions.The goal of this paper is to address these issues and to propose a dimension reduction based solely on the spectral properties of the network structure.In do-ing so, we provide theoretical justifications for the required number of effective dimensions and the nature of the observations.Beyond mere improvement in precision over existing approaches, our method allows the detection of dynamical failures that would be missed altogether with previous reductions.
The paper is structured as follows.We first present the general framework for complex dynamical networks (Sec.II).We then describe a general method to obtain a 1-dimensional reduction for mutualistic networks (Sec.III).We show that the reduction scheme of Gao et al. emerges as an approximation of our general approach when specifically considering random networks.In Sec.IV we develop the cycle reduction, a multidimensional approach useful to reduce heterogeneous and bipartite networks.Next, in Sec.V, we complete the method by including the contribution of subdominant eigenvalues of the adjacency matrix.We finally assess the goodness of these reductions, as a function of the structure, and the nature of the dynamics (Sec.VI).

II. MODEL DEFINITION
The diverse nature of complex systems requires to establish a common ground.In Sec.II A, we regroup dynamical complex networks under a general model that encodes the structure and the dynamics.Then, in Sec.II B, we provide examples of contrasting models of dynamics satisfying the formalism used afterward in the paper to illustrate the dimension reduction methods.

A. General formalism
We consider a complex network of N nodes for which the interactions are encoded in the weighted and directed adjacency matrix W .The element w ij ∈ R of W is interpreted as the strength of the directed interaction from node j to node i.
Each node has an activity x i ∈ R whose evolution is governed by the general equation where F (x i ) and G(x i , x j ) are real-valued functions.For technical reasons that will become clear in the next section, both F (x i ) and G(x i , x j ) are required to have continuous derivatives of second order.The product w ij G(x i , x j ) specifies the type of interactions.If w ij ∂G(x i , x j )/∂x j ≤ 0, the interaction is competitive and the increase of activity of node j tends to decrease the activity of node i.If w ij ∂G(x i , x j )/∂x j ≥ 0, the interaction is mutualistic and therefore node j activity benefits node i.For mixed dynamics of mutualistic and competitive interactions, it is common to fix ∂G(x i , x j )/∂x j ≥ 0 and use negative weights w ij < 0 for competitive interactions.Unless specified, we will only consider mutualistic dynamics with w ij ∂G(x i , x j )/∂x j ≥ 0.
To describe the evolution of the whole system at the macroand the mesoscopic scales, it is convenient to focus on observables.We define an observable as a smooth function mapping the activities x 1 , . . ., x N to a real number.Among all observables, the linear observables, which are functions of the form L(x) = i a i x i ∈ R, are of particular interest.All the observables used throughout the paper belong to that category.A typical example would be the average activity with a i = 1/N .The concept of linear observable will be central to this study in order to assess complex systems through low-dimensional objects.

B. Examples of possible dynamics
A number of dynamical systems satisfy the form of Eq. ( 1).For instance, in computational neuroscience, the Cowan-Wilson model [16] describes the firing-rate activity of a population of neurons as where τ and µ are parameters controlling the steepness of the activation function and the firing-rate threshold, respectively [17].
In biology, the generalized Lotka-Volterra dynamics describes the evolution of the population of species in an ecosystem as where ω is the intrinsic growth rate [18], and x i is the population of individuals of species i.To prevent unbounded growth and account for species migration and the Allee effect, a more complex model of ecological networks has been proposed [14,19]: where all parameters are real-valued, B i accounts for the migration rate, K i > 0 for the ecosystem capacity, and C i > 0 for the minimum abundance for species growth.The parameters D i , E i , H i control the strengths of the interactions between the species.The Michaelis-Menten equation is yet another example [20].It applies to the gene regulatory networks and governs the concentration of substrates as where a, b, c ∈ R are parameters.In social networks, the spreading of a virus or rumors can be described using the Susceptible-Infected-Susceptible model (SIS) [21].In this context, the activity x i ∈ [0, 1] is interpreted as the probability of being infected and evolves according to with γ ≥ 0 as the normalized infection rate.

III. 1-DIMENSIONAL REDUCTION
The systems described by Eqs. ( 2) are N -dimensional as their dynamics are governed by N coupled differential equations.As the number of nodes N grows, the computational cost of solving N coupled equations increases which raises a number of issues [8,22].Moreover, the state of the original system given by the N -dimensional vector x becomes less comprehensive and less insightful, and does not provide insight into the general properties of the solutions.
Hence, we must rely on measures, or observables, to reduce N -dimensional systems to more practical and accessible objects.For instance, the unweighted average activity could be a simple measure on how dissimilar the system state is compared to a specifically chosen state.We also want to make predictions on those measures to anticipate dynamical breakdowns and locate the global state of the system on a standardized bifurcation diagram.However, when solely based on the unweighted average activity, those predictions are often nonrepresentative of the original system [15].
Alternatively, the weighted average activity seems more reliable as we inject additional information on the importance of nodes and has already been proven to be a promising avenue of breakdown predictions [14,15].In the next subsections, we introduce a general procedure to select a weighted average activity and to predict its evolution.

A. Derivation of the reduction formalism
Let us consider a real linear observable R of the activity constructed as a weighted average activity: where a i ∈ R, the i-th component of the column vector a, is a normalized weight so that Indeed, a i will be interpreted as the relative contribution, or centrality, of node i to the observable.When all a i are nonnegative, a is a probability vector.
Therefore, R is a function that takes the instantaneous activity of each node and returns a real number that describes the global state of the network.For instance, for a i = 1/N , R describes the unweighted average activity.Although the average activity is attractive because of its simplicity, it may not be easy to predict its value using only the structure of the network and the nature of the dynamics.Thus, we hypothesize that a should be specific to the structure.
Let us now explain how the weight vector a is constrained by the adjacency matrix W .By taking the time derivative of Eq. (3) and using Eq.(1) (refer to Appendix A for complete derivation), we obtain that the dynamics of R -truncated up to second-order terms O[(x k − R) 2 ] -is given by the 1-dimensional equation where β is a structural parameter given by and K is a N ×N diagonal matrix of diagonal elements w ij , the in-degree of node i.The parameter α can be measured directly on the network as the weighted averaged in-degree, Interestingly, we show, in Appendix A, that the closed form of Eq. ( 5) is satisfied only if a is a normalized eigenvector with eigenvalue α, of the transposed adjacency matrix W T , that is, We now have obtained a single equation [Eq.(5)] that governs the evolution of the weighted average activity of a complex network, and constrained the weight vector a to be adapted to the structure under study.
Clearly, Eq. ( 5) shares similarities with Eq. (1).We can interpret the former as a reduced system of one dynamical node that interacts with itself.The nature of its dynamics is identical to the one from the original system, i.e. specified by F (x i ) and G(x i , x j ), and the coupling is parametrized by α, β.The activity R of the single node describes the weighted average activity of the original network, and is obtained by solving Eq. ( 5) which solutions are solely controlled by the nature of the dynamics and the structural parameters α, β.

B. Choice of a universal weight vector
As briefly discussed, the weight vector a must be a normalized eigenvector of W T .In principle, any eigenvector of W T , except those that satisfy 1 T a = 0, could be used for the dynamical reduction.However, the larger the modulus of α is, the stronger is the influence of the structure on the weighted activity.An eigenvector whose eigenvalue has a low modulus leads to linear observable R that does not properly take into account the network structure, which in turn leads to correction terms O[(x k − R) 2 ] greater than those produced by eigenvectors with a higher modulus.This is the reason behind the choice of a as the eigenvector with the largest eigenvalue modulus: a is the dominant eigenvector.
For an arbitrary weighted adjacency matrix, the dominant eigenvalue and the components of the dominant eigenvector can be complex.In this case, the observable R as well as the structural parameters α and β are complex.The 1dimensional dynamical system of Eq. ( 5) becomes complex too and can be interpreted as a 2-dimensional real dynamical system.
There is however a large class of networks for which the dominant eigenvalue and the components of the dominant eigenvector are all real.For instance, undirected networks and directed networks that have non-negative weights w ij and are strongly connected fall into that class.This is guaranteed by the Perron-Frobenius theorem that states that if the network is strongly connected, i.e. a path exists between each pair of nodes, and that all edge weights satisfy w ij ≥ 0, then the dominant eigenvalue λ D of W T is non-negative λ D ≥ 0, and the dominant eigenvector is elementwise positive [23].Moreover, in practice, the dominant eigenvector can be easily computed using the power method.
The procedure to apply this 1-dimensional dimension reduction is simple.First, we compute the dominant eigenvalue α and the corresponding eigenvector v D of W T .Second, we define the normalized eigenvector a = v D /(1 T v D ), and obtain β ccording to Eq. ( 6).In most cases, we want to determine the weighted average activity at equilibrium R * , determined by solving This is a universal equation in the sense that α and β are independent of the dynamics, controlled by F (x i ) and G(x i , x j ); α and β only depend upon the network structure, encoded in W .The 1-dimensional reduction process has been applied to different dynamics for random uncorrelated networks and led to surprisingly accurate predictions (see Fig. 1).

C. Choice of an approximate weight vector
Recently, Gao et al. [14] have introduced a different 1dimensional reduction for dynamics of the form of Eq. (1).In this section, we show how their reduction is found as a special case of our 1-dimensional reduction when applied to uncorrelated random networks.
Uncorrelated random networks are a family of networks for which the degree distribution can be arbitrary but the probability of connection between two nodes is independent of the presence or absence of any other edge [25].We can easily generate random networks using the configuration model [26].We first sample the nodes in-and out-expected degrees κ in , κ out from an arbitrary degree distribution.Then, we con-  9) while dots are equilibrium states resulting from the evolution of the whole N -dimensional system.For each dynamics, 100 networks have been generated.For each network, the dynamics are integrated to equilibrium and an orange dot is placed at the corresponding point (α, R * ).Next, the network is perturbed by removing an edge and the dynamics is brought back to equilibrium, and a new dot at (α , R * ) is placed.The perturbation step is repeated 50 times for each network.nect node j to node i with probability where m = i κ out i is the expected total number of stubs.If the resulting network is strongly connected, the Perron-Frobenius theorem guarantees that the dominant eigenvector v D of W T will have only non-negative elements.We may then use this dominant eigenvector to construct the observable For networks that satisfy Eq. ( 10), the dominant eigenvector v D of W T can be approximated by the vector of out-degrees k out if [27] It can then be shown that the weights of the reduced system are approximately It results that α measures the average neighbor in-degree and one can show that β = 1.Therefore, as shown by the following equation, R is simply the average neighbor activity: It turns out that this special case is exactly the formalism proposed by Gao et al. [14] with R = x eff and α = β eff , in their notation.It also means that the formalism of Gao et al. is mostly appropriate for random networks [Eq.( 10)] respecting Eq. ( 12).Moreover, a recent work [28] has introduced a corrected eigenvalue approximation for random networks with power-law degree distribution p(k) ∼ k −γ with γ > 5/2.This may further limit the accuracy of Gao et al. approach with respect to our 1-dimensional reduction scheme.

IV. MULTIDIMENSIONAL REDUCTION: DOMINANT EIGENVECTORS
In the 1-dimensional reduction, it has been supposed that the dynamical global state of a network is dominated by the information contained in a single dominant eigenvalue and corresponding eigenvector.Therefore, it also presupposes that other eigenvalues can be safely neglected and do not provide any relevant information about the dynamics on the network.But, if the network admits many eigenvalues of equal and large modulus, it seems plausible to expect that all those eigenvalues are important as well, and should be included in a n-dimensional reduction.In this section, we address this problem by introducing a n-dimensional approach to predict the evolution of n coupled observables.We argue that if the spectrum then n observables should be considered, leading to a ndimensional reduced dynamical system.

A. Cycle reduction
Let us consider n observables R j , 1 ≤ j ≤ n, each being a different linear combination of the activity where a j is a real-value weight vector associated with the observable R j and normalized i [a j ] i = 1.As in the 1dimensional reduction, a j are yet undetermined.While there are several choices for a j that are a priori plausible, we discuss the cycle reduction, which is natural when several dominant eigenvalues are approximately of equal modulus.
Using a similar approach as the 1-dimensional reduction (See Appendix B), one finds that the evolution of the observables is given by where and α j is an observable of the weighted average neighbor in-degree, To satisfy Eqs. ( 16), the weight vectors are constrained by the structure and must transform according to which preserves the positiveness and the required normalization.Moreover, Eq. ( 19) needs to be a periodic application, i.e. a j+n = a j , in order to close the system to n observables.
The initial choice of a 1 is then highly constrained to satisfy this condition.In the following section, we will explain how the weight vectors can be computed using the dominant eigenvectors of W T .In contrast to Eq. ( 5) where a single observable is used, we now have developed a closed n-dimensional system of observables that are coupled by a set of structural parameters {α j , β j } j=1,..,n .

B. Choice of the universal weight vectors
It is yet unclear if one should use a n-dimensional reduction or a 1-dimensional reduction for a certain network structure.By answering this question, we also address how to set the weight vectors of the cycle reduction.
We first recall that the Perron-Frobenius theorem, guarantees that W T has a non-negative dominant eigenvector v D only if W is a connected.Therefore, we can rule out that it is always possible to construct a 1-dimensional reduction relying on the dominant eigenvector.
But the same reasoning also implies that we can always construct a n-dimensional system using the dominant eigenvector v D of W T .One could use a 1 = v D and apply iteratively Eq. ( 19) to obtain the set of weight vectors {a j } 1,...,n , as prescribed.In doing so, the resulting weight vectors would all be identical a 1 = a 2 = ... = a n = v D , as it obviously satisfies both Eq. ( 19) and the periodicity condition a j+n = a j .Hence, we find n identical observables R 1 = R 2 = ... = R n , and the constructed n-dimensional system is no better than the 1-dimensional system.For that reason, a n-dimensional cycle reduction is only advantageous if we can construct a set of distinct weight vectors a 1 = a 2 = ... = a n from the dominant eigenvectors.
The maximum number of significant and distinct observables that we can construct is determined by the periodicity of the transposed adjacency matrix W T .The periodicity n is the number of eigenvalues λ m of modulus equal to the spectral radius r ≥ 0, e.g.|λ m | = r.From the Perron-Frobenius theorem, they must be uniformly distributed on a circle, centered at the origin, in the complex plane.Thus, the m th dominant eigenvalue can be written as λ m = r e 2πim/n , for a given periodicity n.Since λ m is an eigenvalue of W T , it must have an eigenvector v m that satisfies By multiplying both sides by (W T ) (n−1) , we find Therefore, r n is a real-value positive n times degenerated eigenvalue of (W T ) n .Since each eigenvector of W T is also eigenvector of (W T ) n , we can combine those eigenvectors to construct new distinct eigenvectors of (W T ) n with eigenvalue r n and use them as weight vectors.We first construct the weight vector where c m ∈ C are arbitrary coefficients.From Eq. ( 19), we iteratively compute a j from a j−1 .By doing so, we both satisfy the periodic condition a j+n = a j and construct distinct weight vectors.
The reduction only requires to arbitrarily choose c = (c 1 , . . ., c n ) to construct a 1 .We propose to select c by minimizing the scalar product of the first two weight vectors where ).In Appendix C, we give a general and exact solution of c for n = 2.
In summary, the cycle reduction method goes as follows.First, compute a set of n eigenvectors {v m } m=1,...,n of W T whose eigenvalues have a modulus equal to the spectral radius r.Second, obtain {c i } i=1,...,n by solving Eq. ( 21).Third, iteratively construct a i from Eq. ( 19).Finally, compute α i , β i and solve R i at equilibrium from Eqs. (16).
An interesting aspect of this method is that it allows to combine the information of each observable to construct a global observable where π j ∈ R. Since R j = i [a i ] j x j , then the contribution p j of node j to the global observable is We can then tune π i to reach the desired node contributions.For instance, to access the unweighted average activity where A + is the Moore-Penrose pseudo-inverse of matrix A [23].

C. Examples: Star and bipartite networks
We give an example of the cycle reduction for a highly heterogeneous family of networks: star networks.We construct a star network of N nodes where the strength of a directed edge to a periphery nodes is s pc and s cp for edges directed toward the central node.Hence, the adjacency matrix is One finds that W T has two eigenvalues of modulus equal to the spectral radius, λ + = s pc s cp (N − 1) and λ − = − s pc s cp (N − 1).From the previous analysis, this signals that we can construct a 2-dimensional reduction.The associated eigenvectors are We construct the 2-dimensional reduction by combining the eigenvectors v + , v − to minimize |a T 1 a 2 | and to satisfy the normalization 1 T a 1 = 1.As it is easily verified, the following linear combinations fulfill the requirements: where the second vector has been obtained from Eq. (19).
Note that the overlap a T 1 a 2 = 0 exactly.Explicitly, in component form, we have One then finds that Hence, the 2-dimensional reduction reads One may note that R 1 is exactly equal to the activity of the central node and R 2 is the activity of a periphery node.Thus, the 2-dimensional formalism is an exact reduction in this example.This is confirmed in Fig. 2 where simulations and predictions are compared for the Gao et al., the 1-dimensional reduction and the 2-dimensional reduction for the Cowan-Wilson dynamics.Star networks are not the only systems conforming to the 2-dimensional reduction; all bipartite networks also do.Bipartite networks have nodes that can be separated into two groups such that connections only exist between nodes of different groups [26].This network architecture is common in many real systems such as plant-pollinator interactions [29], scientific collaborations [30], and actor-film networks [31].
Bipartite networks exhibit a remarkable and useful spectral property.Each eigenvalue is paired, i.e.
for all j = 1, 2, .., N , assuming that λ 1 ≥ λ 2 ≥ ... ≥ λ N [32].Thus, a bipartite graph contains two eigenvalues λ 1 , λ N of modulus equal to the spectral radius.This suggests that a 2dimensional representation could always be constructed out of these two eigenvalues, following the prescription of Sec.IV A.
We therefore gain a clear understanding of the results of Jiang et al. [15], suggesting that 2-dimensional reductions are better predictors of tipping-points for bipartite networks, compared to 1-dimensional approaches.Moreover, it solves the problem of selecting the right weighted average: The eigenvector-weighted average should always be favored over others.

V. MULTIDIMENSIONAL REDUCTION: INCLUDING SUBDOMINANT EIGENVECTORS
Until now, we have developed a truly direct method to construct n-dimensional reduced systems.Using only the network structure, we can first identify the number of dimensions, i.e. the number of eigenvalues of modulus equal to the spectral radius, and then construct the weight vectors a j to predict the observables R j .
Yet, our method is compelling only if the observables R j are good indicators of the global states of the network, which requires that each region of the network contributes significantly to at least one observable R j .Since we use the dominant eigenvectors, we do not control the contribution of each node.Thus, if the dominant eigenvector assigns negligible weights to some nodes, it may result in an incomplete description of the network.Unfortunately, modular networks fall into this category.Let us introduce the stochastic block model (SBM) to understand the underlying problem of misrepresentation.
The SBM is a generative model of modular networks [33].Nodes are first assigned to modules.Then, we connect a node from module s to a node from module r with probability p rs .This simple method allows us to generate modular random networks.
The spectrum of a SBM network is rather different than the spectrum of a random network (Fig. 4).We first note that we only have a single eigenvalue of modulus equals to the spectral radius, indicating that we should use a 1-dimensional reduction.However, the eigenvalues are distributed in a multimodal distribution with as many dominant and subdominant eigenvalues as they are modules.
For instance, let us consider a network with two communities of equal size N/2.Using the spectral theory of random matrices, we estimate the two dominant eigenvalues and their corresponding eigenvectors .
If p rs /p rr ≈ 0 ∀r = s, the dominant eigenvector v 1 assigns a negligible weight to the nodes in the community r = 1 [34].Thus, if we solely use a 1-dimensional reduction with a = v 1 , the observable R will not take into account the activity of half the network.Fortunately, the second-dominant eigenvector accounts for the remaining nodes.Thus, we must apply twice the 1-dimensional reduction of Sec.III A and construct two uncoupled observables, one with the dominant eigenvector R D = v T 1 x, and one for the subdominant eigenvector R SD = v T 2 x, for which the dynamics follow, To make a global prediction, we simply combine the observables It follows that the global structural parameter is also a linear composition In general, we can construct as many uncoupled observables as the number of modules in the network, by using the eigenvectors associated with the eigenvalues detached from the bulk of the spectrum.We show a numerical example of this method for a network of two communities (Fig 3).The Gao et al. formalism predicts a single bifurcation, which almost coincides with the transition of the densest community.In Fig 3(b), the 1dimensional reduction, using the first dominant eigenvector, predicts the activity accurately.However, as it has been previously discussed, the first dominant eigenvector omits half the network and we only see a single bifurcation.Thus, even if we are highly accurate, we miss characterizing the distinctive multistep bifurcation, a prominent feature of interacting networks [36].Finally, combining observables as in Eq. ( 31), we recover the bifurcations of both modules [Fig 3(d)].

VI. GOODNESS OF REDUCTION
The goodness of the reduction method, i.e. how accurate are the predictions of the low-dimensional representation compared with the observations on the original network, depends on the nature of the dynamics and on the network structure.For instance, some complex patterns of interactions may be less amendable to a low-dimensional formalism, resulting in disparities between the predicted and exact values of the observables.In this section, we explore the impacts of the structure and the dynamics on the goodness of the reduction methods.

A. Impact of the structure
To measure the impact of the structure, we introduce a generative model of networks called generalized preferential attachment model [37].Parameters of this model can be continuously tuned to obtain networks ranging from chain-like networks to star networks, with scale-free networks as an intermediate state.Scale-free systems are an important family of networks, recognizable by their power-law degree distribution p(k) ∼ k −γ [24].Due to their lack of well-defined characteristic scale, it is a priori unclear whether these systems can be efficiently reduced.
The growth process of the generalized preferential attachment goes as follows.We initialize the network with two connected nodes.Then, at each time step t, we add a new node to the network.It is connected to an existing node chosen with probability where ν ∈ R is the exponent of the attachment kernel, s j (t) is the number of connections of node j at time t, and N (t) is the number of nodes at time t.The generative model is solely tuned by the kernel parameter ν ∈ R. It controls the inequalities of the attachment probability.On the one hand, if ν 1, the generated networks are star-like as we always attach new nodes to the richest node.On the other hand, if ν 0, the networks are more chainlike as we always connect to the least connected node [37].The classic preferential attachment model is found for ν = 1.Therefore, for 0 < ν < 1, we observe a continuum of network organizations which gradually become more scale-free as the parameter ν is increased.Examples of networks generated from this model are illustrated in Fig. 5.
From now on, we will distinguish the predicted observable from the reduced system, denoted R(α), and the measured observable R(α) = a T x from the original network.
We have applied the degree-weighted, the 1-dimensional and the 2-cycle reductions to networks generated with the generalized preferential attachment model for ν ∈ [−1, 2].For each network, we have computed the total error ∆ R between the measured activity R * (α) = a T x * at equilibrium on the original network and the predicted activity R * (α) by the reduction system, We have found a transition in the dimension reduction accuracy for all methods at ν = 1, corresponding to the preferential attachment model [ Fig 6].As we enter the star-like region ν > 1, the average error is reaching a plateau to specific values for the 1-dimensional reductions while the 2-dimension reduction becomes highly effective.
We argue that this transition is not dynamics-specific but mostly due to the network architecture.The nature of the error can be interpreted by a careful examination of the generative model.First, negative values of ν tend to homogenize the degree of the nodes.Thus, the more uniform the network is, the easier it is to capture its behavior in a 1-dimension reduction.For positive values of ν, the reduced model tends to favor degree inequalities which are best achieved when the networks are star-like.For 0 < ν < 1, the degree distribution resembles a power law with exponential cutoffs [37].However, a pure power-law distribution is only achieved at precisely ν = 1.Thus, this transition in the degree distribution forces the reduction of the degree-weighted reduction to predict inaccurate observables.Finally, the region ν > 1 is dominated by star-like networks, which has been previously shown to be better represented by the 2-dimensional reduction than any 1-dimensional reduction (See Subsec.IV C).

B. Detection of transitions
Most dynamical systems exhibit activity bifurcations when a certain structural threshold is reached [38].Therefore, the goodness of the reduction should display at least qualitative changes as the structural threshold is crossed.We investigate this kind of prediction using the SIS model.
The SIS model has been thoroughly studied in the last decade [39].In the SIS model, nodes reversibly switch from susceptible to infected states with a certain probability that depends on their neighborhood.We can formulate this dynamics using a mean-field approach, where x i is the probability that node i is infected and γ ≥ 0 is the normalized infection rate.In this model, the average fraction of infected node x = N −1 i x i undergoes a critical transition at a certain threshold γ C .The classical problem in the study of the SIS model consists in estimating the value of γ C above which a significant fraction of the whole system is infected [40].
We will however study a related problem: the parameter γ is fixed and the structure is evolving.Using the dimension reduction procedure, we investigate the critical structural parameter α global , or an equivalent parameter depending on the reduction approach, for which the global state at equilibrium R * undergoes a critical transition characterized by In Fig. 7, we investigate the errors on the position of the critical transition for the degree weighted 1-dimension approach, the eigenvector-weighted 1-dimensional approach, and the 2dimensional cycle reduction.We use a network of N = 60 nodes generated from the generalized preferential attachment model with ν = 1.8.We observe that the two proposed approaches based on dominant eigenvectors are able to accurately predict the critical transition while the degree-weighted approach does not.This behavior is typical for reductions of networks generated from ν ∈ [1, ∞).Our results indicate that even if the 1-dimensional observable fails to predict the true level of activity R, it still predicts with high accuracy the onset of the epidemy.We conclude that the largest eigenvalue is a reliable indicator of the onset for correlated networks.Perhaps, this conclusion is not surprising as it has been previously discovered under a different approach [41].Nonetheless, it supports the proposed reductions as valuable candidates for predicting the onset of critical transitions.

C. Impact of the dynamics
As previously discussed, the goodness of the reduction is highly dependent on the network structure.The other element that impacts the goodness of the reduction is the nature of the dynamics.For instance, linear dynamics such as F (x i ) = x i , G(x i , x j ) = x j lead to exact 1-dimensional reduction.However, typical dynamics are nonlinear and may add significant contributions to the quadratic terms that have been neglected (See Appendix A).We investigate this aspect by looking at two contrasting dynamics: Cowan-Wilson model and Lotka-Volterra model.
First, let us introduce the error ∆ α on the structural parameter that predicts the original network activity.We compare the measured structural parameter α = a T k in on the original network with the structural parameter α that matches, from the 1-dimensional reduction, the measured activity R * = a T x * on the original network.The relative error can be written as Notice that this is different from Eq. ( 34): ∆ α compares the horizontal error in the space (α, R) while ∆ R measures the vertical error.Both errors convey distinct information and are complementary.However, we will use ∆ α since it can be written as a function that depends explicitly on the nature of the dynamics for the 1-dimensional reduction.From Eq. ( 5), we obtain α at the dynamical equilibrium Ṙ = 0.This leads to One can then easily evaluate the error for specific dynamics as in [42].In the following paragraphs, we give two examples of dynamics and compare the errors for the Gao et al. formalism and our 1-dimensional reduction.

Error on the Cowan-Wilson dynamics
The Cowan-Wilson dynamics describes the firing-rate activity of populations of neurons.The evolution of a node activity is given by Eq. (2a), which is transcribed here as where τ > 0. The equilibrium solution x * cannot be found analytically, so it must be evaluated numerically, even for N = 1.However, x * is easily approximated for extreme values of activities.We derive error estimations for extreme regimes x * j µ and x * j µ.In general, from Eq. ( 37), the error can be written as In the limit of high levels of activity x * j µ, the exponential vanishes and one finds that x * ≈ W 1 = k in .Therefore, the error is approximately For the 1-dimensional reduction, a is an eigenvector of W T such that a T W = αa T .Since a is normalized, a T 1 = 1, it follows that a T W 1 = α.Thus, the error ∆ α vanishes in the limit of large activity.The same applies for Gao et al. formalism for high activity.By using a = k in , one finds that α = (k in ) T k out and R * = (k out ) T k in so that ∆ α → 0.
Using a similar procedure for x j µ, one can also show that ∆ α ≈ 0 for both methods.
We conclude that in the extreme regimes of high and low activities, both reduction methods provide a practically exact solution.However, we are more often interested in the hysteresis region where the activity collapses rapidly.Unfortunately, analytic error estimates are lacking in this regime.Still, numerical results and theoretical insights suggest that the proposed 1-dimensional reduction should always be favored over the degree-weighted reduction.This is confirmed below for a more tracktable dynamics.

Error on the Lotka-Volterra dynamics
Let us consider the Lotka-Volterra dynamics governing the evolution of species populations.The N -dimensional system goes as where • denotes an elementwise multiplication.At equilibrium, x * satisfies With Eq. ( 37), we write the expected error as Using the 1-dimensional reduction, Thus, the error depends only on β: Therefore, the difference between the exact value for a T x * and the approximate value derived from the 1-dimensional reduced system is only 1 − 1/β.Given the expression for β [Eq.(A5)], it should be close to β ≈ 1, so the error goes to zero ∆ α ≈ 0. This contrasts with results derived by Tu et al. [42] for the method of Gao et al., where they reported non-vanishing error averages and variances.

VII. CONCLUSION
We have built systematic methods of dimension reduction adapted to different families of networks (random, star-like, bipartite, SBM).The activity of the reduced systems is used as an indicator of the global activity of large networks.Without further restriction than imposing a linear form of the global activity, we have found that the dominant eigenvectors of the adjacency matrix are central to the global states' evolution.Moreover, when considering the cycle reduction, the dimension of the reduced systems is found to correspond to the periodicity of the adjacency matrix.
We have further shown that the proposed reduction of Gao et al. is a special case of the general scheme when applied to uncorrelated random networks.Moreover, the range of applicability of our method extends to modular, heterogeneous and bipartite networks.
Our results suggest, both numerically and theoretically, that the eigenvector-weighted reduction should be preferred over the degree-weighted reduction.The fundamental reason for this should be interpreted as a paradigm shift.Originally, the degree-weighted reduction has been used to approximate the state of a network by the state of the average neighbor.But, the degree of a node is only a local centrality measure since it does not provide any information to whom a node is connected with.In the eigenvector-weighted reduction, the dominant eigenvector provides a more global node centrality as it contains the information on how each node is connected with the rest of the network [26].Therefore, the eigenvector-based reduction brings a new light on the influence of each node on the global states of a network.
On a more practical side, our general method is able to predict the correct number of bifurcation points.This number depends on the dimension of the reduction.As shown in the paper, it is not guaranteed that a 1-dimensional system will capture all the bifurcations that can occur in the original system.In the SBM case, the 2-dimensional reduction reveals additional bifurcation points that are missed altogether by all 1-dimensional reductions.
As a closing remark, observe that our reduction method has been designed to access large dynamical networks through low-dimensional formalisms.More than just a theoretical work, this approach can be used to address concrete problems of real-world systems.For example, it might be used to describe with high accuracy the bifurcation patterns, to identify dynamical vulnerabilities, to suggest intervention strategies to prevent dynamical breakdown, or to classify networks on a standardized diagram.
terms of R only.To do so, we develop each function around the observable: It means that F (x i ) does not impose any constraint on a. Now, for the function G(x i , x j ), we develop around x i = βR and x j = γR: G(x i , x j ) ≈ G(βR, γR) + (x i − βR)G 1 (βR, γR) + (x j − γR)G 2 (βR, γR) where second order terms have been neglected.Letting α = ij a i w ij , we find that i,j a i w ij G(x i , x j ) is given by i,j a i w ij G(x i , x j ) ≈ αG(βR, γR) a i w ij (x j − γR).
The left-hand side is a function of R only if the linear terms cancel out exactly, which is possible if and only if Since R = x T a, we conclude that the last two equations are satisfied for all x ∈ R N only if a is an eigenvector of both matrices K and W T , with corresponding eigenvalues βα and γα.Although we cannot solve these two equations simultaneously in general, we can enforce that at least one equation is satisfied exactly.Choosing a as an eigenvector of W T , we can prove that γ = 1 if the vector a is normalized 1 T a = 1: If W T a = λa and λ = αγ, then λ = 1 T W T a = a T W 1 = a T k in = α, so γ = 1.
We then choose β to best satisfy Eq. (A2a) by minimizing the mean square error (MSE): where the symbol || • || denotes the standard euclidean norm.Basic calculus leads to Note that β * is a ratio of weighted averages where b is normalized 1 T b = 1 and has for elements b i = a 2 i / N i=1 a 2 i .From the construction of b, we deduce that b must be similar to a and β * close to 1, which has been confirmed throughout most of the simulations.

Appendix B: Derivation of the multidimensional cycle formalism
For the cycle reduction, we construct n observables with normalized weights 1 T a j = 1.Using Eq. ( 1), we find that the dynamics of R k is equal to As for the 1-dimensional reduction, one finds that up to the second order of corrections.We then develop G(x i , x j ) around x i = β k R k and x j = γ k R k+1 , which yields Using the same arguments as in the 1-dimensional reduction, one can prove that i,j [a k ] i w ij G(x i , x j ) ≈ α k G(β k R k , γ k R k+1 ), with α k = a T k in , only if the following equations are satisfied simultaneously The second equation is easily satisfied if: with γ k = 1.After n applications of Eq. (B6), we close the system with a n+1 = a 1 , which is the respected if a 1 is eigenvector of (W T ) n .As for the parameter β k , we minimize the MSE and find )

FIG. 1 .
FIG. 1. (Color online) Observable R * = a T x * at equilibrium as a function of the dominant eigenvalue α of W T , for different dynamics on Erdős-Rényi networks of N = 100 nodes and density p = 0.1.(a) Cowan-Wilson dynamics Eq. (2a) with τ = 1, µ = 3, (b) SIS dynamics Eq. (2e) with γ = 1, (c) Mutualistic ecological dynamics Eq. (2c) with Bi = 0.1, Ci = 1, Ki = 5, Di = 6, Ei = 0.9, Hi = 0.1 [24], (d) Michaelis-Menten dynamics Eq. (2d) with a = 1, b = 1, c = 1.Dashed lines are theoretical predictions obtained from Eq. (9) while dots are equilibrium states resulting from the evolution of the whole N -dimensional system.For each dynamics, 100 networks have been generated.For each network, the dynamics are integrated to equilibrium and an orange dot is placed at the corresponding point (α, R * ).Next, the network is perturbed by removing an edge and the dynamics is brought back to equilibrium, and a new dot at (α , R * ) is placed.The perturbation step is repeated 50 times for each network.

FIG. 2 .
FIG. 2. (Color online) (a) Schematisation of the star network of N = 6 nodes where the edge weight toward the core is twice the weight of an edge toward the periphery.(b) Average neighbor activity at equilibrium as a function of its structural parameter using the degree-weighted reduction of Gao et al.[14].(c) Dominant eigenvector weighted average activity at equilibrium R * = v T D x * as a function of the dominant eigenvalue α for the 1-dimensional reduction [Eq.(5)].(d) Average activity at equilibrium obtained by a combination of two observables x * = N −1 [R * 1 + (N − 1)R * 2 ], as a function of the average output degree k out = N −1 [α1+(N −1)α2] computed using the 2-dimensional reduction formalism [Eq.(28b)] .Full lines are results from simulations and dashed lines are theoretical predictions.The network dynamics is the Cowan-Wilson model with τ = 1, µ = 3 [Eq.(2a)].

FIG. 3 .
FIG. 3. (Color online) (a) Schematisation of the undirected planted partition network of two communities of 100 nodes each with indensities pin = 0.4 and pin = 0.7, and out-density pout = 3 × 10 −3 .(b) Average neighbor activity at equilibrium as a function of its structural parameter using the degree-weighted reduction of Gao et al. [14].(c) Dominant eigenvector average activity R = v T D x at equilibrium as a function of the associated eigenvalue α for the 1dimensional reduction [Eq.(9)].(d) Combination of the uncoupled observables Rglobal = (RD + RSD)/2 at equilibrium as a function of the structural parameter αglobal [See Eq. (30b)].Full lines result from simulations and dashed lines are theoretical predictions.The dynamics is the Cowan-Wilson model with τ = 1, µ = 3 [Eq.2a].

FIG. 4 .
FIG. 4. (Color online) Spectral density for (a) random networks of 80 nodes with density p = 0.1, (b) of SBM with two communities of 40 nodes with in-densities p11 = 0.3 and p22 = 0.7 and out-density p12 = p21 = 0.01, and (c) of SBM with four communities of 40 nodes with in-densities p11 = 0.4, p22 = 0.6, p33 = 0.7, and p44 = 0.9 and out-density prs = 0.05 ∀r = s.The orange lines highlight the position of the dominant and subdominant eigenvalues in the expected network.Each density is produced by collecting spectra from 500 network instances.

4 FIG. 5 .
FIG. 5. (Color online) Network instances produced using the generalized preferential attachment model.The model tends to generate chainlike networks for negative values of ν and star networks for positive values.Figure inspired from [35].

FIG. 6 .
FIG. 6. (Color online) (Top) Total error on the prediction of the global state activity at equilibrium for networks obtained from the generalized preferential attachment model [See Eq. (34)].A gray line indicates the classical preferential attachment model at ν = 1 where a critical transition in the reduced descriptions is found.The total error is averaged from 300 networks of N = 200 generated uniformly on the domain ν ∈ [−1, 2].The activity on the network is the SIS model with γ = 1 [See Eq. (2e)].(Bottom) Instances of bifurcation diagrams for the three reduction schemes (columns) and for ν = {−1, 1, 2} (rows).For the 2-dimensional cycle reduction, the x-axis is meant to be the average of the structural parameters: α global = (α1 + α2)/2.The blue dashed lines are predictions from reduced systems and the orange lines are the measured activities on the original networks.The gray regions indicate the absolute errors [See Eq. (34)].

FIG. 7 .
FIG. 7. (Color online) (a) Schematisation of the undirected network generated using the generalized preferential attachment model with ν = 1.8 and N = 60.(b) Average neighbor activity at equilibrium as a function of its structural parameter using the degree-weighted reduction [14].(c) Dominant eigenvector average activity at equilibrium R * = v T D x * as a function of the associated eigenvalue α for the 1-dimensional reduction [Eq.(5)].(d) Solution at equilibrium of the 2-dimensional system R * global = (R * 1 + R * 2 )/2 as a function of the structural parameter αglobal = (α1 + α2)/2.Full lines are results from simulations and dashed lines are theoretical predictions.Dotted lines indicate the position of the transition αglobal.The dynamics is the SIS model with γ = 0.2 [Eq.(2e)].