Spectra of random graphs with community structure and arbitrary degrees

In this section we calculate the spectral density $\rho(z)$ of the centered matrix $\mathbf{X}$, Eq. (4). The spectral density is defined by

where $\lambda_{i}$ is the $i$th eigenvalue of $\mathbf{X}$ and $\delta(z)$ is the Dirac delta. The starting point for our calculation is the well-known Stieltjes–Perron formula, which gives the spectral density directly in terms of the matrix as

where $z-\mathbf{X}$ is shorthand for $z\mathbf{I}-\mathbf{X}$ with $\mathbf{I}$ being the identity.

To calculate the trace, we follow the approach of Bai and Silverstein [24], making use of the result that the $i$th diagonal component of the inverse of a symmetric matrix $\mathbf{B}$ is [14]

where $B_{ii}$ is the $i$th diagonal element of $\mathbf{B}$, $\mathbf{b}_{i}$ is the $i$th column of the matrix, and $\mathbf{B}_{i}$ is the matrix with the $i$th row and column removed. In the limit of large system size, and provided that the degrees of vertices become large as the network does, the distribution of values of $[\mathbf{B}^{-1}]_{ii}$ becomes narrowly peaked about its mean, and one can write the mean value as

If, as in our case, the elements of $\mathbf{B}$ are independent random variables with mean zero, then

where we have made use of $\langle[\mathbf{b}_{i}]_{j}[\mathbf{b}_{i}]_{k}\rangle=\langle[\mathbf{b}_{i}]_{j}\rangle\langle[\mathbf{b}_{i}]_{k}\rangle=0$ when $j\neq k$.

In our particular example we have $\mathbf{B}=z-\mathbf{X}$, which means that

(since $i\neq j$ by definition, the $i$th row having been removed from the matrix), so

where the last equality applies in the limit of large system size (for which it makes a vanishing difference whether we drop the $i$th row and column from the matrix or not, so $\mathbf{B}_{i}$ can be replaced with $z-\mathbf{X}$ for all $i$). Then, noting that $\langle B_{ii}\rangle=z-\langle X_{ii}\rangle=z$, Eq. (10) becomes

Summing this expression over $i$ we then get the trace we were looking for, which we will write in terms of a new function

where we have for convenience defined the vector function

The quantity $g(z)$ (which is just the trace divided by $n$) is called the Stieltjes transform of the matrix $\mathbf{X}$, and it will play a substantial role in the remainder of our calculation.

It remains to calculate the function $\mathbf{h}(z)$, which is now straightforward. Multiplying Eq. (14) by $\mathbf{k}_{i}$ and substituting into (16), we get

The solution for the spectral density involves solving this equation for $\mathbf{h}(z)$, then substituting the answer into Eq. (15) to get the Stieltjes transform $g(z)$. Then the spectral density itself can be calculated from Eq. (8):

Alternatively, we can simplify the calculation somewhat by rewriting Eq. (14) as

then summing over $i$ and dividing by $n$ to get $zg(z)-c\|\mathbf{h}(z)\|^{2}=1$, or

where $c=2m/n$ as previously, which is the average degree of the network, and $\|\mathbf{h}(z)\|$ denotes the vector magnitude of $\mathbf{h}(z)$, i.e., $\mathbf{h}\cdot\mathbf{h}$ (not the complex absolute value). Then the spectral density itself, from Eq. (18), is

If we further suppose that the parameter vectors $\mathbf{k}_{i}$ are drawn independently from some probability distribution $p(\mathbf{k})$, which plays roughly the role played by the degree distribution in other network models, then in the limit of large network size Eq. (17) can be written as

Equations (21) and (22) between them give us our solution for the spectral density. These equations can be regarded as generalizations of the equations for the configuration model given in Ref. [14] and similar equations have also appeared in applications of random matrix methods to other problems [25, 26, 27, 28, 29].

