Approximating the cumulant generating function of triangles in the Erd\"os-R\'enyi random graph

We study the pressure of the"edge-triangle model", which is equivalent to the cumulant generating function of triangles in the Erd\"os-R\'enyi random graph. By analyzing finite graphs of increasing volume, as well as the graphon variational problem in the infinite volume limit, we locate a curve in the parameter space where a one-step replica symmetry breaking transition occurs. Sampling a large graph in the broken symmetry phase is well described by a graphon with a structure very close to the one of an equi-bipartite graph.


Introduction
Sampling random graphs with prescribed macroscopic properties (such as a given density of certain subgraphs) received considerable attention in recent years. From a statistical physics perspective, one can think of two procedures: • the micro-canonical ensemble, where the sampling is performed with a uniform distribution over the set of all graphs that satisfy the macroscopic constraint exactly; • the canonical ensemble, where the sampling is done with respect to a larger set of graphs that satisfies the macroscopic constraint only on average.
We shall discuss here the simplest non trivial case, i.e. the constraint is on the number of edges and the number of triangles in the graph. In the micro-canonical ensemble these numbers are prescribed exactly. The canonical ensemble is instead provided by the so-called edge-triangle model, which is defined by the Boltzamnn-Gibbs distribution in which one tunes the average density of edges and the average density of triangles by varying the corresponding conjugate parameters. The edge-triangle model is in turn the simplest example of the more general exponential random graph class of models, in which one introduces several parameters to control the density of an arbitrary set of subgraphs. Whereas the equivalence in the thermodynamic limit of micro-canonical and canonical ensemble is true for several physical systems of interest (often the system is then studied in the canonical ensemble, that is usually more analytically tractable than the micro-canonical one), for random graphs it has been shown that such equivalence can not be taken for granted. In particular [18] identified a region of values for the densities of edges and triangles where there is ensemble inequivalence, as measured by a positive relative entropy between the micro-canonical and canonical measures. In this paper we address the following problem: How does a large graph look like when it is sampled from the edge-triangle model (i.e. imposing some given average values for the densities of edges and triangle )? Does the sampling from the edge-triangle model give the same result of the sampling with respect to the microcanonical ensemble (i.e. imposing some given exact values of edge and triangle densities)?

The edge-triangle model and the Erdös-Rényi random graph
To define the setting, let us consider a graph with n vertices, that we identify with the elements of the set [n] = {1, 2, 3, . . . , n}. We shall describe the graph using its adjacency matrix x = (x i,j ) i,j∈ [n] , defined as follows: the entry x i,j = 1 if the edge connecting vertex i with vertex j is present, and x i,j = 0 otherwise. Since the graphs considered in this paper are undirected and without loops, the adjacency matrices will be always symmetric, with 0 or 1 entries and zeros on the diagonal. We will denote by X n the set of adjacency matrices of size n, and we also define for later use X = ∪ n≥2 X n the set of adjacency matrices of all sizes. The number of edges and triangles in a graph represented by a matrix x ∈ X n is given by The previous quantities result in random variables if the graph is a random graph, i.e. if it is sampled from the set of 2 ( n 2 ) undirected simple graphs with n vertices according to some probability distribution. The simplest possible distribution is the so-called Erdös-Rényi model, where each pair of vertices is connected with probability p > 0, independently of the other pairs. Thus, in the Erdös-Rényi case the entries of the adjacency matrix, x = (x i,j ) i,j∈ [n] , form a set of independent identically distributed (i.i.d.) Bernoulli variables with P(x i,j = 1) = p. In spite of the simple probabilistic set-up the large deviation principle for the Erdös-Rényi graph is far from simple and has been developed only in recent years [14,16,13,25,12,32,17]. In particular, it has been found that the large deviation function may be non-convex.
The exponential random graph model is devised to enhance or decrease the probability of specific geometric structures in the random graph. Here we define the edge-triangle model that involves only triangles and edges [26]. Let β 1 , β 2 ∈ R, then the probability of a graph with adjacency matrix x ∈ X n in the edge-triangle model is given by: where Z n (β 1 , β 2 ) is the partition function, i.e. the normalizing factor The factors 6 and 2 are conventional in the definition and account for the permutations of 3 vertices of a triangle and the 2 vertices of an edge. The Erdös-Rényi model with paramemeter 0 < p < 1 is embedded in the edge-triangle model, since its distribution: can be obtained from (2) by setting β 1 = h p /2 and β 2 = 0. If β 2 > 0 the probability of finding triangles is enhanced with respect to the Erdös-Rényi case, while it is decreased if β 2 < 0. In the limiting case β 2 → −∞ the edge-triangle model (2) gives zero probability to graphs containing triangles.
A key quantity in the study of the thermodynamic properties is the pressure, that at finite volume is defined as ψ n (β 1 , β 2 ) = 1 n 2 ln Z n (β 1 , β 2 ).
By taking derivatives with respect to the model parameters one computes the averages of the edge density e = 2E n /n 2 and of the triangle density t = 6 T n /n 3 : where · n denotes expectation w.r.t. the measure ν n defined in (2). We shall be interested in the behavior of very large graphs which mathematically is described by the (thermodynamic) limit n → ∞. General convexity arguments imply that the thermodynamic limit is well defined, so that the infinite volume pressure exists and by Lebesgue dominated convergence limits and derivatives can be interchanged so that the relation (6) gives in the thermodynamic limit where e = lim n→∞ e n . We shall work with the parametrization (β 1 , β 2 ) = (h p /2, α/6) where we recall h p = ln p 1−p and 0 < p < 1. In this way the pressure of the edge-triangle model can be read as the cumulant generating function of the number of triangles in the Erdös-Rényi random graph. In other words, defining where · ER n denotes the expectation w.r.t. the measure ν ER n , by a simple computation one can show that Thus, studying the cumulant generating function of the number of triangles in the Erdös-Rényi random graph is equivalent to studying the pressure of the edge-triangle model.

