Spectra of random networks with arbitrary degrees

TL;DR

This paper introduces a message passing approach to efficiently compute the spectral density of large, locally tree-like random networks with arbitrary degree distributions, including bounds on spectral edges.

Contribution

It develops a novel message passing method and an approximation for spectral density in configuration model networks, enabling fast analysis of their spectral properties.

Findings

The approximation performs well for most network densities.

Bounds on spectral band edges help identify structural phase transitions.

The method is applicable to a wide range of degree distributions.

Abstract

We derive a message passing method for computing the spectra of locally tree-like networks and an approximation to it that allows us to compute closed-form expressions or fast numerical approximates for the spectral density of random graphs with arbitrary node degrees -- the so-called configuration model. We find the latter approximation to work well for all but the sparsest of networks. We also derive bounds on the position of the band edges of the spectrum, which are important for identifying structural phase transitions in networks.

Equations

$$\rho (x)=\frac{1}{n}\sum _{i=1}^n\delta (x-{\lambda }_i),$$

$$\delta (x)=\underset{ϵ\to 0^{+}}{\lim }\frac{ϵ/\pi }{x^2+{ϵ}^2}=-\frac{1}{\pi }\underset{ϵ\to 0^{+}}{\lim }\mathrm{missing}Im\frac{1}{x+\mathrm{i}ϵ},$$

$$\rho (x)=-\frac{1}{n\pi }\underset{ϵ\to 0^{+}}{\lim }\mathrm{missing}Im\sum _{i=1}^n\frac{1}{x-{\lambda }_i+\mathrm{i}ϵ}.$$

$$\rho (z)=-\frac{1}{n\pi }\sum _{i=1}^n\frac{1}{z-{\lambda }_i}=-\frac{1}{n\pi }\mathrm{missing}Tr(z\mathbf{I}-\mathbf{A}{)}^{-1},$$

$$\rho (z)=-\frac{1}{n\pi z}\sum _{k=0}^{\infty }\frac{\mathrm{missing}Tr{\mathbf{A}}^k}{z^k}.$$

$$n_{2r}^{uv}$$

$$\displaystyle\hskip 40.00006pt{}\times\delta\biggl{(}r-1,\sum_{i=1}^{m}r_{i}\biggr{)},$$

$$h^{uv}(z)=\sum _{r=1}^{\infty }\frac{n_{2r}^{uv}}{z^{2r}}.$$

$$h^{uv}(z)$$

$$=\frac{1}{z^2}\sum _{m=1}^{\infty }\prod _{i=1}^m\sum _{\begin{array}{c}{w}_i\in N_v\\ w_i\ne u\end{array}}{h}^{v{w}_i}(z)$$

$$\displaystyle={1\over z^{2}}\sum_{m=1}^{\infty}\biggl{[}\sum_{\begin{subarray}{c}w\in N_{v}\\ w\neq u\end{subarray}}h^{vw}(z)\biggr{]}^{m}.$$

$$h^{uv}(z)=\frac{1/z^2}{1-{\sum }_{\begin{array}{c}w\in N_v\\ w\ne u\end{array}}{h}^{vw}(z)}.$$

$$n_{2r}^u$$

$$\displaystyle\hskip 40.00006pt{}\times\delta\biggl{(}r,\sum_{i=1}^{m}r_{i}\biggr{)}.$$

$$g^u(z)$$

$$\displaystyle=\sum_{m=1}^{\infty}\prod_{i=1}^{m}\sum_{v_{i}\in N_{u}}h^{uv_{i}}(z)=\sum_{m=1}^{\infty}\biggl{[}\sum_{v\in N_{u}}h^{uv}(z)\biggr{]}^{m},$$

$$g^u(z)=\frac{1}{1-{\sum }_{v\in N_u}{h}^{uv}(z)}.$$

$$\rho (z)=-\frac{1}{n\pi z}\sum _{u=1}^n{g}^u(z)=-\frac{1}{n\pi z}\sum _{u=1}^n\frac{1}{1-{\sum }_{v\in N_u}{h}^{uv}(z)}.$$

$$z^2{h}^{uv}(z)=1+z^2{h}^{uv}(z)\sum _{\begin{array}{c}w\in N_v\\ w\ne u\end{array}}{h}^{vw}(z),$$

$$z^2{h}_k(z)$$

$$=1+z^2{h}_k(z)\frac{k-1}{nk{p}_k}\sum _{v:k_v=k}\sum _{w\in N_v}{h}^{vw}(z).$$

