Generating directed networks with prescribed Laplacian spectra

Complex real-world phenomena are often modeled as dynamical systems on networks. In many cases of interest, the spectrum of the underlying graph Laplacian sets the system stability and ultimately shapes the matter or information flow. This motivates devising suitable strategies, with rigorous mathematical foundation, to generate Laplacian that possess prescribed spectra. In this paper, we show that a weighted Laplacians can be constructed so as to exactly realize a desired complex spectrum. The method configures as a non trivial generalization of existing recipes which assume the spectra to be real. Applications of the proposed technique to (i) a network of Stuart-Landau oscillators and (ii) to the Kuramoto model are discussed. Synchronization can be enforced by assuming a properly engineered, signed and weighted, adjacency matrix to rule the pattern of pairing interactions.


I. INTRODUCTION
Complex networks play a role of paramount importance for a wide range of problems, of cross-disciplinary breadth. In several cases of interest, networks define the skeleton of pairwise interaction between coupled populations, families of homogeneous constituents anchored to a given node of the collections. The nature of existing paired relationships between mutually entangled populations is encoded in the weight of the links, that bridge adjacent nodes. Cooperative and competitive interactions can, in principle, be accommodated for by allowing for the weights to take positive or negative values. Local and remote non linear couplings shape the ensuing dynamics, and possibly steer the system towards a stationary stable equilibrium which is compatible with the assigned initial condition. The stability of the fixed points can be analyzed by studying the dynamical system in its linearized version. For reaction-diffusion systems defined on networks, the stability of the inspected equilibrium is ultimately dictated by the spectrum of the discrete Laplacian matrix [1][2][3].
The eigenvalues of the Laplacian define, in fact, the support of the dispersion relation, the curve that sets the rate for the exponential growth of the imposed perturbation. More specifically, external disturbances can be, in general, decomposed on the basis formed by the eigenvectors of the Laplacian operator. Each eigenvector defines an independent mode, which senses the web of intricate paths made accessible across the network: the perturbation can eventually develop, or, alternatively, fade away, along the selected direction, depending on the corresponding entry of the dispersion relation, as fixed by its associated eigenvalues.
Stability is an attribute of paramount importance as it relates to resilience, the ability of a given system to oppose to external perturbations that would take it away from the existing equilibrium. Similarly, synchronization, a widespread phenomenon in distributed systems, can be enforced by properly adjusting the spectrum of the matrix which encodes for intertwined pairings. Based on the above, it is therefore essential to devise suitably tailored recipes for generating networks, which display a prescribed Laplacian spectrum, compatible with the stability constraint [13,14].
The problem of recovering a network from a set of assigned eigenvalues has been tackled in the literature both from an algorithmic [4][5][6] and formal [7,8] standpoints. In [8], a procedure is discussed to generate an undirected and weighted graph from its spectrum.
The result extends beyond the well-known theorem of Botti and Merris [9] which states that the reconstruction of non-weighted graphs is, in general, impossible since almost all (non-weighted) trees share their spectrum with another non-isomorphic tree. In [11], a method is proposed to obtain a, directed or undirected, graph whose eigenvalues are constrained to match specific bounds, which ultimately reflect the nodes degrees, as well as the associated weights. In [12], a mathematically rigorous strategy is instead developed to yield weighted graphs which exactly realize any desired spectrum. As discussed in [12], the method translates into an efficient approach to control the dynamics of various archetypal physical systems via suitably designed Laplacian spectra. The results are however limited to undirected Laplacians, characterized by a real spectrum. The purpose of this paper is to expand beyond these lines, by proposing and testing a procedure which enables one to recover a signed Laplacian operator which displays a prescribed complex spectrum. Signed Laplacians are often used in the literature for applications which relate to social contagion, cluster synchronization or repulsive-attractive interactions [19,20]. In engineering, they are often employed in modeling microgrids dynamics [21].
The paper is organized as follows. The first section is devoted to illustrating the devised method, focusing on the mathematical aspects. We then turn to discussing the implementation of the scheme and introducing the sparsification algorithms that are run to cut unessential links. In the subsequent section, we elaborate on the conditions that are to be met to generate a positively weigthed network. This discussion is carried out with reference to a specific setting. Then, we apply the newly introduced technique to the study of an ensemble made of coupled Stuart-Landau oscillators [22][23][24] and to (a simplified version of) the Kuramoto model [27,28]. Finally, in the last section, we sum up the contributions and provide concluding remarks.

