Top Eigenpair Statistics for Weighted Sparse Graphs

TL;DR

This paper introduces a formalism using cavity and replica methods to analytically compute the statistics of the top eigenpair in weighted sparse graphs, validated by numerical simulations.

Contribution

It develops a novel analytical approach to determine top eigenpair statistics in weighted sparse graphs using cavity and replica techniques.

Findings

Analytical expressions match numerical diagonalization results.

Method applies to random regular graphs and sparse Markov transition matrices.

Provides insights into eigenvector component contributions from nodes of different degrees.

Abstract

We develop a formalism to compute the statistics of the top eigenpair of weighted sparse graphs with finite mean connectivity and bounded maximal degree. Framing the problem in terms of optimisation of a quadratic form on the sphere and introducing a fictitious temperature, we employ the cavity and replica methods to find the solution in terms of self-consistent equations for auxiliary probability density functions, which can be solved by population dynamics. This derivation allows us to identify and unpack the individual contributions to the top eigenvector's components coming from nodes of degree $k$. The analytical results are in perfect agreement with numerical diagonalisation of large (weighted) adjacency matrices, and are further cross-checked on the cases of random regular graphs and sparse Markov transition matrices for unbiased random walks.

Figures

Equations

$$J_{ij}=c_{ij}{K}_{ij},$$

$$P(k_i=k)={\mathcal{N}}^{-1}{\mathrm{e}}^{-\bar{c}}{\bar{c}}^k/k!,k=0,\dots ,k_{\max },$$

$$J{\mathbf{v}}_{\alpha }={\lambda }_{\alpha }{\mathbf{v}}_{\alpha },$$

$$N{\lambda }_1=\underset{|\mathbf{v}{|}^2=N}{\min }\left[\hat{H}(\mathbf{v})\right]=\underset{|\mathbf{v}{|}^2=N}{\min }\left[-\frac{1}{2}\left(\mathbf{v},J\mathbf{v}\right)\right],$$

$$-\frac{1}{2}{\lambda }_{1}N\le \hat{H}\left(\mathbf{v}\right)\le -\frac{1}{2}{\lambda }_{N}N,$$

$$P_{\beta ,J}(\mathbf{v})=\frac{1}{Z}\exp \left[\frac{\beta }{2}(\mathbf{v},J\mathbf{v})\right]\delta (|\mathbf{v}{|}^2-N),$$

$$P_{\beta ,J}\left(\mathbf{v}\right)=\frac{1}{Z}\exp \left\{\beta \left[\frac{1}{2}\left(\mathbf{v},J\mathbf{v}\right)-\frac{\lambda }{2}\left(\mathbf{v},\mathbf{v}\right)\right]\right\},$$

$$Z=\int \mathrm{d}\mathbf{v}\exp \left\{\beta \left[\frac{1}{2}\left(\mathbf{v},J\mathbf{v}\right)-\frac{\lambda }{2}\left(\mathbf{v},\mathbf{v}\right)\right]\right\},$$

$$P_i\left(v_i\right)=\frac{1}{Z_i}\exp \left(-\beta \frac{\lambda }{2}{v}_i^2\right)\int \mathrm{d}{\mathbf{v}}_{\partial i}\exp \left(\beta \sum _{j\in \partial i}{J}_{ij}{v}_i{v}_j\right)P^{(i)}\left({\mathbf{v}}_{\partial i}\right),$$

$$P_i\left(v_i\right)=\frac{1}{Z_i}\exp \left(-\beta \frac{\lambda }{2}{v}_i^2\right)\prod _{j\in \partial i}\int \mathrm{d}{v}_j\exp \left(\beta J_{ij}{v}_i{v}_j\right)P_j^{(i)}\left(v_j\right).$$

$$P_j^{(i)}\left(v_j\right)=\frac{1}{Z_j^{(i)}}\exp \left(-\beta \frac{\lambda }{2}{v}_j^2\right)\prod _{\ell \in \partial j\backslash i}\int \mathrm{d}{v}_{\ell }\exp \left(\beta J_{j\ell }{v}_j{v}_{\ell }\right)P_{\ell }^{(j)}\left(v_{\ell }\right),$$

$$P_j^{(i)}\left(v_j\right)=\sqrt{\frac{\beta {\Omega }_j^{(i)}}{2\pi }}\exp \left(-\frac{\beta {H_j^{(i)}}^2}{2{\Omega }_j^{(i)}}\right)\exp \left(-\frac{\beta }{2}{\Omega }_j^{(i)}{v}_j^2+\beta H_j^{(i)}{v}_j\right),$$

$$P_j^{(i)}\left(v_j\right)=\frac{1}{Z_j^{(i)}}\exp \left(-\frac{\beta }{2}\lambda v_j^2\right)\prod _{\ell \in \partial j\backslash i}\exp \left[\frac{\beta }{2}\frac{{\left(J_{j\ell }{v}_j+H_{\ell }^{(j)}\right)}^2}{{\Omega }_{\ell }^{(j)}}\right].$$

$${\Omega }_j^{(i)}$$

$$H_j^{(i)}$$

$$P_i\left(v_i\right)=\frac{1}{Z_i}\exp \left(-\frac{\beta }{2}{\Omega }_i{v}_i^2+\beta H_i{v}_i\right),$$

$${\Omega }_i$$

$$H_i$$

$$P_i(v_i)=\delta \left(v_i-\frac{H_i}{{\Omega }_i}\right),$$

$$q\left(\omega ,h\right)$$

$$={\left(\sum _{i=1}^N{k}_i\right)}^{-1}\sum _{i=1}^N\sum _{j\in \partial i}\delta \left(\omega -{\Omega }_j^{(i)}\right)\delta \left(h-H_j^{(i)}\right),$$

$$q\left(\omega ,h\right)=\sum _{k=1}^{k_{\max }}r\left(k\right)\int \left[\prod _{\ell =1}^{k-1}\mathrm{d}q\left({\omega }_{\ell },h_{\ell }\right)\right]{⟨\delta \left(\omega -\lambda +\sum _{\ell =1}^{k-1}\frac{K_{\ell }^2}{{\omega }_{\ell }}\right)\delta \left(h-\sum _{\ell =1}^{k-1}\frac{h_{\ell }{K}_{\ell }}{{\omega }_{\ell }}\right)⟩}_{{\{K\}}_{k-1}},$$

$$r\left(k\right)=\frac{kp\left(k\right)}{⟨k⟩},$$

$$Q\left(\Omega ,H\right)=\frac{1}{N}\sum _{i=1}^N\delta \left(\Omega -{\Omega }_i\right)\delta \left(H-H_i\right).$$

$$Q\left(\Omega ,H\right)=\sum _{k=0}^{k_{\max }}p\left(k\right)\int \left[\prod _{\ell =1}^k\mathrm{d}q\left({\omega }_{\ell },h_{\ell }\right)\right]{⟨\delta \left(\Omega -\lambda +\sum _{\ell =1}^k\frac{K_{\ell }^2}{{\omega }_{\ell }}\right)\delta \left(H-\sum _{\ell =1}^k\frac{h_{\ell }{K}_{\ell }}{{\omega }_{\ell }}\right)⟩}_{{\{K\}}_k},$$

$$\varrho (u)=⟨\frac{1}{N}\sum _{i=1}^N\delta \left(u-v_1^{(i)}\right)⟩=\int \mathrm{d}\Omega \mathrm{d}HQ\left(\Omega ,H\right)\delta \left(u-\frac{H}{\Omega }\right).$$

$$1=⟨u^2⟩=\int \mathrm{d}\Omega \mathrm{d}HQ\left(\Omega ,H\right)\frac{H^2}{{\Omega }^2}.$$

$$P\left(\left\{J_{ij}\right\}|\left\{k_i\right\}\right)=P\left(\left\{c_{ij}\right\}|\left\{k_i\right\}\right)\prod _{i<j}{\delta }_{K_{ij},K_{ji}}p\left(K_{ij}\right),$$

$$P\left(\left\{c_{ij}\right\}|\left\{k_i\right\}\right)=\frac{1}{\mathcal{M}}\prod _{i<j}{\delta }_{c_{ij},c_{ji}}\left(\frac{c}{N}{\delta }_{c_{ij},1}+\left(1-\frac{c}{N}\right){\delta }_{c_{ij},0}\right)\prod _{i=1}^N{\delta }_{{\sum }_j{c}_{ij},k_i},$$

$$P\left(\left\{c_{ij}\right\}\right)=\prod _{i<j}{\delta }_{c_{ij},c_{ji}}\left(\frac{c}{N}{\delta }_{c_{ij},1}+\left(1-\frac{c}{N}\right){\delta }_{c_{ij},0}\right).$$

$$P\left(\left\{J_{ij}\right\}\right)=P\left(\left\{c_{ij}\right\}\right)\prod_{i<j}\delta_{K_{ij},K_{ji}}p\left(K_{ij}\right)\$$

$${⟨{\lambda }_1⟩}_J=\underset{\beta \to \infty }{\lim }\frac{2}{\beta N}{⟨\ln Z⟩}_J,Z=\int \mathrm{d}\mathbf{v}\exp \left[\frac{\beta }{2}\left(\mathbf{v},J\mathbf{v}\right)\right]\delta \left({\left|\mathbf{v}\right|}^2-N\right),$$

Full text

1. 1 Introduction

2. 2 Formulation of the problem

3. 3 Cavity approach

4. 3.1 Single instance

5. 3.2 Thermodynamic limit $N\to\infty$

6. 4 Replica derivation

7. 4.1 Typical largest eigenvalue

8. 4.1.1 Erdős-Rényi graph: weighted adjacency matrix.

9. 4.1.2 Random regular graph: adjacency matrix.

10. 4.2 Density of the top eigenvector’s components

11. 4.2.1 Erdős-Rényi graph: weighted adjacency matrix.

12. 4.2.2 Random regular graph: adjacency matrix

13. 4.2.3 Large-$c$ limit for weighted adjacency matrices

14. 5 Application: sparse random Markov transition matrices

15. 6 Population dynamics

16. 7 Conclusions

17. A The solution of the single instance self-consistency equations and the non-backtracking operator.