$$\frac{1}{nk{p}_k}\sum _{v:k_v=k}\sum _{w\in N_v}{h}^{vw}(z)\to h(z).$$

$$h_k(z)=\frac{1/z^2}{1-(k-1)h(z)}.$$

$$q_k=\frac{(k+1)p_{k+1}}{c},$$

$$h(z)$$

$$=\frac{1}{z^2}\sum _{k=0}^{\infty }\frac{q_k}{1-kh(z)},$$

$$g_k(z)=\frac{1}{1-kh(z)}.$$

$$\sum _{u=1}^n{g}^u(z)=\sum _{k=0}^{\infty }n{p}_k{g}_k(z)=n\sum _{k=0}^{\infty }\frac{p_k}{1-kh(z)},$$

$$\rho (z)=-\frac{1}{\pi z}\sum _{k=0}^{\infty }\frac{p_k}{1-kh(z)}.$$

$$h(z)=\frac{1}{z^2}{\int }_0^{\infty }\frac{q_k\mathrm{d}k}{1-kh(z)},\rho (z)=-\frac{1}{\pi z}{\int }_0^{\infty }\frac{p_k\mathrm{d}k}{1-kh(z)},$$

$$h(z)=\frac{1/z^2}{1-(c-1)h(z)}.$$

$$h(z)=\frac{1\pm \sqrt{1-4(c-1)/z^2}}{2(c-1)}.$$

$$\rho (z)=\frac{1}{\pi z[ch(z)-1]},$$

$$\rho (x)=(c/2\pi )\frac{\sqrt{4(c-1)-x^2}}{c^2-x^2},$$

Full text

I Introduction

The spectral properties of the adjacency matrix of a network play a central role in the analysis of network structure, for instance in the eigenvector centrality Bonacich87 , in graph partitioning and community detection Fiedler73 ; PSL90 ; Newman06b , in the theory of dynamical systems on networks PG16 , and in the analysis of structural phase transitions such as percolation BBCR10 ; KNZ14 , localization CM11 ; MZN14 , and detectability transitions DKMZ11a ; NN12 .

There exist well-known computer algorithms for calculating the eigenvalue spectrum of matrices, which can be applied directly to the adjacency matrix of a network, but they are numerically demanding, taking time of order $n^{3}$ to calculate all eigenvalues of an $n$-node network. Furthermore, they give only the spectrum of a single network instance, where in many cases we would like to calculate the average spectrum of an entire graph ensemble. In this paper we develop an alternative approach to calculating network spectra, based on message passing methods and focusing particularly on the case of the so-called configuration model—a random graph with given degree sequence and one of the standard models in the theory of complex networks. We first develop a general message passing method, which is exact on arbitrary networks that are free of short loops and which is closely related to previous approaches for calculating graph spectra. Then we further develop an approximation to this method for the particular case of the configuration model in the limit of large network size that allows us to find closed-form expressions for the spectral density in certain cases and to perform numerical calculations in $\mathrm{O}(1)$ time in others.

There has been a significant amount of previous work on the spectra of complex networks. In early work, Farkas et al. FDBV01 computed spectra for a range of networks using standard numerical methods and demonstrated clear deviations of the spectral density from the Wigner form expected for traditional dense random graphs. Goh et al. GKK01a looked at networks with power-law degree distributions (a common feature of many empirical networks), giving both numerical results and analytic bounds for the largest eigenvalues. Dorogovtsev et al. DGMS03 gave an analytic prescription for calculating complete spectra of configuration model networks but were unable to solve the resulting equations. Instead therefore, they develop an approximation, which is similar in some respects to our own though different in both motivation and final form, and which appears to give good results in some cases but is rather inaccurate in others, such as networks with power-law degree distributions. A contrasting approach has been pursued by Semerjian and Cugliandolo SC02b who start from the classic formulation of Edwards and Jones EJ76 of the spectral density in terms of a Gaussian path integral and use a replica-type analysis to derive several different approximations to the spectra of a traditional Erdős–Rényi style random graph. Kühn Kuhn08 applied similar methods to the configuration model, using an analysis reminiscent of the Viana–Bray solution for a dilute spin-glass to derive an approximation to the spectrum that is similar in its accuracy to ours, although quite technically daunting. Rogers et al. RCKT08 , starting again from the Edwards–Jones formulation, derive a message passing method essentially equivalent to the one we use, which they employ as a numerical tool for computing spectra. Chung et al. CLV03 study a slightly different class of model networks, those with given expected degrees and statistically independent edges (sometimes called the Chung–Lu model following earlier work by some of the same authors CL02b ), for which they derive an expression for the single largest eigenvalue, which plays a role for instance in percolation calculations BBCR10 . Using tools from free probability theory, Nadakuditi and Newman NN13 calculated complete spectra for the same class of networks.

