The spectral density of dense random networks and the breakdown of the Wigner semicircle law

A random network or graph is a collection of nodes joined by edges following a probabilistic rule. Random networks are formidable tools to model large assemblies of interacting units, like neurons in the brain, computers and routers in the Internet, or persons forming a friendship network [1]. Motivated by our increasing ability to collect and process vast amounts of empirical data, the theory of random networks has experienced an enormous progress, leading to important insights in physics, biology, and sociology [2]. The implications of the structure of networks to the dynamical processes occurring on them remains a fundamental topic in network theory [3].

Dynamical processes on a random network are to a large extent governed by the spectrum of the corresponding adjacency random matrix. This is a random matrix where each entry equals the strength of the interaction between a pair of nodes. The study of a broad range of problems amounts to linearize a large set of differential equations [4, 5, 6], coupled through a random network, around a stationary state, whose stability and transient dynamics are ultimately determined by the eigenvalue distribution of the adjacency matrix. Examples in this context are the study of the epidemic threshold for the spreading of diseases [7], the synchronization transition in networks of coupled oscillators [8, 9], and the functional stability of large biological systems, such as gene regulatory networks [10, 11], ecosystems [12, 13, 14], and neural networks[15, 16, 17]. The spectrum of the adjacency matrix also contains information about the network structure, since the trace of the adjacency matrix raised to a certain power yields the number of network loops of a given length [18]. Therefore, the problem of how the network architecture influences the spectrum of the adjacency matrix has attracted considerable research efforts [19, 20, 21, 22, 23, 24, 4, 25].

Synthetic models of random networks provide a controlled way to study the role of the network structure on the spectrum of the adjacency matrix. The degree sequence is the most basic tool to characterize the graph structure [1]. The degree of a node counts the number of edges attached to the node, while the degree sequence specifies all network degrees. When a network has an infinitely large number of constituents, it is natural to consider the degree distribution, i.e. the fraction of nodes with a certain degree, instead of dealing with the degree sequence. The configuration model stands out as one of the most fundamental and versatile models of random graphs [26, 27, 1, 28], since it enables the degree distribution to be freely specified, while keeping the pattern of interconnections entirely random. From a practical viewpoint, the configuration model resolves one of the main shortcomings of Erdös- Rényi random graphs [29, 30, 31], namely its Poisson degree distribution, which has little resemblance with the long-tailed degree distributions found in empirical networks [32, 33]. The priceless possibility to fix the degree distribution of a random graph is not only useful to model the structure of empirical networks, but it offers the ideal setting to examine the impact of degree fluctuations on the spectrum of the adjacency matrix.

There has been a significant amount of numeric and analytic work on the spectral properties of adjacency random matrices [34, 19, 35, 20, 21, 36, 37, 38, 22, 24, 23, 39, 40, 4, 25, 41]. A remarkable analytic result is the set of exact and mathematically rigorous equations determining the eigenvalue distribution of the configuration model [36, 37, 38, 42]. These equations form a natural starting point to study the impact of degree fluctuations on network spectra. Unfortunately, aside from a few particular cases [42, 39, 24], these equations can be analytically solved only in the dense limit [36, 37], when the mean degree becomes infinitely large and random networks approach a fully-connected structure.

The analytic solution for the dense limit of the aforementioned equations is at the root of perturbative and non-perturbative approximations for the eigenvalue distribution of large-degree random networks [34, 43, 20]. With the exception of graphs with a power-law degree distribution [19, 35, 21], whose moments are divergent, one expects that the eigenvalue distribution of undirected random networks converges, in the dense limit, to the Wigner semicircle distribution of random matrix theory, reflecting a high level of universality. This expectation has been confirmed numerically [19] and analytically [36, 37] for Erdös-Rényi and regular random graphs. However, the recent work [44] has rigorously shown that the spectrum of the configuration model does depend on the degree distribution in the dense limit. The main theorem in [44] also provides an approach to compute the spectral density for specific degree distributions. In spite of these rigorous results, a more detailed understanding of how the structure of a dense random network influences its spectrum, and the universal status of the Wigner semicircle law, is still lacking.