Ensemble inequivalence
As the pressure ψ n (β 1 , β 2 ) is the crucial quantity in the canonical ensemble, the entropy s n (e, t) defined as is the important quantity in the microcanical ensemble. An heuristic application of the Laplace method would imply ψ(β 1 , β 2 ) = lim n→∞ 1 n 2 ln e n 2 (β1e+β2t−sn(e,t)) de dt = sup i.e. in the thermodynamic limit the canonical pressure can be obtained from the microcanonical entropy by a Legendre transform. One then says that the two ensembles are thermodynamic equivalent if such correspondence also holds in the reversed direction. This is the same as requiring that the microcanonical entropy is strictly convex, i.e. the involution property of the Legendre transform for strictly convex functions. The problem of ensemble inequivalence is, in fact, more general than just thermodynamic inequivalence (see [31] for a recent account). When the correspondence via Legendre transform between pressure and entropy does not hold, the difference between canonical and microcanonical ensemble can then be probed in several ways. In this paper we shall focus on macrostate inequivalence. Namely, we ask if sampling a very large graph uniformly at random from the set of all graphs with given dentities of edges and triangles (e * , t * ) is statistically equivalent to sampling a very large graph from the edge-triangle model with parameter values β 1 (e * , t * ), β 2 (e * , t * ) that are obtained by inverting the relations (8) with e = e * and t = t * . In the thermodynamic limit graphs are described by the notion of graphon, which is discussed below. Thus the equivalence that we consider amounts to ask if the two sampling procedures produce the same graphon.
The sampling from the micro-canonical ensemble has been investigated for instance in [28,29,30,22,3]. It has been found that the structure of graphs drawn from the microcanonical ensemble is very rich and may vary a lot as a function of the number of prescribed edges and triangles. For instance, for a choice (e, t) such that e = 1 2 + with ∈ ( l−2 2l , l−1 2l+2 ) with l ∈ N \ {1} and t on the scallopy curve, the vertex set of a graph drawn from the microcanonical ensemble can be partitioned into subsets ( − 1 of them of the same size and the last of different size). The graph has the form of a complete -partite graph on these pieces, plus some additional edges in the last piece that create no additional triangles [18].
To the best of our knowledge, sampling from the canonical ensemble has been investigated only in a limited region of the parameters (β 1 , β 2 ), see the review of know results in Section 2.1. It is the aim of this paper to conduct a systematic exploration of the full parameter space by means of numerical simulations.