In outline this paper is as follows. In Section II we derive message passing equations for the spectral density which form the basis for subsequent developments. In Section III we derive our approximation to the message passing equations for the case of the configuration model and show how it can be used both to perform fast numerical calculations and, in some cases, to give analytic solutions for the spectral density. We also use it to derive bounds on the position of the edges of the spectrum, which play a central role in the theory of structural phase transitions in networks. In Section IV we give a number of example applications of our methods, showing the accuracy of the approach in some cases, as well as other cases where it breaks down. In Section V we give our conclusions and suggest some directions for future work.

II Spectral density of a locally tree-like network

Our calculations begin with the derivation of a system of message passing equations for computing the spectral density of a locally tree-like network, meaning any network where the local neighborhood of any node, out to any fixed distance, takes the form of a tree with probability one in the limit of large network size. This part of the calculation is similar to developments described previously by Dorogovtsev et al. DGMS03 and Rogers et al. RCKT08 , although our derivation is interesting in its own right because it employs only elementary algebraic methods, where previous approaches have relied on heavier mathematical machinery.

Suppose we are given a single undirected unweighted network of $n$ nodes and asked to calculate its spectral density, which is the function

[TABLE]

where $\delta(x)$ is the Dirac delta function and $\lambda_{1}\ldots\lambda_{n}$ are the $n$ eigenvalues of the adjacency matrix $\mathbf{A}$ of the network—the matrix with binary-valued elements $A_{uv}$ equal to 1 whenever there is an edge between nodes $u,v$ and 0 otherwise. Following a standard line of development we express the delta function as the limit of a Lorenzian (or Cauchy) distribution:

[TABLE]

where the notation $\lim_{\epsilon\to 0^{+}}$ means that the parameter $\epsilon$, which controls the width of the Lorenzian, tends to zero from above. Substituting (2) into (1), we find that

[TABLE]

It will be convenient for subsequent developments to generalize this spectral density into the complex plane, defining $z=x+\mathrm{i}\epsilon$ and

[TABLE]

where $\mathbf{I}$ is the identity. The standard spectral density for the network is the limiting value of the imaginary part of this function as $z$ tends to the real line from above. In fact, in many practical situations it is desirable to retain a small imaginary part for $z$, corresponding to $\epsilon>0$, producing a Lorenzian broadening of the delta-function peaks in the spectral density. For finite networks this gives us a smooth density function $\rho(z)$ rather than a set of spikes, effectively making a kernel density estimate of the spectral density with a Lorenzian kernel.

If we expand the matrix inverse in (4) as a geometric series $(z\mathbf{I}-\mathbf{A})^{-1}=z^{-1}\sum_{k=0}^{\infty}(\mathbf{A}/z)^{k}$ and take the trace term by term, we find that

[TABLE]

The quantity $\mathop{\mathrm{missing}}{Tr}\mathbf{A}^{k}$ is equal to the number of closed walks of length $k$ in the network—paths that start at any node and return there (not necessarily for the first time) exactly $k$ steps later. If we can count such closed walks on our network for all values of $k$ then we can compute the spectral density from Eq. (5). We do this counting using a message passing method.

II.1 Message passing

As we have said, our focus here is on locally tree-like networks, meaning networks in which local neighborhoods are trees, having a vanishing density of short loops. The absence of loops means that any closed walk must necessarily start and end by traversing the same edge—it cannot return to the starting node by any edge other than the one it left by, since in so doing it would complete a loop in the network, of which there are none. Indeed every edge in a closed walk on a locally tree-like network must, for the same reason, be traversed twice, once in each direction, or more generally the same number of times in both directions. This in turn means that all closed walks have an even number of steps.

Let $n^{uv}_{2r}$, with $r$ a positive integer, be the number of closed walks that begin by traversing the edge from $u$ to $v$ and end, after exactly $2r$ steps, by traversing the same edge back again from $v$ to $u$ for the first time. Other edges may be traversed any (even) number of times, but the edge between $u$ and $v$ is traversed only once each way.

The smallest possible value of $r$ in this scenario is 1, for which $n^{uv}_{2}=1$ trivially. For all higher values $r>1$, we can write a self-consistent expression for $n^{uv}_{2r}$ thus:

[TABLE]

