Localization of Laplacian eigenvectors on random networks

In large random networks, each eigenvector of the Laplacian matrix tends to localize on a subset of network nodes having similar numbers of edges, namely, the components of each Laplacian eigenvector take relatively large values only on a particular subset of nodes whose degrees are close. Although this localization property has significant consequences for dynamical processes on random networks, a clear theoretical explanation has not yet been established. Here we analyze the origin of localization of Laplacian eigenvectors on random networks by using a perturbation theory. We clarify how heterogeneity in the node degrees leads to the eigenvector localization and that there exists a clear degree-eigenvalue correspondence, that is, the characteristic degrees of the localized nodes essentially determine the eigenvalues. We show that this theory can account for the localization properties of Laplacian eigenvectors on several classes of random networks, and argue that this localization should occur generally in networks with degree heterogeneity.


Results
Laplacian matrix and its eigenvectors. We consider a network consisting of N nodes. The network topology is specified by a N × N adjacency matrix A, whose element A ij takes a value of 1 if there is an edge between nodes i and j, and 0 otherwise ( =  i j N , 1, 2, , ). We assume that the network is connected (a path exists between arbitrary nodes), and that the connection is non-directed, i.e., = A A ij ji . In this study, we define the Laplacian matrix = L L { } ij of the network as is the degree of the ith node, i.e., the number of edges, and δ i j , is the Kronecker's delta symbol 18,19 . As we will show, it is convenient to sort the node indices {i} in decreasing order of the degree k i , so that inequalities ≥ ≥ ≥  k k k N 1 2 hold. We denote the average degree of the network as Diffusion processes on the network are described by the Laplacian matrix. Suppose that each network node is occupied by some substance X, and denote its concentration on node i as = X x [ ] i i . The change in the concentration by the diffusive transportation of X is described as , where the flux of the substance from node j to node i is proportional to the concentration difference x j − x i (Fick's law).
( ) and the eigenvalue Λ α ( ) of the Laplacian matrix L satisfy the eigenvalue equation is the index of the eigenvector. The eigenvectors can be orthonormalized as 1, 2, , , because L is a real symmetric matrix. The Laplacian matrix L is negative semidefinite, i.e., is satisfied for any vector → =  x x x ( , , ) N 1 6 . Therefore, all Laplacian eigenvalues are non-positive, and only one of them, which corresponds to the uniform eigenvector  N (1, , 1)/ , takes 0 because the network is connected. The eigenvector indices {α} are also sorted in increasing order of the Laplacian eigenvalues so that Λ ≤ Λ ≤ ≤ Λ =  0 Localization of Laplacian eigenvectors. First, let us illustrate the Laplacian eigenvectors for several classes of random networks. In Fig. 1, all eigenvectors (except the uniform eigenvector φ , which has exceptional characteristics and is excluded from the analysis) of the scale-free network generated by the Barabási-Albert preferential attachment algorithm (BA) 20 , classical Erdös-Rényi random network (ER) 21 , and real neural network of C. elegans (CE) 22,23 are displayed in the contour plot, where the horizontal axis is the node index and the vertical axis is the eigenvector index. We show the absolute value φ α i ( ) of the eigenvector components, because each component is statistically symmetric . For the BA network, three typical eigenvectors with different α are also shown for illustration. For the CE network, we focus only on the connectivity and symmetrize the original network, which consists of 277 neurons with directed connections, by defining the adjacency matrix as if there is a edge between nodes i and j. Remarkably, clear diagonal structures are observed in all figures. Because the nodes are sorted by their degrees, this means that, in each eigenvector, only the nodes sharing similar degrees take large vector components, while other nodes have very small components. Indeed, for each eigenvector shown in Fig. 1, the mean difference in the degrees of the localized nodes where φ > .
is 1.33 for the BA network, 2.01 for the ER network, or 2.10 for the CE network. Each of these numbers is much smaller than the entire range of the degree distribution in each network, k 1 − k N , which is 106 for the BA, 25 for the ER, or 75 for the CE.
Moreover, the visible diagonal structures indicate that the characteristic degree of each localized subset linearly correlates with the eigenvalue index, which is also sorted by their eigenvalues. Thus, clear degree-eigenvalue correlation exists in the Laplacian eigenvectors (See Fig. 2(c-f)). It is also notable that the patterns of the localization are qualitatively different among the networks. In the BA network [ Fig. 1(a)], the localization is stronger near the hubs, i.e., the nodes with large degrees (e.g., α = 50 in panel (b)), while comparatively weak at the peripheries, i.e., the nodes with small degrees (e.g., α = 350 in panel (b)). In contrast to the BA network, in the ER [ Fig. 1(c)] and CE networks [ Fig. 1(d)], the localization is stronger both at hubs and peripheries, and weaker at the intermediate nodes. In McGraw and Menzinger 13 , the level of localization has been quantified by using the inverse participation ratio, i.e., () 2 2 , a standard quantity used in the analysis of Anderson localization.

