Exactly Solvable Random Graph Ensemble with Extensively Many Short Cycles

TL;DR

This paper introduces a new ensemble of 2-regular random graphs with tunable short cycle distributions, analyzing their phase transitions, spectral properties, and equivalence of different statistical formulations.

Contribution

It presents an exactly solvable model of random graphs with many short cycles, including a detailed phase diagram and spectral density analysis.

Findings

Identifies a second order phase transition between connected and disconnected phases.

Derives the spectral density with discrete eigenvalues for short cycles.

Shows equivalence of canonical and grand canonical ensembles in the thermodynamic limit.

Abstract

We introduce and analyse ensembles of 2-regular random graphs with a tuneable distribution of short cycles. The phenomenology of these graphs depends critically on the scaling of the ensembles' control parameters relative to the number of nodes. A phase diagram is presented, showing a second order phase transition from a connected to a disconnected phase. We study both the canonical formulation, where the size is large but fixed, and the grand canonical formulation, where the size is sampled from a discrete distribution, and show their equivalence in the thermodynamical limit. We also compute analytically the spectral density, which consists of a discrete set of isolated eigenvalues, representing short cycles, and a continuous part, representing cycles of diverging size.

Equations

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

$$Z_N(\mathbf{\alpha })$$

$$m_{\ell }$$

$${\varphi }_N(\mathbf{\alpha })$$

$$m_{\ell }$$

$$Z_N(\mathbf{\alpha })=\sum _{\mathbf{n}}\text{D}(\mathbf{n}){\text{rme}}^{{\sum }_{\ell =3}^K\ell {\alpha }_{\ell }{n}_{\ell }},$$

$$\text{D}(\mathbf{n})$$

$$Z_N(\mathbf{\alpha })$$

$${\varphi }_N(\mathbf{\alpha })$$

$$f_N(\omega ,\mathbf{\alpha })$$

$$\varphi (\mathbf{\alpha })=\underset{N\to \infty }{\lim }{\mathrm{extr}}_{\omega }{f}_N(\omega ,\mathbf{\alpha }).$$

$${\alpha }_{\ell }$$

$$\phi (\tilde{\mathbf{\alpha }})$$

$$1$$

$$m_{\ell }=\frac{1}{2}{\text{rme}}^{({\tilde{\alpha }}_{\ell }-\text{rmi}{\omega }_N)\ell }.$$

$$m_{\infty }$$

$$1$$

$$m_{\ell }$$

$$m_{\infty }$$

$$m_{\infty }=0:$$

$$\sum _{\ell =3}^K{\text{rme}}^{\ell {\tilde{\alpha }}_{\ell }}=2.$$

$$m_{\ell }=\frac{1}{2}{\text{rme}}^{\ell {\tilde{\alpha }}_{\ell }},m_{\infty }=1-\frac{1}{2}\sum _{\ell =3}^K{\text{rme}}^{\ell {\tilde{\alpha }}_{\ell }},$$

$$S_N$$

$$\varrho (\mu |\mathbf{A})$$

$$\varrho (\mu )=\underset{N\to \infty }{\lim }\sum _{\mathbf{A}\in {\mathcal{G}}_N}p(\mathbf{A})\varrho (\mu |\mathbf{A}).$$

$${\varrho }_{\ell }(\mu )$$

$${\varrho }_{\infty }(\mu )$$

$$\varrho (\mu |\mathbf{A})$$

$$\varrho (\mu )$$

$$Q(\mathbf{\alpha })=\sum _{N=1}^{\infty }{w}_N{Z}_N(\mathbf{\alpha }),$$

$$Q(\mathbf{\alpha })$$

Full text

1 Introduction

