Imaginary replica analysis of loopy regular random graphs

TL;DR

This paper develops an analytical method using imaginary replica indices to analyze the eigenvalue spectra of regular loopy graphs under spectral constraints, revealing phase transitions and validating results with simulations.

Contribution

Introduces a novel imaginary replica approach for spectral analysis of loopy graphs, enabling the study of phase transitions in constrained random graph ensembles.

Findings

Derived expressions for eigenvalue spectra in weak loop interactions

Identified shattering phase transitions in graph structures

Validated analytical predictions with MCMC simulations

Abstract

We present an analytical approach for describing spectrally constrained maximum entropy ensembles of finitely connected regular loopy graphs, valid in the regime of weak loop-loop interactions. We derive an expression for the leading two orders of the expected eigenvalue spectrum, through the use of infinitely many replica indices taking imaginary values. We apply the method to models in which the spectral constraint reduces to a soft constraint on the number of triangles, which exhibit `shattering' transitions to phases with extensively many disconnected cliques, to models with controlled numbers of triangles and squares, and to models where the spectral constraint reduces to a count of the number of adjacency matrix eigenvalues in a given interval. Our predictions are supported by MCMC simulations based on edge swaps with nontrivial acceptance probabilities.

Equations

$$\text{Tr}({\mathbf{A}}^{\ell })$$

$$\varrho (\mu |\mathbf{A})=\frac{1}{N}\sum _{i=1}^N\delta (\mu -{\mu }_i(\mathbf{A})),$$

$$p(\mathbf{A})$$

$$Z[\hat{\varrho }]$$

$$\varrho (\mu )$$

$$p(\mathbf{A})$$

$$\displaystyle p_{\rm ER}(\mathbf{A})=\Big{(}\frac{q}{N}\Big{)}^{\!\frac{Nq}{2}}\Big{(}\frac{N-q}{N}\Big{)}^{\!\frac{N}{2}(N-1-q)}~{}~{}~{}{\rm for~{}all~{}}\mathbf{A}\!\in\!G~{}~{}{\rm with}~{}~{}\prod_{i=1}^{N}\delta_{q,\sum_{j}A_{ij}}=1.$$

$$\varphi [\hat{\varrho }]$$

$${\text{rme}}^{N\int \text{rmd}\mu \hat{\varrho }(\mu )\varrho (\mu |\mathbf{A})}$$

$$Z({\mu }_{\epsilon }|\mathbf{A})$$

$$\displaystyle\prod_{\mu}Z({\mu_{\varepsilon}}|\mathbf{A})^{n_{\mu}}\overline{Z({\mu_{\varepsilon}}|\mathbf{A})}^{m_{\mu}}=\int\!\rmd[\boldsymbol{\Phi},\boldsymbol{\Psi}]\exp\Big{(}\!-\frac{\rmi}{2}\sum_{{\mu,\alpha_{\mu}}}\sum_{ij}\phi^{i}_{\mu,\alpha_{\mu}}\phi^{j}_{\mu,\alpha_{\mu}}(A_{ij}\!-\!{\mu_{\varepsilon}}\delta_{ij})\Big{)}$$

$$\displaystyle\hskip 71.13188pt\times\exp\Big{(}\frac{\rmi}{2}\sum_{\mu,\beta_{\mu}}\sum_{ij}\psi^{i}_{\mu,\beta_{\mu}}\psi^{j}_{\mu,\beta_{\mu}}(A_{ij}\!-\overline{\mu_{\varepsilon}}\delta_{ij})\Big{)},$$

$$\displaystyle\boldsymbol{\varphi}=\left(\begin{array}[]{c}\{\phi_{\mu,\alpha_{\mu}}\}\\ \{\psi_{\mu,\beta_{\mu}}\}\end{array}\right),$$

$$\mathbf{\phi }\cdot {\mathbf{\phi }}^{\prime }=\sum _{\mu ,{\alpha }_{\mu }}{\varphi }_{\mu ,{\alpha }_{\mu }}{\varphi }_{\mu ,{\alpha }_{\mu }}^{\prime }-\sum _{\mu ,{\beta }_{\mu }}{\psi }_{\mu ,{\beta }_{\mu }}{\psi }_{\mu ,{\beta }_{\mu }}^{\prime },$$

$$\displaystyle\boldsymbol{\mu}=\left(\begin{array}[]{cccccc}{\mu_{1}}_{\varepsilon}\mathbf{I}_{n_{\mu_{1}}}&\mathbf{0}&&\dots&&\mathbf{0}\\ \mathbf{0}&\ddots&\mathbf{0}&&&\\ &\mathbf{0}&{\mu_{M}}_{\varepsilon}\mathbf{I}_{n_{\mu_{M}}}&\mathbf{0}&&\vdots\\ \vdots&&\mathbf{0}&\overline{{\mu_{1}}_{\varepsilon}}\mathbf{I}_{m_{\mu_{1}}}&\mathbf{0}&\\ &&&\mathbf{0}&\ddots&\mathbf{0}\\ \mathbf{0}&&\dots&&\mathbf{0}&\overline{{\mu_{M}}_{\varepsilon}}\mathbf{I}_{m_{\mu_{M}}}\end{array}\right).~{}$$

$$\nu (\mathbf{\phi })={\text{rme}}^{\frac{1}{2}\text{rmi}\mathbf{\phi }\cdot \mathbf{\mu }\mathbf{\phi }},\lim =\underset{\Delta \to 0}{\lim }\underset{n_{\mu }\to \text{rmi}\frac{\Delta }{\pi }{\hat{\varrho }}^{\prime }(\mu )}{\lim }\underset{m_{\mu }\to -n_{\mu }}{\lim }.$$

$$\varphi [\hat{\varrho }]$$

$$P(\mathbf{\phi },\omega )=\frac{1}{N}\sum _{i=1}^N\delta (\mathbf{\phi }-{\mathbf{\phi }}^i)\delta (\omega -{\omega }_i)$$

$$1=\int \mathcal{D}P\mathcal{D}\hat{P}{\text{rme}}^{N\text{rmi}\int \text{rmd}\mathbf{\phi }\text{rmd}\omega P(\mathbf{\phi },\omega )\hat{P}(\mathbf{\phi },\omega )-\text{rmi}\int \text{rmd}\mathbf{\phi }\text{rmd}\omega \hat{P}(\mathbf{\phi },\omega ){\sum }_i\delta (\mathbf{\phi }-{\mathbf{\phi }}^i)\delta (\omega -{\omega }_i)}.$$

$$\varphi [\hat{\varrho }]$$

$$S[P,\hat{P}]$$

$$\varphi [\hat{\varrho }]$$

$$W_0(\mathbf{\phi })$$

$$Z_q$$

$$T(\mathbf{\phi },{\mathbf{\phi }}^{\prime })$$

$$r[W_0(\mathbf{\phi })]$$

$$\mathrm{Tr}({\mathbf{T}}^{\ell })=(q-1{)}^{\ell }[\mathrm{Tr}({\mathbf{M}}^{\ell })-1]+(-1{)}^{\ell }$$

$$\varphi [\hat{\varrho }]$$

$$W_0(\mathbf{\phi })$$

$$Z(\mathbf{X})$$

$$\displaystyle\mathbf{X}=\left(\begin{array}[]{cccccc}x\left(\mu_{1}\right)\mathbf{I}_{n_{\mu_{1}}}&\mathbf{0}&&\dots&&\mathbf{0}\\ \mathbf{0}&\ddots&\mathbf{0}&&&\\ &\mathbf{0}&x\left(\mu_{M}\right)\mathbf{I}_{n_{\mu_{M}}}&\mathbf{0}&&\vdots\\ \vdots&&\mathbf{0}&\overline{x\left(\mu_{1}\right)}\mathbf{I}_{m_{\mu_{1}}}&\mathbf{0}&\\ &&&\mathbf{0}&\ddots&\mathbf{0}\\ \mathbf{0}&&\dots&&\mathbf{0}&\overline{x\left(\mu_{M}\right)}\mathbf{I}_{m_{\mu_{M}}}\end{array}\right).$$

$$\int \text{rmd}{\mathbf{\phi }}^{\prime }{\text{rme}}^{-\text{rmi}\mathbf{\phi }\cdot {\mathbf{\phi }}^{\prime }-\frac{1}{2}\text{rmi}{\mathbf{\phi }}^{\prime }\cdot \mathbf{X}{\mathbf{\phi }}^{\prime }}=Z(\mathbf{X}){\text{rme}}^{\frac{1}{2}\text{rmi}\mathbf{\phi }\cdot {\mathbf{X}}^{-1}\mathbf{\phi }}.$$

Full text

1 Introduction

A huge amount of scientific work has been devoted in recent decades to the study of random graphs, motivated partly by their interesting mathematical properties [1] and partly because of their frequent and fruitful use in the modelling of complex systems [2]. The most famous random graph ensembles are probably the Erdös-Renyi model (ER) [3, 4, 5] and the configuration model (CM) [6]. These are popular because models of interacting systems defined on their typical graph instances can often be solved analytically, and because sampling graphs from these ensembles is easy. Both properties derive mainly from the fact that ER and CM graphs are typically locally tree like, which enables the application of many relatively simple mathematical and numerical approaches. In fact, nearly all mathematical and computational techniques currently available for analysing processes on large graphs (cavity methods, belief propagation and other message passing algorithms, generating functional analysis, conventional replica methods, etc) rely explicitly or implicitly on being able to neglect the presence of short loops, or on being able to treat such loops as perturbations of a fundamentally tree-like architecture.

Ironically, their built-in locally tree like topology makes both ER and CM ensembles, and their numerous tree-like variations such as scale-free network models, also rather unsuitable as models for real networks. Simplicity comes at a nontrivial cost. Real networks are usually not tree like, but have many short loops, i.e. closed non backtracking paths [7]. Moreover, their significant ‘loopiness’ cannot be ignored, since it has profound functional implications111For instance, extensive loopiness is the main topological difference between the trivially solved ferromagnets on Bethe lattices or on random regular graphs with degree six, and the as yet unsolved three-dimensional Ising model. . Hence there is a clear need for simple loopy graph ensembles, and for new mathematical tools with which to analyse them.

Studying random graphs is also important in the area of null models for real observed networks. Null models of networks are unbiased random graph ensembles describing graphs that share with the networks the values of a given set of observables, i.e. maximum entropy ensembles [8, 9] subject to specific topological constraints. The constraints can be imposed as ‘soft’ conditions, where observables are matched on average, or as ‘hard’ conditions, where observables are matched by every graph in the ensemble. Realistic null models for graphs must be sparse and loopy. The simplest approach would be to control the number of edges and the number of triangles via soft constraints, leading to e.g. exponential ensemble of [10]. Unfortunately, upon varying the control parameters this ensemble was found to switch between very weak clustering and dominance by dense graphs [11, 12, 13]. See also [14]. The alternative is to impose both the number of single edges and of edges in a triangle around a node via hard constraints [15], but these graphs do not lend themselves as easily to tractable solutions, see e.g. [16, 17, 18].

In the present paper we focus only maximum entropy ensembles of $q$-regular graphs, with soft-constrained adjacency matrix eigenvalue spectra to introduce ‘loopiness’. We use the imaginary replica approach to obtain analytical results on the dependence of the spectrum on the functional Lagrange parameter, following [19], and we build on recent studies such as [20, 21, 22, 23, 24, 25]. The functional nature of the constraint in our ensemble allows us to bias the number of eigenvalues in all infinitesimal intervals simultaneously, while recovering the average spectral density through a simple functional derivative. We compute analytically in leading two orders in the system size the expected spectrum, a calculation found to take its most natural form using Chebyshev polynomials, and is valid in the regime where the loops are still sufficiently rare to prevent loop-loop interactions from becoming relevant. The simplest nontrivial members of our family of spectrally constrained maximum entropy ensembles are those where the spectral constraint reduces to a constraint on the expected number of triangles (equivalently, on the average clustering coefficient), with a single Lagrange parameter $\alpha$. We show that these models always exhibit transitions into a phase where the graphs shatter into an extensive number of disconnected cliques of $q+1$ nodes, analogous to the condensation transition in the model in [10]. We are able to present a description of the ensemble for all values of $\alpha$ for large enough $N$. We also show how our general theory can be applied to compute spectra for other ensembles with more complicated Lagrange parameters.

Our paper is organised as follows. In section 2 we give the relevant definitions of our spectrally constrained ensembles. Section 3 is devoted to the derivation of the generating function for the adjacency matrix spectrum. The theory is then applied in section 4, first to existing models (recovering known results as a test), followed by applications to other loopy graph ensembles. Technical details are often delegated to Appendices, to improve the flow of the paper. We end with a discussion of the results obtained, and an outlook on future work.

2 Definitions

We study ensembles of simple nondirected $N$-node regular graphs with degree $q$. Each graph is defined by its symmetric $N\times N$ adjacency matrix $\mathbf{A}$, with entries $A_{ij}=1$ if nodes $i$ and $j$ are connected, and $A_{ij}=0$ otherwise ($A_{ii}=0$ for all $i$). A path of length $\ell$ on a graph is a sequence of $\ell+1$ pairwise connected nodes. A closed path starts and ends in the same node. Loops are closed paths without repetition of nodes, except for the first and the last node. The formula for the total number of closed paths of a given length in a graph $\mathbf{A}$ is, apart from a simple overcounting factor,

[TABLE]

It follows from the relation $\Tr(\mathbf{A}^{\ell})=N\int\!\rmd\mu\,\varrho(\mu|\mathbf{A})\mu^{\ell}$. that controlling the numbers of closed paths of all lengths $\ell$ in random graphs is equivalent to controlling the moments of the spectral density $\varrho(\mu|\mathbf{A})$ of $\mathbf{A}$,

[TABLE]

where $\mu_{i}(\mathbf{A})$ is the $i$-th eigenvalue of $\mathbf{A}$.

In exponential spectrally constrained ensembles [19], the graph probabilities $p(\mathbf{A})$ on the set $G$ of simple nondirected $N$-node graphs are determined by maximising the Shannon entropy $S[p]=\sum_{\mathbf{A}\in G}p(\mathbf{A})\log p(\mathbf{A})$, subject to prescribed values of all degrees and a prescribed expectation value of the spectral density. For $q$-regular random graphs this gives

[TABLE]

Here $\hat{\varrho}(\mu)$ is a functional Lagrange multiplier. By construction, (3) defines the most unbiased ensemble of $q$-regular nondirected graphs with a prescribed adjacency matrix spectrum. We write averages over (3) as $\langle f(\mathbf{A})\rangle=\sum_{\mathbf{A}\in G}p(\mathbf{A})f(\mathbf{A})$. The expected eigenvalue density $\varrho(\mu)=\langle\varrho(\mu|\mathbf{A})\rangle$ can be obtained from a generating function $\phi[\hat{\varrho}]$:

[TABLE]

Our main interest is in finding an analytical expression for the expected density $\varrho(\mu)$ in terms of the functional Lagrange parameter $\hat{\varrho}(\mu)$.

For the simple choice $\hat{\varrho}(\mu)=\alpha\mu^{3}$ we recover from (3) the model of [26], in which the number of loops of length three (i.e. of triangles) is constrained:

[TABLE]

The number of triangles $n_{\triangle}(\mathbf{A})$ differs from the trace only by an overcounting factor, viz. $n_{\triangle}(\mathbf{A})=\frac{1}{6}\Tr(\mathbf{A}^{3})$. The generating function associated with (6) is $\phi(\alpha)=N^{-1}\log Z(\alpha)$. From $m(\alpha)=\partial\phi(\alpha)/\partial\alpha=N^{-1}\langle{\rm Tr}(\mathbf{A}^{3})\rangle$ follows the average clustering coefficient $\langle C(\mathbf{A})\rangle=m(\alpha)/q(q-1)$. For $q=2$ one can compute $\phi(\alpha)$ using standard combinatorics [26], for both $\alpha={\mathcal{O}\left(N^{0}\right)}$ and $\alpha={\mathcal{O}\left(\log N\right)}$. For arbitrary $q$, the expected loop density $m(\alpha)$ always vanishes for $N\to\infty$ if $\alpha={\mathcal{O}\left(1\right)}$. As in [26], one could rescale $\alpha$ with a factor $\log N$. While this could give an asymptotic theory with finite loop densities, we will find that it would not be a useful network model for applications.

3 Evaluation of the generating function

3.1 Imaginary replica approach

We note that all $q$-regular graphs of size $N$ have identical probabilities in an Erdös-Renyi (ER) ensemble with average degree $q$,

[TABLE]

Hence we can rewrite (5) as an average over this ER ensemble:

[TABLE]

Upon using the Edwards-Jones formula [27] for $\varrho(\mu|\mathbf{A})$, writing the integral over eigenvalues as $\int\!\rmd\mu\ldots=\lim_{\Delta\to 0}\Delta\sum_{\mu}\ldots$, and after some modest manipulations, the key quantity in this expression can be written as follows, with infinitely many imaginary replicas (two for each eigenvalue$\mu$, $n_{\mu}$ and $m_{\mu}$):

[TABLE]

Here ${\mu_{\varepsilon}}=\mu+\rmi\varepsilon$ and

[TABLE]

One initially takes $n_{\mu},m_{\mu}\in\mathbb{N}$, in order to perform the calculation, followed by analytic continuation to the relevant imaginary values. In its above form, (9) appeared first in [28, 19], but similar formulae involving limits of replica dimensions to non-zero values have been introduced in different contexts, in particular when counting the number of eigenvalues in certain intervals for random matrix ensembles, see e.g. [20, 21, 22, 29, 23, 24, 25]. We can combine the integrals in (9) as follows:

[TABLE]

where d[\boldsymbol{\Phi},\boldsymbol{\Psi}]\!=\!\prod_{i=1}^{N}\prod_{\mu}\big{[}\big{(}\prod_{\alpha_{\mu}=1}^{n_{\mu}}\!\rmd\boldsymbol{\phi}^{i}_{\mu,\alpha_{\mu}}\big{)}\big{(}\prod_{\beta_{\mu}=1}^{m_{\mu}}\!\rmd\psi^{i}_{\mu,\beta_{\mu}}\big{)}\big{]}. To simplify our notation we introduce the vector $\boldsymbol{\varphi}\in\mathbb{R}^{d}$, where $d=\sum_{\mu}n_{\mu}\!+\!\sum_{\mu}m_{\mu}$, with entries

[TABLE]

the dot product

[TABLE]

and the following $d\!\times\!d$ diagonal matrix $\boldsymbol{\mu}$, in which $\mu\in\{\mu_{1},\ldots,\mu_{M}\}$ and $\mathbf{I}_{n}$ denotes the $n$ dimensional identity matrix:

[TABLE]

Here $M$ denotes the number of $\mu$-values in the discretized eigenvalue integral, so $M\to\infty$ when we take the limit $\Delta\to 0$. We finally introduce the two shorthands

[TABLE]

The above definitions, together with the integral form of the Kronecker delta, $\delta_{nm}=(2\pi)^{-1}\!\int_{-\pi}^{\pi}\!\rmd\omega~{}\exp[\rmi\omega(n-m)]$, enable us to compute the generating function (8) following [19], which in turn is reminiscent of previous spectral calculations for sparse random graphs [30, 31]. Upon dropping the irrelevant constant in (8) we get

[TABLE]

Since we intend to compute finite size spectrum fluctuations, in expanding the logarithm for large $N$ we keep both the ${\mathcal{O}\left(N\right)}$ and ${\mathcal{O}\left(1\right)}$ terms in the exponent:

[TABLE]

We now introduce the order parameter

[TABLE]

We enforce it by inserting the following functional integral, obtained by writing delta functions for each $(\boldsymbol{\varphi},\omega)$ in integral representation, and with the usual path integral measure $\mathscr{D}P=\prod_{\boldsymbol{\phi}}\prod_{\omega}[\rmd P(\boldsymbol{\phi},\omega)\sqrt{N\Delta_{\varphi}\Delta_{\omega}/2\pi)}]$ (where $\Delta_{\varphi},\Delta_{\omega}\to 0$):

[TABLE]

The result is:

[TABLE]

It was shown in [31] how this type of integral can be reduced to an integral over a single functional variable. In A we work out the details, leading to

[TABLE]

in which $U_{1}(\boldsymbol{\varphi},\boldsymbol{\varphi}^{\prime})=\rme^{-\rmi\boldsymbol{\varphi}\cdot\boldsymbol{\varphi}^{\prime}}$, and the function $W_{0}(\boldsymbol{\varphi})$ is to be solved from

[TABLE]

and

[TABLE]

Finally, following [32] we may use the following identity222The proof of this interesting identify follows directly from the operator properties $\mathbf{M}\mathbf{B}=\mathbf{B}\mathbf{M}=\mathbf{B}^{2}={1\!{\rm I}}$, where $M(\boldsymbol{\varphi},\boldsymbol{\varphi}^{\prime})=r[W_{0}(\boldsymbol{\varphi})]U_{1}(\boldsymbol{\varphi},\boldsymbol{\varphi}^{\prime})$ and $B(\boldsymbol{\varphi},\boldsymbol{\varphi}^{\prime})=W_{0}(\boldsymbol{\varphi})\int\!\rmd\boldsymbol{\psi}~{}U_{1}(\boldsymbol{\varphi}^{\prime},\boldsymbol{\psi})W_{0}(\boldsymbol{\psi})$.:

[TABLE]

with $M(\boldsymbol{\varphi},\boldsymbol{\varphi}^{\prime})=r[W_{0}(\boldsymbol{\varphi})]U_{1}(\boldsymbol{\varphi},\boldsymbol{\varphi}^{\prime})$, to simplify (30) modulo an additive constant to

[TABLE]

Expression (36) was originally presented in [31], and we have indeed chosen our notation at the start deliberately to emphasize and exploit the similarity. However, although identical in structure, the present formula (36) differs from the one in [31] in the underlying definitions of the fields $\boldsymbol{\varphi}$, the dot product $\boldsymbol{\varphi}\cdot\boldsymbol{\varphi}^{\prime}$ and the function $\nu(\boldsymbol{\varphi})$, which here all involve the full eigenvalue spectrum and describe our present controlled non-uniform measures over the space of graphs.

3.2 Replica symmetric solution

In order to continue, we assume that the order parameter $W_{0}(\boldsymbol{\varphi})$ is replica symmetric (RS), i.e. invariant under all permutations of all replicas (noting that in the present problem we have a separate replica index for each eigenvalue $\mu$, and that $W_{0}(\boldsymbol{\varphi})$ is not a normalized distribution). We write $W_{0}(\boldsymbol{\varphi})$ as a superposition of zero mean complex Gaussian distributions:

[TABLE]

where $\mathbf{X}\in\mathbb{R}^{d\times d}$ is a diagonal matrix with the following structure:

[TABLE]

Expression (38) is indeed invariant under all permutations of replica indices with any fixed value of $\mu$, that is $\{\phi_{\mu,1},\ldots,\phi_{\mu,n_{\mu}}\}$ and $\{\psi_{\mu,1},\ldots,\psi_{\mu,m_{\mu}}\}$. The new RS order parameter is the distribution $W(\mathbf{X})$, where each $\mathbf{X}$ is specified by $M$ complex numbers $x(\mu)$. Hence the integration in (15) is over the real and imaginary part of each $x(\mu)$, so $\rmd\mathbf{X}=\prod_{\mu}\rmd\,[x(\mu)]\rmd\,\textrm{Im}[x(\mu)]$. For (37) to be well defined, we must restrict all $x(\mu)$ to have $\textrm{Im}\;x(\mu)<0$. We may assume that $\int\!\rmd\mathbf{X}~{}W(\mathbf{X})=1$, since possible non-normalization of $W_{0}$ is reflected in the inclusion in (38) of a constant $\mathcal{C}$.

We note that for the present definition (15) of the dot product, the following identity is still valid.

[TABLE]

Insertion of our RS ansatz (38) into the full order parameter equation (31) shows, using the above identity, that the RS ansatz indeed gives a solution of (31), provided the RS order parameter satisfies

[TABLE]

We similarly derive $Z_{q}=C^{q}\mathcal{Z}_{q}$, giving in combination with (49):

[TABLE]

The above RS order parameter equations (47,48) have one specific simple solution, namely the delta distribution $W(\mathbf{X})=\delta(\mathbf{X}-\mathbf{X}^{\star})$, in which the entries of $\mathbf{X}^{\star}$ satisfy

[TABLE]

Of the two possible solutions of this equation we must choose the one with $\textrm{Im}\;x(\mu)<0$:

[TABLE]

For this special RS solution we have

[TABLE]

and hence also $\mathcal{C}^{2}=Z(\mathbf{X}^{\star})/Z(\mathbf{X}^{\star}\!-(\mathbf{X}^{\star})^{-1})$. The kernel $\mathbf{M}$ which appears in the generating function (36) will now have the following entries:

[TABLE]

In B we work out the traces ${\rm Tr}(\mathbf{M}^{\ell})$ for the RS solution, and find

[TABLE]

Where $Z({\mu_{\varepsilon}}|\mathbf{A}^{\star}_{\ell,\mu})$ denotes the original complex Gaussian integral defined in (10), and $\mathbf{A}^{\star}_{\ell,\mu}$ is now the $\ell\times\ell$ adjacency matrix of a loop of length $\ell$ in the presence of a complex field acting on the diagonal, of value $(2\!-\!q)/x^{\star}(\mu)$:

[TABLE]

Substituting (50) and (58) into expression (36) for the generating function, folowed by using formulae (54,3.2) for the constants $\mathcal{Z}_{q-1}$ and $\mathcal{Z}_{q}$, then gives

[TABLE]

3.3 Imaginary replica limits

At this stage we can safely take the three limits defined in (23), where first for each discretized eigenvalue $\mu$ of the adjacency matrix the replica dimensions $n_{\mu}$ and $m_{\mu}$ take specific imaginary values, followed by the limit $\Delta\to 0$ that converts discretized eigenvalues of adjacency matrices into continuous ones. The objects in (59) affected by these limits are all of the following form, with $\hat{\varrho}^{\prime}(\mu)=\rmd\hat{\varrho}(\mu)/\rmd\mu$:

[TABLE]

In particular, application to $f(\mu)=[2\pi/\rmi x^{\star}(\mu)]^{\frac{1}{2}}$ and to $f(\mu)=Z({\mu_{\varepsilon}}|\mathbf{A}^{\star}_{\ell,\mu})$ gives

[TABLE]

and

[TABLE]

This gives us

[TABLE]

Since $\delta\phi[\hat{\varrho}]/\delta\hat{\varrho}(\mu)=\langle\varrho(\mu|\mathbf{A})\rangle$, the first line of (63) is the generator of the asymptotic spectrum in the limit $N\to\infty$, whereas the second line will give us the ${\mathcal{O}}(N^{-1})$ finite size corrections to the spectrum. In C we show that, upon taking the limit $\varepsilon\to 0$, the factor inside the curly brackets in the first line indeed works out to be exactly the Kesten-McKay law (KM) [33, 34] for random regular graphs:

[TABLE]

This shows that, for regular graphs, the deformation of the measure in the ensemble (6) does not alter the resulting spectrum in leading order, but in sub-leading order ${\mathcal{O}}(N^{-1})$. In regular graphs, the Lagrange parameter $\hat{\varrho}(\mu)$ apparently needs to be rescaled further with $N$ to induce a spectrum that in leading order differs from (64), similar to what was found in [26].

Having simplified the first line of (63) to $\int\!\rmd\mu~{}\hat{\varrho}(\mu)\varrho_{0}(\mu)$, we now work out further the exponent in the second line of (63). First, in C we show that

[TABLE]

We can evaluate the second term in the exponent using the eigenvalues $\lambda_{k}=2\cos(2\pi k/\ell)$ of the adjacency matrix $A_{ij}=\delta_{i,j+1}+\delta_{i,j-1}~{}({\rm mod}~{}\ell)$ of a length-$\ell$ loop, and the identity $\rmd x^{\star}(\mu)/\rmd\mu=-\rmi x^{\star}(\mu)/\sqrt{4(q\!-\!1)\!-\!\mu^{2}}$:

[TABLE]

With the above simplifications we can write both the leading two orders in $N$ of the generating function $\phi[\hat{\varrho}]$ and of the resulting average spectrum $\varrho(\mu)=\delta\phi[\hat{\varrho}]/\delta\hat{\varrho}(\mu)$ for our ensemble (3), for Lagrange parameters $\hat{\varrho}(\mu)={\mathcal{O}}(1)$, in the following transparent form, which represents one of the main results of this paper:

[TABLE]

The sum over $\ell$ need not always be convergent, but this has been shown not to necessarily pose problems [35, 32, 31]. The coefficient $(q\!-\!1)^{\ell}/2\ell$ in $\phi[\hat{\varrho}]$ is exactly the asymptotic expected number of loops of finite length $\ell$ inside a random regular graph [36]. It is amazing that this can be recovered with a replica calculation, and gives an intuitive interpretation of the series in (68): $\ell[g_{\ell}(\mu)\!-\!h(\mu)]$ is the correction to the Kersten-McKay spectrum formula due to the appearance of a loop of length $\ell$, as discussed in [32, 31]. In our present ensemble (3), the number of loops of a given length are given by the usual number found in random regular graphs, multiplied by a factor that depends on the spectral Lagrange parameter $\hat{\varrho}(\mu)$. Given that the effect of the loops is additive in (68), we must expect that these spectral corrections come from isolated and well separated loops in the graph. We wish to point out that, while a cavity approach could account for the presence of loops, it would not be able to provide information on their average number in an ensemble such as (3). The imaginary replica approach presented here, in contrast, has simultaneously provided for the ensemble (3) both the spectrum formula and the expected number of loops.

3.4 Remaining integrals over eigenvalues

In our present theory we have an as yet arbitrary functional Lagrange parameter $\hat{\varrho}(\mu)$ which controls the dependence of the graph probabilities on their expected spectra. In (68) we still have integrals over $\hat{\varrho}(\mu)$, of the form:

[TABLE]

While expressions (65,66) for $h(\mu)$ and $g_{\ell}(\mu)$ will turn out useful in establishing links with previous research in a subsequent section, here we will continue the further evaluation of ${\cal J}_{\ell}[\hat{\varrho}]$ using the earlier forms

[TABLE]

These give

[TABLE]

We now use $x^{\star}(\mu)+{x^{\star}(\mu)}^{-1}=-\mu-(q\!-\!2)/x^{\star}(\mu)$, which follows directly from (51):

[TABLE]

In this expression we recognize the generating function of the Chebyshev polynomials $T_{n}(t)$ [37]. These are defined for $t\in[-1,1]$, and can be written in explicit form as

[TABLE]

They obey the orthogonality relation

[TABLE]

as well as

[TABLE]

The first five Chebyshev polynomials are [37]:

[TABLE]

For the evaluation of (73), in particular, we may apply the generating function identity

[TABLE]

to the choices $x=1/x^{\star}(\mu)$ and $t=-\cos(2\pi k/\ell)$, in order to obtain

[TABLE]

We next use $\rmd x^{\star}(\mu)/\rmd\mu=-\rmi x^{\star}(\mu)/\sqrt{4(q\!-\!1)\!-\!\mu^{2}}$ and the short-hands

[TABLE]

This results in

[TABLE]

If $\mu^{2}<4(q\!-\!1)$ one has $|x^{\star}(\mu)|^{2}=q-1$ and $\sqrt{4(q\!-\!1)\!-\!\mu^{2}}\in[0,\infty)$. For $\mu^{2}>4(q\!-\!1)$, on the other hand, we have $x^{\star}(\mu)\in{\rm I\!R}$ and $\sqrt{4(q\!-\!1)\!-\!\mu^{2}}$ is purely imaginary. We also note that for $\mu^{2}<4(q\!-\!1)$ we may write

[TABLE]

In combination these properties allow us to simplify ${\cal I}_{\ell}[\hat{\varrho}]$ to

[TABLE]

The above result shows that the Chebyshev polynomials form the natural basis in terms of which to express the functional Lagrange parameter $\hat{\varrho}(\mu)$, and is inserted into our spectrum formula (68) to give

[TABLE]

It is instructive to work out ${\cal J}_{\ell}[\hat{\varrho}]$ for $\ell\geq 3$ and some simple choices of $\hat{\varrho}(\mu)$:

• •

$\hat{\varrho}(\mu)=\mu^{3}$:

This choice corresponds to random regular graphs in which the number of triangles is controlled. We use $t^{3}=\frac{1}{4}T_{3}(t)+\frac{3}{4}T_{1}(t)$ and the orthogonality relation (75):

[TABLE]

• •

$\hat{\varrho}(\mu)=\mu^{4}$:

This choice corresponds to random regular graphs in which the number of squares is controlled. We use $t^{4}=\frac{1}{8}T_{4}(t)+\frac{1}{4}T_{2}(t)+\frac{3}{8}T_{0}(t)$ and the orthogonality (75):

[TABLE]

4 Applications of the general theory

4.1 Recovering previous results as a test

Upon making the trivial choice $\hat{\varrho}(\mu)=0$ we return to the conventional ensembles with uniform probabilities, and our equations (67,68) recover the natural spectrum fluctuations of random regular graphs, as previously studied in detail in [38] and with the traditional replica method (where $n\to 0$) in [31]:

[TABLE]

with

[TABLE]

This series was summed in [31], and we can connect the result of the summation, in the notation of [31], directly to the theory developed in the present paper as follows:

[TABLE]

in which now

[TABLE]

Here $h(\mu)$ and $x^{\star}(\mu)$ are given in (65) and (52), respectively.

As a second test we can make the special choices $q=2$ and $\hat{\varrho}(\mu)=\alpha\mu^{3}$, resulting in the ensemble that was studied in [26] via direct combinatorics, i.e. without the replica method. This particular model represents the simplest solvable non-uniform random graph ensemble with tuneability of the frequency of short loops. First, by setting $q=2$ our general results (65,66,68) simplify greatly. We now find that

[TABLE]

Upon inserting also $\hat{\varrho}(\mu)=\alpha\mu^{3}$ into (68) we need the values of

[TABLE]

With this the ensemble spectrum becomes

[TABLE]

and for the triangle density $m(\alpha)=\int\!\rmd\mu~{}\varrho(\mu)\mu^{3}$ we obtain

[TABLE]

These results are indeed identical to those derived combinatorially in [26].

4.2 Triangularly constrained regular graph ensemble with arbitrary degree

We proceed to apply the general theory developed in the previous section to the graph ensemble (6) with controlled numbers of triangles, i.e. with $\hat{\varrho}(\mu)=\alpha\mu^{3}$, but now for arbitrary values of the degree $q$ where the direct combinatorial approach of [26] is no longer possible. We can start directly by inserting (87) into (86), upon adding the control parameter $\alpha$, giving

[TABLE]

Here $\varrho_{1}(\mu)$ is the function (90) that already appeared in the spectrum of the non-deformed ensembles of [31], and for which we can use the resummation (91). The impact of controlling the graph probabilities (6) with a nonzero Lagrange parameter $\hat{\varrho}(\mu)=\alpha\mu^{3}$ is fully concentrated in $\tilde{\varrho}_{1}(\mu)$. We next insert our earlier expressions for $g_{\ell}(\mu)$ and $h(\mu)$ into (101) and simplify the result where possible:

[TABLE]

The results of testing this prediction against numerical simulations are shown in Figure 1, and reveal excellent agreement. In the simulations we sampled numerically from (6) with an edge swap based Markov Chain Monte Carlo algorithm (MCMC) with nontrivial acceptance probabilities. Edge swaps are accepted or rejected depending on the change in the number of loops and on the change in the possible number of possible swaps. This corrects for entropic effects, see e.g. [39] or [9]. Since we work with a system of finite size $N$, our predictions refer to the average eigenvalue density, not to the density of individual graph instances. The error bars in Figure 1 are computed following 10 different initialization seeds of the MCMC algorithm, consisting of distinct regular graphs sampled uniformly with a configuration model algorithm. Following each initialization, 20 samples were taken, separated in algorrithmic time by $\sim 10^{3}$ accepted MCMC swaps per link in the graph.

We can also calculate the expected triangle density $m(\alpha)=\int\!\rmd\mu~{}\varrho(\mu)\mu^{3}$ for the ensemble (6). It is easier to do this by integrating over (86) rather than via (100,LABEL:eq:rho1tilde), although both routes give the same result:

[TABLE]

This formula gives very accurate results for $\alpha$ values up to a certain point, defined as $\alpha_{1}(N)$ in the next section. This can be seen very clearly in figures 2 and 4, where we test its predictions for ensembles (6) with $q=3$. Since $m(\alpha)$ represents an ensemble average, we compare (103) against the average loop density over multiple graphs drawn from the ensemble, $\hat{m}(\alpha)=M^{-1}\sum_{m=1}^{M}\Tr(\mathbf{A}_{m}^{3})$. Figures 2 and 4 show averages and standard deviations of the estimator $\hat{m}(\alpha)$, for 50 different small samples sampled during MCMC simulations, separated in algorithmic time by $\sim 10^{3}$ swaps per link (in order to ensure independence of the $M$ samples).

We have developed a theory that quantifies the ${\mathcal{O}\left(1/N\right)}$ effects on the eigenvalue spectrum of probability deformations in ensembles of the general family (3), in which loops can be induced via the functional Lagrange parameter $\hat{\varrho}(\mu)$, and we applied our results to a specific member (6) of this family. As mentioned before, the theory implicitly assumes that loops inside the graph are far away from each other. As the control parameter $\alpha$ in (6) is increased for fixed system size $N$, we must therefore expect the behaviour of the ensemble to start deviating from the predictions (100,LABEL:eq:rho1tilde,103) as soon as the loops start to interact. This can indeed be seen in figure 2. As $\alpha$ increases the clustering coefficient starts deviating from (103), which is shown as a dashed line. Mathematically, one can explain the deviations from (103) as the emergence of higher order corrections to the saddle point approximation, ${\mathcal{O}\left(N^{-\gamma}\right)}$ with $\gamma>1$, that were not incorporated into the replica calculation. These would account for the presence of loops that are not isolated from each other. The accuracy of (100,LABEL:eq:rho1tilde) and (103) suggests that calculating higher order corrections in the replica calculation would improve our predictions, but this would require of course a much more complicated calculation. In the next section we will explore what happens as we keep increasing $\alpha$ beyond the validity of our replica calculation.

4.3 Phases of the ensemble and the shattering transition

We will now give a qualitative picture of the behaviour of the ensemble (6) for all values of $\alpha\in[0,\infty)$. We will focus on $q\geq 3$, since the case $q=2$ was already covered in [26]. In MCMC simulations one observes three distinct regimes, which are not phases in a rigorous thermodynamic sense, but size dependent ranges of $\alpha$ values that exhibit qualitatively different phenomenology:

• •

Small $\alpha$: connected phase

The triangle promoting probability bias in the ensemble introduces isolated loops embedded in the giant component. Here the analysis of the previous section should apply, as is confirmed in figures 2 and 4 for different values of $q$. Indeed one observes only small deviations from (103), as one approaches the next phase.

• •

Intermediate $\alpha$: partially connected phase

Edges can now be part of more than one triangle, and the graphs contain an increasing number of cliques of $q+\!1$ nodes, denoted by $K_{q+1}$. The triangle density and hence the clustering coefficient grow considerably faster than in the previous phase, increasingly so for larger degrees $q$.

• •

Large $\alpha$: disconnected phase

Here the graphs break down completely into large collections of those cliques that had started to appear in the previous phase. The resulting configurations correspond to $q$ regular graphs with the maximum possible number of triangles. In analogy with physics, we call these ground states.

We label the transition points between the phases $\alpha_{1}(N)$ and $\alpha_{2}(N)$, see figure 2. Both $\alpha_{1}(N)$ and $\alpha_{2}(N)$ grow logarithmically with $N$. We refer to the transition from connected to partially connected as the shattering transition to highlight its topological nature. In figure 3 we show typical graphs sampled via MCMC in the three phases.

In order to complement the previous distinctions with quantitative estimates, we will next give an alternative derivation of (103) that incorporates higher order effects. We can always write averages over (6) in terms of averages over random regular graphs (RRG), described by the ensemble $p_{0}(\mathbf{A})={\cal N}_{q}^{-1}\prod_{i\leq N}\delta_{q,\sum_{j}A_{ij}}$ with uniform probabilities. We will write average over the unbiased random regular graph ensemble $p_{0}(\mathbf{A})$ as $\langle\ldots\rangle_{\rm RRG}$. In particular,

[TABLE]

We can therefore use some of the rigorous results from random graph theory [40] established for unbiased ensembles. For instance, a standard result on RRGs [36] concerns the asymptotic distribution of isolated triangles:

[TABLE]

in which $\triangle(\mathbf{A})$ is the number of isolated triangles in graph $\mathbf{A}$. To understand the shattering transition, we also need to know the statistics of cliques $K_{q+1}$ in RRGs. We do this by splitting the triangle count in two contributions, one from the isolated triangles inside the giant component and the other from the triangles in the cliques $K_{q+1}$. We denote the numbers of each type by $T(\mathbf{A})$ and $K(\mathbf{A})$ respectively.333We hereby disregard as insufficiently relevant those triangles that are neither isolated, nor in cliques. This step, which is motivated by what is observed in the MCMC simulations, will be validated a posteriori by the accuracy of the resulting prediction. Thus

[TABLE]

We now assume that also $K(\mathbf{A})$ follows a Poisson distribution, and that for large $N$ the joint distribution of $T$ and $K$ factorizes asymptotically. This assumption is to be confirmed a posteriori in MCMC simulations. Now

[TABLE]

Here $\lambda_{T}=\frac{1}{6}(q\!-\!1)^{3}$, and we can obtain the parameter $\lambda_{K}$ for the Poisson distribution for the cliques $K_{q+1}$ from [41]. Details are given in D. Upon approximating $N!/(N\!-\!q\!-\!1)!\approx N^{q+1}$ we obtain

[TABLE]

We can now proceed with the evaluation of (104):

[TABLE]

From this directly follows $m(\alpha)=\partial\phi(\alpha)/\partial\alpha$:

[TABLE]

The first term coincides with formula (103) from the previous section. The second term in (110) represents the impact of cliques, and, according to figure 4, accounts for most of the deviations from (103). In spite of our approximation of only accounting for isolated loops and isolated cliques, the resulting description is seen to give very good agreement with simulations for the whole range of $\alpha$ values.

As one might expect, the MCMC sampling algorithm slows down as it approaches the ground state. While we will not carry out a detailed dynamical analysis, we will mention the MCMC process slows down considerably precisely in the partially connected phase. To obtain good (i.e. sufficiently independent) samples even close to the ground state, we increased the number of accepted swaps per link in between samples beyond $\alpha_{1}(N)$ to values in the range of $10^{4}-10^{5}$ accepted swaps per link. We also increased considerably the waiting time before the first sample to $\sim 10^{8}$ accepted swaps per link, to allow the system to escape from possible metastable states.

Expression (110) has a clear interpretation: the expected number of subgraphs of the types $T$ and $K$ are boosted independently when increasing $\alpha$, each with an exponential factor in accordance with the model (6). This already provides an explanation for the phases described previously. If we denote by $m_{T}(\alpha)$ and $m_{K}(\alpha)$ the first and second term of (110), then we can describe the phases in terms of the relation between these two terms.

• •

Small $\alpha$: connected phase

Here we have $m_{T}(\alpha)\gg m_{K}(\alpha)$. Even though cliques $K_{q+1}$ may be present, their probability is too small to be relevant.

• •

Intermediate $\alpha$: partially connected phase

Here $m_{K}(\alpha)$ becomes significant. We may define the onset $\alpha_{1}(N)$ of the partially connected phase to be the point where $m_{K}(\alpha)=\eta m_{T}(\alpha)$ for some finite $\eta\in(0,1)$. Here we chose $\eta=1/10$, which was found appropriate in the ranges $q\leq 6$ and $N\leq 2000$. The shattering transition point is then given by

[TABLE]

• •

Large $\alpha$: disconnected phase

Here the contribution from $m_{K}(\alpha)$ dominates, and the whole graph is made of disconnected cliques. The critical value $\alpha_{2}(N)$ marking the start of this phase is defined by the instance where $m(\alpha_{2}(N))=q(q\!-\!1)$, i.e. where the maximum possible density of loops is achieved. We can replace $m(\alpha_{2}(N))$ by $m_{K}(\alpha_{2}(N))$ since the contribution from disconnected triangles is now very small, giving

[TABLE]

Even though $m_{K}(\alpha)$ in (110) has a higher power of $1/N$ compared to $m_{K}(\alpha)$, the prefactor of $\alpha$ for $m_{K}(\alpha)$ is higher in the exponential, viz. $(q\!+\!1)q(q\!-\!2)>6$ for $q\geq 3$. This means that for any $N$, the clique contribution $m_{K}(\alpha)$ will always grow faster with $\alpha$ than $m_{T}(\alpha)$, which implies it will always surpass it for large enough $\alpha$. We do not claim that graphs will only be made either of isolated loops in a giant component, but only that knowing the behaviour of the two quantities $m_{T}(\alpha)$ and $m_{K}(\alpha)$ appears sufficient to understand the overall behaviour of the loopy ensemble (6).

Let us briefly discuss the potential practical utility of (6) in light of the previous results. Our ensemble is the maximally unbiased random graph ensemble over regular graphs that satisfies the condition of having a particular clustering coefficient $C$. In order to achieve one’s desired value of $C$ it is only necessary to find the appropriate $\alpha(C)$ by solving $C=C(\alpha)=m(\alpha)/q(q\!-\!1)$ using (110). However, if one’s interest is in using (6) as a null model for real networks with link clustering, the presence of cliques is undesirable. If we aim to generate graphs with a single component and a nontrivial number of loops, we need to stay in the connected phase. Moreover, in this phase we have a very accurate control of $m(\alpha)$ and the spectral density through (100). We conclude that (6) can be a useful null model when $C\in(0,C(\alpha_{1}(N)))$. In that clustering range we can simply take

[TABLE]

While the shattering transition occurs at $\alpha=\alpha_{1}(N)$, the finite size nature of the problem makes it possible that some cliques appear already somewhat earlier. However, it is clear from Figure 2 that an upper bound to the level of clustering achievable without cliques is given by $C_{\rm max}=C_{T}(\alpha_{2}(N))$, the contribution to clustering from disconnected triangles at the transition point to the disconnected phase. Values $C\geq C_{\rm max}$ are impossible to achieve in the ensemble (6) without triangles appearing outside the giant component, and additional constraints would have to be introduced into the model to prevent the formation of isolated cliques. This dependence on $N$ of the hard upper bound, shown also in Figure 5, is somewhat unexpected. In particular,

[TABLE]

Hence for very large sizes, $N\gg 1$, even very small clustering coefficients are not accessible in the connected regime. This can be understood intuitively as an entropic effect. For fixed value of $C(\alpha)$ and large enough $N$, there are simply many more graphs that achieve $C(\mathbf{A})=C(\alpha)$ through cliques than graphs that achieve it through loops embedded in the giant component.

4.4 Other ensembles

So far we have focused specifically on the ensemble (6) as the simplest nontrivial instance of the more general family (5), suitable for testing limits and for developing further our intuition for the phenomenology of ‘loopy’ random graph ensembles. However, we have the more general results (85,86), applicable to any functional Lagrange parameter $\hat{\varrho}[\mu]$, with the key integral expressed as an expansion in Chebyshev polynomials. We will now turn to other choices for $\hat{\varrho}[\mu]$.

Our first choice is $\hat{\varrho}[\mu]=\alpha\mu^{3}+\beta\mu^{4}$, which generalizes the ensemble (6) in that we now control the number of closed paths of both length 3 and length 4:

[TABLE]

The calculations for this ensemble (115) are very similar to those carried out for (6), which allows us to be brief. We have already computed the relevant integrals in (87,88), and we can therefore immediately proceed to the spectral density:

[TABLE]

in which the functions $g_{\ell}(\mu)$ and $h(\mu)$ are given in (70,71). A comparisons of the predicted spectrum (117) with measurents in MCMC simulations, for $N=2000$ and $\alpha=\beta=0.2$, is shown in the left panel of Figure 6. The MCMC algorithm used was similar to the one described before, but now they also require monitoring the evolution of $\Tr(\mathbf{A}^{4})$ (in addition to $\Tr(\mathbf{A}^{3})$), as both appear in the move acceptance probabilities. In each run 100 samples were generated from each initial seed, after a burn-in (waiting time) of $\sim 10^{3}$ swaps per link. Error bars give the standard deviation corresponding to fluctuations between 10 different initial seeds, so that a total of 1000 graphs were averaged. As in the previous case, we recover the results from [26] when setting $q=2$. As one would expect, (117) is only valid in the vicinity of $(\alpha,\beta)=(0,0)$, to avoid the emergence of extensively many small fully connected $q$-regular bipartite graphlets, which maximize the number of $4-$loops around a node.

Our second alternative choice for $\hat{\varrho}[\mu]$ is the following block function, which introduces a bias in the graph probabilities depending on the number of eigenvalues inside the interval $[-1,1]$:

[TABLE]

Now we have $N\int\!\rmd\mu~{}\hat{\varrho}(\mu)\varrho(\mu|\mathbf{A})=\alpha\mathcal{I}(\mathbf{A}|[-1,1])$, where $\mathcal{I}(\mathbf{A}|[-1,1])$ denotes the number of eigenvalues of $\mathbf{A}$ inside the interval $[-1,1]$. In contrast to powers of $\mu$, understanding intuitively the topological effects of the choice (118) is not straightforward, notwithstanding the clear nontrivial effect on the observed spectrum. In this case we have ${\cal J}_{\ell}[\hat{\varrho}]=\alpha\int\!\rmd\mu~{}\theta[1\!-\!|\mu|][g_{\ell}(\mu)\!-\!h(\mu)]\not=0$ for all $\ell$, so we introduce a (sufficiently large) cutoff $L$ in the summation of (86). Since $|{\cal J}_{\ell}[\hat{\varrho}]|\ll 1$ we set this integral to zero for $\ell>L$, as was done previously in [32], leaving the truncation

[TABLE]

in which the integrals can be worked out in more explicit form, as we did for the previous cases. In figure 6 (right panel) we compare the prediction (119) with results from numerical MCMC samples, and observe a good agreement. We point out that generating graph samples from the spectrally constrained ensemble (5, 118) numerically is considerably more computationally expensive than the for the previous models. Here, instead of the number of triangles or squares, the number of eigenvalues inside the interval $[-1,1]$ has to be monitored. This requires that the full set of eigenvalues of the graph $\mathbf{A}$ has to be calculated after each edge swap, which necessitated parallel execution in multi-core computers, to reduce CPU time to a few weeks. We seen in Figure 6 that the deviations from $\varrho_{0}(\mu)$ are quite small, nevertheless they are nontrivial and are predicted accurately. To measure spectra at this level of detail, we averaged over $10^{4}$ graphs, separated during the MCMC process by $\sim 1$ swaps per link. Error bars are obtained by splitting this data set in groups of $10$.

The above results are similar in form to the ones derived for weighted graphs in [25], the main difference is that in [25] a second set of replicas with the traditional limit $n\to 0$ is introduced to get the spectrum. It is interesting to see that with the functional formalism (5) both the observable $\int\!\rmd\mu~{}\hat{\varrho}(\mu)\varrho(\mu)$ and the spectrum $\varrho(\mu)$ itself are calculated at the same time.

5 Discussion

In this paper we have extended and applied an analytic approach for describing constrained maximum entropy ensembles of finitely connected random loopy graphs of large but finite size. We focused on regular random graphs with soft constrained adjacency matrix eigenvalue spectra. We were able to develop a general theory describing the ${\mathcal{O}\left(1/N\right)}$ contributions to the expected eigenvalue spectrum, through the use of an infinite number of replica indices taking values in the imaginary axis [19], and building on techniques from earlier studies such as [30, 31, 32, 25].

The simplest nontrivial spectrally constrained ensembles are those in which the spectral constraint reduces to a soft constraint on the number of triangles. We quantified the behaviour of such systems, which following [26] we have come to regard as the ‘harmonic oscillators’ of loopy graph ensembles, and showed how they allow for fine tuning of their average clustering coefficients. A limitation on their use as null models for regular graphs with nontrivial clustering is that there is a maximum achievable clustering coefficient, whose value depends on the size of the graph, beyond which the ensemble undergoes a transition into a new phase, where high clustering levels are achieved by the graph fracturing into extensively many disconnected cliques. We presented a precise analytic estimate for an upper bound on the maximum clustering coefficient that is achievable without fracturing of the graph. We also showed how the general theory applies to other spectrally constrained ensembles, such as those where both the number of triangles and the number of squares are controlled, and to ensembles where the spectral constraint reduces to a count of the number of adjacency matrix eigenvalues in a given interval. We carried out numerical simulations via MCMC processes based on edge swaps with nontrivial acceptance probabilities, which are themselves generally nontrivial in view of the need to recompute eigenvalue spectra after each accepted move. In comparing triangle counts and spectra, we found excellent agreement between the theoretical predictions and the MCMC measurements in all cases, provided we remain in the parameter regime where higher orders in $N$ of the generating function are not yet important.

The most natural generalization of the presently studied family of models would be to extend the imaginary replica approach to sparse graphs with an arbitrary degree distribution $p(k)$. Preliminary numerical simulations show that these non-regular graph ensembles behave in a very similar way to what has been observed for regular graphs. In addition it would be interesting to explore further the possibility of controlling short loops in finitely connected graphs without this being realized microscopically by such graphs fracturing into extensively many disconnected graphlets, even at high loop densities. This would seem to require more complicated choices of the functional Lagrange parameter $\hat{\varrho}(\mu)$ than the ones studied so far, possibly including choices that scale differently with $N$. Both these directions for further research are now being explored by the authors.

Acknowledgement

FAL gratefully acknowledges financial support through a scholarship from Conacyt (Mexico). The authors are grateful to Alexander Mozeika and Diego Vidal Cruzprieto for very valuable discussions.