Perturbation analysis of the Laplacian matrix.
To analyze the origin of this intriguing localization property, we apply the perturbation theory 24,25 to the Laplacian matrix. A similar perturbation approach was used by Kim and Motter to analyze the Laplacian eigenvalues of scale-free networks 26 .
In the present problem, the Laplacian matrix L has two types of elements of distinct orders. The diagnoal elements δ −k i i i , are order 〈k〉, while the non-diagional elements, which take 0 or 1, are (1)  . By introducing an expansion parameter ε = − k 1 , we can rewrite the Laplacian matrix as L = L 0 + εL 1 , whose elements δ = − L k ij i i j 0, , and = L kA ij ij 1, are of the same order,  k ( ). When the network is sufficiently dense, i.e.,  k 1, the expansion parameter ε is small, and it is expected that the perturbation theory yields reasonable approximation of the Laplacian eigenvectors. For convenience, we employ the bra-ket notation to denote the Laplacian eigenvector, i.e., φ α → = α ( ) , and drop the summation symbol as . Expanding the Laplacian eigenvectors |α〉 and eigenvalues Λ α ( ) in series of ε as α α ε α ε α = + + + , and substituting into the eigenvalue equation (1), the following set of equations is obtained up to  ε ( ) 2 : Let us first consider the unperturbed system (2). One can easily find that the eigenvectors α 0 and eigenvalues Λ α 0 ( ) are given exactly as , and the corresponding eigenvalue Λ α 0 ( ) is simply equal to the negative of the characteristic node degree α k . Thus, strictly localized eigenvectors are obtained at the zeroth-order, where the network topology is completely ignored and the non-diagonal elements are assumed to be vanishingly small. Note that L 0 is no longer a Laplacian matrix , which always vanishes, the zeroth-order eigenvalue Λ N 0 ( ) is not zero but equal to −k N , i.e., the smallest degree of the network.
In order to analyze the localization property of the Laplacian eigenvectors, we should consider the higher-order perturbation terms and, in particular, the fact that networks generally possess multiple nodes with the same degrees. From the zeroth-order solution (5), this indicates that the zeroth-order eigenvectors are degenerate. Therefore, we need to employ the degenerate perturbation theory 24,25 . From Eqs (2-4), we can compute the approximate eigenvectors and eigenvalues by the first-and second-order degenerate perturbation theory, respectively. The complete derivation of the perturbation corrections for a general class of degenerate systems has been reported, e.g., in ref. 25 . See Methods and Supplementary Information for details. Accuracy of the perturbation approximation is also discussed in the Methods.
Approximate eigenvectors. We now apply the perturbation theory to the networks used in Fig. 1, and demonstrate that it can predict the localization properties of random networks. The results are shown in Figs 2-6. Figure 2(a-c) show the node degrees of the BA, ER, and CE networks as functions of the node index. Because the node indices are sorted in decreasing order of degrees, the curves monotonically decreases with the node index. Figure 2(d-f) show scatter plots of degree-eigenvalue pairs, 1, , ), for the three networks. We can observe that the data points approximately lie along the diagonal line in each figure, implying that the eigenvalues and node degrees are closely correlated in these networks. Such correlation between eigenvalues and degrees has also been reported in a preceding study 27 . Now, Fig. 2(g-i) show the zeroth-and second-order approximations of the Laplacian eigenvalues. We can see that the zeroth-order result already provides a good approximation to the true Laplacian eigenvalues for the BA and CE networks. For the ER network, the zeroth-order eigenvalues somewhat deviate from the true eigenvalues, but higher-order approximation gives closer values. Thus, the perturbation theory accounts for the Laplacian eigenvalues of these networks, which indicates that they are essentially determined by the node degrees. Figure 3 displays the approximated eigenvectors as a function of the eigenvector index α and the node index i, similarly to Fig. 1. As can be seen, the first-order predictions are in good qualitative agreement with the true Laplacian eigenvectors shown in Fig. 1. The diagonal structures indicating localization of eigenvectors are well reproduced for all networks. Moreover, the different patterns of localization among the networks are correctly reproduced. That is, the localization is stronger at hubs and weaker at peripheries in the BA network, while it is strong both at hubs and peripheries and weak at the intermediate nodes in the ER and CE networks.
The slightly broadened localization patterns along the diagonal line can be interpreted as follows. If some nodes in the network share the same degree, the corresponding zeroth-order eigenvectors are degenerate. The eigenvectors in the degenerate subspace are mutually mixed, yielding block-diagonal structures of various sizes in the contour plot. The size of each block-diagonal component is equal to the number of the degenerate eigenvectors, i.e., the number of nodes sharing the same degree. Therefore, the degree distribution of the network determines the pattern of localization.
Indeed, in the BA network, the degrees obey a scale-free distribution, where only a small number of nodes have large degrees (hubs) and the majority of the nodes have small degrees (peripheries) [ Fig. 3(d)]. Correspondingly, the degeneracy of the degrees, i.e., the zeroth-order eigenvalues, is small for the hubs and large for the peripheries. Therefore, the localization is stronger at the hubs than that at the peripheries because less eigenvectors are involved.
In contrast, the ER network has a binomial distribution of the degrees [ Fig. 3(e)]. The majority of the nodes belong to the intermediate degrees, and the hubs and peripheries are composed of relatively small numbers of nodes. This leads to stronger localization at both hubs and peripheries, and weaker localization at intermediate nodes, in contrast to the BA case. Similarly, in the CE network, the degrees obey a binomial-like distribution [ Fig. 3(f)], so the localization pattern is also similar to that of the ER network.
In Indeed, quantitative node-wise comparison of the true and approximate vectors yields considerable discrepancy. The fact that the essential localization property of the vectors is qualitatively reproduced in Fig. 3 suggests that the true and approximate eigenvectors share similar characteristics when they are averaged over degenerate nodes and eigenvectors. In order to evaluate the performance of approximation while excluding the effect of degeneracy, we construct reduced degree-wise vectors from the true node-wise Laplacian eigenvectors, where vector components of the original node-wise Laplacian vectors are averaged over degenerate nodes having the same degree and over the degenerate eigenmodes having the same zero-th order eigenvalues as follows: Here, the degree index k runs from the minimum degree to the maximum degree of the network nodes, the reduced eigenvalue index β runs from the minimum to maximum of the zero-th order eigenvalues, N k is the number of degenerate nodes with degree k, and β N is the number of degenerate eigenmodes with the zero-th order eigenvalue − β k , respectively. We then calculate the correlation coefficient σ between the true and approximate eigenvectors, defined by is the reduced degree-wise vector obtained from the approximate eigenvectors. For comparison, we also generate N independent random eigenvectors whose components are randomly drawn from a uniform distribution over [−1, 1], normalize the vectors so that their norms become equal to 1, and calculated their correlations similarly. We exclude the uniform eigenvector φ → =  N (1, 1, , 1)/ N ( ) from the analysis, which is exceptional and cannot be predicted by the perturbation theory. Figure 4(e-g) show the correlation coefficient σ between the true and approximate eigenvectors with respect to the reduced eigenvalue index β. As can be seen in the figures, the correlation coefficient σ takes large values for non-degenerate eigenmodes, which are much larger than the correlation coefficient for random vectors and thus indicate similarity between the true and approximate vectors. It can be seen that non-degenerate vectors show larger correlations than degenerate vectors. It can be also observed that σ is higher near hubs for the BA, and near hubs and peripheries for the ER and CE. For other eigenvectors, σ can be as small as those of random vectors, indicating that the perturbation theory does not predict some of the eigenvectors well.
We stress that, although the correlation coefficients can be small for some of the eigenvectors, that is, the perturbation approximation does not predict them quantitatively, essential qualitative properties of the true eigenvectors such as the localizing nodes and the degree of localization are still reproduced well. This is because such properties are mainly determined by the degree of the nodes in the same degenerate block, which share statistically similar connectivities to the rest of the network. See Methods for the discussion on the accuracy of the perturbation approximation.
Thus, the perturbation theory can account for the eigenvector localization and degree-eigenvalue correspondence reasonably well. It reveals how degree heterogeneity and degeneracy lead to the eigenvector localization and the degree-eigenvalue correspondence. Furthermore, it clarifies why the representation of the Laplacian eigenvectors in Fig. 1, with respect to the eigenvector index α and the node index i both sorted in decreasing order of the eigenvalues and degrees, yields the clearly visible localized structures.
Our perturbation analysis also explains how similar the degrees of the nodes should be in order that the Laplacian eigenvector localizes on these nodes. From Eqs (9-11), we observe that the difference in the zero-th order eigenvalues in the denominator, Λ − Λ , which is equal to the difference in the node degrees from Eq. (5), should be sufficiently smaller than the maximal range of the eigenvalues, Λ − Λ N 0 ( ) 0 (1) , which is equal to − k k max m in , in order to give a dominant contribution to the first-order correction to the eigenvector. (If this is not the case, all nodes in the network will give similar contributions to the perturbation correction and localization will not be observed). Thus, the degrees of the nodes should be similar in the sense that their difference is much smaller than the range of the entire degree distribution, − k k N 1 , as we mentioned in the introduction. Weighted and directed networks. Although we have so far presented the results only for non-directed and non-weighted networks, our analysis can straightforwardly be extended to directed and weighted networks. For directed networks, the adjacency matrix A is generally asymmetric. The weight of the edge from node j to node i is specified by the element W ij of the weight matrix W. The Laplacian matrix of such a network is defined as  . Thus, one can follow the perturbation analysis as described above with the generalized Laplacian matrix L.
As an illustrative example, Fig. 5 compares the true Laplacan eigenvectors and the result of the perturbation approximation of the real asymmetric neural network of C. elegans, which we used in Figs 1-4 after symmetrization. Note that the Laplacian matrix is now asymmetric and the elements of its eigenvectors can take complex values. We focus only on the localization pattern of the Laplacian eigenvectors and plot the absolute value of each vector component. As can be seen in the figure, the approximate eigenvectors can reproduce the localization pattern of the true eigenvectors qualitatively well.