Models with pairwise interacting elements are ubiquitous in physics and are sufficient to capture the phenomenology of many systems, ranging from condensed matter via biology to the social sciences and informatics. The properties of the network of interactions strongly affects the properties of the system under study, and hence the analysis of networks is central in modern physics. Testing statistically whether a specific network property influences the dynamics of a system requires sampling networks where such a property is controllable, and for which the probability measure is known. This is why random graph ensembles, especially maximum entropy ones from which it is possible to sample networks systematically with controlled properties, have gained increasing popularity. They range from degree configuration models [1, 2, 3, 4, 5, 6], i.e. where the number of connections or nodes are fixed, to more complicated models where a block structure is given to the nodes in the network, as in stochastic block models [7], or other ensembles where clustering, i.e. the tendency of nodes with common neighbours to be connected, is enhanced [8, 9, 10, 11, 12, 13, 14, 15]. However, controlling analytically and numerically second or higher-order properties of networks, i.e. node properties that not only depend on first neighbours, such as the density of cycles, is still a great mathematical and analytical challenge, whose range of applications continues to grow [16, 10, 17, 18, 19, 20, 21, 22].

So far, nearly all analytical results obtained for random graph ensembles rely on the assumption of the absence of short cycles, the tree-like approximation, and we have analytical solutions only for random graphs were clustering is absent or too weak (or improbable) to be relevant. One of the first random graph ensembles in literature to include short cycles was [23], where a term depending on the number of 3-cycles of the graph was included as a modification to the well known Erdös-Rényi model (ER). This was done in order to encourage this connection transitivity in the graph. However, as was found in simulations [23] and in a more rigorous way in [8, 11], unless the graph is particularly small, this approach does not allow for a tuneable number of triangles. Depending on the values and the scaling of the parameters, the model of [23] either stays in a phase very close to the ER model, with a very slight increase in triangles, or it collapses to a condensed phase, where the complete clique has probability one. This abrupt transition was found to be a generic feature of exponential random graph models. As was shown in [24, 25], this phenomenon will be observed not only in two-parameter models like the Strauss model, but in any exponential graph ensemble that is biased such as to induce a finite number of subgraph densities.

The natural way to prevent clique formation in the condensed phase is to study random graph ensembles with hard degree constraints. Here all graphs have exactly the same degree distribution, and this distribution is chosen such that the complete clique is not an allowed state. However, this constraint makes analytical solution intractable, leaving numerical sampling from the ensemble as the only route for investigation. Examples are the Poisonnian graphs studied numerically in [26], where it was found that a triangle bias induced finite size graphs to break down into small clusters to maximize the triangle density. Regular graphs with triangle bias were numerically explored in [6], and showed similar phenomenology. However, both Poissonnian and regular graph enembles with triangle bias have so far resisted analytical solution.

In this paper we introduce an exactly solvable ensemble of 2-regular random graphs, with an exponential measure that controls the presence of short cycles up to any finite length. The imposition of 2-regularity removes the possibility of a complete clique forming, and forces the graph instead to be partitioned into a set of disconnected cycles of different lengths. This makes the ensemble analytically solvable and perfectly tuneable. The model displays a second-order transition, from a phase dominated by extensively long cycles, to a phase where only (extensively many) cycles of short lengths are present.

In section 2 we introduce and solve the model, in its canonical formulation. In section 3 we describe analytically the phases of the ensemble and the critical hyper-surface in the space of parameters; from this result we also compute analytically the spectral density of the ensemble. In section 4 we demonstrate the equivalence of the canonical and grand canonical formulations of the model, and in section 5 we show the agreement of the analytical predictions with numerical experiments. In a final discussion section we summarize the results and delineate the future directions for this research, which are twofold. The first is to relax the 2-regularity constraint of the ensemble, in order to make it more directly comparable to realistic networks. The second is to understand better recent analytical approaches to random graph ensembles that involve constraints on the number of closed paths of all lengths, which is equivalent to constraining random graphs via their spectra [27].

2 Definitions

We define a random graph ensemble over the set of undirected simple regular graphs of degree 2, which we denote by $\mathcal{G}_{N}$. Any graph in $\mathcal{G}_{N}$ is necessarily a set of disjoint cycles. The probability assigned to each graph $\mathbf{A}\in\mathcal{G}_{N}$ is chosen proportional to the exponential of a weighted sum of the number of triangles, squares, pentagons, $\ldots$ , K-cycles present in $\mathbf{A}$. We refer to this as biasing with respect of the number of short cycles. Thus