where $N_{v}$ denotes the set of neighbor nodes of $v$ and $\delta(i,j)$ is the Kronecker delta. To break it down, this expression says that a walk starting and ending along the edge $uv$ makes some number $m$ of subsequent excursions from node $v$ each of which has a first step to one of the neighbors of $v$ other than $u$, that the total number of such walks is the product of the numbers of distinct excursions, and that the individual lengths $2r_{1}\ldots 2r_{m}$ of these excursions necessarily sum to $2r-2$.

To solve this system of equations for the numbers $n^{uv}_{2r}$, we define the useful quantity

[TABLE]

Substituting for $n^{uv}_{2r}$ from Eq. (6) into Eq. (7) and performing the sum over $r$, we get

[TABLE]

The remaining sum over $m$ is a simple geometric series, which can be completed to give

[TABLE]

This is our fundamental message passing equation. We can think of $h^{uv}(z)$ as a message passed from node $v$ to its neighbor $u$, whose value can be computed from the values of the messages received by $v$ from its other neighbors $w$.

If one can compute the values of the messages $h^{uv}(z)$ for any given value of $z$, one can compute the spectral density itself as follows. By analogy with Eq. (6), the number $n^{u}_{2r}$ of closed walks of length $2r$ that start and end at node $u$ can be written as

[TABLE]

Then we define

[TABLE]

and hence

[TABLE]

The spectral density, Eq. (5), can now be written in terms of this quantity as

[TABLE]

Between them, Eqs. (9) and (13) give us our prescription for calculating the spectral density. Note that the variable $z$ enters Eq. (9) only as $z^{2}$, which means that the spectral density will always be symmetric about the origin.

The message passing equations can be used as a numerical tool for computing network spectra, as Rogers et al. RCKT08 do with the equivalent equations they derive. One simply chooses initial values of the messages (for instance at random) and iterates Eq. (9) repeatedly until convergence is achieved. Direct solution of the message passing equations, however, will not be our primary goal here.

III Degree-based approximation

The developments presented so far provide an elementary derivation of a message passing method for calculating network spectra. Though interesting, however, the method derived performs essentially the same calculation as the previously proposed method of Rogers et al. RCKT08 and in this sense does not add much to our toolkit. In this section, however, we go further, focusing specifically on the configuration model and introducing an approximation to the message passing algorithm that reduces its accuracy very little while making it enormously faster and, in some cases, allowing us to compute analytic solutions for the spectral density.

The configuration model MR95 ; NSW01 , a random graph model with arbitrary node degrees, is one of the most fundamental of network models. To generate a configuration model network one fixes the degree of each node separately, then connects nodes at random while respecting the degrees. The configuration model is widely used both for understanding network structure in general and as a starting point for further calculations of network properties and processes.

The configuration model generates networks that are locally tree-like and hence the message passing approach can be applied to them. We do this in Fig. 1 for an example network in which nodes have just two distinct degrees, 5 and 10, with equal probability. The figure shows a scatter plot in the complex plane of the resulting values of the messages $h^{uv}(z)$ at the point $z=3+\mathrm{i}\epsilon$. As the figure shows, the values in this case form two distinct, compact, non-overlapping clouds, which correspond to the degrees of the nodes $v$ “sending” the messages: the right-most cloud corresponds to nodes of degree 5 while the left-most one corresponds to nodes of degree 10. This simple observation suggests a possible approximation to the message passing equations, in which we approximate each message with the mean or centroid value of the cloud to which it belongs (shown by the plus symbols in the figure), in effect assuming that the messages are a function of degree only. This approximation, as we will see, turns out to give remarkably accurate estimates of the spectral density in many cases.

Beyond numerical results like those in Fig. 1, there is some precedent for a degree-based approximation of this kind to be found in the previous literature on complex networks. In studies of epidemic models, for instance, which can be written in the form of a message passing process KN10a , Pastor-Satorras and Vespignani PV01a ; PV01b found a similar degree-based approximation to work well. The approximation also bears some conceptual similarity to the “effective medium approximation” introduced by Semerjian and Cugliandolo SC02b for the Erdős–Rényi random graph (building on previous work by Biroli and Monasson BM99 ) and by Dorogovtsev et al. DGMS03 for the configuration model, although the details of the functional forms are different in both cases.

So consider a configuration model with given node degrees such that the fraction of nodes with degree $k$ is $p_{k}$. Our approximation consists of replacing each message $h^{uv}(z)$ by the mean value $h_{k}(z)$ of messages sent from nodes with the same degree $k$ as node $v$. Rearranging our message passing equation, Eq. (9), as