As an example of the methods of the previous section, consider a network of $n$ vertices with two communities of $\tfrac{1}{2}n$ vertices each. Let the first group consist of vertices $1\ldots\tfrac{1}{2}n$ and the second of vertices $\tfrac{1}{2}n+1\ldots n$. Vertices in the first group will have parameter vector $\mathbf{k}_{i}=(\kappa_{i},\theta)$ and those in the second group will have $\mathbf{k}_{i}=(\kappa_{i-n/2},-\theta)$, where the quantities $\kappa_{i}$ and $\theta$ are positive constants that we choose and $\kappa_{i}\geq\theta$ for all $i$, to ensure that the expected values $\langle A_{ij}\rangle=\mathbf{k}_{i}\!\cdot\!\mathbf{k}_{j}/2m$ of the adjacency matrix elements are non-negative.

This particular parametrization is attractive for a number of reasons. First, it already takes the form of the rank-2 eigenvector decomposition of Eq. (5), which simplifies the our calculations—the two (unnormalized) eigenvectors are the $n$-element vectors $(\boldsymbol{\kappa},\boldsymbol{\kappa})$ and $(1,1,\ldots,-1,-1,\dots)$ where $\boldsymbol{\kappa}$ is the $(\tfrac{1}{2}n)$-element vector with elements $\kappa_{1},\ldots,\kappa_{n/2}$. Also the expected degrees take a particularly simple form. The expected degree of vertex $i$ for $i\leq\tfrac{1}{2}n$ is

But, applying Eq. (1), we have $m=\sum_{j=1}^{n/2}\kappa_{j}$ and hence the expected degree of vertex $i$ is simply $\kappa_{i}$. By a similar calculation it can easily be shown that for $i>\tfrac{1}{2}n$ the expected degree is $\kappa_{i-n/2}$, and the average degree in the whole network is

The parameter $\theta$ also has a simple interpretation in this model: it controls the strength of the community structure. For instance, when $\theta=0$ vertices in the two communities are equivalent and there is no community structure at all.

To calculate the spectrum for this model, we substitute the values of $\mathbf{k}_{i}$ into Eq. (22) to get equations for the two components of the vector function $\mathbf{h}(z)$ thus:

where $p(\kappa)$ is the probability distribution of the quantities $\kappa_{i}$. Equation (26) has the trivial solution $h_{2}(z)=0$, so the two equations simplify to a single one:

and then

which is independent of the parameter $\theta$. These results are identical to those for the corresponding quantities in the ordinary configuration model with no community structure and expected degree distribution $p(\kappa)$, as derived in Ref. [14], and hence we expect the spectrum of the centered adjacency matrix to be the same for the current model as it is for the configuration model with the same distribution of expected degrees.

To give a simple example application, suppose that there are only two different values of $\kappa$. Half the vertices in each community have a value $\kappa_{1}$ and the other half $\kappa_{2}$. Then $p(\kappa)=\tfrac{1}{2}[\delta(\kappa-\kappa_{1})+\delta(\kappa-\kappa_{2})]$, where $\delta(x)$ is the Dirac delta, and $c=\tfrac{1}{2}(\kappa_{1}+\kappa_{2})$. With this choice

which can be rearranged to give the cubic equation:

which can be solved exactly for $h_{1}(z)$ and hence we can derive an exact expression for the spectral density. The expression itself is cumbersome (like the solutions of most cubic equations), but Fig. 1 shows an example for the choice $\kappa_{1}=60$, $\kappa_{2}=120$, along with numerical results for the spectrum of a single random realization of the model. As the figure shows, the two agree well. (The histogram in the left-hand part of the figure represents the spectrum of the centered matrix. The two outlying eigenvalues that appear to the right belong to the full, non-centered adjacency matrix and are calculated in the following section.)

Note also that in the special case where $\kappa_{1}=\kappa_{2}=c$, so that $\kappa$ is constant over all vertices, Eq. (29) simplifies further to

which is a quadratic equation with solutions