[TABLE]

Here $n_{\ell}(\mathbf{A})$ denotes the number of length-$\ell$ cycles, i.e. closed paths of length $\ell$ without backtracking and without over-counting, and $\boldsymbol{\alpha}=(\alpha_{3},\dots,\alpha_{K})\in{\rm I\!R}^{K-2}$ is a vector of control parameters. Note that isolated nodes ($\ell=1$) and dimers ($\ell=2$) cannot occur due to the degree constraint. The factors $\ell$ in (1) are included for later convenience. We are effectively biasing with respect to the total number of $\ell$-cycles starting at a given node through the introduction of the field $\alpha_{\ell}$.

The partition function $Z_{N}(\boldsymbol{\alpha})$ is given by

[TABLE]

Expression (1) defines a maximum entropy random graph ensemble with respect to the $K-2$ observables $n_{\ell}(\mathbf{A})$, whose ensemble averages are controlled by varying the parameters $\boldsymbol{\alpha}$. We choose $K$ to be a fixed number for all values of $N$. This exponential form is a particular version of the one presented in equation (1.1) of [24]. It is an ensemble where we are interested in controlling the expected values of a finite number of graph observables.

The average fraction of the $N$ nodes that will be found in an $\ell$-cycle is given by

[TABLE]

where $\langle f(\mathbf{A})\rangle=\sum_{\mathbf{A}}p(\mathbf{A})f(\mathbf{A})$. Following the statistical mechanics route, we define a generating function $\phi_{N}(\boldsymbol{\alpha})$:

[TABLE]

The main quantities of interest (3) for our graph ensemble (1) can be computed from (4) via

[TABLE]

The generator $\phi_{N}(\boldsymbol{\alpha})$ is minus the free energy density, apart from a factor $N!$ which reflects (topologically irrelevant) node label permutations. Including this factor will ensure that the limit $\phi(\boldsymbol{\alpha})=\lim_{N\to\infty}\phi_{N}(\boldsymbol{\alpha})$ exists.

3 Analytical solution

3.1 Summation over graphs

To evaluate the partition function (2) we need to perform a sum over graphs. Such sums are usually not analytically tractable, especially when the ensemble definition involves cycles, as is the case in (1). Here we are able to perform the summation by rewriting it as

[TABLE]

with $\boldsymbol{n}=(n_{3},\ldots,n_{N})\in{\rm I\!N}^{N-2}$. This decomposition reflects the fact that, in the particular case of $\mathcal{G}_{N}$, we are fortunate that each graph has to be a collection of cycles, and can therefore be identified fully by a sequence $\boldsymbol{n}=(n_{3},\dots,n_{N})$ that specifies the number of cycles of each possible length up to $N$, and a labelling of the nodes. The sum over graphs is then performed by summing over all possible sequences $\boldsymbol{n}$, keeping track of the multiplicity of each sequence via an associated density of states $\textit{D}(\boldsymbol{n})$:

[TABLE]

Apart from the condition $N=\sum_{\ell=3}^{N}\ell n_{\ell}$, this density is proportional to $N!$ but corrected for over-counting due to the indistinguishability of different length-$\ell$ cycles, giving a divisor $n_{\ell}!$, and due to the different ways one can number the nodes in each $\ell$-cycle without altering the graph ($\ell$ cyclic permutations, plus $\ell$ anti-cyclic permutations), giving a further divisor $(2\ell)^{n_{\ell}}$. Using the integral form of the Kronecker delta $\delta_{nm}=\int_{-\pi}^{\pi}(\rmd\omega/2\pi)\rme^{\rmi\omega(n-m)}$, we can thus write the partition function as

[TABLE]

From this, in combination with (4), we infer that

[TABLE]

with

[TABLE]

An exact expression for (9), valid for any finite $N$, would require to perform the integral in it. Instead, we proceed in the usual way as in statistical physics. We look at the thermodynamic limit, focusing then on $\phi(\boldsymbol{\alpha})=\lim_{N\to\infty}\phi_{N}(\boldsymbol{\alpha})$. This will allow us to calculate the asymptotic expressions for (3), which should differ from the finite size values by ${\cal O}(1/N)$ corrections.