EIGENVALUES.
Consider a network made of Ω nodes and denote by A the (weighted) adjacency matrix, where structural information is encoded. More precisely, the element A ij is different from zero when a directed link exists from j to i. The entries of the matrix A are real numbers and their signs reflect the specificity of the interaction at play: negative signs stand for inhibitory (or antagonistic) couplings, while positive entries point to excitatory (or cooperative) interaction. From the adjacency matrix, one can define its associated Laplacian operator. This is the matrix L, whose elements are represents the natural extension of the concept of (incoming) connectivity to the case of a weighted network and δ ij denotes the Kronecker delta.
We shall here discuss a procedure to generate a Laplacian matrix, which displays a prescribed set of eigenvalues. As anticipated above, we will focus in particular on directed Laplacians, which yields, in general, complex spectra. Concretely, we begin by introducing a collection of Ω = 2N + 1 complex quantities defined as 1 (1) The first 2N elements come in complex conjugate pairs and we set in particular Λ i = Λ * i+N , ∀ i = 1, ..., N , where (·) * stands for complex conjugate. The aim of this section is to develop a rigorous procedure to construct a directed (and weighted) graph G with 2N + 1 vertices, whose associated Laplacian has {Λ i } for eigenvalues. Recall that {Λ i } contains the null element, since this latter is, by definition, a Laplacian's eigenvalue.
The procedure that we are going to detail in what follows exploits the eigenvalue decomposition of the Laplacian matrix. To this end we will seek to introduce a proper eigenvector 1 We shall assume all the eigenvalues but 0 to be complex numbers. Let us remark that the method developed readily adapts to the case where also real eigenvalues are present. In this case, the eigenvectors associated to real eigenvalues are generated according to the prescriptions of [12]. Eigenvectors linked to complex eigenvalues are instead assigned following the procedure outlined here. basis such that is a Laplacian. In Eq. (2), D is a diagonal matrix where the sought Laplacian eigenvalues are stored. More specifically, D ii = Λ i−1 for i = 2, .., 2N + 1 and D 11 = Λ 2N +1 = 0. The problem is hence traced back to constructing V , whose columns are the right eigenvectors of L. We also recall that rows of the inverse matrix V −1 are the left eigenvectors of L. As outlined in [15], Laplacian (right and left) eigenvectors should satisfy a set of conditions: 1. The columns of V , which refer to complex conjugate eigenvalues, must be complex conjugate too.
2. The same condition holds for the rows of V −1 .
3. Moreover, the columns of V (resp. the rows of V −1 ) corresponding to eigenvalues different from 0 should sum up to zero.
4. Finally, the right eigenvector relative to the null eigenvalue should be uniform (i.e. display identical components).
In light of the above, we put forward for V the following structure: where, i stands for the imaginary unit, U is an invertible N × N matrix having real entries, the vector u = (1 . . . 1) T has dimension N × 1 and the vector v is defined as The first column of V is hence a uniform vector, corresponding to the eigenvector associated to the null eigenvalue. By construction, every other column sums up to zero, that is (4) holds. In the following, we will write D jj = α j + iβ j , which, in turn, implies D j+N,j+N = α j −iβ j , for j = 2, . . . , N +1. Here, α j and β j are real quantities and respectively denote the real and imaginary parts of the j-th eigenvalue. To proceed further, one needs to determine the inverse of V .
To achieve this goal we begin by considering a generic matrix W , which satisfies the general constraints that are in place for V −1 . In formulae: Note that the jth and (N + j)th rows of W , for j = 1, .., N + 1, are complex conjugated, as required. Moreover, summing all the elements of each row (but the first) yields zero, a condition that the inverse of V should meet, as anticipated above. Building on these premises, we shall here determine the unknown S, w and d so as to match the identity W V = I, where I stands for the (2N + 1) × (2N + 1) identity matrix. This implies, in turn, that W ≡ V −1 due to the uniqueness of the inverse matrix.
A straightforward manipulation yields the following conditions for, respectively, d and where use has been made of the identity u T u = N and where The quantity d is completely specified by the first of Eqs. (7) and solely depends on c and N , the size of the system. By making use of the identities (4) and (6), one can progress in the analysis of the second and third conditions (7) to eventually get: where: It is, therefore, immediate to conclude that The analysis can be pushed further to relate S to matrices U and E. The calculation, detailed in Appendix A, yields The expression for S * can be immediately obtained by taking the complex conjugate of the above equation. In conclusion, the matrix W defined in (5) is the inverse matrix of V , provided that d and S are respectively assigned as specified above.
Clearly, matrix L defined in (2) has the desired spectrum (1). We should, however, prove that L is a Laplacian. This amounts to showing that L is a real matrix, whose columns sum up to zero. The proof is given hereafter.
Proposition. Matrix L is real.
From Eqs. (3) and (5), one can readily compute the elements of L via matrix products and, taking advantage of the block structures of V and W , prove that L is real. ℜ(·) is introduced to represent the real part of (·). The generic element L st can be written as: due to the diagonal structure of the matrix D and recalling that D 11 = 0. By making use of the specific form of V and W , one gets: since αβ + α * β * = 2ℜ(αβ), for any complex numbers α and β. One can thus conclude that L, as generated by the above procedure, is real.
Proposition. Each column of L sums up to zero.
Because of the diagonal structure of D: Then: since (i) D 11 = 0 and (ii) the components of all eigenvectors corresponding to non-null eigenvalues, sum up to zero. Notice that this result can also be proven by observing that the uniform vector d1 is the left eigenvector corresponding to the null eigenvalue, that is From (23), it follows that i L ij = 0, for every j.
Proposition. L is balanced.
We can also show that j L ij = 0 i.e. that the sum of all the elements of any given row i returns zero. According to (2), the first column of V is the right eigenvector corresponding to eigenvalue 0, namely This implies, in turn, j L ij = 0 ∀i, which ends the proof. The Laplacian is hence balanced, as the sums on the rows and on the (corresponding) columns return the same result.
From the generated Laplacian operator, one can readily calculate the adjacency matrix of the underlying network. In general, for any assigned spectrum, the recovered adjacency matrix is fully connected, meaning that there exists links connecting each pair of nodes.
Notice that links are weighted and signed. The weights can be small or have a modest impact on the spectrum of the associated Laplacian. This motivates the implementation of a dedicated sparsification procedure, which seeks to remove unessential links, in terms of their reflection on the ensuing Laplacian spectrum. The next section is devoted to elaborating along these lines.