[TABLE]

making the replacement $h^{uv}(z)\to h_{k}(z)$, and averaging over all edges $(u,v)$ where $v$ has degree $k$, of which there are $nkp_{k}$, we get

[TABLE]

In a configuration model network, however, the degrees of adjacent nodes $v,w$ are uncorrelated, so the average of $h^{vw}(z)$ over many different $w$ is independent of $k$ and, in the limit of large network size, simply equal to the average message in the network as a whole, which we will denote $h(z)$:

[TABLE]

Thus we have $z^{2}h_{k}(z)=1+z^{2}(k-1)h_{k}(z)h(z)$ or

[TABLE]

Furthermore, the average message $h(z)$ is itself equal to the average of $h_{k}(z)$ over all degrees $k$, but here we must be careful. The node $v$ appearing in the message $h^{uv}(z)$ is by definition reached by following an edge from node $u$, and the degrees of nodes reached by following an edge are drawn not from the overall degree distribution $p_{k}$ of the network as a whole, but from the modified distribution $kp_{k}/c$ where $c=\sum_{k}kp_{k}$ is the average degree NSW01 . Commonly this is expressed in terms of the so-called excess degree distribution:

[TABLE]

which is the probability distribution of the number of edges attached to the node other than the one we followed to reach it. In terms of this quantity, the average message is given by

[TABLE]

where we have made the change of variables $k\to k+1$ in the final equality.

We can use this result to calculate the spectral density itself by making a similar degree-based approximation to the quantities $g^{u}(z)$ appearing in Eq. (12). We assume that $g^{u}(z)$ is well approximated by the mean $g_{k}(z)$ of $g^{u}(z)$ over all nodes $u$ with degree $k$. Rearranging Eq. (12), making this approximation, and averaging again over nodes of degree $k$, we find that $g_{k}(z)=1+kg_{k}(z)h(z)$, or

[TABLE]

Summing over all nodes, we then get

[TABLE]

and substituting this result into Eq. (13) we find that, within this approximation, the spectral density is given by

[TABLE]

Between them, Eqs. (19) and (22) now give us a complete formula for calculating the spectral density. Notice that, by contrast with the original message passing method, these equations do not depend on the precise form of the network, or even on its size—they are a function only of the degree distribution. Moreover, the average over messages becomes more and more accurate as network size grows, so that in effect Eqs. (19) and (22) give us (an approximation to) the spectrum of the configuration model in the limit of large size.

Equation (19) can be solved for $h(z)$ for any given value of $z$ by simple iteration, starting from a suitable initial value, such as $h(z)=0$. Then we substitute the result into Eq. (22) to get $\rho(z)$. Note that for each value $z$ at which we want to calculate the spectral density we need now iterate only one equation, Eq. (19), in contrast to the full message passing method of Eq. (9), which involves iterating equations equal in number to twice the number of edges in the network, which can be thousands or millions in real-world situations. Equation (19) still contains a sum over degrees $k$, but this sum has only as many terms as there are distinct degrees in the network, which in most cases is far smaller than the number of edges. Indeed, as we have said the equations given here apply in the limit of infinite size, for which the direct iteration of Eq. (9) would of course be impossible.

As an example application of our approach consider Fig. 2, which shows the spectrum of a network drawn from the same model as Fig. 1, in which nodes have just two distinct degrees, 5 and 10, with equal probability. The histogram in the figure shows the results of a direct numerical diagonalization of the adjacency matrix of a single such network with $10\,000$ nodes, while the solid curve shows the spectral density calculated from the full message passing method, Eq. (13). The dashed curve shows the spectrum calculated using our degree-based approximation, Eqs. (19) and (22). As we can see, the agreement between all three calculations is excellent, and in particular there is barely any perceptible difference between the full and approximate versions of the message passing calculation. To the extent that the two differ, it is mostly in the center of the figure in the region close to $x=0$, a pattern that we will see repeated in other examples.

Note moreover that, while the full message passing calculation took several minutes of computer time to complete, the approximate calculation took only a few seconds, essentially all of which was spent on iterating Eq. (19). A faster calculation still might be possible by using a more efficient method of solution than simple iteration, such as Newton’s method.

As an aside, we note that essentially the same method of calculation can be applied to the Chung–Lu model, the model in which the edges are independent random variables and only the expected degrees of the nodes are fixed, not the actual degrees CL02b . This model also produces locally tree-like networks, meaning that the message passing method is applicable, and furthermore one can usefully approximate the message passing equations by assuming that messages originating at all nodes with the same expected degree take the same value. Following a similar argument to that for the configuration model above, we then arrive at the following equations for the spectral density:

[TABLE]

where $k$ now represents not the actual degree of a node but its expected degree, which can take any non-negative real value.

Interestingly, after some translation of notation these are exactly the same equations that were derived previously for the Chung–Lu model using methods from free probability theory NN13 . Superficially, they look identical to Eqs. (19) and (22) for the configuration model, apart from the replacement of the sums by integrals, but the similarity is deceptive. The definition of the excess degree distribution $q_{k}$ for the Chung–Lu model is different from that for the configuration model. (For the Chung–Lu model the correct definition is $q_{k}=kp_{k}/c$, compared with $q_{k}=(k+1)p_{k+1}/c$ for the configuration model.) This can make a substantial difference to the shape of the spectrum and if one simply applies the solution for the Chung–Lu model to the configuration model the results are quite poor, particularly for sparse networks. Nonetheless, the degree-based approximation does give a more straightforward derivation of the equations for the Chung–Lu model, requiring less advanced techniques than the free probability approach. (The reverse procedure does not seem to work, however: there is no obvious way to use free probability theory to derive the equations for the configuration model. The sticking point is that the adjacency matrix and the degree sequence are not “free” with respect to one another in the free probability sense.)

III.1 Analytic solutions

In some cases Eq. (19) allows a closed-form solution for $h(z)$, depending on the form of the degree distribution. As a simple example consider a regular network, meaning one in which all nodes have the same degree $c$, for which Eq. (19) takes the form

[TABLE]

This can be rearranged into the form of a quadratic equation $(c-1)h^{2}-h+1/z^{2}=0$, with solutions

[TABLE]

Then

[TABLE]

and, after taking the imaginary part and going to the real line, we find the spectral density to be

[TABLE]

which recovers the standard Kesten–McKay distribution for a random regular graph McKay81 . Thus in this case our approximation is not an approximation at all: the Kesten–McKay distribution gives the exact spectral density of a random regular graph in the limit of large size. The calculation is exact because in a regular graph the neighborhood of every node has the same network structure, so all messages are in fact exactly equal.

For a slightly more complex example consider again a network of the kind in Fig. 2, with nodes of two different degrees, which we will denote $a$ and $b$. For such a network Eq. (19) becomes

[TABLE]

which gives a cubic equation for $h(z)$, which is solvable in closed form though the solution is complicated and we will not reproduce it here. Applied, for instance, to the example network in Fig. 2, it gives essentially the same result as our numerical solution for the same system.

More broadly, we can define a moment generating function $\mu_{p}(x)$ for a distribution $p_{k}$ to be the power series in $x$ whose coefficients are the moments $\left\langle k^{r}\right\rangle_{p}$ of $p_{k}$ thus:

[TABLE]

In terms of such generating functions, the fundamental equation (19) can be written

[TABLE]

and Eq. (22) becomes

[TABLE]

If one knows the moment generating functions for a particular degree distribution, one can use these equations to derive the spectral density.

III.2 Position of the band edges

One issue of particular interest is the position of the edges of the band of nonzero spectral density in the spectrum of a network. The upper band edge plays a role for instance in locating localization transitions in networks MZN14 and the position of the so-called detectability threshold for community detection DKMZ11a ; NN12 .

To understand why there is a finite band at all, and where its edges fall, consider Fig. 3, which sketches a graphical solution to the fundamental equation (19) for $h(z)$. We rearrange the equation in the form

[TABLE]

and then separately plot the left and right sides of the equation as a function of $h$. Where the two cross are our solutions for $h(z)$.

The right-hand side of the equation, represented by the rising curves in Fig. 3, has simple poles at $h=1/k$ for all nonzero values $k$ of the excess degree that occur in the network. In the figure we assume that values $k=0\ldots 4$ occur, so that there are poles at $h=1$, $\frac{1}{2}$, $\frac{1}{3}$, and $\frac{1}{4}$. Meanwhile, the left-hand side $z^{2}h$ takes the form of a simple straight line through the origin with slope $z^{2}$. Lines are sketched in the figure for three different values of $z^{2}$. The solutions of (32) are marked with dots.