The limit $N\to\infty$ of (9) can now be obtained by evaluating the integral over $\omega$ in (8) via steepest descent:

[TABLE]

The extremum is found by solving $\partial f(\omega,\boldsymbol{\alpha})/\partial\omega=0$.

3.2 Scaling with $N$ of control parameters

We observe that for finite $\{\alpha_{\ell}\}$ our model cannot exhibit nonzero cycle densities $m_{\ell}$ in the infinite size limit, since the $\boldsymbol{\alpha}$-dependent term in (10) vanishes for $N\to\infty$. We are therefore led to redefining the parameters $\boldsymbol{\alpha}$ with a size dependent shift,

[TABLE]

where $\tilde{\alpha}_{\ell}={\cal O}(1)$. An intuitive explanation for this scaling is presented in section 4. We denote the vector of shifted ${\cal O}(1)$ control parameters by $\tilde{\boldsymbol{\alpha}}=(\tilde{\alpha}_{3},\dots,\tilde{\alpha}_{K})$, and we define $\phi_{N}(\boldsymbol{\alpha})=\varphi_{N}(\tilde{\boldsymbol{\alpha}})$. This implies that for $N\to\infty$ we will have $m_{\ell}=\partial\varphi(\tilde{\boldsymbol{\alpha}})/\partial\tilde{\alpha}_{\ell}$, in which now

[TABLE]

Differentiation of this latter expression reveals that the value $\omega_{N}$ at the extremum is to be solved from

[TABLE]

and that the asymptotic values of the observables $m_{\ell}$ are subsequently given by

[TABLE]

This last identity, in combination with (14), prompts us to introduce $m_{\infty}=1-\sum_{\ell\leq K}m_{\ell}\in[0,1]$, which gives the fraction of the nodes that are not in cycles of length $K$ or less. It is for $N\to\infty$ apparently given by

[TABLE]

It follows from (15), that the physical saddle point $\omega$, after contour deformation, must be purely imaginary. We switch accordingly to the new variable $x=\rme^{-\rmi\omega}\in{\rm I\!R}_{0}^{+}$, in terms of which our equations become:

[TABLE]

3.3 Phase phenomenology of the ensemble

We will now demonstrate that the solutions to the coupled equations (17,18) give rise to two phases of our graph ensemble. A disconnected phase is characterized by the fact that all nodes are typically in cycles of length $K$ or less, so $m_{\infty}=0$. A second phase, the connected phase, is characterized by finding a finite fraction of the nodes in longer cycles, so here $m_{\infty}>0$. The transition separating the phases is marked by bifurcation of $m_{\infty}>0$ solutions.

If $\lim_{N\to\infty}x_{N}=x<1$, the second term of (17) vanishes for $N\to\infty$, and we immediately obtain $m_{\infty}=0$. Hence we are in the disconnected phase, and here the asymptotic observables $m_{\ell}$ are simply found by solving

[TABLE]

The condition $x<1$ for this solution to exist will be met for large values of $\{\tilde{\alpha}_{\ell}\}$. Upon reducing the control parameters $\{\tilde{\alpha}_{\ell}\}$, the value of $x$ will increase, and a transition to the connected phase occurs exactly when $x=1$. This happens at the critical manifold in the $K\!-\!2$ dimensional parameter space, defined by validity of

[TABLE]

To confirm the equations of the connected phase, we need to investigate how the solution $x_{N}$ of (17) scales with $N$ as we approach $x=1$. Substituting $x_{N}=1-\xi/N$, expanding (17) in $N$, and taking the limit $N\to\infty$ gives

[TABLE]

and the link between $\xi$ and $m_{\infty}$ is $m_{\infty}=(1\!-\!\rme^{-\xi})/2\xi$.