18. B Exact replica calculation for the largest eigenvalue for any bounded degree distribution $p(k)$.

1 Introduction

The largest eigenvalue and the associated top eigenvector of a $N\times N$ matrix $J$ play a very important role in many applications. In multivariate data analysis and Principal Component Analysis, the top eigenpair of the covariance matrix provides information about the most relevant correlations hidden in the dataset [1, 2]. These extremal questions also arise in connection with synchronisation problems on networks [3], percolation problems [4], linear stability of coupled ODEs [5], financial stability [6] and several other problems in physics and chemistry, connected to the applications of Perron’s theorem [7]. Also in the realm of quantum mechanics, the search for the ground state of a complicated Hamiltonian essentially amounts to solving the top eigenpair problem for a differential operator [8]. The top eigenpair is also relevant in signal reconstruction problems employing algorithms based on the spectral method [9]. In the context of graph theory, the eigenvectors of both adjacency and Laplacian matrices are employed to solve combinatorial optimisation problems, such as graph 3-colouring [10] and to develop clustering and cutting techniques [11, 12, 13]. In particular, the top eigenvector of graphs is intimately related to the “ranking” of the nodes of the network [14]. Indeed, beyond the natural notion of ranking of a node given by its degree, the relevance of a node can be estimated from how “important” its neighbours are. The vector expressing the importance of each node is exactly the top eigenvector of the network adjacency matrix. Google PageRank algorithm works in a similar way [15, 16]: the PageRanks vector is indeed the top eigenvector of a large Markov transition matrix between web pages.

When the matrix $J$ is random and symmetric with i.i.d. entries, analytical results on the statistics of the top eigenpair date back to the classical work by Füredi and Komlós [17]: the largest eigenvalue of such matrices follows a Gaussian distribution with finite variance, provided that the moments of the distribution of the entries do not scale with the matrix size. This result directly relates to the largest eigenvalue of Erdős-Rényi (E-R) [18] adjacency matrices in the case when the probability $p$ for two nodes to be connected does not scale with the matrix size $N$. This result has been then extended by Janson [19] in the case when $p$ is large. However, in our analysis we will be mostly dealing with the sparse case, i.e. when $p=c/N$, with $c$ being the constant mean degree of nodes (or equivalently, the mean number of nonzero elements per row of the corresponding adjacency matrix). In this sparse regime, Krivelevich and Sudakov [20] proved a theorem stating that for any constant $c$ the largest eigenvalue of Erdős-Rényi graph diverges slowly with $N$ as $\sqrt{\log{N}/\log{\log{N}}}$. To ensure that the largest eigenvalue remains $\sim\mathcal{O}(1)$, the nodes with very large degree must be pruned (see [21]).

The characterisation of eigenvectors properties has proved to be much harder and is generally a less explored area of random matrix theory. Excluding the cases of i) invariant ensembles, where eigenvector components follow the celebrated Porter-Thomas distribution [22, 23], ii) dense non-Hermitian matrices (see for instance the seminal works of Chalker and Mehlig [24] along with results about correlations between eigenvectors [25, 26] and some more recent applications [27, 28, 29, 30]) and iii) perturbed matrices [31, 32, 33, 34, 35], systematic results are scarcer for sparse Hermitian matrices, especially in the limit of high sparsity. Indeed, although Gaussian statistics and delocalisation of eigenvectors are known properties of adjacency matrices of Erdős-Rényi and random regular graphs in the case where the mean degree $c=c(N)$ diverges with $N$ [36, 37, 38], very few results are available for the high sparsity regime, i.e. with fixed $c$. In this limit, numerical studies have shown that most of the eigenvectors of a random regular graph follow a Gaussian distribution [39], as well as almost-eigenvectors [40], whereas Erdős-Rényi eigenvectors are localised especially for low values of $c$.

The statistics of the first eigenvector components for very sparse symmetric random matrices was first considered in the seminal works by Kabashima and collaborators [41, 42, 43], which constitute the starting point of our analysis. The focus there is on specific classes of real sparse random matrices, i.e. when the matrix connectivity is either a random regular graph or a mixture of multiple degrees, and the nonzero elements are drawn from a Bernoulli distribution. More precisely, in [41] and [43] the cavity method was employed for the top eigenpair problem, while in [42] the replica formalism was instead adopted to study the same problem in the thermodynamic limit, recovering cavity results. Our aim is to analyse and develop both the cavity and replica formalisms they pioneered even further, and to present them in a unified way that looks - at least to our eyes - more transparent.

We will be implementing a Statistical Mechanics formulation of the top eigenpair problem, using both the cavity (Section 3) and replica (Section 4) methods - borrowed from the standard arsenal of disordered systems physics - as main solving tools.

The replica method, widely used in the physics of spin glasses [44], was first introduced in the context of random matrices by Edwards and Jones [45] to compute the average spectral density of random matrices defined in terms of the joint probability density function (pdf) of their entries. Building on this approach, Bray and Rodgers in their seminal paper [46] were able to express the spectral density of Erdős-Rényi adjacency matrices as the solution of a (nearly intractable) integral equation. Therefore, asymptotic analyses for large average connectivities [46], and approximation schemes such as the single defect approximation (SDA) and the effective medium approximation (EMA) [47, 48] were first developed as a way around this hindrance. An alternative approach was pursued in [49] (see also [50]): starting from Bray-Rodgers replica-symmetric setup [46], the functional order parameters of the theory are expressed as continuous superpositions of Gaussians with fluctuating variances, as suggested by earlier solutions of models for finitely coordinated harmonically coupled systems [51]. This formulation gives rise to non-linear integral equations for the probability densities of such variances, which can be efficiently solved by a population dynamics algorithm. Our paper will follow a similar approach in Section 4.

The cavity method [52], also known as Bethe-Peierls or belief-propagation method, was introduced in the context of disordered systems and sparse random matrices as a more intuitive and straightforward alternative to replicas: the two methods are known to provide the same results for the spectral density of graphs [53], even though a general, first-principle proof of their equivalence does not seem to be currently available. A rigorous proof of the correctness of cavity method and the tree-like approximation for finitely coordinated graphs is given in [54]. One of the advantages of the cavity method is that it allows one to solve the spectral problem for very large single instances of sparse random graphs, as done in [55]. Both the replica and cavity approaches in [49] and [55] retrieve known results such as the Kesten-McKay law for the spectra of random regular graphs [56, 57], the Marčenko-Pastur law and the Wigner’s semicircle law respectively for sparse covariance matrices and for Erdős-Rényi adjacency matrices in the large mean degree limit. Both approaches have also been employed to characterise the spectral density of sparse Markov matrices [58, 59] and graphs with modular [60] and small-world [61] structure and with topological constraints [62]. The localisation transition for sparse symmetric matrices was studied in [63]. The two methods have also been extended to the study of the spectral density of sparse non-Hermitian matrices [64, 65], whereas eigenvalue outliers have been considered in [66]; for an excellent review, see [67]. The spectral properties of the Hashimoto non-backtracking operator - arising in the cavity solution (see A for details) have been investigated in [68, 69, 70]. In this paper, we propose a “grand canonical” cavity derivation that differs in the details from [41] (see Section 3). We also provide a detailed analysis of the single-instance recursion equations, showing that their convergence is strictly related to the spectral properties of a modified non-backtracking operator associated with the single-instance matrix. At the same time, building on the insights coming from the replica treatment, we are able to better understand the behaviour of the stochastic recursions that provide the solution of the top eigenpair problem in the thermodynamic limit. Furthermore, the population dynamics algorithm employed to solve these recursions allows us to characterise the distributions of the cavity fields in the thermodynamic limit and identify the individual contributions of nodes of different degrees $k$ to the top eigenvector’s entries.

The plan of the paper is as follows. In Section 2, we will formulate the problem and provide the main starting points. In Section 3, we will describe the cavity approach to the problem, first for the single instance case (in 3.1), and then in the thermodynamic limit (in 3.2). In Section 4, we formulate the replica approach to the same problem, first focussing on the largest eigenvalue problem (in 4.1) and then on the density of top eigenvector’s components (in 4.2). For both problems, we take the weighted Erdős-Rényi and random regular graphs as representative examples. In Section 5 we build on our previous results to complete the picture for Markov transition matrices on a random graph structure. In Section 6, we provide the details of the population dynamics algorithm, and in Section 7 we offer a summary and outlook for future research. In A, we provide a detailed discussion of the single-instance cavity approach and associated non-backtracking operator. In B, we offer a detailed replica derivation of the typical location of the largest eigenvalue for sparse graphs characterised by a generic degree distribution $p(k)$.

2 Formulation of the problem

We consider a sparse random $N\times N$ symmetric matrix $J=\left(J_{ij}\right)$, with real i.i.d. entries. The matrix entries are defined as

[TABLE]

where the $c_{ij}\in\{0,1\}$ constitute the connectivity matrix, i.e. the adjacency matrix of the underlying graph, and the $K_{ij}$ encode bond weights. We will typically consider the case of Poissonian highly sparse connectivity - where the node degrees $k_{i}$ (or equivalently the number of nonzero elements per row of $J$) fluctuate according to a bounded Poisson distribution

[TABLE]

with the mean degree a finite constant $c\equiv\left\langle k\right\rangle$ and $\mathcal{N}=\sum_{k=0}^{k_{\mathrm{max}}}\mathrm{e}^{-\bar{c}}\bar{c}^{k}/k!$ to ensure normalisation. The bond weights $K_{ij}$ will be i.i.d. random variables drawn from a parent pdf $p(K)$ with bounded support. This setting is sufficient to ensure that the largest eigenvalue $\lambda_{1}$ of $J$ will remain of $\mathcal{O}(1)$ for $N\to\infty$.

The spectral theorem ensures that $J$ can be diagonalised via an orthonormal basis of eigenvectors $\bm{v}_{\alpha}$ with corresponding real eigenvalues $\lambda_{\alpha}$ ($\alpha=1,\ldots,N$),

[TABLE]

for each eigenpair $\alpha=1,\ldots,N$. We assume that there is no eigenvalue degeneracy, and that they are sorted $\lambda_{1}>\lambda_{2}>\ldots>\lambda_{N}$.