Generically, Eq. (32) gives us a polynomial equation of degree $m+1$, where $m$ is the number of distinct nonzero values of the excess degree, or equivalently the number of poles in the right-hand side of the equation. This in turn means the equation has $m+1$ solutions for $h(z)$. As we can see from the geometry of Fig. 3 there are always at least $m-1$ real solutions, one in each of the intervals between the $m$ poles. Real solutions, however, cannot give us a nonzero spectral density—we need complex $h(z)$ to get a nonzero density when we take the imaginary part of Eq. (22). Thus we focus on the remaining two solutions, which can be either real or complex, with the band edge corresponding to the point at which complex solutions first emerge.

There are three different possible forms of the solution, corresponding to the three different values of $z^{2}$ depicted in Fig. 3. If $z^{2}$ is sufficiently large (the steepest line in the figure), the line $z^{2}h$ intersects the first segment of the curve representing the right-hand side of Eq. (32), giving us two real solutions for $h(z)$. If this happens then all $m+1$ solutions for $h$ are real, there are no complex solutions at all, and the spectral density is zero. In this regime we are outside the spectral band. Conversely, if $z^{2}$ is sufficiently small, then the line of $z^{2}h$ does not intersect the first segment, there are only $m-1$ real solutions, and the two remaining solutions must necessarily be complex, placing us inside the band. Between these two regimes lies the borderline case—the band edge—represented by the dashed diagonal line in the figure, which is precisely tangent to the first segment of the curve.

While it is difficult, or sometimes impossible, to solve exactly for the roots of polynomial equations, we can derive useful bounds on the position of the band edge by inspecting the geometry of Fig. 3. Note that the curve representing the right-hand side of Eq. (32) intercepts the vertical axis at $y=\sum_{k}q_{k}=1$, since $q_{k}$ is a properly normalized probability distribution. Thus the tangent point must fall at $y>1$. At the same time the horizontal coordinate of the tangent point must satisfy $h<1/K$, where $K$ is the largest value of the excess degree $k$ in the network (which is 1 less than the largest value of the ordinary degree). Thus the critical slope of the tangent line at the band edge satisfies

[TABLE]

meaning that the upper band edge falls at a point $z\geq\sqrt{K}$ and the lower one falls at $z\leq-\sqrt{K}$. For the network in Fig. 2, for example, which has nodes of degree 5 and 10 only, the largest excess degree is $K=9$, and hence the upper band edge satisfies $z\geq 3$. From an inspection of Fig. 2 this appears to be correct—the band edge looks to fall at around $z=5.5$. It is not a very good bound, though it is interesting nonetheless, since it implies that in a network in which node degrees are unbounded there will be no upper edge to the eigenvalue spectrum: the band edge diverges as $K$ diverges. In a network with an exponential or power-law degree distribution, for instance, there will be no upper limit to the spectral band.

We can derive a better bound on the position of the band edge by first computing a lower bound on the right-hand side of (32) in the region $h<1/K$ thus:

[TABLE]

and the value of $z^{2}$ at the band edge is always greater than the slope of the tangent line to this curve. The tangent falls at the point where there is a double root of

[TABLE]

which is equivalent to the quadratic equation

[TABLE]

This has a double root when its discriminant vanishes, which gives us another quadratic equation $z^{4}-2(1+q_{K})Kz^{2}+(1-q_{K})^{2}K^{2}=0$ for $z^{2}$, whose solution now gives us an improved bound on the position of the band edge:

[TABLE]

Taking the example of our network with equal fractions of nodes of degree 5 and 10 again, so that $K=9$ and $q_{K}=\frac{2}{3}$, we find that $z\geq 3(1+\sqrt{2/3})\simeq 5.449$.

We can also derive an upper bound on the position of the band edge by rearranging Eq. (32) into the form

[TABLE]

For any real $h<1/K$ this gives an upper bound on the value of $z^{2}$ at the band edge, with values of $h$ closer to the true double root giving better bounds. We use the value of $h$ at the double root of the approximation, Eq. (34), which is $h=1/[K(1+\sqrt{q_{K}})]$. For our network with nodes of degree 5 and 10 this yields an upper bound of $5.609$ on the position of band edge, meaning that the band edge falls in the relatively narrow interval $5.449\leq z\leq 5.609$. This interval is shown on Fig. 2 by the dashed vertical lines, and appears to agree well with both the exact position of the band edge from the full message passing calculation and the approximate position derived from the degree-based approximation.

IV Examples

In this section we look at some applications of our methods to example networks, illustrating the advantages and limitations of the approach.