It turns out that all cycles of finite length $L>K$ will always have vanishing density for $N\to\infty$. This can be seen simply by replacing $K\to L$ in the previous analysis, but with $\alpha_{\ell}=0$ for all $K<\ell\leq L$. The newly added control parameters with $\ell>K$ will give $\tilde{\alpha}_{\ell}=-\ell^{-1}\log N$, and hence $m_{\ell}=\lim_{N\to\infty}\frac{1}{2}x^{\ell}\rme^{\ell\tilde{\alpha}_{\ell}}=\lim_{N\to\infty}\frac{1}{2}x^{\ell}/N=0$, in both phases. We knew that in the disconnected phase all nodes will typically be in the controlled short cycles of length $K$ or less. We may now conclude that, in the connected phase, those nodes that are not in the controlled short cycles (the fraction $m_{\infty}>0$) will typically be found in cycles of diverging length.

The ensemble’s Shannon entropy [28] is given by

[TABLE]

Since $\sum_{\ell=3}^{K}(m_{\ell}/\ell)\leq\frac{1}{3}\sum_{\ell=3}^{K}m_{\ell}\leq\frac{1}{3}$, the leading order will for large $N$ always scale as $N\log(N)$, and be bounded according to $S_{N}\geq\frac{2}{3}N\log(N)+{\cal O}(N)$, but with a reduced prefactor if we increase the fraction of nodes in short cycles. The lower bound is achieved in the disconnected phase, when $m_{3}=1$ and $m_{\ell>3}=0$.

3.4 Spectral densities of adjacency matrices

A graph can be represented uniquely by its adjacency matrix $\{A_{ij}\}$, where $A_{ij}\in\{0,1\}$, and $A_{ij}=1$ if and only if there is a link from $j$ to $i$. The set $\mathcal{G}_{N}$ contains only simple nondirected graphs, so our adjacency matrices are symmetric and with zero diagonal elements. The eigenvalue density of the adjacency matrix of a graph $\mathbf{A}$,

[TABLE]

contains valuable information on the statistics of cycles in the graph. Here the sum runs over the set of (real) eigenvalues $\{\mu_{i}(\mathbf{A})\}_{i=1,\dots,N}$ of $\mathbf{A}$, taking into account multiplicities. For instance, the number of closed paths in $\mathbf{A}$ is proportional to $\int\!\rmd\mu~{}\varrho(\mu|\mathbf{A})\mu^{\ell}$. Our main quantity of interest will be the expected density, averaged over the ensemble probabilities (1), in the infinite size limit,

[TABLE]

The adjacency matrix of a graph that consists of a single cycle of length $\ell$ has the Toeplitz form, and is therefore diagonalized trivially, leading to the density

[TABLE]

If the cycle length $\ell$ diverges, this density becomes continuous (in a distributional sense), see e.g. [29],

[TABLE]

The set of eigenvalues for each $\mathbf{A}\in\mathcal{G}_{N}$ will just be the union of all the sets of eigenvalues of the disjoint cycles of which it is composed, taking multiplicities into account:

[TABLE]

Upon averaging over the ensemble, using (3) and our earlier observation that for $N\to\infty$ the fraction of nodes in cycles of finite length $L>K$ vanishes, we immediately obtain the asymptotic ensemble-averaged spectrum corresponding to (1), expressed in terms of (26) and (27):

[TABLE]

Since we are working with regular graphs, we can immediately recover the spectrum of the Laplacian operator ($\mathbf{L}=2\mathbf{I}-\mathbf{A}$) by the change of variable $\mu\to 2-\lambda$.

4 Grand Canonical approach

Within the canonical approach one finds that, if $N$ is sufficiently large, graphs generated randomly from (1) will all display the same values of the main intensive quantities, such as the fraction of $\ell$-cycles (modulo finite size fluctuations). We expect a similar claim to hold if we sample randomly both the graphs and the number $N$ of nodes, i.e. if we work with grand canonical graph ensembles. The grand partition function of our ensemble with weights $w_{N}=\rme^{-\mu N}/N!$ (where $\mu>0$) is given by

[TABLE]