The goal of this work is to set up a formalism based on the statistical mechanics of disordered systems to find:

• •

The average (or typical value) $\langle\lambda_{1}\rangle_{J}$ of the largest eigenvalue $\lambda_{1}$.

• •

The density \varrho(u)=\Big{\langle}\frac{1}{N}\sum_{i=1}^{N}\delta(u-v_{1}^{(i)})\Big{\rangle}_{J} of the top eigenvector’s components, $\bm{v}_{1}=(v_{1}^{(1)},\ldots,v_{1}^{(N)})$ ,

where the average $\langle\cdot\rangle_{J}$ is taken over the distribution of the matrix $J$.

The problem can be formulated as the optimisation problem of a quadratic function $\hat{H}(\bm{v})$, according to which $\bm{v}_{1}$ is the vector normalized to $N$ that realises the condition

[TABLE]

as dictated by the Courant-Fischer definition of eigenvectors. The round brackets $\left(\cdot,\cdot\right)$ indicate the dot product between vectors in $\mathbb{R}^{N}$. It is easy to show that $\hat{H}\left(\bm{v}\right)$ is bounded

[TABLE]

and attains its minimum when computed on the top eigenvector.

For a fixed matrix $J$, the minimum in (4) can be computed by introducing a fictitious canonical ensemble of $N$-dimensional vectors $\bm{v}$ at inverse temperature $\beta$, whose Gibbs-Boltzmann distribution reads

[TABLE]

where the delta function enforces normalisation. Clearly, in the low temperature limit $\beta\to\infty$, only one ’state’ remains populated, which corresponds to $\bm{v}=\bm{v}_{1}$, the top eigenvector of the matrix $J$. The hard normalisation constraint can also be relaxed for our purposes, and replaced with a soft, ”grand canonical” version

[TABLE]

where $\lambda$ is an auxiliary Lagrange multiplier. The two versions above are expected to provide the same physical results in the limits $\beta,N\to\infty$, as we explicitly demonstrate by using (7) for our cavity treatment in Section 3, and (6) as a starting point of our replica calculation in Section 4.

3 Cavity approach

In what follows, we will use a cavity method formulation for the top eigenpair problem which is deeply rooted in the statistical mechanics approach to disordered systems. Our formulation provides equations for the statistics of the top eigenpair that are fully equivalent to those found earlier by Kabashima et al. in [41]. Our treatment, however, brings more neatly to the surface a few subtleties related to the solution of self-consistency equations and their range of applicability, this way providing a more transparent derivation.

The central idea of the cavity method [52] consists in computing observables related to a given node, relying on some information concerning its neighbourhood when the node of interest is removed from the network. It is useful every time the underlying graph has a finite connectivity structure: its predictions become exact for trees and approximately exact for tree-like structures (where loops are negligible) such as graphs in the high sparsity regime.

3.1 Single instance

Consider for the time being a single instance of the random matrix $J$. Starting from the soft-constraint distribution (7), whose partition function is

[TABLE]

it is trivial to notice that the condition $\lambda>\lambda_{1}$ is necessary to ensure convergence for all $\beta$.

The marginal distribution of the component $v_{i}$, obtained by integrating out all other components in (7), and using the sparsity condition $J_{ij}=0$ if $j\notin\partial i$ (where $\partial i$ denotes the immediate neighbourhood of $i$) is

[TABLE]

where $P^{(i)}(\bm{v}_{\partial i})$ is the joint distribution of the components pertaining to the immediate neighbourhood of $i$, $\partial i$, when the node $i$ has been removed. Indeed, all the components outside $\partial i$ can be integrated out without difficulty, and the resulting constant term can be just reabsorbed in the normalisation constant. $P^{(i)}(\bm{v}_{\partial i})$ is also known as cavity probability distribution.

Adopting now a tree-like approximation, which is accurate for very sparse graphs, all nodes $j$ in $\partial i$ are connected with each other only through $i$ (see Fig. 1), therefore they get disconnected when the node $i$ is removed from the network: this implies that the integral appearing in (9) factorises as

[TABLE]

In the same way, a similar expression can be derived for the marginal cavity distribution $P_{j}^{(i)}(v_{j})$ now appearing in (10). Iterating the reasoning as before, and further removing the node $j\in\partial i$ in the network in which the node $i$ had already been removed, one can write

[TABLE]

where the symbol $\partial j\backslash i$ denotes the neighbourhood of $j$ excluding $i$.

Equation (11) has now become a self-consistent equation for the cavity probability distributions, which can be solved by a Gaussian ansatz for $P_{j}^{(i)}\left(v_{j}\right)$

[TABLE]

where the parameters $\Omega_{i}^{(i)}$ and $H_{j}^{(i)}$ are called cavity fields. This ansatz is chosen to obtain a solution $\bm{v}$ whose components are not peaked at zero in the $\beta\rightarrow\infty$ limit. Inserting the Gaussian ansatz (12) in (11) and performing the resulting Gaussian integrals, one obtains

[TABLE]

Comparing the coefficients of the same powers of $v_{j}$ between (12) and (13), we obtain the following two self-consistent relations which define the cavity fields $\Omega_{i}^{(i)}$ and $H_{j}^{(i)}$

[TABLE]

These equations have been obtained before in [41].

The Gaussian ansatz (12) can then be inserted in (10), resulting in a Gaussian distribution for the single-site marginals

[TABLE]

where the $N$ coefficients $\Omega_{i}$ and $H_{i}$ are given by the following equations

[TABLE]

Here, $\Omega_{j}^{(i)}$ and $H_{j}^{(i)}$ are the fixed-point solutions of (14) and (15).

In the limit $\beta\to\infty$, the marginal distribution (16) converges to

[TABLE]

from which one concludes that the components of the top eigenvector of the fixed matrix $J$ (a single instance of the ensemble) must be given by $v_{1}^{(i)}=H_{i}/\Omega_{i}$, where $H_{i}$ and $\Omega_{i}$ are the values obtained from (17) and (18), after the fixed-points of the recursions (14) and (15) have been obtained.

A detailed discussion on how to solve the above recursions in practice and on the role of the (yet unspecified) multiplier $\lambda$ is deferred to A. Although this derivation only relies on the tree-like approximation for the local connectivity and is arguably very easy and intuitive, it is not particularly interesting as it stands: the complexity of the cavity algorithm for a single instance is actually higher than a high-precision, direct diagonalisation of the matrix $J$, therefore it is of little practical use per se. It is, however, a conceptually necessary ingredient to discuss infinite-size matrices, as we do in the next subsection.

3.2 Thermodynamic limit $N\to\infty$

In an infinitely large network, it is no longer possible to keep track of an infinite number of cavity fields. Following [41], we consider first the joint probability density that the cavity fields of type $\Omega_{j}^{(i)}$ and $H_{j}^{(i)}$ take up values around $\omega$ and $h$

[TABLE]

where $N$ is now large but finite. This is a properly normalised pdf: indeed, we can associate two cavity fields $\Omega_{j}^{(i)}$ and $H_{j}^{(i)}$ to any link $(i,j)$ of the network. Since every node $i$ is the source of $k_{i}$ links, their total number is given by $\sum_{i=1}^{N}k_{i}$.

Next, one may appeal to the single-instance update rules given by (14) and (15) to characterise the above distribution self-consistently, as is done in [41]. It should be stressed that in an infinitely large network links can only be distinguished by the degree of the node they are pointing to. Thus, for a given edge $(i,j)$ pointing to a node $j$ of degree $k$, the values $\omega$ and $h$ of the pair of cavity fields $\Omega_{j}^{(i)}$ and $H_{j}^{(i)}$ living on this edge are determined respectively by the $k-1$ values $\left\{\omega_{\ell}\right\}$ and $\left\{h_{\ell}\right\}$ of the cavity fields $\Omega_{\ell}^{(j)}$ and $H_{\ell}^{(j)}$ living on each of the edges connecting $j$ with its neighbours $\ell\in\partial j\backslash i$. In an infinite system, these values can be thought of as $k-1$ independent realisations of the random variables of types $\Omega_{j}^{(i)}$ and $H_{j}^{(i)}$, drawn from their joint pdf $q(\omega,h)$. The entries of $J$ that appear in the single instance recursions (14) and (15) are replaced by a set $\left\{K_{\ell}\right\}_{k-1}$ of $k-1$ independent realisations of the random variable $K$, each distributed according to the pdf $p(K)$ of bond weights. The full distribution $q(\omega,h)$ is then obtained by weighing each edge contribution with the probability $r(k)$ of having a random link pointing to a node of degree $k$ and summing up over all possible degrees up to $k_{\mathrm{max}}$, leading to the self-consistency equation

[TABLE]

where $\mathrm{d}q\left(\omega_{\ell},h_{\ell}\right)\equiv\mathrm{d}\omega_{\ell}\mathrm{d}h_{\ell}q\left(\omega_{\ell},h_{\ell}\right)$, and the average $\langle\cdot\rangle_{{\{K\}}_{k-1}}$ is taken over $k-1$ independent realisations of the random variable $K$. We recall that

[TABLE]

where $p\left(k\right)$ is the probability of having a node of degree $k$ and $\left\langle k\right\rangle=\sum_{k}kp\left(k\right)$ [71]. The sum in (21) starts from $k=1$ since we should not be concerned with isolated nodes.

Eq. (21) is generally solved via a population dynamics algorithm (see Section 6 for details). In some exceptional cases, such as for adjacency matrices of random regular graphs, it can be solved analytically (see discussion in sections 4.1.2 and 4.2.2 below).

In a similar fashion, the joint pdf of the coefficients $\Omega_{i}$ and $H_{i}$ can be expressed as

[TABLE]

In this case, there is a pair of marginal coefficients $\Omega_{i}$ and $H_{i}$ living on each node. Since in the infinite size limit the nodes can only be distinguished by their degree, following the same line of reasoning that led to (21), the joint pdf of the random variables of the type $\Omega_{i}$ and $H_{i}$ in the thermodynamic limit can be written as

[TABLE]