Figure 4 shows spectra for four networks generated from configuration models with different degree distributions. Panel (a) shows the spectrum of a network with a Poisson degree distribution with mean degree $c=5$—effectively a standard Erdős–Rényi style random graph in which each possible edge is present with the same probability $c/n$. Since the elements of the adjacency matrix in such a network are independent identically-distributed random variables, one might expect the spectrum to follow the standard Wigner semicircle distribution. However, as shown previously by many authors EJ76 ; FDBV01 ; SC02b ; Kuhn08 , this is not the case when, as here, the graph is very sparse, meaning $c$ is small. The deviation from the Wigner law is clear in Fig. 4a, with the spectrum having a distinctly non-semicircular shape with a peak at $x=0$. The histogram in the figure shows the distribution of eigenvalues calculated by direct diagonalization of a single instance of the model with $n=10\,000$ nodes, while the solid curve shows the spectrum calculated using the full message passing method of Section II.1. The dashed curve shows the results of the degree-based approximation introduced here and, as we can see, it works well in this case, being barely distinguishable from the full calculation.

The Poisson distribution, however, is not a good approximation to the degree distributions of most real-world networks, which are typically strongly right-skewed BA99b ; ASBS00 . Many networks are observed to have degree distributions that approximately follow a power law or Pareto distribution. Figure 4b shows the spectrum for such a network, which displays characteristic long tails with no clear band edge, as noted previously for instance in GKK01a , where it is shown that the leading eigenvalue of a power-law network scales as $n^{1/4}$ with system size, and is as a result unbounded in the limit of large $n$. Again the degree-based approach does a good job of approximating the spectrum of the network, although there are now clear deviations visible for small absolute values of $x$ (i.e., close to the origin). In particular, notice that, by contrast with panel (a), the peak at $x=0$ is not well reproduced by our approximation (but is captured by the full message passing calculation).

In general, while we find that the degree-based approximation does well in many cases, it is weakest when the network is particularly sparse. In panel (c) of the figure, for example, we show results for a network with nodes of degree 2, 3, and 4 only, in equal numbers. This network is now very sparse—the average degree is only 3—and differences between the true spectrum (solid line) and approximation (dashed line) are becoming visible.

Panel (d) shows a particularly extreme case, of a network with nodes of degree 1, 2, and 3 only for which the spectrum becomes quite ill-behaved with notable spikes and other irregularities. While these are once again reproduced faithfully by the full message passing calculation, our approximation fails to capture them, and moreover does a relatively poor job of the overall shape of the spectrum. The spikes in the spectrum are due to the presence of numerous small components in the network. The peaks at $\pm 1$ are due to components of size 2, for example. Our approximation fails to pick these out because it does not distinguish between messages in small components and those in the giant component. The method does, however, still capture the basic shape of the spectrum further from the origin, and gets the positions of the band edges approximately correctly.

V Conclusions

In this paper we have derived a set of message passing equations that allow one to calculate the spectrum of the adjacency matrix of an arbitrary, tree-like network. For the particular case of the configuration model, we approximate these equations by assuming that the messages are a function of node degree only, which gives a much simplified form for the spectral density that applies in the limit of large network size. Test applications of these methods on a range of networks shows that the approximation works well for all but the sparsest of networks, and where it does show deviations from the true spectral density it is typically in the region close to the origin and away from the band edges. We have also used our approximate equations as a starting point for deriving bounds on the positions of the edges of the spectrum, and in particular the upper edge of the continuous spectral band, which plays a role in determining the locations of structural phase transitions in networks.

It is interesting to ask if there are better approximations we could make to the message passing equations than the one introduced here, and there are a number of possibilities. One could imagine an approximation where the clouds of points appearing in Fig. 1 are represented by more than one different value: a large cloud could be divided into two parts, each of which is approximated by its own individual centroid. Another possibility is that we could attempt to represent the messages for the lowest degree nodes more accurately in some way. Most of the error in our approximation is in the low-degree nodes, since these are the ones for which the sum of incoming messages will show the largest statistical fluctuation. It is for this reason that the approximation works poorly for very sparse networks, such as the one in Fig. 4d. One could imagine, for instance, explicitly summing over the possible degrees of the neighbors of a low-degree node, producing a kind of two-step approximation for these nodes that would presumably be more accurate than the one-step approach we currently employ. These possibilities, however, we leave for future work.

Acknowledgements.

This work was funded in part by the NSF under grants DMS–1407207 and DMS–1710848 (MEJN) and by ONR grant N00014-15-1-2141, DARPA Young Faculty Award D14AP00086, and ARO MURI grants W911NF-11-1-0391 and 2015-05174-05 (RRN).