and hence the spectral density is

where we take the negative square root in Eq. (32) to get a positive density. Equation (33) has the form of the classic semicircle distribution for random matrices. This model is equivalent to the standard stochastic block model and (33) agrees with the expression for the spectral density derived for that model by other means in Ref. [13].

So far we have derived the spectral density of the centered adjacency matrix $\mathbf{X}=\mathbf{A}-\langle\mathbf{A}\rangle$. We can use the results of these calculations to compute the spectrum of the full adjacency matrix by generalizing the method used in [14], as follows.

Using Eq. (5) we can write the adjacency matrix as

Let us first consider the effect of adding just one of the terms in the sum to the centered matrix $\mathbf{X}$, calculating the spectrum of the matrix $\mathbf{X}+\alpha_{1}\mathbf{u}_{1}\mathbf{u}_{1}^{T}$. Let $\mathbf{v}$ be an eigenvector of this matrix with eigenvalue $z$:

Rearranging this equation we have $\alpha_{1}\mathbf{u}_{1}\mathbf{u}_{1}^{T}\mathbf{v}=(z-\mathbf{X})\mathbf{v}$ and, multiplying by $\mathbf{u}_{1}^{T}(z-\mathbf{X})^{-1}$, we find

Note that the vector $\mathbf{v}$ has canceled out of the equation, leaving us with an equation in $z$ alone. The solutions for $z$ of this equation give us the eigenvalues of the matrix $\mathbf{X}+\alpha_{1}\mathbf{u}_{1}\mathbf{u}_{1}^{T}$.

Expanding the vector $\mathbf{u}_{1}$ as a linear combination of the eigenvectors $\mathbf{x}_{i}$ of the matrix $\mathbf{X}$, the equation can also be written in the form

where $\lambda_{i}$ are the eigenvalues of $\mathbf{X}$. Figure 2 shows a graphical representation of the solution of this equation for the eigenvalues $z$. The left-hand side of the equation, represented by the solid curves, has simple poles at $z=\lambda_{i}$ for all $i$. The right-hand side, represented by the horizontal dashed line, is constant. Where the two intercept, represented by the dots, are the solutions for $z$. From the geometry of the figure we can see that the values of $z$ must fall between consecutive values of $\lambda_{i}$—we say that the $z$’s and $\lambda$’s are interlaced. If we number the eigenvalues $\lambda_{i}$ in order from largest to smallest so that $\lambda_{1}\geq\lambda_{2}\geq\ldots\geq\lambda_{n}$, and similarly for the $n$ solutions $z_{i}$ to Eq. (37), then $z_{1}\geq\lambda_{1}\geq z_{2}\geq\lambda_{2}\geq\ldots\geq z_{n}\geq\lambda_{n}$. In the limit of large system size, as the $\lambda_{i}$ become more and more closely spaced in the spectrum of the matrix, this interlacing places tighter and tighter bounds on the values of $z_{i}$, and asymptotically we have $z_{i}=\lambda_{i}$ and the spectral density of $\mathbf{X}+\alpha_{1}\mathbf{u}_{1}\mathbf{u}_{1}^{T}$ is the same as that of $\mathbf{X}$ alone.

There is one exception, however, in the highest-lying eigenvalue $z_{1}$, which is bounded below by $\lambda_{1}$ but unbounded above, meaning it need not be equal to $\lambda_{1}$ and may lie outside the band of values occupied by the spectrum of the matrix $\mathbf{X}$. To calculate this eigenvalue we observe that, the matrix $\mathbf{X}$ being random, its eigenvectors $\mathbf{x}_{i}$ are also random and hence $\mathbf{x}_{i}^{T}\mathbf{u}_{1}$ is a zero-mean random variable with variance $1/n$. Taking the average of Eq. (37) over the ensemble of networks, the numerator on the left-hand side gives simply a factor of $1/n$ and we have

The solution to this equation gives us the value of $z_{1}$.