where $p(k)$ is the degree distribution. Here, $q\left(\omega_{\ell},h_{\ell}\right)$ is the fixed-point distribution of cavity fields, i.e. the solution of the self-consistency equation (21), which should therefore be solved beforehand.

The distribution of the top eigenvector’s components in the thermodynamic limit is then obtained in terms of the pdf $Q\left(\Omega,H\right)$ in (24), exploiting the analogy with the single-instance case in (19), and reads

[TABLE]

Both equations (21) and (24) still depend on the parameter $\lambda$: it must be fixed taking into account the normalisation of the top eigenvector. This condition amounts to requiring that

[TABLE]

Crucially, the value of $\lambda$ for which the above normalisation condition is satisfied turns out to be exactly equal to the typical largest eigenvalue, $\lambda\equiv\left\langle\lambda_{1}\right\rangle_{J}$. Indeed, for every $\lambda>\left\langle\lambda_{1}\right\rangle_{J}$, the distribution of the $h$’s shrinks to a delta peak located at zero, whereas for $\lambda<\left\langle\lambda_{1}\right\rangle_{J}$, negative values of the $\omega$’s start to appear while the $h$’s grow without bounds in the self-consistency solution of (21). This is not surprising, since $\lambda<\left\langle\lambda_{1}\right\rangle_{J}$ is precisely the range of values for $\lambda$ that makes the Gibbs-Boltzmann distribution (7) not normalisable.

As a final remark on the cavity solution, the equations (21) and (25) will match respectively (69) and (111) obtained via the replica method in Section 4 below.

The discussion above has the advantage of leading rather quickly to the results (24) and (25). It is, however, instructive to reconsider this problem from the point of view of the replica approach, which provides a lengthier but rather systematic procedure, and arrives at the very same equations while departing from very different premises. Both approaches (cavity or replicas) present different advantages and drawbacks - especially if seen through the prism of full mathematical rigour - and it is therefore of interest to compare them back to back. For the sake of clarity, we will keep the two pathways (typical largest eigenvalue vs. density of top eigenvector’s components) clearly separate until the point where we realise that the same self-consistency equation governs the statistics of both quantities.

4 Replica derivation

In this section, we evaluate the average location of the largest eigenvalue and the density of top eigenvectors’ components within the replica framework. The starting point of our analysis is the formalism pioneered in [42]. However, our derivation is not confined to specific connectivity distributions of the matrix entries as in [42], and thus provides a rather general and robust methodology that can be applied to any graph with finite mean connectivity and bounded maximal degree. We also make a quite transparent and convincing case for the equivalence between the cavity and replica methods in these problems. Moreover, as we did for the cavity approach, we thoroughly discuss bounds on the values of parameters that guarantee a converging solution.

4.1 Typical largest eigenvalue

Consider again a $N\times N$ symmetric matrix $J_{ij}=c_{ij}K_{ij}$. The joint distribution of the matrix entries is

[TABLE]

where, in the framework of the configuration model [61], the distribution $P\left(\left\{c_{ij}\right\}\middle|\left\{k_{i}\right\}\right)$ of connectivities $\left\{c_{ij}\right\}$ compatible with a given degree sequence $\left\{k_{i}\right\}$ is given by

[TABLE]

and the pdf $p\left(K_{ij}\right)$ of bond weights (with compact support and upper edge $\zeta$) can be kept unspecified until the very end.

It has been shown in many works [49, 59] that a convenient shortcut for the calculation consists in replacing the “microcanonical” Eq. (28) with the standard Erdős-Rényi connectivity distribution

[TABLE]

Although Eq. (29) technically gives rise to an unbounded Poisson degree distribution with mean $c$ – and therefore a largest eigenvalue whose location typically grows with $N$ [20] – the final results (e.g. Eq. (68)) can be easily adjusted and extended to cover any degree distribution $p(k)$ with finite mean and bounded largest degree. For simplicity, we will therefore consider the distribution of the matrix entries to be simply

[TABLE]

at the outset, where $P\left(\left\{c_{ij}\right\}\right)$ is given by (29). Once the Erdős-Rényi Poissonian degree distribution has appeared in the formulae, it will be straightforward to replace it with the actual finite-mean degree distribution of interest (for instance, the truncated Poisson distribution (2)). In B, we will however provide a first-principle derivation for sparse graphs with a generic degree distribution $p(k)$, without relying on any shortcut.

The average of the largest eigenvalue can be computed as the formal limit

[TABLE]

in terms of the quenched free energy of the model defined in (6).

The average over $J$ is computed using the replica trick as follows

[TABLE]

where $n$ is initially taken as an integer, and then analytically continued to real values in the vicinity of $n=0$.

The replicated partition function is

[TABLE]

Taking the average w.r.t the joint distribution (30) of matrix entries yields [49]

[TABLE]

where $\left\langle\cdot\right\rangle_{K}$ denotes averaging over the single-variable pdf $p\left(K\right)$ characterising the i.i.d. bond weights $K_{ij}$.

We also employ a Fourier representation of the Dirac delta enforcing the normalisation constraints

[TABLE]

The replicated partition function thus becomes

[TABLE]

In order to decouple sites, we introduce the functional order parameter

[TABLE]

where the symbol $\vec{v}$ denotes a $n$-dimensional vector in replica space. We enforce its definition using the integral identity

[TABLE]

In terms of this order parameter and its conjugate, the replicated partition function can be written as

[TABLE]

The multiple integral in the last line above factorises into $N$ identical copies of the same $n$-dimensional integral, and can thus be written as

[TABLE]

where $\mathrm{Log}$ denotes the principal branch of the complex logarithm.

Therefore, the replicated partition function takes a form amenable to a saddle point evaluation for large $N$ (where we assume we can safely exchange the limits $n\to 0$ and $N\to\infty$)

[TABLE]

where

[TABLE]

and

[TABLE]

The stationarity of the action $S_{n}$ w.r.t. variations of $\varphi$ and $\hat{\varphi}$ requires that the order parameter at the saddle point $\varphi^{\star}$ and its conjugate $\hat{\varphi}^{\star}$ satisfy the following coupled equations

[TABLE]

which have to be solved together with the stationarity conditions w.r.t each component $\lambda_{\bar{a}}$ of $\vec{\lambda}$

[TABLE]

The equations (47) and (48) bear a striking resemblance with the saddle-point equations leading to the spectral density of Erdős-Rényi random graphs [46, 49], except for the fact that the “Hamiltonian” of our problem is real-valued and includes the inverse temperature $\beta$. Following [49], we will now search for replica-symmetric solutions written in the form of superpositions of uncountably infinite Gaussians with a non-zero mean. This ansatz will be preserving permutational symmetry between replicas, but (at odds with the choice in [49]) not the rotational invariance in the space of replicas111A rotationally invariant ansatz would not produce a physically meaningful result for this problem.:

[TABLE]

where

[TABLE]

To justify the procedure above, on one hand the replica symmetric ansatz has been known for quite a while to lead to the correct results for the spectral problem of sparse random matrices [45, 46, 49, 72]. On the other hand, it is known that expressing the order parameter as a superposition of Gaussian pdfs provides the correct solution for harmonically coupled system [51].

In (51) and (52), $\pi$ and $\hat{\pi}$ are normalised joint pdfs of the parameters appearing in the Gaussian distributions, while $\hat{c}$ is introduced taking into account that $\mathrm{i}\hat{\varphi}(\vec{v})$ needs not be normalised. The advantage of writing an ansatz in this form is that - once inserted into (47) and (48) - it makes it possible to perform explicitly the $\vec{v}$-integrals, eventually leading to simpler coupled equations for $\pi$ and $\hat{\pi}$, as detailed below. The convergence of the $\vec{v}$-integrals will also impose the following conditions on $\omega$ and $\hat{\omega}$: $\omega>\hat{\omega}$ and $\omega>\zeta$ (where $\zeta$ is the upper edge of the support of the pdf $p(K)$ of bond weights).

As a further remark, the different signs in front of $\omega$ and $\hat{\omega}$ in (51) and (52) are picked with an eye towards performing the subsequent $\vec{v}$-integrals: since $\mathrm{i}\hat{\varphi}^{\star}(\vec{v})$ is not a pdf, $\hat{\omega}$ being positive is not problematic.

Rewriting the action in terms of $\pi$ and $\hat{\pi}$, after performing the $\vec{v}$-integrations, and extracting the leading $n\to 0$ contribution yields

[TABLE]

where we have introduced the shorthands

[TABLE]

and $\{\mathrm{d}\hat{\pi}\}_{s}=\prod_{\ell=1}^{s}\mathrm{d}\hat{\omega}_{\ell}\mathrm{d}\hat{h}_{\ell}\hat{\pi}(\hat{\omega}_{\ell},\hat{h}_{\ell})$, along with $\{\hat{\omega}\}_{s}=\sum_{\ell=1}^{s}\hat{\omega}_{\ell}$ and $\{\hat{h}\}_{s}=\sum_{\ell=1}^{s}\hat{h}_{\ell}$. The symbol $p_{\hat{c}}(s)$ denotes a Poissonian degree distribution $p_{\hat{c}}(s)=\hat{c}^{s}\mathrm{e}^{-\hat{c}}/s!$ with mean $\hat{c}$, which naturally arises in the calculation. We note that the $\mathcal{O}(1)$ terms in $S_{1}$ and $S_{4}$ cancel, so that $S_{n}=\mathcal{O}(n)$ as expected.

The full action in terms of $\pi$ and $\hat{\pi}$ now reads

[TABLE]

The stationarity condition w.r.t $\lambda$ entails

[TABLE]

where the average $\langle\cdot\rangle_{\bar{P}}$ is taken with respect to the Gaussian measure

[TABLE]

More explicitly, (60) reads

[TABLE]

We note that the $\beta$-dependent term vanishes as $\beta\rightarrow\infty$.

The stationarity condition with respect to variations of $\pi$, $\frac{\delta S}{\delta\pi}=0$, entails the condition

[TABLE]

where $\gamma$ is a Lagrange multiplier introduced to enforce the normalisation of $\pi$. Given the definition of $Z^{(2)}_{\beta}$, (63) is equivalent to

[TABLE]

The condition that (64) must hold for all $\omega$ and $h$ can be translated into

[TABLE]