In this work we study the dense limit of the eigenvalue distribution of undirected networks drawn from the configuration model. We analyze four examples of degree distributions in which all moments are finite and the variances scale differently with the average degree, and we show that the dense limit of the eigenvalue distribution is determined by the degree fluctuations. In particular, we derive an equation yielding the eigenvalue distribution of dense random graphs with an exponential degree distribution [28], from which we unveil the existence of a logarithmic divergence in the spectral density and the absence of sharp spectral edges, i.e., the eigenvalue distribution of exponential random graphs is supported on the entire real line. These are remarkable differences with respect to the Wigner semicircle distribution of random matrix theory [45]. We also discuss how the analytic equations for the spectral density can be derived in two different ways: from the dense limit of the exact resolvent equations in [36, 37, 42], and from the main theorem proved in [44]. Based on an exact calculation of the fourth moment of the eigenvalue distribution, we obtain a sufficient condition for the breakdown of the Wigner semicircle law. This condition is given only in terms of the variance of the degree distribution. We finish the paper by proposing a classification of the different universal behaviours of the eigenvalue distribution of the configuration model in the dense limit, based on an exponent characterizing how the variance of the degree distribution diverges for increasing mean degree.

Our paper is organized as follows. In the next section we define the adjacency matrix of random networks and the spectral density. Section III introduces the configuration model and defines the degree distributions studied in this work. In section IV, we present numeric and analytic results for the dense limit of the spectral density for each example of degree distribution. Section IV explains how the analytic equations for the spectral density are obtained from the dense limit of the resolvent equations, while in section V we derive these equations from the main theorem in reference [44]. The condition for the breakdown of the Wigner semicircle law and the classification of the different universal behaviours are discussed in section VI. We summarize our results and conclude in section VII.

Let $N$ be the total number of network nodes. The set of binary random variables $\{C_{ij}\}_{i,j=1,\dots,N}$ specifies the network structure: if $C_{ij}=1$, there is an undirected link $i\leftrightarrow j$ between nodes $i$ and $j$, while $C_{ij}=0$ means this link is absent. We also associate a random weight $J_{ij}\in\mathbb{R}$ to each edge $i\leftrightarrow j$, which accounts for the interaction strength between nodes $i$ and $j$. We consider undirected simple random networks that have no self-edges nor multi-edges, in which $C_{ij}=C_{ji}$, $J_{ij}=J_{ji}$, and $C_{ii}=0$.

The degree $K_{i}=\sum_{j=1}^{N}C_{ij}$ of node $i$ is a random variable that counts the number of nodes attached to $i$, while the degree sequence $K_{1},\dots,K_{N}$ provides global information on the fluctuations of the network connectivity. In the limit $N\rightarrow\infty$, it is more convenient to work with the degree distribution

where $\delta$ is the Kronecker symbol. The average degree $c$ and the variance $\sigma^{2}$ of $p_{k}$ read

The degree distribution $p_{k}$ is one of the primary quantities to characterize the structure of random networks in the limit $N\rightarrow\infty$.

We will study the eigenvalue distribution of the $N\times N$ symmetric adjacency matrix $\bm{A}$, with elements defined as

The adjacency matrix fully encodes the network structure, along with the coupling strengths between adjacent nodes. The empirical spectral measure of $\bm{A}$ is given by

where $\lambda_{1},\dots,\lambda_{N}$ are the (real) eigenvalues of $\bm{A}$. By introducing the $N\times N$ resolvent matrix

with $z=\lambda-i\eta$ on the lower half complex plane, the empirical spectral measure follows from the diagonal elements of $\bm{G}(z)$

In the limit $N\rightarrow\infty$, the empirical mean of ${\rm Im}G_{ii}(z)$ normally converges to its ensemble averaged value, obtained from the distribution of $\bm{A}$, which implies that $\rho(\lambda)=\lim_{N\rightarrow\infty}\rho_{N}(\lambda)$ is well-defined. Here we will study the spectral density $\rho(\lambda)$ when the average degree grows to infinity, hence the limit $c\rightarrow\infty$ is performed after $N\rightarrow\infty$. The scaling of the elements $A_{ij}$ with the average degree $c$ in Eq. (2) ensures the spectral density has a finite variance when $c\rightarrow\infty$.