III. EXAMPLES AND SPARSIFICATION.
In this section we discuss a sparsification procedure, which aims at a posteriori simplifying the structure of the recovered network. To this end, we begin by generating a network following the strategy outlined in the preceding section, and which yields an assigned spectrum for the associated Laplacian. The Laplacian spectrum that we seek to recover consists of Ω = 2N + 1 complex entries, the eigenvalues, which are here confined in the left portion of the complex plane, by setting ℜ(Λ j ) = α j < 0 for j ≥ 2, see blue crosses in Fig. 1

(a).
This choice is somehow arbitrary, and ultimately amounts to enforce stability into a linear system of the form: where x i is the i-th entry of the Ω-dimensional state vector x. In the final part of the paper, we will turn to considering more complex scenarios where the stability of the examined dynamics is also influenced by local reaction terms.
The network that we obtain, following the scheme outlined in the preceding sections yielding a Laplacian with the prescribed spectrum, is in general fully connected. In other words, a weighted link exists between any pair of nodes. The weights of the link can be, in principle, very small and, as such, bear a modest imprint on the ensuing Laplacian spectrum. Motivated by this observation, we perform an a posteriori sparsification of the obtained network: this aims to identifying and then removing the finite subset of links that appear to have a modest impact on the eigenvalues of the associated Laplacian.
The first sparsification procedure that we have considered, aims at removing unessential links while confining the spectrum of the Laplacian operator within a bounded region of the complex plane. More precisely, we focus on the links which display a weight in the range (-σ, σ), where σ is a small, arbitrarily chosen, cut off. All links whose weight is smaller that σ in absolute value are selected, in a random order. The selected link is removed and the modified Laplacian spectrum computed. Denote byΛ j , for j = 2, ..., 2N + 1, the Laplacian eigenvalues obtained upon removal of the link. The change to the network arrangement becomes for j = 2, ..., N . Here ℑ(·) stands for the imaginary part of (·) and δ is an arbitrary threshold which quantifies the amount of perturbation that is deemed acceptable for the problem at hand. As a further condition, we check that ℜ(Λ j ) < 0 for j = 2, ..., 2N + 1, which, in turn, corresponds to preserving the stability of the linear system (25). Clearly, the order of extraction of the links, which are candidate to be trimmed, matters. Different realizations of the procedure of progressive sparsification illustrated above might hence result in distinct final outcomes. In Fig. 1(a), the eigenvalues obtained after the sparsification algorithm are plotted (red circles) for one choice of the cutoff δ. The sparsity pattern of the adjacency matrix obtained at the end of the above procedure is displayed in panel (b) of Fig. 1.
In panels (c) and (d) of Fig. 1 we plot, with an appropriate color code, the entries of the adjacency matrices, before and after the sparsification. Only weights which are significantly different from zero (see annexed colorbars) are displayed. As appreciated by visual inspection, the skeleton of the network is not altered by the applied sparsification. To monitor how the eigenvalues get redistributed within the bounded domain to which they belong, we introduce the following indicators: The quantity I x measures the dispersion along the imaginary axis, by weighting the squared distance of each eigenvalue from the horizontal axis. Conversely, I y reflects the scattering of the eigenvalues about their mean, in the direction of the real axis. In Fig. 2, I x and I y , normalized to their respective values obtained before application of the sparsification algorithm, are shown against N , an indicator of the size of the generated networks.
The sparsification procedure shrinks the eigenvalues in the x-direction, while the opposite tendency is observed for the distribution along the y-direction.
The second sparsification method implements a more stringent constraint. Just like before, we select the links with weights in the range (-σ, σ), where σ acts as a small threshold amount. Unlike with the former case, we now eliminate the selected link only if the change produced on the modulus of each of the N eigenvalues is smaller than δ, namely if |Λ j −Λ j | < δ, for j = 2, ..., 2N + 1. In Fig. 3, the eigenvalues obtained after the sparsification algorithm are plotted (red circles) for two choices of the cutoff δ. The number of links that can be effectively removed grows with δ, the size of the allowed perturbation, as clearly demonstrated in Fig. 4.
Summing up, we have developed and tested a procedure to generate a network which returns an associated Laplacian matrix with a prescribed complex spectrum. The weighted network obtained following the above procedure is, in general, fully connected. Dedicated sparsification strategies can be applied to remove the links which carry a small weight, and bear a modest imprint on the ensuing Laplacian spectrum. In the following, we will consider a specific setting of the aforementioned generation scheme, which makes it possible for the Laplacian elements to be computed analytically.