where $c=\hat{c}$ to enforce normalization of $\hat{\pi}$.

Similarly, the stationarity condition with respect to variations of $\hat{\pi}$, $\frac{\delta S}{\delta\hat{\pi}}=0,$ produces the condition

[TABLE]

where $\hat{\gamma}$ is the Lagrange multiplier enforcing the normalisation of $\hat{\pi}$. We can then conclude that the saddle-point pdf $\pi$ must satisfy

[TABLE]

Inserting (65) into (67) yields, after simple algebra

[TABLE]

where the brackets $\left\langle\cdot\right\rangle_{\{K\}_{s-1}}$ denote averaging with respect to a collection of $s-1$ i.i.d. random variables $K_{\ell}$, each drawn from the bond weight pdf $p(K)$.

We recall at this point that the replica derivation started under the simplifying assumption that the connectivity distribution was that of a standard Erdős-Rényi graph (see (30)). This implies that the degree distribution $p_{c}(s)$ - naturally appearing in (68) - is a Poisson distribution with unbounded support. However, Eq. (68) remains formally valid for any degree distribution $p_{c}(s)$ with finite mean $c$. In our case, it is then necessary to consider (2) and manually correct222Obviously, the “truncated” Eq. (69) would have been obtained anyway without any shortcuts, had we started from the exact connectivity distribution (28) at the outset. This is explicitly shown in B. (68) to account for the existence of a maximal degree, therefore yielding

[TABLE]

where $r(s)$ is the link-degree distribution (22). Note that (69) is formally identical to the self-consistent equation (21) found for the cavity field pdf, after the identification $\pi(\omega,h)\equiv q(\omega,h)$.

The constant term $\lambda\equiv\mathrm{i}\lambda^{\star}$ – which turns out to be real-valued – needs to be tuned so as to enforce (62) for $\beta\to\infty$, which reads (trading $\hat{\pi}$ for $\pi$)

[TABLE]

where – to avoid introducing more cumbersome notations – $p_{c}(s)$ now indicates the actual bounded degree distribution (2).

Surprisingly, even though the cavity and replica methods depart from completely different assumptions, they converge towards the same result: this has been already shown in [53] for the spectral problem in the Erdős-Rényi case.

A few remarks are in order:

• •

For the action to converge, we have obtained the following conditions $\omega>\zeta$, $\omega>\hat{\omega}$ and $\lambda\equiv\mathrm{i}\lambda^{\star}>\{\hat{\omega}\}_{s}$, where $\zeta$ is the upper bound of the support of the bond weights $p(K)$.

• •

Thanks to the structure of $\hat{\pi}$ (65), the entire action can be just expressed in term of $\pi$ (68) and $\lambda$ (70).

• •

The value of $\lambda\equiv\mathrm{i}\lambda^{\star}$ is real, and corresponds to the typical value of the largest eigenvalue $\langle\lambda_{1}\rangle_{J}$, as will be shown in subsection 4.1.1. This is of course again compatible with the cavity results.

• •

In Eq. (69), the contribution corresponding to $s=1$ in the sum gives rise to the term $\delta(\omega-\lambda)$ on the right hand side. Therefore, we expect to see a pronounced peak at the location of $\lambda=\langle\lambda_{1}\rangle_{J}$ in the plot of the marginal pdf $\pi(\omega)=\int\mathrm{d}h~{}\pi(\omega,h)$, once the contributions coming from nodes of different degrees are “unpacked”. This is confirmed in Fig. 4 below.

• •

Both the cavity and replica approaches can be safely extended to non-Poissonian degree distributions as well, as long as the mean connectivity $c$ remains finite as $N\rightarrow\infty$, thus considerably enlarging the class of models for which the equivalence between cavity and replicas holds true.

4.1.1 Erdős-Rényi graph: weighted adjacency matrix.

We proceed here with the case of a weighted adjacency matrix of sparse Erdős-Rényi graphs, with bounded maximal degree and bond weights drawn from the pdf $p(K)$. The pure $\{0,1\}$-adjacency matrix case is recovered considering $p(K)=\delta(K-1)$. Given the distributions (69) and (65) at stationarity and recalling (53), the $\mathcal{O}(n)$ terms of the action $S_{n}$ in (59) - keeping only the leading $\beta\rightarrow\infty$ term - are expressed as:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Multiplying and dividing the integrand of (74) by $\lambda-\sum_{\ell=1}^{s}\hat{\omega}_{\ell}$, and using (62) (for $\beta\to\infty$), we get

[TABLE]

Multiplying the second term by $1=\int\mathrm{d}\omega\mathrm{d}h\delta\left(\omega-(\lambda-\{\hat{\omega}\}_{s-1})\right)\delta(h-\{\hat{h}\}_{s-1})$, and using (67), we obtain (after some manipulations)

[TABLE]

Summing up all terms, the action at the saddle point reads

[TABLE]

which would imply from (32) for the average of the largest eigenvalue the formula

[TABLE]

However, we were able to numerically show that at the saddle point

[TABLE]

implying that

[TABLE]

as expected from the corresponding cavity calculation. The identity (79) can be more easily checked numerically once expressed in the alternative way

[TABLE]

which has the additional advantage of showing explicitly that $\lambda\equiv\mathrm{i}\lambda^{\star}$ is indeed a real-valued quantity.

The bottom panels in Fig. 2 show the marginal distributions $\pi(\omega)=\int\mathrm{d}h~{}\pi(\omega,h)$ and $\pi(h)=\int\mathrm{d}\omega~{}\pi(\omega,h)$ for the case of a pure Erdős-Rényi $\{0,1\}$-adjacency matrix, for which $p(K)=\delta(K-1)$. Figure 4 instead shows $\pi(\omega)$ and $\pi(h)$ for the case of a weighted Erdős-Rényi adjacency matrix, with a uniform bond pdf $p(K)=1/2$ for $K\in(1,3)$. In Fig. 3, we plot the behaviour of the typical largest eigenvalue as the maximum degree $k_{\mathrm{max}}$ is varied.

4.1.2 Random regular graph: adjacency matrix.

We now consider the simpler and analytically tractable case of the random regular graph (RRG). A RRG with connectivity $c$ is characterized by the property that every node has exactly $c$ neighbours, or equivalently every row of its $\{0,1\}$-adjacency matrix has exactly $c$ nonzero entries. This implies that the largest eigenvalue of such matrix is $\left\langle\lambda_{1}\right\rangle_{J}=\lambda=c$ (deterministically), and its corresponding eigenvector has all identical components $\bm{v}_{1}=(1,1,...,1)^{T}$.

In this case, the Poissonian degree distribution featuring in (68) can be safely replaced by $\delta_{s,c}$. Furthermore, if we consider the pure adjacency matrix case (i.e. with $p(K)=\delta(K-1)$), (68) and (70) become

[TABLE]

which can be exactly solved by the ansatz

[TABLE]

leading to the following equations for the parameters $\bar{\omega},\bar{h}$ and $\lambda$

[TABLE]

Eq. (86) entails that $\bar{\omega}=c-1$. Then, inserting this value in (85), we find $\lambda=c$.

The value of $\bar{h}$ can then be found exploiting the normalization condition (87), yielding $\bar{h}=c-2$.

The action at the saddle-point reads then

[TABLE]

and therefore, the typical largest eigenvalue is

[TABLE]

equal to $\lambda$ as expected.

4.2 Density of the top eigenvector’s components

In this statistical mechanics framework, the quantity

[TABLE]

is defined such that in the limit $\beta\rightarrow\infty$ it gives the density of the top eigenvector components for a given $N\times N$ sparse symmetric random matrix $J$. The simple angle brackets $\left\langle...\right\rangle$ stands for thermal averaging, i.e. with respect to the Gibbs-Boltzmann distribution (6) of the system

[TABLE]

Defining an auxiliary partition function as

[TABLE]

where $\delta_{\epsilon}$ is a smooth regulariser of the delta function, the quantity (90) can be formally expressed as

[TABLE]

Averaging now over the matrix ensemble

[TABLE]

and sending $\beta\to\infty$ at the very end, the density of the top eigenvector’s components is eventually given by the remarkable formula

[TABLE]

To compute the average of the logarithm of the auxiliary partition function $Z^{(\beta)}_{\epsilon}(t,J;u)$, we will employ the replica trick once again

[TABLE]

We can anticipate that the replicated partition function will take the form

[TABLE]

where $\varphi$ and $\hat{\varphi}$ are functional order parameters. In a saddle point approximation for large $N$

[TABLE]

where the starred objects satisfy self-consistency equations in which $t$ can be safely set to zero: indeed, the partial derivative $\frac{\partial}{\partial t}$ in (95) only acts on terms containing any explicit dependence on $t$, and not through any other indirect functional dependence. Inserting (98) into (96), and assuming that

[TABLE]

as $n\to 0$ (in a replica-symmetric setting), the final expression for the average density of top eigenvector’s components for $N\to\infty$ reduces to

[TABLE]

where $(\cdot)^{\prime}$ stands for differentiation with respect to $t$.

Interestingly, we will find that the stationarity conditions defining $\varphi^{\star}$, $\hat{\varphi}^{\star}$ and $\lambda^{\star}$ at the saddle point for $t=0$ are just identical to those found in the replica calculation for the largest eigenvalue. The explicit $n$-dependence of the action $S^{(\beta)}_{n}\left(\varphi^{\star},\hat{\varphi}^{\star},\vec{\lambda}^{\star};t,\epsilon;u\right)$ is extracted by representing the order parameters $\varphi$ and $\hat{\varphi}$ as an infinite superposition of Gaussians, as previously done for the leading eigenvalue calculation.

In the next subsections, we will apply this formalism to weighted Erdős-Rényi and random regular adjacency matrices.

4.2.1 Erdős-Rényi graph: weighted adjacency matrix.

The average replicated partition function becomes

[TABLE]

in complete analogy with (36).

By introducing the functional order parameter

[TABLE]

via the usual integral identity

[TABLE]

the replicated partition function can be once again cast in a form that allows for a saddle point approximation

[TABLE]

where the action $S^{(\beta)}_{n}\left[\varphi,\hat{\varphi},\vec{\lambda};t,\epsilon;u\right]$ is the sum of four terms