We study random networks with arbitrary degree distributions, known as the configuration model of networks [26, 27, 28, 1], where $p_{k}$ is specified at the outset. A single instance of the adjacency matrix of the configuration model is created as follows. First, the degrees $K_{1},\dots,K_{N}$ are drawn independently from a common distribution $p_{k}$. After assigning $K_{i}$ stubs of edges to each node $i$ ($i=1,\dots,N$), a pair of stubs is uniformly chosen at random and then connected to form an edge. This last step is repeated on the remainder stubs until there are no stubs left, and the outcome is a particular matching of the stubs with the prescribed random degrees $K_{1},\dots,K_{N}$. We do not allow for the existence of self-edges and multi-edges in the random network. We set $A_{ij}=J_{ij}/\sqrt{c}$ if there is an edge connecting nodes $i$ and $j$, and $A_{ij}=0$ otherwise. The coupling strengths $\{J_{ij}\}_{i,j=1,\dots,N}$ are i.i.d. random variables drawn from a common distribution $P_{J}$. We refer to [1] for other properties and subtleties of the configuration model. In the limit $N\rightarrow\infty$, the ensemble of adjacency random matrices $\bm{A}$ constructed from this procedure is specified by the degree distribution $p_{k}$ and the distribution of weights $P_{J}$. This is probably the simplest network model that allows to clearly exploit the influence of degree fluctuations on the spectral properties of $\bm{A}$.

We will present results for regular, Poisson, exponential, and Borel degree distributions. The analytic expression for $p_{k}$ in each case, together with the variance $\sigma^{2}$, is displayed in table 1. The properties of the configuration model with Poisson and exponential degree distributions have been extensively discussed in [28, 1]. The degree distributions in table 1 obey the following properties: (i) $p_{k}$ is parametrized solely in terms of the mean degree $c$; (ii) all moments of $p_{k}$ are finite for $c<\infty$; (iii) for sufficiently large $c$, the variances obey $\sigma_{\rm reg}^{2}<\sigma_{\rm pois}^{2}<\sigma_{\rm exp}^{2}<\sigma_{\rm bor}^{2}$. These four examples of degree distributions are very convenient to study, in a controllable way, the effect of degree fluctuations on random networks as we increase $c$. Indeed, we can compare for $c\gg 1$ the spectral density of networks with the same average degree, but with increasing variances $\sigma^{2}$. We do not consider here random networks with power-law degree distributions [32, 33], since in this case higher-order moments of $p_{k}$ diverge already for finite $c$.

Although the Borel degree distribution is not commonly employed in the study of random networks, the rate of its exponential decay is smaller than the one of the exponential degree distribution with the same $c$, which makes the Borel distribution well-suited to our analysis. The Borel distribution, introduced in the context of queuing theory [46, 47], also appears as the distribution of the total progeny in a Galton-Watson branching process with Poisson distributed degrees [48]. Figure 1 illustrates the Borel degree distribution for different average degrees. For fixed $k\ll c$, the Borel distribution $p_{k}$ attains a finite limit as $c\rightarrow\infty$, in contrast to the Poisson and the exponential degree distribution. For $k=O(c)$, the Borel distribution $p_{k}$ is proportional to $1/\sqrt{c}$ for $c\gg 1$. Therefore, Borel networks with large $c$ contain a finite fraction of nodes with $c$-independent degrees and a smaller fraction of nodes with degrees proportional to $c$. Note also that $p_{k}$ in figure 1 is highly skewed, its mode is located at $k=1$, and the tail of $p_{k}$ decays as $\ln p_{k}\sim\left(\frac{1}{c}+\ln\left(\frac{c-1}{c}\right)\right)k$ for $k\gg 1$. All network models considered here have exponentially decaying degree distributions.