Main results and paper organization
We investigate the sampling from the canonical ensemble by studying the pressure of the edge-triangle model, or equivalently the cumulant generating function µ p (α) of triangles in the Erdös-Rényi model with parameter p. We perform numerical simulations for finite graphs and compare them to the variational formulation describing the infinite volume. We shall collect multiple evidences that the structure of graphs in the canonical ensemble has only two possibilities: it is either the constant graphon describing the Erdös-Rényi graph (i.e. independent edges, yet with a modified parameter for the probability of edges accounting for the imposed number of triangles) or it is the graphon describing the 1-step replica symmetric breaking solution (generalizing the bipartite random graphs that is know to be the exact solution for α → −∞). By means of different numerical analysis ("cloning" method for a direct measurement and "gradient" method for the solution of the pressure variational problem) we shall identify a curve α c (p) in the plane (p, α) separating these two regimes called, respectively, the replica symmetric phase and the replica symmetry broken phase. We do not have a proof that replica symmetry broken phase α < α c (p) is entirely described by the 1-step replica symmetric breaking solution. However, in contrast to the microcanonical sampling, our numerical analysis suggests that in the description of the canonical sampling no higher level of replica symmetry breaking is required. As a consequence of the numerical analysis the value α c (p) may be identified as a bona-fide critical value for a 1-step replica symmetry breaking transition.
The paper is structured as follows: • in section 2 we review the variational formulation of the pressure. We recall the results that are known in the literature for the solution of the variational problem and describe the 1-step replica symmetry breaking solution.
• In section 3 we present the numerical analysis of the finite volume pressure based on the 'cloning method', which is a population dynamics algorithm.
• In section 4 we solve a discretized version of the variational problem in the infinite volume by a gradient projection method. Here we do not fix a-priori a specific structure for the optimal graphon.
• In section 5 we solve the variational problem restricted to a specific class of graphons, those corresponding to the 1-step replica symmetry breaking solution.
We will argue that the results of sections 4 and 5 coincide (within numeral accuracy) and they are well approximated by the finite volume direct measurements of section 3.