[TABLE]

where for simplicity we omit the full dependence on variables on the right hand side. The first three contributions are identical to those appearing in the largest eigenvalue calculation (see (43), (44) and (45)). The explicit $t$ and $u$ dependence is confined to the fourth contribution,

[TABLE]

The saddle point equations for $\varphi^{\star},\hat{\varphi}^{\star}$ (where we can safely set $t=0$) are then identical to those (see (47) and (48)) appearing in the calculation for the average largest eigenvalue. Therefore we can follow the same strategy as before, and represent $\varphi^{\star}$ and $\hat{\varphi}^{\star}$ as uncountably infinite superposition of Gaussians, whose parameters fluctuate according to joint pdfs $\pi$ and $\hat{\pi}$ as in (51) and (52). Such joint pdfs satisfy the very same self-consistency equations as in (68) and (65) and for these reasons we can use the same labels as before. The only difference with respect to the previous case is in the extra $t$-derivative that we have to take from the contribution $S_{4}(\hat{\varphi}^{\star},\lambda;t,\epsilon;u)$.

Inserting the ansatz

[TABLE]

into (106), and expanding $\mathrm{e}^{\mathrm{i}\hat{\varphi}^{\star}(\vec{v})}=\sum_{s\geq 0}\frac{(\mathrm{i}\hat{\varphi}^{\star}(\vec{v}))^{s}}{s!}$, we obtain (in the limit $n\to 0$)

[TABLE]

Therefore, we can isolate the function $s_{\beta}(t,\epsilon;u)$ in (99) as

[TABLE]

in view of the identifications $\hat{c}=c$ and $\mathrm{i}\lambda^{\star}\equiv\lambda$ as before. Taking the $t$-derivative and setting $t$ and $\epsilon$ to zero, we get

[TABLE]

Taking the $\beta\to\infty$ limit as in (100), we eventually find

[TABLE]

Expressing everything in terms of the $\pi$-distribution, indicating with $p_{c}(s)$ the actual degree distribution (2) and truncating the series at the largest degree $k_{\mathrm{max}}$ (as we did in previous sections), we eventually obtain

[TABLE]

where $\pi(\omega,h)$ satisfies the self-consistent equation (69) (to be solved via population dynamics), supplemented with the normalisation condition (70). Once again, the brackets $\left\langle\cdot\right\rangle_{\{K\}_{s}}$ denote averaging w.r.t to a collection of $s$ i.i.d random variables $K_{\ell}$, each drawn from the bond weight pdf $p(K)$.

Eq. (111) essentially recovers Eq. (25) found with the cavity method. As a general remark, it is worth noticing that the $\beta$-dependent distribution $\varrho_{\beta}(u)$ had already arisen naturally in the eigenvalue calculation when evaluating the stationarity conditions with respect to $\lambda$. In fact, the distribution in (61) is exactly identical to $\varrho_{\beta}(u)$. Moreover, in the cavity formalism, $\varrho_{\beta}\left(u\right)$ is closely related to the single-site marginal of a single instance (16).

We remark once again that – in analogy with the typical largest eigenvalue calculation – the validity of Eq. (111) is not restricted to a truncated Poisson degree distribution (2). It actually provides the density of the top eigenvector’s components for the weighted adjacency matrix of any configuration model with finite connectivity and bounded maximal degree as a weighted superposition of delta functions, one for each degree of the graph. It is then natural to identify the quantity $\frac{\sum_{\ell=1}^{s}\frac{h_{\ell}K_{\ell}}{\omega_{\ell}}}{\lambda-\sum_{\ell=1}^{s}\frac{K_{\ell}^{2}}{\omega_{\ell}}}$ as the contribution to the density coming from nodes of degree $s$.

The $s=0$ contribution from isolated nodes indeed gives rise to the sharp peak at $u=0$. The $\varrho(u)$ of a Erdős-Rényi $\{0,1\}$-adjacency matrix is shown in Figure 2 (top panels), whereas the case of weighted Erdős-Rényi adjacency matrices is shown in Figure 5.

4.2.2 Random regular graph: adjacency matrix

In this case, building on subsection 4.1.2 and recalling that $p_{c}(s)=\delta_{s,c}$ and $p(K)=\delta(K-1)$, the ratio in (111) simply becomes $c(c-2)/[c(c-1)-c]=1$, entailing

[TABLE]

as expected.

4.2.3 Large-$c$ limit for weighted adjacency matrices

We consider now the large $c$-limit of Erdős-Rényi graphs (more generally, any configuration model graph for which $\frac{\sigma_{k}^{2}}{\langle k\rangle^{2}}=\frac{\langle k^{2}\rangle-\langle k\rangle^{2}}{\langle k\rangle^{2}}\to 0$ as $\langle k\rangle=c\to\infty$). A meaningful large-$c$ limit is obtained for Eq. (68) by rescaling each instance of the bond random weights as $K_{ij}=\mathcal{J}_{ij}/\sqrt{c}$, leading to

[TABLE]

In the $c\gg 1$ limit, the $s$-sum in Eq. (113) is restricted to $s=c\pm\mathcal{O}(\sigma_{k})$ (with $\sigma_{k}=\sqrt{c}$ for Erdős-Rényi graphs), so that the argument appearing in the first $\delta$-function on the r.h.s of this equation can be evaluated by appeal to the Law of Large Number (LLN). This entails that

[TABLE]

is non-fluctuating, hence the self-consistency equation demands that

[TABLE]

with (by the LLN)

[TABLE]

Specializing to $\langle\mathcal{J}^{2}\rangle_{\mathcal{J}}=1$, we see that

[TABLE]

which requires $\lambda\geq 2$ to have real positive $\bar{\omega}$.

Similarly, the argument of the second $\delta$-function on the r.h.s of (113) exhibits a scaling that allows one to conclude (for $\langle\mathcal{J}_{\ell}\rangle_{\mathcal{J}}=0$) that

[TABLE]

by appeal to the Central Limit Theorem. The variance follows using independence of the $\{h_{\ell}\}$ and $\{\mathcal{J}_{\ell}\}$

[TABLE]

This equation allows a finite variance if and only if $\bar{\omega}^{2}=1$, which requires $\lambda=\pm 2$, i.e. that $\lambda$ – the most probable location of the largest eigenvalue – is at the edge of the Wigner semi-circle (and we require the positive solution).

To obtain the distribution $\varrho(u)$ of eigenvector components, it is instructive and more direct to look back at the cavity equations (24), (25) and (26). After the rescaling $K_{\ell}=\mathcal{J}_{\ell}/\sqrt{c}$ and in the large $c$-limit, it is easy to see from (24) that $\Omega=\bar{\omega}$ and that $H$ is a sum of Gaussians, and thus itself Gaussian, of variance $\sigma_{h}^{2}/\bar{\omega}^{2}\equiv\sigma_{h}^{2}$ by (118). It then follows from the normalisation condition (26) that $\sigma_{h}^{2}=1$, so eventually

[TABLE]

Looking now at the variable $\eta=u^{2}$, and noting that positive and negative $u$ give rise to the same $\eta$, one obtains by the simple transformation of pdf’s

[TABLE]

which is the standard form of Porter-Thomas distribution for real-valued (invariant) random matrices (see [23], Eq. (9.10)).

5 Application: sparse random Markov transition matrices

In this section, we cross-check the formalism with an ensemble of transition matrices $W$ for discrete Markov chains in an $N$-dimensional state space. The evolution equation for the probability vector $\bm{p}(t)$ is given by

[TABLE]

The transition matrix $W$ is such that $W_{ij}\geq 0\ \forall(i,j)$ and $\sum_{i}W_{ij}=1\ \forall j$. For an irreducible chain, the top right eigenvector of the matrix $W$ corresponding to the Perron-Frobenius eigenvalue $\lambda_{1}=1$ represents the unique equilibrium distribution, i.e. $\bm{v}_{1}=\bm{p}^{\mathrm{eq}}$. The matrix $W$ is in general not symmetric: however, if the Markov process satisfies a detailed balance condition, i.e. $W_{ij}p_{j}^{\mathrm{eq}}=W_{ji}p_{i}^{\mathrm{eq}}$, it can be symmetrised via a similarity transformation, yielding

[TABLE]

The symmetrised matrix $W^{S}$ will be the target of our analysis: even though it is not itself a Markov matrix since the columns normalisation constraint is lost, in view of the detailed balance condition $W^{S}$ has the same (real) spectrum of $W$, and its top eigenvector $\bm{v}_{1}$ is given in terms of the top right eigenvector of $W$, $\bm{p}_{\mathrm{eq}}$, as

[TABLE]

We will consider the case of an unbiased random walk: the matrix $W$ is then defined as

[TABLE]

where $c_{ij}$ represents the connectivity matrix and $k_{j}=\sum_{i}c_{ij}$ is the degree of the node $j$. In this case, the top right eigenvector of $W$ is proportional to the vector expressing the degree sequence: for our purposes, we choose the inverse of the mean degree as proportionality constant, i.e. $p_{i}^{\mathrm{eq}}=k_{i}/\left\langle k\right\rangle$. The symmetrised matrix $W^{S}$ is expressed as

[TABLE]

with its top eigenvector being $v_{1}^{(i)}=\sqrt{k_{i}/\left\langle k\right\rangle}$. Therefore, we expect that

[TABLE]

where $p(k)$ is the degree distribution of the connectivity matrix $\{c_{ij}\}$.

In order to avoid isolated nodes and isolated clusters of nodes, we consider a shifted Poissonian degree distribution with $k_{\mathrm{min}}=2$, i.e.

[TABLE]

with mean degree $\left\langle k\right\rangle=c+2$.

The single-instance cavity treatment starts from the Gibbs-Boltzmann distribution

[TABLE]

which, after the change of variable $\tilde{v_{i}}=v_{i}/\sqrt{k_{i}}$, becomes

[TABLE]

It is convenient to frame and solve the problem in terms of the vector $\tilde{\bm{v}}$, since in this case the matrix involved in the analysis is just the standard $\{0,1\}$-adjacency matrix of the underlying graph, as in [58, 59]. The cavity single-instance equations for this problem read