with $Z_{N}(\boldsymbol{\alpha})$ defined in (2). The divisor $N!$ in $w_{N}$ will simplify our calculation, without losing the benefits of the thermodynamic limit (since we will find that for $\mu\to 0$ the expected system size still diverges). Direct calculation of $Q(\boldsymbol{\alpha})$ now circumvents the integration over $\omega$:

[TABLE]

where we used $\sum_{\ell>0}x^{\ell}/\ell=-\log(1\!-\!x)$. From $Q(\boldsymbol{\alpha})$ we obtain, in turn, the grand potential $\Omega(\boldsymbol{\alpha})=-\log Q(\boldsymbol{\alpha})$:

[TABLE]

Its partial derivatives with respect to $\mu$ and $\boldsymbol{\alpha}$ yield the average system size, via $\langle N\rangle=\partial\Omega(\boldsymbol{\alpha})/\partial\mu$, and the average number of length-$\ell$ cycles (for $\ell=3,\ldots,K$), via $\langle n_{\ell}(\mathbf{A})\rangle=-\ell^{-1}\partial\Omega(\boldsymbol{\alpha})/\partial\alpha_{\ell}$. We thereby find that

[TABLE]

and

[TABLE]

Clearly, $\langle N\rangle$ diverges for $\mu\to 0$, which gives our thermodynamic limit. In this limit we can then work out for $\ell\in\{3,\ldots,K\}$ the ratios

[TABLE]

Similar to the canonical case, any $\mu$-independent $\boldsymbol{\alpha}$ will asymptotically always yield a vanishing fraction of nodes in cycles of length $\ell\leq K$. It is clear from (34) that without a re-parametrization, the expected value of $\ell$-cycles only increases exponentially with $\alpha_{\ell}$ in the thermodynamic limit. We need to re-parametrize in such a way that the expected number of $\ell$-cycles increases as the expected system size increases. The re-parametrization required to obtain a non-trivial thermodynamic limit is $\alpha_{\ell}=\tilde{\alpha}_{\ell}+\ell^{-1}\log\langle N\rangle$. Upon following this prescription, we then reproduce the canonical result

[TABLE]

and our expression (33) for $\langle N\rangle$ now becomes

[TABLE]

The re-parametrization of $\boldsymbol{\alpha}$ now depends on $\boldsymbol{\alpha}$ itself, via $\langle N\rangle$, and has to be consistent with a nonnegative value for (37), i.e. with $\frac{1}{2}\sum_{\ell=3}^{K}\rme^{(\tilde{\alpha}_{\ell}-\mu)\ell}\leq 1$. Expression (36) gives us the physical interpretation $\sum_{\ell=3}^{K}\ell\langle n_{\ell}\rangle/\langle N\rangle\leq 1$. In the limit $\mu\downarrow 0$ the condition becomes $\sum_{\ell=3}^{K}\rme^{\ell\tilde{\alpha}_{\ell}}\leq 2$. In the case of inequality we have $\sum_{\ell=3}^{K}\ell\langle n_{\ell}\rangle/\langle N\rangle<1$, so we are in the connected phase. The case of equality reproduces our earlier phase transition condition (21) and we enter the disconnected phase; here the thermodynamic limit is reached already for nonzero $\mu$, and we can again recover our canonical equations, with $\exp(-\mu)$ now playing the role of the canonical order parameter $x$.

5 Numerical simulations

Calculating $Z_{N}(\boldsymbol{\alpha})$ by numerical enumeration for nontrivial values of $N$ is not a realistic option, since the size of the set $\mathcal{G}_{N}$ grows super-exponentially with $N$. Instead, to test our theoretical predictions we have sampled graphs from the ensemble (1) using the Markov Chain Monte Carlo (MCMC) method described in e.g. [30] or [6]. Starting from an arbitrary 2-regular $N$-node graph, this stochastic process is based on executing repeated (degree-preserving) edge swap moves with appropriate nontrivial move acceptance probabilities, constructed such that the Markov chain’s equilibrium distribution is the target measure (1). In each simulation experiment, the MCMC process was first run for $10^{5}$ to $10^{6}$ accepted moves per link, and equilibration was confirmed by measuring the Hamming distance between the instantaneous and the initial state. After this randomization stage, the instantaneous state $\mathbf{A}$ arrived at by the chain was defined to be our graph sample. We have limited our simulations to ensembles with $K=3$ and $K=4$. The degree of equilibration achieved by the MCMC during a run of $10^{5}$ accepted moves per link is illustrated in Figure 1, where we show typical evolution curves of the order parameter $m_{3}$ during the stochastic process.