The configuration model of networks has a key property: in the limit $N\rightarrow\infty$, the set of nodes in the neighborhood of a given node, drawn uniformly from the network, is arranged in a tree-like structure [42]. Nevertheless, the topology of random networks is fundamentally different from a Cayley tree [49], in the sense that boundary nodes are absent from the configuration model and cycles do survive in the limit $N\rightarrow\infty$, but their average length scales as $\ln N$ for $N\gg 1$. Roughly speaking, a random network can be seen as a Cayley tree with fluctuating degrees, wrapped onto itself.

In the limit $N\rightarrow\infty$, the spectral density is obtained from

where $d^{2}g\equiv d{\rm Re}g\,d{\rm Im}g$, and $\mathcal{P}(g)$ is the joint distribution of the real and imaginary parts of $G_{ii}(z)$. The symbol $\int_{{\rm Im}g\geq 0}$ represents an integral over $g\in\mathbb{C}$ restricted to the upper half complex plane. The distribution of the resolvent follows from the solution of the coupled equations [36, 37]

where $\langle\dots\rangle_{J_{1},\dots,J_{L}}$ denotes the average over the independent interaction strengths $J_{1},\dots,J_{L}$, distributed according to $P_{J}$. The quantity $\mathcal{P}_{\rm cav}(g)$ is the joint distribution of the real and imaginary parts of the resolvent elements $G_{ii}^{(j)}(z)$ on the cavity graph [36, 50], defined as the graph where an arbitrary node $j$ and all its edges have been removed. Equations (7) and (8) have been derived using both the cavity [36, 50] and the replica method [37] of disordered systems, based on the locally tree-like structure of the configuration model. Equations (6-8) have been rigorously proven in [42] and they are exact in the limit $N\rightarrow\infty$, constituting an interesting point of departure to study the dense limit $c\rightarrow\infty$ of the spectral density $\rho(\lambda)$.

Let $F(G)$ be an arbitrary function of the complex random variable $G$ distributed as $\mathcal{P}_{\rm cav}(g)$. By defining the average of $F(G)$

we can use Eq. (8) and rewrite the above expression as

where we introduced the distribution of $Q^{\prime}\equiv\frac{1}{c}\sum_{\ell=1}^{K-1}J_{\ell}^{2}G_{\ell}$

with the random variable $K$ denoting the degree of an uniformly chosen node. Following an analogous procedure, we can also rewrite the spectral density in terms of the distribution $\mathcal{W}(q)$ of the random variable $Q\equiv\frac{1}{c}\sum_{\ell=1}^{K}J_{\ell}^{2}G_{\ell}$

defined as

We note that $\mathcal{W}(q)$ and $\mathcal{W}_{\rm cav}(q)$ are distributions of sums of independent complex random variables containing a random and large number of terms. For the examples of $p_{k}$ in table 1, one can show that, in the limit $c\rightarrow\infty$, the number of summands in $Q$ and $Q^{\prime}$ diverges. Thus, instead of working with the distributions of the resolvent, it is more convenient to extract the $c\rightarrow\infty$ limit of $\mathcal{W}(q)$ and $\mathcal{W}_{\rm cav}(q)$. Let us introduce the characteristic functions $\varphi(p,t)$ and $\varphi_{\rm cav}(p,t)$ of, respectively, $\mathcal{W}(q)$ and $\mathcal{W}_{\rm cav}(q)$

with

In order to study the dense limit $c\rightarrow\infty$, we expand $S_{c}(p,t)$ in powers of $1/c$, keeping in mind that the moments of $\mathcal{P}_{\rm cav}(g)$ depend on $c$. The leading term

should yield the dense limit $c\rightarrow\infty$ of the spectral density. Note that we assumed $\langle G\rangle$ converges to a well-defined limit $\langle G\rangle_{\infty}$ when $c\rightarrow\infty$. In order to proceed further, we specify the degree distribution $p_{k}$ in Eqs. (13) and (14).

In this section we state the main theorem of reference [44] and explain how the resolvent equations for $c\rightarrow\infty$, obtained in the previous section, are recovered from this theorem. We use the notation of the previous sections and we discuss the theorem in a less technical language, which helps to establish the link between [44] and the results presented in this paper in a more straightforward way.