This then gives us the complete spectrum for the matrix $\mathbf{X}+\alpha_{1}\mathbf{u}_{1}\mathbf{u}_{1}^{T}$. It consists of a continuous spectral band with spectral density equal to that of the matrix $\mathbf{X}$ alone, which is calculated from Eq. (21), plus a single eigenvalue outside the band whose value is the solution for $z$ of $g(z)=1/\alpha_{1}$.

We could have made the same argument about any single term $\alpha_{r}\mathbf{u}_{r}\mathbf{u}_{r}^{T}$ appearing in Eq. (34) and derived the corresponding result that the continuous spectral band is unchanged from the centered matrix but there can be an outlying eigenvalue $z_{r}$ given by

The calculation of the spectrum of the full adjacency matrix requires that we consider all terms in Eq. (34) simultaneously, but in practice it turns out that it is enough to consider them one by one using Eq. (39). The argument for this is in two parts as follows.

1. 1.

   We have shown that the spectral density of the continuous band in the spectrum of the matrix $\mathbf{X}+\alpha_{1}\mathbf{u}_{1}\mathbf{u}_{1}^{T}$ is the same as that for the matrix $\mathbf{X}$ alone, and there is one additional outlying eigenvalue, which we denote $z_{1}$. Now we can add another term $\alpha_{2}\mathbf{u}_{2}\mathbf{u}_{2}^{T}$ and repeat our argument for the matrix $\mathbf{X}+\alpha_{1}\mathbf{u}_{1}\mathbf{u}_{1}^{T}+\alpha_{2}\mathbf{u}_{2}\mathbf{u}_{2}^{T}$, finding the equivalent of Eq. (37) to be

   $$\sum_{i=2}^{n}{(\mathbf{x}_{i}^{T}\mathbf{u}_{2})^{2}\over z^{\prime}-z_{i}}+{(\mathbf{x}_{1}^{T}\mathbf{u}_{2})^{2}\over z^{\prime}-z_{1}}={1\over\alpha_{2}},$$

   (40)

   where $z^{\prime}$ is the eigenvalue of the new matrix and $z_{i}$ are the solutions of (37). As before, this implies there is an interlacing condition and that the spectral density of the perturbed matrix is the same within the spectral band as that for the unperturbed matrix. We can repeat this argument as often as we like and thus demonstrate that the shape of the spectral band never changes, so long as the number of perturbations (which is also the rank of $\langle\mathbf{A}\rangle$) is small compared to the size of the network, i.e, $q\ll n$.

2. 2.

   This argument pins down all but the top two eigenvalues of $\mathbf{X}+\alpha_{1}\mathbf{u}_{1}\mathbf{u}_{1}^{T}+\alpha_{2}\mathbf{u}_{2}\mathbf{u}_{2}^{T}$. These two we can calculate by a variant of our previous argument. We average Eq. (40) over the ensemble, noting again that $\langle(\mathbf{x}_{i}^{T}\mathbf{u}_{2})^{2}\rangle=1/n$ and find that

   $${1\over n}\sum_{i=2}^{n}{1\over z^{\prime}-z_{i}}+{1/n\over z^{\prime}-z_{1}}={1\over\alpha_{2}}.$$

   (41)

   For large $n$ the first sum is once again equal to the Stieltjes transform $g(z)$ and hence the top two eigenvalues are solutions for $z^{\prime}$ of

   $$g(z^{\prime})+{1/n\over z^{\prime}-z_{1}}={1\over\alpha_{2}}.$$

   (42)

   But $g(z)$ and $\alpha_{2}$ are of order 1, while the term $n^{-1}/(z^{\prime}-z_{1})$ is of order $1/n$ and hence can in most circumstances be neglected, giving $g(z^{\prime})=1/\alpha_{2}$, which recovers Eq. (39). The only time this term cannot be neglected is when $z^{\prime}$ is within a distance of order $1/n$ from $z_{1}$, in which case we have a simple pole in the left-hand side of the equation as $z^{\prime}$ approaches $z_{1}$. Thus the left-hand side has the form sketched in Fig. 3, following $g(z)$ closely for most values of $z$, but diverging suddenly when very close to $z_{1}$. Equation (42) then has two solutions, as indicated by the dots in the figure, one given by $g(z)=1/\alpha_{2}$ and one that is asymptotically equal to $z_{1}$, which is the solution of $g(z)=1/\alpha_{1}$.