For $K=3$ we have just one control parameter $\tilde{\alpha}_{3}$, and the order parameter is the fraction $m_{3}$ of nodes in triangles. The theory claims that, for large $N$, the graphs from our ensemble will be collections of triangles and large rings. The key equations (20,21,22) reduce to the following predictions, with $\tilde{\alpha}_{c}=\frac{1}{3}\log(2)\approx 0.23105\dots$:

[TABLE]

Numerical simulations with sizes $N=1000$ and $N=5000$ show excellent agreement with these predictions, as shown in Figure 2, both in terms of the values of $m_{3}$ and in terms of the location of the transition.

For $K=4$ we have two control parameters, $\tilde{\alpha}_{3}$ and $\tilde{\alpha}_{4}$, and the theory claims that for large $N$ the graphs from our ensemble will now be collections of triangles, squares and large rings. Here the key equations (20,21,22) predict that the transition line in parameter space is given by $\rme^{3\tilde{\alpha}_{3}}+\rme^{4\tilde{\alpha}_{4}}=2$, and that the fractions $m_{3}$ and $m_{4}$ of nodes found in triangles and squares, respectively, are solved (together with the auxiliary order parameter $x$, in the disconnected phase) from:

[TABLE]

Figure 3 (left panel) shows the resulting predicted phase diagram in the $(\tilde{\alpha}_{3},\tilde{\alpha}_{4})$ plane. The mapping $(\tilde{\alpha}_{3},\tilde{\alpha}_{4})\mapsto(m_{3},m_{4})$ will map the lower region of the phase diagram (the connected phase) to the interior of the triangle $m_{3}+m_{4}<1$ in the right panel of Figure 3. The upper region of the phase diagram on the left (the disconnected phase), including the critical line, will be mapped to the line $m_{3}+m_{4}=1$ in the right panel. To test also these predictions against numerical simulations, we have chosen multiple points $(\tilde{\alpha}_{3},\tilde{\alpha}_{4})$ in both regions of the phase diagram, grouped such that the predicted values of $m_{\infty}=1-m_{3}-m_{4}$ were always in the set $\{0.25,0.5,0.75\}$. The prediction would therefore be that in the $(m_{3},m_{4})$ plane these groups of points should be found on the lines $m_{3}+m_{4}=1-m_{\infty}$. Upon measuring the fractions $m_{3}$ and $m_{4}$ via MCMC in the corresponding graph ensembles, these predictions are once more validated convincingly. See Figure 3.

Graphs sampled from our ensemble with $K=4$ do indeed typically consist of controlled numbers of triangles and squares, and a long ring. Figure 4 shows examples of such graphs, obtained via MCMC, together with the eigenvalue spectra of their adjacency matrices (obtained by numerical diagonalization). Also the observed spectra agree with the corresponding theoretical predictions (29).

6 Discussion

In this paper we presented an analytical solution for an exponential random graph ensemble with a controllable density of short cycles. Whereas one would normally not expect such non-treelike graph ensembles to be solvable, here this is possible as a consequence of imposing a local degree constraint of strict 2-regularity. We found a second order phase transition, which separates a connected phase with large and small cycles from a disconnected phase where the graphs are typically formed only of extensively many short cycles. The short cycles appear in controlled proportions, for which we found analytical expressions in terms of the ensemble’s parameters. We also derived an analytical expression for the critical submanifold in the phase diagram, and for the expected eigenvalue spectrum of the graphs’ adjacency matrices.

We analysed both the canonical and the grand canonical formulation of the ensemble. The canonical version was solved via steepest descent integration. In the grand canonical version one avoids steepest descent integration, but (as always) the chemical potential takes over the role of the steepest descent integration variable of the canonical version. In the thermodynamic limit, the canonical and grand canonical routes result in identical equations. These equations are found to give highly accurate predictions already for modest graph sizes, such as $N=1000$, as we confirmed in numerical simulations.