[TABLE]

whereas the equations for the single-site marginal coefficients read

[TABLE]

In the thermodynamic limit $N\to\infty$, the equations (129) and (130) lead to

[TABLE]

in complete analogy with (21).

Similarly, equations (131) and (132) lead to

[TABLE]

entailing that the density of the top eigenvector’s components in the space of vectors $\tilde{\bm{v}}$ is given by

[TABLE]

which follows from the general theory.

As before, (133) and (134) are efficiently solved via a population dynamics algorithm: as expected, the convergence is attained for $\lambda=1$, i.e. in correspondence of the largest eigenvalue of $W^{S}$. Running the simulations, we find that the distribution $\varrho(\tilde{u})$ converges to a delta peak centered at a finite real positive value: this behaviour agrees perfectly with the theoretical predictions, because it precisely implies that $\varrho(u)$ must be given by (125). Indeed, the two quantities are related via the aforementioned change of variables, $u\leftarrow\tilde{u}\sqrt{k}$, and the constant value the variables $\tilde{u}$ converge to corresponds to $1/\sqrt{\left\langle k\right\rangle}$, once the normalisation is fixed according to (26). In Fig. 6, we compare the density of the top eigenvector’s components for sparse Markov matrices (representing the transition matrices of unbiased random walks) with numerical diagonalisation.

As a concluding remark, we notice that the same route can be followed to characterise the top eigenpair statistics of the so-called tilted Markov transition matrix [73] appearing in the context of rare events for random walks on networks [74]. This will be discussed in a separate publication.

6 Population dynamics

The population dynamics algorithm employed to solve (69) is deeply rooted in the statistical mechanics of spin glasses [75, 76]. In our context, it can be summarised as follows.

Two coupled populations with $N_{P}$ members each $\left\{\left(\omega_{i},h_{i}\right)\right\}_{1\leq i\leq N_{P}}$ are randomly initialised, taking into account that $\omega_{i}>\zeta$, where $\zeta$ is the upper edge of the support of the pdf $p(K)$.

For any suitable value of $\mathrm{i}\lambda^{*}\equiv\lambda\in\mathbb{R}$, the following steps are iterated until stable populations are obtained:

Generate a random $s\sim\frac{s}{c}p_{c}\left(s\right)$, where $c=\left\langle s\right\rangle$ 2. 2.

Generate $s-1$ i.i.d random variables $K_{\ell}$ from the bond weights pdf $p(K)$ 3. 3.

Select $s-1$ pairs $\left(\omega_{\ell},h_{\ell}\right)$ from the population at random; compute

[TABLE]

and replace a randomly selected pair $\left(\omega_{j},h_{j}\right)$ where $j=1,...,N_{P}$ with the pair $\left(\omega^{(new)},h^{(new)}\right)$. 4. 4.

Return to (i).

Convergence is assessed by looking at the first moments of the vector formed by the $N_{P}$ samples. A sweep is completed when all the $N_{P}$ pairs $(\omega_{j},h_{j})$ of the population have been updated at least once according to the steps above.

The procedure to evaluate (24) (or alternatively (25)) is almost identical, except for the details concerning the $s$-sampling. Starting from two coupled populations with $N_{P}$ members $\{(\Omega_{i},H_{i})\}_{1\leq i\leq N_{P}}$, the following steps are iterated:

Generate a random $s\sim p_{c}\left(s\right)$, where $c=\left\langle s\right\rangle$ 2. 2.

Generate $s$ i.i.d random variables $K_{\ell}$ from the bond weights pdf $p(K)$ 3. 3.

Select $s$ pairs $\left(\omega_{\ell},h_{\ell}\right)$ from the population $\left\{\left(\omega_{i},h_{i}\right)\right\}_{1\leq i\leq N_{P}}$ at random; compute

[TABLE] 4. 4.

Replace a randomly selected pair $\left(\Omega_{j},H_{j}\right)$ where $j=1,...,N_{P}$ with the pair $\left(\Omega^{(new)},H^{(new)}\right)$, which is then a new sample from $Q(\Omega,H)$. It can be used via Eq. (25) to create $u^{(new)}=H^{(new)}/\Omega^{(new)}$ as a new sample from $\varrho(u)$. 5. 5.

Return to (i).

The value of the parameter $\lambda$ controls the convergence of the algorithm: indeed, the convergence to a non-trivial distribution is achieved only when $\lambda$ is equal to the typical largest eigenvalue $\left\langle\lambda_{1}\right\rangle_{J}$, as prescribed by the theory: for any $\lambda>\left\langle\lambda_{1}\right\rangle_{J}$, the variables of type $h$ will shrink to zero, whereas for $\lambda<\left\langle\lambda_{1}\right\rangle_{J}$ they will blow up in norm. Hence, the value $\lambda=\left\langle\lambda_{1}\right\rangle_{J}$ is the only value for which the normalisation condition (70) (or equivalently (26)) can be satisfied, in complete agreement with the replica predictions.

In view of the expected behaviour described above, we will initially start from a large value of $\lambda$, which is then progressively decreased until convergence is achieved. A suitable starting value is given by the largest degree $k_{\mathrm{max}}$ that appears in the connectivity distribution. The value of $k_{\mathrm{max}}$ is fixed in such a way that $p_{c}(k_{\mathrm{max}})N_{P}\geq 1$: only if this condition is met, the value $k_{\mathrm{max}}$ appears at least once in the degree array that is created to sample from $p_{c}(k)$. Because of this choice, the largest degree depends on the limits of the machine that is used to run the population dynamics algorithm: by using a population size $N_{P}=10^{6}$ and a parameter $\bar{c}=4$ in (2), we are able to reach $k_{\mathrm{max}}=16$. Thus, the normalisation constant $\mathcal{N}$ in (2) is not very different from $1$ and $\bar{c}\simeq c=\langle s\rangle$, making the truncation of the Poisson distribution - for all practical purposes - ineffective.

Once $\lambda$ has been set to the only value ($=\left\langle\lambda_{1}\right\rangle_{J}$) for which a non-trivial finite normalisation can be found, the value of such normalisation can be adjusted by properly rescaling the $h$’s. Such rescaling is always allowed due to the linear nature of the recursion that governs their update. This recursion will be discussed in detail in A.

The population dynamics algorithm can also be employed to evaluate numerically the integral in (81). The integral has the following structure:

[TABLE]

where $f$ is a generic function of the cavity fields and $K$. Once the correct value of $\lambda=\left\langle\lambda_{1}\right\rangle_{J}$ has been found, a number $N_{\mathrm{eq}}$ of equilibration sweeps is performed, following the protocol illustrated above.

After equilibration, a variable $F=0$ is initialised. Then for $j=1,\ldots,N_{\mathrm{meas}}$:

Perform a sweep 2. 2.

Pick $(\omega,h)$ and $(\omega^{\prime},h^{\prime})$ at random, generate $K\sim p(K)$ and compute $f\left(\omega,h,\omega^{\prime},h^{\prime},K\right)$. 3. 3.

Update $F$: $F=F+f\left(\omega,h,\omega^{\prime},h^{\prime},K\right)$.

The value of the integral (140) is approximated by invoking the law of large numbers, as

[TABLE]

where the typical fluctuation is of the order of $1/\sqrt{N_{P}N_{\mathrm{meas}}}$.

7 Conclusions

In summary, we have further developed a formalism - pioneered by Kabashima and collaborators - to compute the statistics of the largest eigenvalue and of the corresponding top eigenvector for some ensembles of sparse symmetric matrices, i.e. (weighted) adjacency matrices of graphs with finite mean connectivity. The top eigenpair problem can be recast as the optimisation of a quadratic Hamiltonian on the sphere: introducing the associated Gibbs-Boltzmann distribution and a fictitious inverse temperature $\beta$, the top eigenvector represents the ground state of the system, which is attained in the limit $\beta\to\infty$. In order to extract this limit, we have employed two methods, cavity and replicas, both borrowed from the Statistical Mechanics approach to disordered systems. We first analysed the case of a single-instance matrix within a “grand canonical” cavity framework. The single-instance cavity method leads fairly quickly to superficially appealing recursion equations, however it has the obvious drawback of enlarging - and not reducing - the complexity of the problem: indeed, it turns a $N$-dimensional problem involving the single matrix $J$ into an $Nc$ dimensional problem - where $c=\left\langle k\right\rangle>1$ is the mean degree - involving the non-backtracking operator $B$, as detailed in A.

However, the cavity single-instance recursions constitute an essential ingredient to arrive at the equations (21), (25) and (26) for the associated joint probability densities of the auxiliary fields of type $\Omega$ and $H$ that characterise the typical largest eigenvalue and the statistic of the top eigenvector in the thermodynamic limit $N\to\infty$. Moreover, the exact same equations (see (69), (111) and (70)) are found via the completely alternative replica derivation, entailing that the two methods are equivalent in the thermodynamic limit. Within the population dynamics algorithm employed to solve the stochastic recursion (21) (or equivalently (69)), we are able to identify the typical largest eigenvalue as the parameter controlling the convergence of the algorithm, and unpack the contributions coming to nodes of different degrees to the average density of the top eigenvector’s components. The simulations show excellent agreement of the theory with the direct diagonalisation of large matrices. As a further cross-check of the formalism, we computed the average density of the top eigenvector’s components of sparse Markov matrices representing unbiased random walks on a sparse network under the detailed balance condition, thus retrieving the expected relation between such components and the node degrees of the underlying network.

The authors acknowledge funding by the Engineering and Physical Sciences Research Council (EPSRC) through the Centre for Doctoral Training in Cross Disciplinary Approaches to Non-Equilibrium Systems (CANES, Grant Nr. EP/L015854/1).

Appendix A

The solution of the single instance self-consistency equations and the non-backtracking operator.

The set of self-consistency equations (14) and (15) for the cavity fields, supplemented with (17) and (18) for the coefficients of the marginal distributions, constitutes the full solution of the top eigenpair problem for a single instance of a sparse matrix. Even in this case, the convergence of the update equations (14) and (15) is dictated by the value of the parameter $\lambda$, which once again is related to the possibility to normalise the resulting top eigenvector.