We can repeat this argument as many times as we like to demonstrate that the outlying eigenvalues are just the $q$ solutions of Eq. (39) for each value $r=1\ldots q$. Thus our final solution for the complete spectrum of the adjacency matrix has two parts: a continuous spectral band, given by Eqs. (21) and (22), and $q$ outlying eigenvalues, given by the solutions of Eq. (39), with $g(z)$ given by Eq. (20).

Let us return to the examples of Section III.2 and apply the methods above to the calculation of their outlying eigenvalues. Recall that we looked at networks with two communities and chose parameter vectors $\mathbf{k}_{i}=(\kappa_{i},\theta)$ for vertices in the first community and $\mathbf{k}_{i}=(\kappa_{i-n/2},-\theta)$ for those in the second. For such networks the vector function $\mathbf{h}(z)$ reduces to a single scalar function $h_{1}(z)$ that satisfies Eq. (27). At the same time, Eq. (20) tells us that for this model $zg(z)=1+ch_{1}^{2}(z)$ and hence from Eq. (39) the positions of the outlying eigenvalues are solutions of

for $r=2\ldots q$. Locating the outliers is thus a matter of solving (27) for $h_{1}$, substituting the result into (43), and then solving for $z$.

Consider, for instance, the choice we made in Section III.2, where there were just two values of $\kappa$, denoted $\kappa_{1}$ and $\kappa_{2}$, with half the vertices in each community taking each value. Then $h_{1}$ obeys the cubic equation (30), which can be solved exactly, and hence we can calculate the position of the outliers. Figure 1 shows the results for the choice $\kappa_{1}=60$, $\kappa_{2}=120$, $\theta=50$, along with numerical results for the same parameter values. As the figure shows, analytic and numerical calculations again agree well—so well, in fact, that the difference between them is quite difficult to make out on the plot.

We also looked in Section III.2 at the simple case where $\kappa=c$ for all vertices, so that they all have the same expected degree, in which case the model becomes equivalent to the standard stochastic block model and the continuous spectral band takes the classic semicircle form of Eq. (33). For this model we have $\alpha_{1}=c$ and $\alpha_{2}=\theta^{2}/c$. Using Eq. (32) for $h_{1}(z)$ and solving (43) for $z$, we then find the top two eigenvalues of the adjacency matrix to be

which agrees with the results given previously for the stochastic block model in Ref. [13].

One of the primary uses of network spectra is for the detection of community structure [5, 13]. As we have seen, the number of eigenvalues above the edge of the spectral band is equal to the number of communities in the network, and hence the observation of these eigenvalues can be taken as evidence of the presence of communities and their number as an empirical measure of the number of communities. The identity of the communities themselves—which vertices belong to which community—can be deduced, at least approximately, by looking at the elements of the eigenvectors [5].

However, as shown previously in [13] for the simplest two-community block model, the position of the leading eigenvalues varies as one varies the strength of community structure, and for sufficiently low (but still nonzero) strength an eigenvalue may meet the edge of the spectral band and hence become invisible in the spectrum, meaning it can no longer be used as evidence of the presence of community structure. Moreover, as also shown in [13], the elements of the corresponding eigenvector become uncorrelated with group membership at this point, so that any algorithm which identifies communities by examining the eigenvector elements will fail. The point where this happens, at least in the simple two-community model, coincides with the known “detectability threshold” for community structure, at which it is believed all algorithms for community detection must fail [17, 18, 19].