IV. FOCUSING ON THE SPECIAL CASE U = qI
In the previous sections we described a general method to generate a Laplacian matrix which displays a designated spectrum. The method assumes a generic matrix U , which can be randomly assigned. In the following, we will focus on the specific case where U is proportional to the identity matrix and progress with the analytic characterization of the obtained Laplacian. As we shall argue in the following, working in this framework allows us to derive a set of closed conditions for constraining the weights of the underlying network to strictly positive values. To proceed in this direction we set: where I stands for the identity matrix and q is scalar.  A straightforward calculation returns the following expression for matrix S: while w is a uniform vector with identical entries equal to N a + ((N − 1)a − b)i and the quantities d, a, b are specified by The diagonal elements satisfy: where use has been made of the identity ℜ(D kk ) = α k . The first row and column are given by: while the remaining elements read: The Laplacian matrix obtained with the procedure illustrated above has both positive and negative entries. Signed Laplacians are often used in consensus problems, where nega-tive weights model antagonistic interactions. In other contexts, when e.g. the Laplacian is stemming from diffusive interactions, non diagonal entries are constrained to positive values. In the following, we will provide a set of necessary conditions for the assigned spectrum to eventually yield a Laplacian with positive extra diagonal elements, L ij > 0 for i = j.
Clearly, L ii < 0, as summing on the rows should return zero. The underlying network, therefore, displays positive weights, as its adjacency matrix is basically obtained from the Laplacian matrix by replacing the diagonal elements with zeros.
Further, we will set ℜ(Λ j ) = α j < 0 for j ≥ 2, an assumption which corresponds to dealing with a stable linear system of the type given in Eq. (25). This requirement immediately yields L 11 < 0, as it follows from relation (33). We will also operate in the setting analyzed above, i.e. assuming U = qI. The obtained expressions for the Laplacian elements allow to recast the sought conditions on their signs as: where the inequalities hold for t = 2, . . . , N + 1.
The above conditions can be simplified, after some algebraic manipulations, to return and The interested reader can access the detailed steps in Appendix C.
Assume that the assigned Laplacian spectrum matches the above conditions, while having α k < 0 for k > 2. Then the Laplacian matrix obtained with the above procedure , with U = qI, displays positive non diagonal entries.
Let us explore the consequences of conditions (47)  In the following section, we will apply the method of Laplacian generation here discussed to the study of two prototypical examples of dynamical systems on networks.