Regular lattices.
Our argument suggests that the block-diagonal components representing the eigenvector localization can also be more or less observed even when the network is not random, but is formed by nodes with non-identical degrees. In order to verify this statement, we here consider a regular lattice network as shown in Fig. 6(a), which is composed of three types of nodes with different degrees. Specifically, one third of the nodes have degree k = 6, another one third have degree k = 4, and the rest have degree k = 2.
The zeroth-order unperturbed result is shown in Fig. 6(b). All eigenvectors degenerate into three classes at this stage, corresponding to the characteristic degrees k = 6, k = 4, and k = 2. Higher-order perturbations solve this degeneracy and mix the eigenvectors in each subset into three blocks, corresponding to k = 6, k = 4, and k = 2. Thus, at the first-order perturbation, the Laplacian eigenvectors show block-diagonal structures in the contour plot as shown in Fig. 6(c).
The above prediction is in good agreement with the true Laplacian eigenvectors obtained by direct numerical calculation, shown in Fig. 6(d). Thus, the degree heterogeneity generally leads to localized eigenvectors even in regular lattice networks. This result also suggests that the degree heterogeneity is the origin of the localization property of the Laplacian eigenvectors.

Discussion
The Laplacian eigenvectors of networks with degree heterogeneity generally exhibit localization on the subset of nodes with close degrees and the Laplacian eigenvalues show clear degree-eigenvalue correspondence. We have proposed a simple explanation for the localization property of Laplacian eigenvectors using the degenerate perturbation theory. It clarifies how degree heterogeneity and degeneracy lead to the eigenvector localization and degree-eigenvalue correspondence. We analyzed three kinds of random networks with different statistical properties, and confirmed that our approach can reasonably account for the true Laplacian eigenvectors. We have also shown that the analysis can straightforwardly be extended to directed and weighted networks.
Our results show that the node degrees in heterogeneous networks correspond to the wavenumbers in regular lattices. Therefore, the degree of the node can play an essential role as the "natural coordinate" in describing the dynamics or patterns on networks with heterogeneous degree distributions, because they determine the eigenvalue and the subset of nodes that participate in the eigenvector. We conjecture that this partly accounts for why various network dynamics, plotted with respect to node degrees, often exhibit ordered patterns and provide us with physical interpretations.
In this study, we used the C. elegans neuronal network to illustrate the generality of the eigenvector localization and did not analyze its particular functional structures. Though extraction of functional structures from the C. elegans neuronal network is beyond the scope of the present study, we can observe a possible sign of such structures from the eigenvalues of the C. elegans network shown in Fig. 2(f); that is, there exists a small cluster of eigenvalues (3 ≤ α ≤ 12) separated from other eigenvalues, and correspondingly a tiny block structure exists on top of the diagonal structure in the Laplacian eigenvectors in Fig. 1(d). Such detailed structures of the Laplacian eigenvectors could reflect some functional structure of the neuronal network of C. elegans. In analyzing such detailed structures, broad degree heterogeneity, which yields the diagonal localized structure, might be disturbing, and spectral clustering methods based on the normalized Laplacian matrix 28 , which removes the effect of degree heterogeneity, could provide more useful information of the network.
Finally, we note that the localization property of eigenvectors is not restricted to network Laplacian matrices. We have recently formulated the advection equation for random networks 29 and found that the eigenvectors of the advection matrix are also localized on a subset of nodes. This localization can also be accounted for by a similar perturbation analysis of the network, and the homogenization process on the network due to advection can be clearly visualized when plotted with respect to the node degrees. Further investigation of eigenvector localization on networks will provide us with insights into complex dynamics on networks.