Note that (15) is a linear recursion driven by the operator $B$, whose elements can be defined as

[TABLE]

$B$ is an example of non-backtracking operator, first introduced by Hashimoto in [77]. For a given graph, the Hashimoto non-backtracking operator $\tilde{B}$ in its original form counts the number of paths from a node $i$ to a node $\ell$ passing through a third node $j$, for every choice of these three different nodes. It is defined as

[TABLE]

In our case, if the absolute value of the largest eigenvalue of the modified non-backtracking operator $B$ is greater than $1$, the absolute values of the cavity fields $H_{j}^{(i)}$’s will blow up, whereas if it smaller than $1$, they will shrink to zero. Therefore, $\lambda$ must be tuned appropriately in (14) to prevent the linear recursion (15) from landing on a trivial solution. Indeed, when $\lambda$ is “too large”, the $\Omega_{j}^{(i)}$’s will be large too, resulting in a largest eigenvalue of $B$ with magnitude smaller than $1$. This would suggest to progressively decrease $\lambda$ from a large value down to its lower bound $\lambda_{1}$, necessary to ensure that the optimisation problem is well-defined. In other words, the largest eigenvalue of the operator $B$ must be exactly $1$ for the $H_{j}^{(i)}$’s to have a finite norm. This will happen only when $\lambda=\lambda_{1}$.

Collecting the $H_{j}^{(i)}$’s in a $2M=\sum_{i=1}^{N}k_{i}$ dimensional vector, Eq. (15) can be rewritten as a linear vector iteration driven by $B$ as

[TABLE]

where the entries $B_{(i,j),(k,\ell)}$ are defined in (142). Relabelling with a new, single index $a$ any pair of connected indices $(i,j)$, (144) reads

[TABLE]

which can interpreted as a vector linear iteration,

[TABLE]

with the index $t$ labelling each iteration.

Starting from a certain initial condition $\underline{H}_{0}$, the solution of (15) is obtained after successive iterations according (146) until $\underline{H}_{t}$ stabilises. The stability can be assessed by looking at the norm of the vector $\underline{H}_{t}$. After a suitable number of iterations $t$, expanding the initial condition vector in the basis $\left\{\bm{b}_{i}\right\}$ formed by the right eigenvectors of $B$, the leading contribution is expressed in terms of its top eigenpair

[TABLE]

where the contributions coming from the other eigenpairs $\left\{\bm{b}_{i},\gamma_{i}\right\}$ are exponentially suppressed, all the other eigenvalues of $B$ being smaller than $\gamma_{1}$.

The ratio $\eta_{t}$ of the norms of two successive iterations approaches a constant value $\eta^{\star}$ as $t\rightarrow\infty$, corresponding to the absolute value of largest eigenvalue of $B$,

[TABLE]

Thus, the convergence of (15) is attained when the value of $\eta^{\star}=\left|\gamma_{1}\right|$ reaches $1$ as $\lambda$ approaches $\lambda_{1}$ from above. We again recall that $\lambda=\lambda_{1}$ is the smallest possible value such that the cavity partition function (8) is well defined, and so the actual value $\lambda_{1}$ can be found by asymptotic extrapolation. Figure 7 shows an example of this procedure.

We remark that the procedure above holds only if the largest eigenvalue of $B$ is real: if it is complex, there will be a pair of complex conjugate first eigenvalues, i.e. those with the largest norm, which dictate the asymptotic behavior of (146). In this case, the bi-orthonormal basis of left and right eigenvectors must be taken into account

[TABLE]

where the coefficients $c_{1}(0)$ and $c_{2}(0)$ are in general complex. Therefore, the quantity $\eta_{t}$ does not approach a steady limit for large $t$ in this case, and oscillations arise. In fact, it can be shown that

[TABLE]

where

[TABLE]

Here, ($|c_{1}|,|c_{2}|)$ and $(\phi_{1},\phi_{2})$ are the moduli and phases of the complex coefficients $c_{1}(0)$ and $c_{2}(0)$, $\eta$ is the ratio of the radial part of the vectors $\underline{H}_{t}$ and $\underline{H}_{t-1}$, $\alpha$ and $\theta$ are respectively the modulus and the phase of the dot product between the right (and left) eigenvector $\bm{b}_{1}$ (respectively $\bm{b}_{1}^{\star}$) with itself, and $\varrho$ and $\phi$ are the modulus and phase of the pair of the complex eigenvalues with the largest norm.

In this case, the recursion (146) does not converge to a single limit, and the cavity formalism does not lead to an acceptable solution. Therefore, the strongest limitation of the single instance cavity method is that the largest eigenvalue $\gamma_{1}$ of the non-backtracking operator $B$ associated to the matrix $J$ must be real. This restriction unfortunately rules out a variety of interesting sparse matrix ensembles.

Appendix B

Exact replica calculation for the largest eigenvalue for any bounded degree distribution $p(k)$.

In this appendix, we show how to get the typical largest eigenvalue with the replica method without any shortcut in the calculation. We will thus employ the distributions (27) and (28) to perform the averaging w.r.t the matrix ensemble. We recall that the parameter $c$ appearing in (28) stands for the actual mean of the bounded degree distribution of interest, which may in general differ from the parameter $\bar{c}$ of the truncated Poisson distribution (see (2)). They tend to coincide only if $k_{\mathrm{max}}$ is large. This procedure is general and holds for any graph within the configurational model with degree sequence originated by a finite-mean degree distribution $p(k)$.

Following the same reasoning in Section 4.1, the replicated partition function is given by (33). Taking the average w.r.t the joint distribution (28) of matrix entries yields [61]

[TABLE]

where $\left\langle\cdot\right\rangle_{K}$ denotes averaging over the single-variable pdf $p\left(K\right)$ characterising the i.i.d. bond weights $K_{ij}$. A Fourier representation of the Kronecker deltas expressing the degree constraints in (28) has been employed. As in Section 4.1, we also employ a Fourier representation of the Dirac delta enforcing the normalisation constraint. The replicated partition function thus becomes

[TABLE]

In order to decouple sites, we introduce the functional order parameter

[TABLE]

where the symbol $\vec{v}$ denotes a $n$-dimensional vector in replica space. We then consider its integrated version [61]

[TABLE]

and enforce the latter definition using the integral identity

[TABLE]

In terms of the integrated order parameter (156) and its conjugate, the replicated partition function can be written as

[TABLE]

The multiple integral in the last line above is the product of $N$ $n$-dimensional integrals, each related to a degree $k_{i}$. It can be written as

[TABLE]

where $\mathrm{Log}$ denotes the principal branch of the complex logarithm, and

[TABLE]

Each of the $\phi_{i}$ integrals can be performed by rewriting the last exponential factor as a power series, viz.

[TABLE]

with $i=1,\ldots,N$. Thus, by invoking the Law of Large Numbers, the single site integral $I$ (159) can be expressed as

[TABLE]

where we have used

[TABLE]

where $p(k)$ is the actual degree distribution of the graph.

As in Section 4.1, the replicated partition function takes a form amenable to a saddle point evaluation for large $N$

[TABLE]

where

[TABLE]

The terms $S_{1},S_{2}$ and $S_{3}$ are equal to those found in Section 4.1, respectively (43), (44) and (45), whereas

[TABLE]

As in Section 4.1, we then search for replica-symmetric saddle-point solutions written in the form of superpositions of uncountably infinite Gaussians with a non-zero mean,

[TABLE]

where

[TABLE]

and – with a modest amount of foresight – we use the same notation as before for the distributions $\pi$ and $\hat{\pi}$. The $\varrho_{0}$ and $\hat{\varrho}_{0}$ are determined such that the distributions $\pi(\omega,h)$ and $\hat{\pi}(\hat{\omega},\hat{h})$ are normalised. The $\varrho_{0}$ in (168) is needed since $\varrho^{\star}(\vec{v})$ is the saddle-point expression of the integrated order parameter.

Rewriting the action in terms of $\pi$ and $\hat{\pi}$, after performing the $\vec{v}$-integrations, and extracting the leading $n\to 0$ contribution yields

[TABLE]

with

[TABLE]

where we have taken into account that $k_{\mathrm{min}}=0$ and we have introduced the shorthands

[TABLE]

and $\{\mathrm{d}\hat{\pi}\}_{s}=\prod_{\ell=1}^{s}\mathrm{d}\hat{\omega}_{\ell}\mathrm{d}\hat{h}_{\ell}\hat{\pi}(\hat{\omega}_{\ell},\hat{h}_{\ell})$, along with $\{\hat{\omega}\}_{s}=\sum_{\ell=1}^{s}\hat{\omega}_{\ell}$ and $\{\hat{h}\}_{s}=\sum_{\ell=1}^{s}\hat{h}_{\ell}$.

We note that the action contains $\mathcal{O}(1)$ and $\mathcal{O}(n)$ terms as $n\to 0$: the $\mathcal{O}(1)$ terms are cancelled by the $\mathcal{O}(1)$ terms arising from the evaluation of the normalisation constant $\mathcal{M}$ at the saddle-point. Indeed, by following a very similar reasoning as in (153), we find that

[TABLE]

We then introduce in (177) the scalar order parameter

[TABLE]

via the integral representation

[TABLE]

By using the same argument as in (162), the normalisation constant $\mathcal{M}$ can be written in a form amenable to a saddle-point evaluation,

[TABLE]

The stationarity conditions for $S_{\mathcal{M}}$ are

[TABLE]

and

[TABLE]

entailing that

[TABLE]

The two conditions above exhibit a gauge invariance[61]. Once the same gauge has been chosen for the saddle-point solution of $\mathcal{M}$ and the $\mathcal{O}(1)$ terms of the action (59) in the numerator, they cancel out so that the action (59) is $\mathcal{O}(n)$ as expected.

Thus, taking into account the cancellation coming from (183) and (184), the action terms in (59) read exactly as those found in Section 4.1, thus proving that the “shortcut” derivation in 4.1 is perfectly legitimate. According to the present derivation, the degree distribution $p(k)$ appearing in the single-site term $S_{4}$ is already the true degree distribution of the graph, and does not require any a posteriori correction.

—————–