V. SELECTED APPLICATIONS
In this Section we consider two different models of interacting oscillators, defined on a network. In both cases, the coupling between individual oscillators is implemented via a discrete Laplacian operator, which reflects the specific network arrangement. It will be shown that a suitable network arrangement can be a priori established, building on the procedure illustrated above, so as to make the inspected systems stable against external perturbations.

A. Coupled Stuart-Landau oscillators
Consider an ensemble made of 2N + 1 nonlinear oscillators and denote by W i their associated complex amplitude. We assume the oscillators to be mutually coupled via a diffusive-like interaction which is mathematically modeled by a discrete Laplacian operator.
Each oscillator obeys a complex Stuart-Landau equation. The dynamics of the system can be cast in the form: where c 1 and c 2 are real parameters. The index j runs from 1 to 2N + 1, the total number of oscillators. Here, K is a suitable parameter setting the coupling strength. Without loss of generality, in what follows it is assumed that K = 1.
A ij is the generic entry of the directed and weighted adjacency matrix A and k i = j A ij .
The system admits a homogeneous limit cycle solution in the form W LC (t) = e −ic 2 t . To characterize the stability of the cycle, one can introduce a non homogeneous perturbation in polar coordinates as: By linearizing around the limit cycle solution (ρ i (t) = 0, θ i (t) = 0), one gets: To proceed further, expand the perturbations ρ j and θ j on the Laplacian eigenvectors basis, that is By inserting this expansion in (51) and using the relation for α = 1, . . . , 2N + 1, we obtain a condition formally equivalent to the expression of the continuous dispersion relation If λ Re = ℜ(λ max ) is positive for some Λ (α) , the perturbation grows exponentially in time, and the initial homogeneous state proves unstable. Conversely, if λ Re = ℜ(λ max ) < 0, for every Λ (α) , the perturbation gets re-absorbed and the system converges back to the fully synchronized state. The condition λ Re < 0 can be further processed analytically, as discussed in [15]. In particular, it can be shown that the latter condition is fulfilled, if the Laplacian eigenvalues fall in a specific portion of the parameter plane, which reflects the choice made for the reaction parameters c 1 , c 2 and K. The region of interest is the one enclosed between the two solid lines, displayed in Fig. 6(a) for the specific selection of the parameters. The blue symbols depicted in Fig. 6(a) are randomly generated so as to fall in the region of the complex plane where stability holds. They represent the spectrum of the Laplacian that we seek to recover following the method illustrated above. In Fig. 6(b), λ Re is plotted against −Λ Re = −ℜ(Λ), confirming the stability of the homogeneous solution.
We now proceed by generating a Laplacian matrix, which is constructed so as to yield the spectrum depicted in Fig. 6(b). From this, we compute the corresponding adjacency matrix A and use it to define the interactions between coupled oscillators, as follows from Eqs. (49). We then integrate numerically the governing equations, assuming the initial state to be a perturbation of the homogenous synchronized equilibrium. As expected, the perturbation fades away and the system regains its unperturbed, fully synchronized, equilibrium.
As a second example we set to study the Kuramoto model. Consider a system made of 2N + 1 oscillators, denote by θ i the phase of the i-th oscillator, and ω i its natural frequency.
The oscillators evolve as dictated by the following system of 2N + 1 coupled differential equations:θ Here, A ij stands for the entries of the adjacency matrix A which sets the interactions between pairs of oscillators. The matrix is, in principle, weighed, and may display positive and negative entries as reflecting the specific interaction (excitatory or inhibitory) being at play.
As an additional assumption, we will here focus in the simplified setting where ω i = ω ∀i. We can then introduce the new variable ψ i = θ i − ωt, and write the governing equation in the equivalent form:ψ A homogeneous solution always exists with ψ i = Ψ ∀i, and for any constant Ψ ∈ [0, 2π), as it can be immediately checked by substitution. To assess the stability of the solution, one sets ψ i = Ψ + δ i , and expands (56) at the leading order in the δ i . In this way, one gets: where L ij = A ij − k i δ ij is the Laplacian operator which stems from A ij . The stability of the simplified Kuramoto model here considered is controlled by a linear system of the type introduced in (25), with the obvious replacement of x i with δ i . The system proves hence stable if the (non trivial) eigenvalues of the Laplacian operator display negative real parts. Our aim, here, is to generate a Laplacian (and therefore a matrix of binary weighted connections among oscillators) which warrants the stability of the system. To this end we assign the eigenvalues (which appear in conjugate pairs) to belong to the negative portion of the complex plane, see Fig. 7(a). The null eigenvalue is clearly included into the spectrum.
Running the procedure discussed in the first part of the paper, we obtain the corresponding Laplacian and compute the associated adjacency matrix. The Kuramoto model (56) is then integrated numerically by assuming the recovered expression for A. As predicted, the system is stable to external perturbations as one can clearly appreciate by inspection of Fig. 7(b).
which proves the results.
Appendix B: About the computation of L ij The aim of this section is to detail the computations needed to obtain explicitly the entries of the Laplace matrix L under the assumption U = qI starting thus from Eqs. (16) and (18).
Let us begin by computing the diagonal elements of L, namely L ii for i = 1, .., N , . For i = 1, one gets: From (4) and (6) we obtain for k = 2, . . . , N + 1: Then where use has been made of the identity ℜ(D kk ) = α k .
Focus now on the conditions for β t . We assume that the following condition holds: and we set to explore its consequences. Eq. (C4) yields: For N > 1, conditions (C3) maps therefore in the following equivalent system: Remark that: due to the inequality which holds for each N > 1. Hence, system (C6) takes the form: Then, (C9) has no solutions, under the working hypothesis that we have put forward to deriving it. In fact: and summing on every t we get which, in turn, implies 4N 2 4N 2 −1 < 1, a condition that is obviously never met. We now go back to revise ansatz (C4), and consider the alternative scenario: Following a path analogous to the one discussed above, we get: and