The theorem concerns the spectral density $\rho(\lambda)$ of the adjacency matrix $\bm{A}$ (see Eq. (2)) of the configuration model in the dense limit. Let $\nu(\kappa)$ be the probability density of the rescaled degrees $\frac{K_{i}}{c}$ ($i=1,\dots,N$) in the limit $c\rightarrow\infty$

The main result of [44] is the following theorem:

Let $\{K_{i}\}_{i=1}^{N}$ be i.i.d random variables forming a degree sequence of the configuration model such that $c=o(N)$ and $J_{ij}=1$. In the limit $N\rightarrow\infty$, if the distribution of $K_{i}/c$ converges to a limit distribution $\nu$ with a finite second moment, then $\rho(\lambda)=\lim_{N\rightarrow\infty}\rho_{N}(\lambda)$ is given by $\nu\boxtimes\rho_{\rm w}$, where $\rho_{\rm w}$ is the Wigner semicircle law.

Thus, the spectral density $\rho(\lambda)$ of the configuration model with $c=o(N)$ is given by the free multiplicative convolution $\boxtimes$ between the probability density $\nu$ of rescaled degrees and the Wigner semicircle law $\rho_{\rm w}$. The free multiplicative convolution $\nu\boxtimes\rho_{\rm w}$ between the probability measure associated to $\rho_{\rm w}$ and the probability measure corresponding to $\nu$ is the unique probability measure associated to $\rho$ such that its $S$-transform $S_{\rho}(\omega)$ ($\omega\in\mathbb{C}$) fulfills the convolution property $S_{\rho}(\omega)=S_{\nu}(\omega)S_{\rm\rho_{w}}(\omega)$, where $S_{\nu}(\omega)$ and $S_{\rm\rho_{w}}(\omega)$ are the $S$-transforms of $\nu$ and $\rho_{\rm w}$, respectively.

The free multiplicative convolution and the $S$-transform were introduced by Voiculescu to deal with the multiplication of free non-commuting random variables [58]. The $S$-transform finds applications in random matrix theory, since it is an important analytic tool to determine the spectral properties of products of large and independent random matrices from the spectra of the individual matrices [59, 60].

For unbounded measures on $\mathbb{R}_{+}$ and for measures with zero mean and all moments finite, such as the Wigner semicircle distribution, the $S$-transform $S_{\varphi}(z)$ ($z\in\mathbb{C}$) of the corresponding probability density $\varphi$ is obtained from [61, 62]

where

The function $\omega_{\varphi}^{-1}(z)$ is the functional inverse of $\omega_{\varphi}(z)$ with respect to composition and $\mathcal{C}_{\varphi}(z)$ is the Cauchy transform of $\varphi$

where $z$ lies outside the support of $\varphi(\lambda)$. Instead of working directly with the inverse $\omega_{\varphi}^{-1}(z)$ in Eq. (33), it is more convenient to express the $S$-transform in terms of the Cauchy transform [63, 60]. Since $\omega_{\varphi}^{-1}\left(\omega_{\varphi}(z)\right)=z$, we can rewrite Eq. (33) as follows

Hence, given the $S$-transform of a probability density $\varphi$, we are able to uniquely determine $\varphi$.

In the context of random matrix theory, $\mathcal{C}_{\rho}(z)$ is the trace of the resolvent matrix

which, after setting $z=\lambda-i\eta$, gives access to the spectral density $\rho$ through Eq. (5). Combining Eqs. (5), (34), and (37), we can also express the spectral density in terms of $\omega_{\rho}(z)$

The above theorem holds for networks where $c$ scales slowly with $N$ such that $\lim_{N\rightarrow\infty}\frac{c}{N}=0$. In the previous section, we have derived results by performing the limit $c\rightarrow\infty$ after $N\rightarrow\infty$, which also implies that $\lim_{N\rightarrow\infty}\frac{c}{N}=0$. Consequently, with the exception of dense Borel networks for which $\nu(\kappa)$ has a divergent second moment (see table 1), the results for $\rho(\lambda)$ derived in the previous section should be recovered from the above theorem.