The parameter $K$ represents the largest cycle length that is controlled in our model. For $K=3$ one controls only the number of triangles, and our ensemble becomes similar to that of Strauss [23, 11, 8] with average degree two. The remaining difference is that in the Strauss model the average degree is imposed implicitly via an overall ‘soft’ constraint, while in the present model all degree values are imposed as local ‘hard’ constraints. Due to this difference, the degeneration of the Strauss model to a phase where the complete clique has probability one (so the number of triangles can no longer be tuned) is avoided in the present ensemble. The complete clique is simply no longer an allowed configuration, and hence the number of triangles becomes fully tuneable, if the model parameters scale appropriately with the system size. In addition, in [25] it is shown that the ’soft’ version of our model would have a phase diagram reminiscent of ours. In both cases the sign of a linear combination of functions of the parameters determines the phase of the ensemble. However, in the ’soft’ case of [25] there is a transition from an almost ER-like phase to a clique, while our model exhibits tuneability of the densities in both phases.

As we mentioned before, the generalization of the present model to other degree distributions has been studied for the Poissonian and the $q$-regular cases. Numerical explorations for 3-regular versions with $K=4$ have been reported in [6], and show phenomenology similar to that found here for $q=2$. In particular, one again observes a disconnected phase for large values of $\alpha_{3}$ and $\alpha_{4}$. One could have thought that the phenomenology of our model, in particular the emergence of a large number of small clusters, is specific to the simplifications induced by the 2-regularity condition, but this is not the case. To emphasize this fact, we redid simulations in the same fashion as in Figure 4, but now for 3-regular as opposed to 2-regular graphs, where analytical solution along the lines followed for $q=2$ is no longer feasible. The results are shown in Figure 5. It is clear that as the bias towards increased numbers of triangles and/or squares is increased, the graph breaks down into small regular graphlets that maximize the cycle density per node. For Poissonian graphs, simulations revealed in [26] that, upon boosting triangles, they also break down into small graphlets, similar to the present model. Also the effect of boosting triangles on the spectrum of the adjacency matrix was similar to what we observed here. However, since the small graphlets that appear in Poissonian models are in general different from isolated triangles and from each other, the associated eigenvalues are described by a different distribution. Similarly, the spectrum of the large component changes in a more complicated way than just by scaling down. Yet, overall we find a similar phenomenology. In fact, our present analysis predicts that the parameters in [26] would need a scaling with $N$, in order for the transition not to be a finite size effect but to persist in the $N\to\infty$ limit.

We could also combine our present model with the Erdös-Rènyi ensemble, to produce connected random graphs with a varying number of short cycles. Again, while the phenomenology of such variations could be explored via simulations, it is not clear how one would be able to obtain analytical solutions without the benefit of locally tree-like topology.

In our view, the main merit of the present model is that its analytical solution helps us understand more complicated ‘loopy’ graph ensembles. We are aware that the analytical route taken in this case is surely impossible for other models. Nevertheless, it provides an explicit analytical solution that reproduces the main features of non-treelike random graph ensembles with hard degree constraints. It helps us understand phenomenology that had so far only been studied numerically. It can also serve as a benchmark model against which more general solution strategies for non-treelike random graphs can be tested, such as [27], which deals with spectrally constrained maximum entropy graph ensembles. The moments of a graph’s spectral density are related to its numbers of cycles, via the traces of powers of the adjacency matrix. In fact, the present model is a special case of the family of ensembles studied in [27], from which it can be obtained by choosing 2-regular degrees and an appropriate polynomial functional Lagrange parameter. The analytical and numerical results of this paper suggest that, to obtain phase transitions, the functional Lagrange parameters in spectrally constrained maximum entropy graph ensembles [27] may need to have a specific scaling with the system size.

Acknowledgement

FAL gratefully acknowledges financial support through a scholarship from Conacyt (Mexico).