Methods
Degenerate perturbation theory. Following a standard argument from quantum mechanics, we classify each eigenvector into the following three types according to its degeneracy: (A) non-degenerate, (B) degeneration that is solved at the first order, and (C) otherwise. We compute the approximate eigenvectors and eigenvalues by the first-and second-order perturbation theory, respectively. For each case, the perturbation corrections are given as follows (see Supplementary Information for the derivation).
For type (A), the first order correction is The first order correction to the eigenvalue vanishes, i.e., 1, , ). Namely, they satisfy the secular eigenvalue equation gives the first-order correction to the Laplacian eigenvalue. The first-order correction to the eigenvector is where the summation symbol with β α ≠ indicates that the index β runs over all eigenvectors except for the degenerate ones, i.e., α α  , , m 1 . The second order correction to the eigenvalue is given by For type (C), suppose that a subset of the eigenvalues α α , , n 1 ( ≤ n m) is still degenerate at the first order perturbation. In this case, the zeroth-order eigenvector is further redefined as α α where c i j , are given by the eigenvectors of the matrix W defined as ( ) . The first-order correction to the eigenvector is given by i . If some eigenvectors are still degenerate even at the second-order perturbation, one can further introduce new zeroth-order eigenvectors so that the degeneracy is solved at the higher-order perturbation. However, for simplicity, we do not consider the higher-order perturbations. Therefore, the eigenvectors are not completely determined in this case.
Accuracy of the perturbation approximation. As is well known, it is generally difficult to prove the convergence of the perturbation series. Accuracy of the perturbation approximation may roughly be assessed by looking at the ratio β α Λ − Λ α β L / 0 1 0 0 ( ) 0 ( ) for each pair of non-degenerate eigenmodes α and β. If this ratio is sufficiently small, the contribution of the zeroth-order eigenvector β 0 to the first-order correction α 1 will be accurately evaluated.
In the present case, β α = βα L A 0 1 0 is the element of the adjacency matrix and takes either 1 or 0, while is the difference between the characteristic degrees of the corresponding eigenvectors and is greater than 1 (See Methods). Therefore, if = βα A 1 and k β is close to k α , the above ratio may not be small, namely, the contribution from β 0 to the first-order correction α 1 may be inaccurate at the nodes with degree k β .
We note, however, that the perturbation theory can still qualitatively account for the localization property in such cases. In Fig. 1, the eigenvector α is localized at the diagonal nodes whose indices satisfy α  i and whose degrees are close to k α , because the node indices {i} are sorted so that ≥ ≥ ≥  k k k N 1 2 . Our main aim is to explain that α is almost vanishing at the non-diagonal nodes with α −  i 1. The degree k β of such non-diagonal nodes are generally far from k α for networks with degree heterogeneity, so that the above ratio would generally be small and the first-order correction α 1 would reliably be obtained for such nodes. Thus, the perturbation theory can account for why the eigenvector takes tiny components at non-diagonal nodes even if they can give inaccurate results for diagonal nodes.
Indeed, as explained in the Results, the theory can reproduce the true eigenvalues very accurately and account for the localization property qualitatively well for all three networks shown in Fig. 1. Furthermore, it can even predict precise localizing patterns quantitatively well for some of the eigenvectors.