Variational formulation 2.1 Review of known results
The theory of graph limits [24,23,8,9,10] relies on the notion of graphon, which describes a random graph in the limit n → ∞. A graphon is defined as a bounded Borel measurable function f : The idea behind this definition is a mapping of a graph to the unitary square: intuitively the interval [0, 1] represents a continuum of vertices and f (x, y) is associated to the probability of connecting with an edge two vertices x and y. For example, the graphon describing the Erdös-Rényi random graph with parameter p is the constant function f identically equal to p. The set of graphons is denoted by W. On this set an equivalence relation is introduced according to which f, g ∈ W are equivalent if there exists a bijection σ : [0, 1] → [0, 1] such that σ and σ −1 are Borel measurable and preserve the Lebesgue measure, such that g(x, y) = f (σ(x), σ(y)). The set of equivalence classes is denoted by W, while f is the equivalence class containing f ∈ W.
In reviewing the known results in the thermodynamic limit, we follow here [15] (for a more general overview of large deviations for random graph see [12]). The thermodynamic limit of the pressure of the edge-triangle model is given by [15,Theorem 3 where t(f ), the density of triangles in f , is and the entropic term is Here f is a representative element of the equivalence classf and I p (u) for u ∈ [0, 1] denotes the Bernoulli relative entropy Then, from (10) and (13) we get with When the set of maximizers only consists of constant functions we say that we are in the replica symmetric phase. Conversely, if the elements of the maximizing set are non-constant functions, then we say that we are in the replica symmetry breaking phase. The infinite-dimensional variational problem (17), involving the non-linear functional H, has been analytically solved only in a region of the parameter values. In particular [15,Theorem 6.2] proves that for all 0 < p < 1 the system is in the replica symmetric phase for α > −2. The variational problem is then reduced to a scalar one, see [15,Theorem 4 where u * (α) is the optimizer that solves the fixed-point equation: From this, one infers that for α > −2 a graph sampled from the edge-triangle model in the limit n → ∞ will look like an Erdös-Rényi random graph with parameter u * (α), i.e. edges are independent from each others and present with a probability u * (α) (we refer to [5,15] for the precise statement). As for the region −∞ < α ≤ −2 the solution of the variational problem (18) is unknown. The Euler-Lagrange equation giving the stationarity condition are given by the following equation, which is the generalization of (20) to the case of non constant functions [15, Theorem 6.1]: It has been proved that for α small enough, a graph sampled from the edge-triangle model no longer resembles an Erdös-Rényi graph. In particular [15,Theorem 6.3] show that for α small enough (and for any value of p) the functional H(f ) is not maximized at any constant function. This result is based on the fact [15, Theorem 7.1] that one actually proves that in the limit α → −∞, the solution of the variational problem (17) is provided by the so-called equi-bipartite graphon defined as and in this limit one has lim For later use, we remark that using the constant graphon Thus, in the limit α → −∞ the replica symmetric solution is wrong by a factor 2!

A conjecture
The graphon (22) is the infinite volume correspondent of the equi-bipartite graph. In the latter the n vertices are partitioned into two disjoint sets of equal size with no edges connecting two vertices belonging to the same set. Clearly, in such a graph triangles do not exist (more generally bipartite graphs do not contain odd cycles [7,Theorem 4]). Analogously, the density of triangles in the equi-bipartite graphon (22) vanishes since t(g) = 0. Thus, as expected, the edge-triangle model in the limit α → −∞ is free of triangles.
Guided by the results of our numerical analysis that we illustrate below, we conjecture that for any fixed 0 < p < 1, there exists a unique finite and negative value α c (p) such that, crossing α c (p) from above, the optimizer of H defined in (18) switches from the constant function f α , identically equal to the solution of the fixed-point equation (20), to the graphon: where p 1 (α) and p 1 (α) are functions taking value in (0, 1) that satisfy the following conditions: We observe that lim α→−∞ g α = g, with g the graphon describing the equi-bipartite graph defined in (22). The rationale behind our conjecture is that the structure of (25) represents the simplest geometry that may emerge from the breaking of the homogeneous graphon. For the resemblance of the structure of the overlap matrix in the 1RSB solution of spin-glasses, we call this graphon the 1-step replica symmetry breaking solution.

Direct measurements for finite graphs
In this section we compute numerically the cumulant generating function of the number of triangles of a finite graph of size n. This expectation if substantially affected by events that, although rare, give a large contribution to the average defining the generating function. A standard tool for estimating the probability of rare events in the Erdös-Rényi random graph is the importance sampling technique, see for instance [6]. Here we follow an approach based on population dynamics, called 'cloning' [20,19,21,27,11,1,2]. In this section we adapt the method to a purely geometric problem, by introducing a dynamics for the graph construction.

Implementing cloning
The cloning algorithm is obtained by tilting a Monte Carlo dynamics that samples from a target distribution. In our case we would like to compute expectations w.r.t. the law of the Erdös-Rényi graph (4). We first describe the dynamical process generating the Erdös-Rényi graph and then we recall how the tilted dynamic arises in the cloning algorithm.
In the Erdös-Rényi graph each each edge is present, independently, with probability p ∈ (0, 1). To generate a graph of size n, we consider the Markov chain {X t , t ∈ N}, taking values on the set X of adjacency matrices, defined as follows. We label the n vertices in an arbitrary order. We start by selecting the first two vertices and we connect them with probability p, thereby obtaining a graph of size two. Then we select the third vertex and try to connect it to the first two vertices independently with probability p, thus obtaining a graph of size three. This procedure is repeated until the graph of size n is formed: each time a new vertex is selected, it is connected independently with probability p to each of those already visited. We stipulate that the discrete time step of this process corresponds to the attempt of adding a single new edge. Since the evolution from the graph of size i to the graph of size i + 1 requires i attempts, then the evolution starting form size two and leading to a graph of size n will require N n = n i=2 i = n 2 − 1 steps. We now introduce the tilted dynamics. Denoting by P the transition matrix of the Markov chain described above, i.e. P (x, y) = P(X t+1 = y|X t = x) for x, y ∈ X , the cumulant generating function of triangles reads where ν 0 is the initial distribution (i.e. Bern(p)). We rewrite the number of triangles for the graph of size n as where ∆T (x t , x t+1 ) = T (x t+1 ) − T (x t ) is the increment of triangles between two consecutive steps. In this way we obtain: The average in the previous equation can be also computed according to a different dynamics. Indeed, introducing the quantity and the stochastic matrix p α with elements P α (x, y) := P (x, y)e α n ∆T (x,y) we can write: As observed in [19] this representation of the average suggests a population dynamics scheme that, starting from a bunch of M initial individuals (clones) with distribution ν 0 (·), make them evolve according to the transition kernel P α (·, ·) and reproduce according to the rate k α (·). Denoting by M Nn the size of the corresponding population of clones at the final time N n , we get: The procedure is further illustrated with the example of size n = 3 in the Appendix.

Numerical results
We present here the results of the cloning algorithm for the cumulant generating function of the triangles. We fix a population of M = 7000 clones and a value of p = 0.4, and we vary the variable α and the graph size n. We denote by µ Cl n (α) the cumulant generating function returned by the cloning algorithm for the graph of size n with p = 0.4. Similarly, we shorthand µ RS (α) := µ RS 0.4 (α) and we denote the asymptotic valuesμ :=μ 0.4 andμ RS :=μ RS 0.4 . The outcomes of our simulations show that: • for large values of α, the algorithm reproduces the replica symmetric solution (i.e. µ Cl n → µ RS as n increases); • for small values of α, the algorithm provides a cumulant generating function which is independent of α within statistical fluctuations. This "constant" function approximates better and better the value of the exact solution at α = −∞ (i.e. µ Cl n →μ as n increases); • for all values of α and n, the curves furnished by the cloning algorithm are always above the curve that is obtained by taking the maximum between the replica symmetric solution and the exact solution at α = −∞ (i.e. µ Cl n ≥ max{µ RS ,μ}). The discrepancy reduces and n increases.
We observe that the intersection between the replica symmetric solution and the exact solution at α = −∞ occurs at a valueα −110, thus a substantially small value where the rate of reproduction of triangles in the cloning algorithm is very small. Despite this "extreme" situation, we see strong evidence that the  replica symmetry solution does not hold for small α's, and for values α ≤ −110 the solution furnished by the equi-bipartite graphon is a very good approximation for the curve returned by the algorithm.
We discuss our results below. Figure 1 explores the region around α = 0. Figure 1(a) shows with dots the results of the cloning algorithm for n = 3, n = 4, n = 20, n = 110 and α ∈ [−3, 2]. The picture also reports, with orange dashed-dotted and maroon dotted lines, the exact curves µ n,0.4 for n = 3 and n = 4, that can be computed explicitly. The curve µ RS p (α), obtained by solving numerically the scalar variational problem (19) characterizing the replica symmetric regime for p = 0.4, is also displayed (pink continuous line). We observe that for n = 3 and n = 4 the cloning algorithm perfectly reproduces the exact result and, as long as n grows, the curves for increasing n's settle on the curve µ RS (α).
A further check on the behavior of the algorithm is provided by Fig.1(b). It shows, as a function of α, the average density of edges e n (α) Cl in the population of clones, that is e n (α) . This quantity is an estimator of the probability p of finding an edge connecting any two vertices of the graphs produced by the cloning. By effect of the tilting, p is different from the original value p = 0.4 and is a function of α, as expected. Since in replica symmetric regime the edge-triangle model converges (in the limit n → ∞) to an Erdös-Rényi model with parameter u * (α), see Section 2, we have reported in Fig.1(b) the optimizer u * (α). The fair agreement of the values e n (α) Cl on the curve of u * (α), gives a further evidence that the cloning algorithm is well working in this regime of α and reproduces the expected behavior of the edge-triangle model. The same conclusion can be drawn form Fig.1(c) that shows the estimate of the average density of triangles provided by cloning, i.e. t n (α) and the density of triangles in the symmetric replica regime, i.e. (u * (α)) 3 .
To check if the algorithm is able to reproduce the change occurring around α −110, we have run simulations in a strongly negative region of α. We plot in Fig.2 the results of the cloning method µ Cl n (α) (dots) together with the replica symmetric cumulant generating function µ RS (α). In the same picture we have also reported the asymptotic valuesμ = −0.255 andμ RS −0.510. The left panel of Fig.2 displays the cloning data for several n values, showing the converge towards a limiting profile as n is increased. Figure 2(b) exhibits the results for the largest graph size n = 110.

Numerical solution of the variational problem
The cloning algorithm implicitly solves the variational problem (17) by producing a population of graphs that approximate the optimizing graphon. Unfortunately, it is hard to scrutinize the structure of the graphon from the adjacency matrix of the graphs, when their size is large. Thus, in order to study the graphon in the replica symmetry broken phase, we solve (17) via a numerical discretization.

Discretization
The spatial discretization of graphon can be obtained considering the set of m × m symmetric matrices {f i,j } i,j with elements in [0, 1]. However, since the graphon that solves (17) does not take the values 0 and 1 [15, Theorem 6.3], we restrict our set of matrices to the following set: where ε > 0 is a small parameter that bounds the functions away from the singularity of the logarithm in the discretization of H(f ), that we define as follows: We remark that, as it happens in the continuous setting, H m (f ) enjoys a symmetry property. Indeed, given f, g ∈ Γ m,ε we have H m (g) = H m (f ) if there exists a permutation σ over {1, . . . , m} such that g i,j = f σ(i),σ(j) . In this case we call g and f equivalent. We solve numerically the discretization of (17) by applying the Gradient Projection (GP) method with both a constant and variable steplength, the latter chosen according to the Barzilai-Borwein rules [4].

Numerical results
Being interested in the structure of the optimizer in the replica symmetry breaking region, we have solved (38) for a set of p values ranging from 0.2 to 0.6 and α below −2, varying it with unitary step. For each value of p and α we have started the iterations of the GP method from a set of 12 initial conditions, using a grid of size m = 40 and setting ε = 10 −4 . Our main findings are: a) the presence of only two different geometrical structures of maximizers: the constant ones and chessboard-like one (see Fig.3). The values of the constant maximizers turn out to be equal, within the numerical approximation, to the solution of the fixed-point equation (20); b) there exists a critical value α c (p) such that when α > α c (p) the optimizer of H m (f ) is constant whereas when α < α c (p) it assumes a chessboard-like structure.
We observe that there are different chessboard-like structures that can be reached starting from different initial conditions, as shown in Fig.3. All of them are equivalent to the discretization of the 1-RSB graphon, that by abuse of notation we also denote by g α : With p fixed, and varying α, we sought the critical value α c (p) by comparing the values of H m computed at the homogeneous solutions f α and at the chessboard-like ones g α . We found that the values of f α coincide, within our approximation and for all α, with u * (α) (the solution of the replica symmetric fixed-point equation (20)), while the values p 1 (α) and p 2 (α) occurring in g α are close to 0 and p, see Tab.1 for large negative α values. This implies that g α is close to the equi-bipartite graphon g defined in (22). The evidence that the phase transition occurs at the critical value α c (p) is given by observing . An example is given in Table 1 where, for the case p = 0.4, the change of the optimizer shows that the position of the critical value α c (p) is between −109 and −110 (we recall that α varies with unitary step).  Table 1: Numerical optimizers of (38) for p = 0.4 and α = −109, −110. We report the solution u * (α) to equation (20), the constant solution f α and the two values p 1 (α) and p 2 (α) taken by the chessboard solution g α . The last two columns show the values of H 40 (f ) that give evidence of the phase transition between −109 and −110.
Denoting by h α the global optimizer found by the GP method, in Fig.4 we represent H 40 (h α ) as a function of α for p ∈ {0.2, 0.4, 0.6}. In the replica symmetric region, i.e for α > α c (p), H 40 (h α ) is expected to approximate the solution µ RS (α) = α (u * (α)) 3 3 − I p (u * (α)) . The overlapping of the two curves above α = −455, α = −110, α = −69 and α = −47 is shown in Fig.4. Below such thresholds, that we identify as approximations of critical value α c (p), the optimal solution h α of H 40 switches from the constant one f α (approximating the fixed-point equation (20)) to a chessboard function g α . The function H 40 (h α ) is nearly flat for α < α c (p) and very close to the asymptotic value H(g) = 1 2 ln(1 − p) of µ p (α), see (23). of Section 4, here we assume that the breaking of the homogeneous graphon gives rise to a solution with the chessboard structure, see Fig.3, p1,p2 (x, y) = that we call the generalized 1-step replica symmetry breaking solution. We introduce the set R = {g (a) p1,p2 (x, y)| a, p 1 , p 2 ∈ [0, 1]} ⊂ W. This set contains the constant graphon (obtained by setting p 1 = p 2 for any a or, equivalently, by setting either a = 0 or a = 1) and, for a = 1 2 , the graphon (25), that we claim to be the solution of (17) for α < α c (p). Also, the limiting graphon (22) is contained in R.
We thus restrict the infinite dimensional problem (17) to the finite dimensional one obtained by restricting to the set R: Being interested in the phase transition, we consider this problem for α ≤ −2. Using this approach, we aim at locating the critical value α c (p) denoting the transition that we claim to occur in R between the homogeneous and the 1-step replica breaking solution. Obviously we can state that α c (p) ≤ α c (p), but we do not have any argument to assert that the transition in R is the same occurring in the whole space W and, as a consequence, that α c (p) coincides with α c (p). However, some evidence for the equality of the two critical points will be obtained in this section. The function to be maximized in (41) can be written as follows: since the density of triangles (14) in the generalized 1-step replica symmetry breaking graphons (40) is This function, that satisfies the symmetry F α (p 1 , p 2 , a) = F α (p 1 , p 2 , 1−a), can be defined by continuity up to the boundaries of [0, 1] 3 by setting I p (0) = I p (1) = 0. Therefore F α (p 1 , p 2 , a) attains its maximum on [0, 1] 3 . Due to the fact that the optimizer of (17) is bounded away from 0 and 1 [15], we assume that the coordinates (p 1 , p 2 ) of the maximum point lie in the interior of [0, 1] 2 and thus satisfies stationarity condition ∇F α (p 1 , p 2 , a) = 0, given by the following system of equuations: with p 1 , p 2 , a ∈ [0, 1] and α ≤ −2.
It is simple to check that the point (p * 1 (α), p * 2 (α), a * ) with p * 1 (α) = p * 2 (α) = u * (α) and any a * ∈ [0, 1] is a solution to equation (43) for any α. Indeed, when p 1 = p 2 equations (a) and (b) are equal and coincide with the fixed-point equation (20). From the numerical solution of system (43) as α is varied, we got evidence of the existence of a thresholdα c (p) above which (u * (α), u * (α), a * ) is the maximum (this holds for any a , being such parameter meaningless in the homogeneous case). Crossingα c (p) from above, a non-homogeneous solution with a * = 1/2 and p * 1 (α) ≈ 0, p * 2 (α) ≈ p arises; it turns out to be the maximum of the problem for α <α c (p).
We give more details on the procedure we followed in order to analyze the solutions of (43). First we observe that there are three special values of a, namely a = 0, 1 2 , 1. As we said above, a = 0, 1 correspond to the constant graphon that can be obtained also as a special case ( with p 1 = p 2 ) of a = 1 2 . Thus since the cases a = 0, 1 can be absorbed in a We solved numerically the previous equation for several values of p and α. The function G α (p 2 ) turns out to be a concave function whose maximun, that is also the solution to eq.(44), coincides with the solution that has the solution p 1 = u * (α). Indeed, the equation αx 2 = ln x(1−p) p(1−x) is equivalent to the fixed-point equation (20) that has an unique solution for α < 0.
• Case a = 1 2 . In this case the function to be maximized takes the much simpler form F α (p 1 , p 2 , 1 2 ) = α 12 (p 3 1 + 3p 1 p 2 2 ) − 1 2 (I p (p 1 ) + I p (p 2 )) with the stationarity condition: We have computed the numerical solution (p * 1 (α), p * 2 (α)) for several values of p, as in Fig.6. The left panels show that crossing a critical value α c (p) from above the transition between the constant graphon (p * 1 (α) = p * 2 (α))) and a 1-step replica symmetry breaking regime takes place. In particular, at α c (p) the solution jumps towards the point (0, p). The central and right columns of Fig.6 represent the solution below the critical value showing that, in the limit α → −∞, the values p * 1 (α) and p * 2 (α) converge to the limits 0 and p, respectively, as conjectured in (26) and (27). A further evidence of this transition is given in Fig.7 in which the curves F α p * 1 (α), p * 2 (α), 1 2 and µ RS (α) are displayed. Above the critical value α c (p) the two curves overlap thus revealing the replica symmetric phase whereas they separate below (replica breaking regime). For α <α c (p), the right panel Fig.7 shows that F α p * 1 (α), p * 2 (α), 1 2 is very close to the the asymptotic value (23) (the discrepancy vanishes as α → −∞). Figure 8 represents the density of triangles and edges corresponding to the maximizer g The former quantity is given in (42) while the latter is The jump discontinuity shown in Fig.8, located atα c (p), makes evident the existence of the transition that separates the replica symmetric phase, in which both quantities decrease for decreasing α, from the replica breaking phase. In this phase the density of edges gets close to the value p 2 , i.e. the density of the limiting equipartite graphon (22), while the density of triangle jumps towards a value close to zero, being zero the density of triangles of the same graphon. Figure 9 shows the dependence of the critical valueα c (p) from p. The data are in good agreement with a fit |α c (p)| ∼ p −2 .
We conclude our analysis by Fig.10, which displays the replica symmetric pressure (19), the 1-step replica symmety breaking pressure (41) and the pressure obtained from the discretized variational problem. We can observe that the three curves perfectly match above the critical thresholdα c (p), whereas in the subcritical phase the pressure obtained from the solution of the discretized problem agrees with the 1RSB pressure and is strictly larger than the replica symmetric pressure. Figure 10 shows that α c (p) is a singularity of α → F α (p * i (α), p * b (α), 1 2 ), at which the left and right derivatives are different.
We compute again (46) by applying the dynamics described in Sec.3 to a family of M clones. The three steps required to construct the edges of the graph G 3 are represented in Fig.11. The leafs of the tree represent the occupation variables of the three possible edges, x 1,2 , x 1,3 , x 2,3 . In the first step the edge connecting two vertices, say 1 and 2, of each clone is added with probability p. Thus, in the clone population about pM graphs have the edge (1, 2) and (1 − p)M graphs have not this edge, see the first level in Fig.11 . In the second step the edge (1, 3) is added, still with probability p, leading to four possible values for the pair (x 1,2 , x 1,3 ). Thus, the expected numbers of types are M p 2 , M p(1 − p), M p(1 − p), M (1 − p) 2 , see level 2 in Fig.11. Let us observe that in the first two steps, since ∆T (x, y) = 0, the original and tilted transition probabilities coincide: P α (x, y) = P (x, y), see (32). The situation changes in the last step. Indeed, the configuration of edges x = 11 (which means that x 1,2 = 1 and x 1,3 = 1) may evolve to y = 111, with ∆T (x, y) = 1, or to y = 110, in which case ∆T (x, y) = 0. For all the other configurations x we have ∆T (x, y) = 0. Thus, see the definition (32) of the tilted probability: being P (11, 111) = p, P (11, 110) = 1 − p and k α (11) = pe α 3 + 1 − p. Then, the probability of the paths connecting the root φ to the leafs 111, respectively 110, will be, : (48) The average in (33) can be computed as the sum over the paths from the root to the leafs (equivalently, as the sum over the leafs). Recalling that k α (11) = pe α 3 + 1 − p and observing also that k α (x) = 1 if x = 11, from (33) we have: µ 3,p (α) = 1 3 ln   P(φ, 111)k α (11) + P(φ, 110)k α (11) + x / ∈{111,110}