We expect qualitatively similar behavior in the present model as well. Consider the Stieltjes transform $g(z)$ defined in Eq. (15). Inside the spectral band the transform is complex by definition—see from Eq. (18). Above the band it is real and monotonically decreasing in $z$, as we can see by evaluating the trace in the basis in which $\mathbf{X}$ is diagonal:

where $\lambda_{i}$ are the eigenvalues of $\mathbf{X}$ as previously. Above the band, where $z>\lambda_{i}$ for all $i$, every term in this sum is monotonically decreasing, and hence so is $g(z)$. This implies via Eq. (39) that larger values of $\alpha_{r}$ give larger eigenvalues and that the largest real value $g_{\textrm{max}}$ of the Stieltjes transform occurs exactly at the band edge. Moreover, as shown in Ref. [14], the edge of the band is marked generically by a square-root singularity in the spectral density, which implies that $g_{\textrm{max}}$ is finite—see Fig. 3 for a sketch of the function. Thus when we make the community structure in the network weaker, meaning we decrease the values of the $\alpha_{r}$, we also decrease the outlying eigenvalues of the adjacency matrix and eventually the lowest of those eigenvalues will meet the edge of the band and disappear at the point where $1/\alpha_{r}=g_{\textrm{max}}$. If we continue to weaken the structure, more eigenvalues will disappear, in order—smallest first, then second smallest, and so forth.

Thus we expect there to be a succession of detectability transitions in the network, $q-1$ of them in all, where $q$ again is the number of communities. At the first of these transitions the $q$th largest eigenvalue will meet the band edge and disappear, meaning there will only be $q-1$ outlying eigenvalues left and hence there will be observational evidence of only $q-1$ communities in the network, even if in fact we know there to be $q$. At the next transition the number will decrease further to $q-2$, and so forth. One thus loses the ability to detect community structure in stages, one community at a time. Final evidence of any structure at all disappears at the point where the second largest eigenvalue meets the band edge.

Consider, for instance, the example network from Section III.2 again, in which there are two groups with parameter vectors of the form $(\kappa_{i},\pm\theta)$, where the parameters $\kappa_{i}$ control the expected degrees and $\theta$ controls the strength of the community structure. As before, let us study the case where the $\kappa_{i}$ take just two different values with equal probability, so that $h_{1}$ satisfies the cubic equation (30) (and $h_{2}=0$). Then we can calculate the maximal real value of $g(z)$ as follows.

Like $g(z)$, the function $h_{1}(z)$ is real outside the continuous spectral band but complex inside it, as one can see from Eq. (28). The band edge is thus the point at which the solution of the cubic equation becomes complex, which is given by the zero of the discriminant of the cubic. Take, for example, the case where $\kappa_{1}=\kappa$ and $\kappa_{2}=2\kappa$ for some constant $\kappa$. Then, employing the standard formula, the discriminant of (30) is

This is zero when, and hence the band edge falls at, $z=\sqrt{x\kappa}$, where $x\simeq 7.058$ is the sole real solution of the cubic equation $27x^{3}-216x^{2}+252x-512=0$. Substituting into Eq. (30), we then find that the value of $h_{1}$ at the band edge is $y/\sqrt{\kappa}$ where $y=0.723$ is the smallest real solution of the cubic equation $2y^{3}-3\sqrt{x}y^{2}+(x+\tfrac{4}{3})y-\sqrt{x}=0$. Then, using Eq. (20) and the fact that the average degree is $c=\frac{3}{2}\kappa$, the value of $g(z)$ at the band edge is

In this case there is only one parameter $\alpha_{r}$ with $r\geq 2$, which is $\alpha_{2}=\theta^{2}/c$. Hence there is a single threshold at which we lose the ability to detect communities, falling at

or

If $\theta$ is smaller than this value then spectral methods will fail to detect the communities in the network. We have checked this behavior numerically and find indeed that spectral community detection fails at approximately this point.