Ground States for Exponential Random Graphs

TL;DR

This paper introduces a perturbative approach to estimate normalization constants in exponential random graph models and provides evidence of phase transition phenomena in the edge-triangle model.

Contribution

It presents a novel perturbative method for normalization estimation and demonstrates discontinuity in model behavior at critical parameter values.

Findings

Perturbative method effectively estimates normalization constants.

Evidence of phase transition in the edge-triangle model.

Discontinuity observed along critical directions.

Abstract

We propose a perturbative method to estimate the normalization constant in exponential random graph models as the weighting parameters approach infinity. As an application, we give evidence of discontinuity in natural parametrization along the critical directions of the edge-triangle model.

Figures

Tables

Table 1. Table 1: Asymptotic comparison for perturbation around neighboring Turán graphons, where β = ( 10 , − 7.5 ) 𝛽 10 7.5 \beta=(10,-7.5) diverges along the critical direction o 1 = ( 1 , − 0.75 ) subscript 𝑜 1 1 0.75 o_{1}=(1,-0.75) . a 𝑎 a , b 𝑏 b and ψ β superscript 𝜓 𝛽 \psi^{\beta} are calculated from ( 32 ) and ( 49 ) with first order perturbation; a opt subscript 𝑎 opt a_{\text{opt}} , b opt subscript 𝑏 opt b_{\text{opt}} and ψ opt β subscript superscript 𝜓 𝛽 opt \psi^{\beta}_{\text{opt}} are based on numerical optimization for ( 45 ). The bipartite feature is associated with a bigger limiting normalization constant and is favored, matching the asymptotic predictions of Theorem III.2 .

	$a$	$b$	${\psi }^{\beta }$	$a_{\text{opt}}$	$b_{\text{opt}}$	${\psi }_{\text{opt}}^{\beta }$
bipartite	0.999999998	$0.0759$	5.0197	almost 1	0.069	5.019
tripartite	$0.9933$	$0.0000454$	5.0022	0.9943	0.000064	5.0021

Equations

$$t(H,G)=\frac{|\text{hom}(H,G)|}{|V(G){|}^{|V(H)|}}.$$

$$T^{\beta }(G_n):=\beta \cdot T(G_n)=\sum _{i=1}^d{\beta }_{i}t(H_i,G_n),$$

$$Z_n^{\beta }=\sum _{G_n\in {\mathcal{G}}_n}\exp \left(n^2{T}^{\beta }(G_n)\right).$$

$$\mathcal{W}=\{f:[0,1{]}^2\to [0,1]\text{ and }f\text{ symmetric}\}.$$

$$f^{G_{n}}(x,y)=\left\{\begin{array}[]{ll}1,&\hbox{if $(\lceil nx\rceil,\lceil ny\rceil)$ is an edge in $G_{n}$,}\\ 0,&\hbox{otherwise,}\\ \end{array}\right.\quad(x,y)\in[0,1]^{2}$$

$$d_{\square }(f,h)=\underset{S,T\subseteq [0,1]}{\sup }\left|{\int }_{S\times T}\left(f(x,y)-h(x,y)\right)dxdy\right|.$$

$${\delta }_{\square }(\tilde{f},\tilde{h})=\underset{{\sigma }_1,{\sigma }_2}{\inf }{d}_{\square }(f_{{\sigma }_1},h_{{\sigma }_2}),$$

$$t(H,f)={\int }_{[0,1{]}^k}\prod _{\{i,j\}\in E(H)}f(x_i,x_j)d{x}_1\cdots d{x}_k,$$

$$\underset{n\to \infty }{\lim }t(H,G(n,f))=t(H,f).$$

$${\mathbb{P}}_n^{\beta }(G_n)=\exp \left(n^2(T^{\beta }(G_n))-{\psi }_n^{\beta }\right),$$

$${\psi }_n^{\beta }=\frac{1}{n^2}\log Z_n^{\beta }.$$

$$I(u)=\frac{1}{2}u\log u+\frac{1}{2}(1-u)\log (1-u),$$

$$I(\tilde{f})={\int }_{[0,1{]}^2}I(f(x,y))dxdy,$$

$${\delta }_{\square }({\tilde{f}}^{G_n},{\mathcal{F}}^{\beta })\to 0\; in\; probability\; asn\to \infty .$$

$${\Delta }_{H,r,s}f(x_r,x_s):={\int }_{[0,1{]}^{|V(H)\text{}\{r,s\}|}}\prod _{\begin{array}{c}(r^{\prime },s^{\prime })\in E(H)\\ (r^{\prime },s^{\prime })\ne (r,s)\end{array}}f(x_{r^{\prime }},x_{s^{\prime }})\prod _{\begin{array}{c}v\in V(H)\\ v\ne r,s\end{array}}d{x}_v.$$

$${\Delta }_{H}f(x,y):=\sum _{(r,s)\in E(H)}{\Delta }_{H,r,s}f(x,y),$$

$${\Delta }_{H_i}f=a_i{\chi }_A+b_i{\chi }_B,$$

$$t(H_i,h)=t(H_i,f)+{\int }_{[0,1{]}^2}g(x,y){\Delta }_{H_i}f(x,y)dxdy+\text{remainder terms}.$$

$${\int }_{[0,1{]}^2}g(x,y){\Delta }_{H_i}f(x,y)dxdy$$

$$a=\frac{1}{|A|}{\int }_{A}h(x,y)dxdy,b=\frac{1}{|B|}{\int }_{B}h(x,y)dxdy.$$

$${\int }_{[0,1{]}^2}g(x,y){\Delta }_{H_i}f(x,y)dxdy=a_i(a-1)|A|+b_{i}b|B|.$$

$${\int }_{[0,1{]}^4}|g(x,y)||g(u,v)|dxdydudv={\left({\int }_{[0,1{]}^2}|g(x,y)|dxdy\right)}^2=||g|{|}_1^2,$$

$$|g(x,y)|=\left\{\begin{array}[]{ll}1,&\hbox{if $x>c$ and $y>c$,}\\ 0,&\hbox{otherwise,}\\ \end{array}\right.$$

$$||h-f|{|}_1={\int }_A(1-h(x,y))dxdy+{\int }_{B}h(x,y)dxdy=(1-a)|A|+b|B|=||X_h-f|{|}_1,$$

$$t(H_i,X_h)=t(H_i,f)+{\int }_{[0,1{]}^2}{g}^{\prime }(x,y){\Delta }_{H_i}f(x,y)dxdy+\text{remainder terms}.$$

$${\int }_{[0,1{]}^2}{g}^{\prime }(x,y){\Delta }_{H_i}f(x,y)dxdy$$

$${\int }_{[0,1{]}^3}|g^{\prime }(x,y)||g^{\prime }(x,z)|dxdydz\le {\left(\max \{1-a,b\}\right)}^2,$$

$${\psi }^{\beta }(h)=r{T}^u(f)+r\left(\sum _{i=1}^d{a}_i{u}_i\right)(a-1)|A|+r\left(\sum _{i=1}^d{b}_i{u}_i\right)b|B|-I(h)+\text{remainder terms}.$$

$${\psi }^{\beta }(X_h)=r{T}^u(f)+r\left(\sum _{i=1}^d{a}_i{u}_i\right)(a-1)|A|+r\left(\sum _{i=1}^d{b}_i{u}_i\right)b|B|-I(X_h)+\text{remainder terms}.$$

$$I(X_h)=I(a{\chi }_A+b{\chi }_B)=I(a)|A|+I(b)|B|.$$

Full text

exponential random graphs; perturbation analysis; phase transitions; critical directions

I Introduction

Over the last decades, the availability of network data on typically very large scales has created the impetus for the development of new theories and methods for modeling and describing the properties of large networks. The introduction of exponential random graphs has aided in this pursuit, as they are able to capture a wide variety of common network tendencies by representing a complex global structure through a set of tractable local features. From the point of view of extremal combinatorics and statistical mechanics, investigations have been focused on the variational principle of the limiting normalization constant, concentration of the limiting probability distribution, phase transitions, and asymptotic structures. See for example Chatterjee and Varadhan CV2011 , Chatterjee and Diaconis CD2013 , Radin and Yin RY , Lubetzky and Zhao LZ1 ; LZ2 , Radin and Sadun RS1 ; RS2 , Radin et al. RRS , Kenyon et al. KRRS1 , Yin Yin1 , Kenyon and Yin KY , Aristoff and Zhu AZ2 , and Chatterjee and Dembo CD2 . The main techniques used in these papers are variants of statistical physics, but the elegant theory of graph limits as developed by Lovász and coauthors (V.T. Sós, B. Szegedy, C. Borgs, J. Chayes, K. Vesztergombi, …) BCLSV1 ; BCLSV2 ; BCLSV3 ; Lov ; LS , also plays an important role in the interdisciplinary inquiry. Building on earlier work of Aldous Aldous1 and Hoover Hoover , the graph limit theory connects sequences of graphs to a unified graphon space equipped with a cut metric. Though the theory itself is tailored to dense graphs, parallel theories for sparse graphs are likewise emerging. See Benjamini and Schramm BS , Aldous and Steele AS , Aldous and Lyons AL , and Lyons Lyons where the notion of local weak convergence is discussed and the recent works of Borgs et al. BCCZ1 ; BCCZ2 that are making progress towards enriching the existing $L^{\infty}$ theory of dense graph limits by developing a limiting object for sparse graph sequences based on $L^{p}$ graphons.

In this paper, we study the “standard” family of exponential graph models in the asymptotic regime as the exponential parameters approach infinity. As the model name indicates, we are associating exponential weights to graphical ensembles. For each $n$ let ${\mathcal{G}}_{n}$ be the ensemble of simple graphs on $n$ vertices and let ${\mathcal{G}}=\cup_{n}\mathcal{G}_{n}$ be the collection on all simple graphs. The exponential weights are defined in terms of subgraph densities. For any $H\in{\mathcal{G}}$ the homomorphism density of $H$ in a graph $G\in{\mathcal{G}}$ is defined as the probability that a random map on the vertex set of $H$ into the vertex set of $G$, $V(H)\to V(G)$ is edge preserving. We write the homomorphism density as

[TABLE]

Now, for all $n$ we define a probability distribution on ${\mathcal{G}}_{n}$ in terms of the homomorphism density. Let $H_{1},..,H_{d}$ be a given selection of simple graphs, where $H_{1}$ is an edge: $H_{1}=K_{2}$. Let $T:{\mathcal{G}}_{n}\to[0,1]^{d}$, where the components of $T$ are homomorphism densities $\{T(G_{n})\}_{i}=t(H_{i},G_{n})$. Given $\beta\in{\mathbb{R}}^{d}$, define the functional on $G_{n}\in\mathcal{G}_{n}$:

[TABLE]

and weight $G_{n}\in{\mathcal{G}}_{n}$ by $e^{n^{2}T^{\beta}(G_{n})}$. The normalization constant for the ensemble ${\mathcal{G}}_{n}$ is then given by the partition function,

[TABLE]

The terminology for the partition function is borrowed from thermodynamics. In this context, $|\beta|$ is the inverse temperature. Renormalizing $T^{\beta}$ in (2) and fixing $\beta/|\beta|$, one obtains a Hamiltonian $H(G_{n})=-\frac{n^{2}}{|\beta|}T^{\beta}(G_{n})$. As $n\to\infty$ the major contribution to the partition function concentrates around the thermal states (labeled by $f^{\beta}$), in a manner analogous to standard thermodynamic models (see Section I.2). These thermal states are well understood in the large temperature regime $|\beta|\rightarrow 0$. For any selection of subgraphs $H_{1},..,H_{d}$ and $|\beta|$ sufficiently small, the associated thermal state lies in the replica symmetric phase CD2013 , i.e., $f^{\beta}\equiv\rho$. On the other hand, the ground states, which are defined as the limit of the thermal states $f^{\beta}$ as $|\beta|\to\infty$, are not so simple. In some cases the ground state is known to be in the replica symmetric phase CD2013 , while in other cases the ground state concentrates around a simple graph in ${\mathcal{G}}$ YRF .

Our motivation for this paper comes from the edge-triangle model, obtained by setting $d=2$ and $H_{2}$ a triangle: $H_{2}=K_{3}$. It was shown in YRF that there are countably many critical directions of $\beta$, along which the ground state of the model is chosen from finitely many simple graphs with some unknown distribution, and our goal is to develop a mechanism that determines which of these simple graphs is the proper ground state.

I.1 Graphon topology

The thermal states $f^{\beta}$ belong to the space of graph functions “graphons” which may be understood as generalizations of graphs. The set of graphs $\mathcal{G}$ may be embedded into the space of graphons which consist of symmetric measurable functions from $[0,1]^{2}$ into $[0,1]$,

[TABLE]

For any $n$ and a graph $G_{n}\in\mathcal{G}_{n}$, the graphon representation is the function

[TABLE]

where the interval $[0,1]$ may be intuitively thought of as a ‘continuum’ of graph vertices. The distance between graphons is given in terms of the “cut distance”, defined for $f,h\in{\mathcal{W}}$, as

[TABLE]

However, a nontrivial difficulty arrises from the arbitrary labeling of vertices as they are embedded in $[0,1]$. Thus we introduce the equivalence $f\sim f_{\sigma}$ where $f_{\sigma}(x,y)=f(\sigma x,\sigma y)$, where $\sigma:[0,1]\to[0,1]$ is any measure preserving bijection. We write the quotient space of graphons under the equivalence $\sim$ as $\widetilde{\mathcal{W}}$, and the equivalance class under $\sim$ of $f\in{\mathcal{W}}$ as $\widetilde{f}$. Incorporating the equivalence relation $\sim$ yields a distance

[TABLE]

where the infimum ranges over all measure preserving bijections $\sigma_{1}$ and $\sigma_{2}$, making $(\widetilde{\mathcal{W}},\delta_{\Box})$ a compact metric space (see Section 9.3 of Lovász Lov ). With some abuse of notation we also refer to $\delta_{\Box}$ as the “cut distance”.

All graphons arise as the limit of some sequence of graphs. Given a graphon $f\in{\mathcal{W}}$ one may construct a graph $G_{n}=G(n,f)$ by selecting iid points $x_{1},..,x_{n}$ uniformly from $[0,1]$ which represent the vertices of $G_{n}$, and then connect vertices $i,j$ with probability $f(x_{i},x_{j})$. In this context the expected subgraph density is given by the graphon homomorphism density

[TABLE]

which generalizes (1) and is continuous in the metric $\delta_{\Box}$. Indeed, for any graph $H$, the subgraph homomorphism density of random graphs selected by the above construction converges almost surely to the graphon homomorphism density,

[TABLE]

I.2 Exponential random graphs

As discussed above, we will define measures on ${\mathcal{G}}_{n}$ in terms of subgraph densities. For $G_{n}\in{\mathcal{G}}_{n}$ define the probability

[TABLE]

where we have introduced the normalization constant (free energy density),

[TABLE]

By replacing the subgraph homomorphism density (1) with graphon homomorphism density (8) in (2), $T^{\beta}$ extends naturally to $(\widetilde{\mathcal{W}},\delta_{\Box})$. Since $T^{\beta}$ is continuous and bounded on the compact set $\widetilde{\mathcal{W}}$, there is a nonempty compact subset ${\mathcal{K}}^{\beta}$ of $\widetilde{\mathcal{W}}$ so that $T^{\beta}$ is maximized on ${\mathcal{K}}^{\beta}$. Take $I:[0,1]\to{\mathbb{R}}$ as

[TABLE]

and then extend the domain of $I$ to $\widetilde{\mathcal{W}}$ by

[TABLE]

where $f$ is any representative element of $\tilde{f}$. It follows from Lemma 2.1 in Chaterjee and Varadhan CV2011 that $\psi^{\beta}(\widetilde{f}):=T^{\beta}(\widetilde{f})-I(\widetilde{f})$ is well defined on $\widetilde{\mathcal{W}}$ and upper semi-continuous under the cut metric $\delta_{\Box}$. Let ${\mathcal{F}}^{\beta}$ be the subset of $\widetilde{\mathcal{W}}$ where $\psi^{\beta}$ is maximized. Then like ${\mathcal{K}}^{\beta}$, ${\mathcal{F}}^{\beta}$ is a nonempty compact subset of $\widetilde{\mathcal{W}}$.

For the purpose of this paper, two theorems from Chatterjee and Diaconis CD2013 (both based on the large deviation result in CV2011 ) merit some special attention. Together they connect the occurrence of a phase transition in the exponential random graph model with the solution of a certain maximization problem. The first theorem (Theorem 3.1 in CD2013 ) states that for any $\beta$, the limiting normalization constant $\lim_{n\to\infty}\psi_{n}^{\beta}$ of the exponential random graph always exists and is equal to $\psi^{\beta}$. The second theorem (Theorem 3.2 in CD2013 ) states that in the large $n$ limit, the quotient image $\tilde{f}^{G_{n}}$ of a random graph $G_{n}$ drawn from (10) must lie close to ${\mathcal{F}}^{\beta}$ with high probability,

[TABLE]

Since the limiting normalization constant $\psi^{\beta}$ is the generating function for the limiting expectations of other random variables on the graph space such as expectations and correlations of homomorphism densities, a phase transition occurs when $\psi^{\beta}$ is non-analytic or when ${\mathcal{F}}^{\beta}$ is not a singleton set. Although it is difficult to evaluate $\psi^{\beta}$ and determine the maximizing set ${\mathcal{F}}^{\beta}$ for most $\beta$, we will derive an efficient method to approximate ${\mathcal{F}}^{\beta}$ and estimate $\psi^{\beta}$ for $\beta$ sufficiently far from the origin.

II The approximation scheme

Take a finite simple graph $H$ with vertex set $V(H)$ and edge set $E(H)$. Consider a graphon $f\in{\mathcal{W}}$. For each $(r,s)\in E(H)$ and each pair of points $x_{r},x_{s}\in[0,1]$, define

[TABLE]

The above definition appears rather complicated, but in essence may be identified with the homomorphism density (8), if we remove edge $(r,s)$ from $H$ and do not integrate out the associated vertices $x_{r}$ and $x_{s}$. For $(x,y)\in[0,1]$, define

[TABLE]

which corresponds to the total homomorphism density generated after all possible ways of removing one edge from $H$. We give some examples to illustrate this idea. When $H$ is an edge, $\Delta_{H}f(x,y)\equiv 1$. When $H$ is a triangle, by symmetry, $\Delta_{H}f(x,y)=3\int_{[0,1]}f(x,z)f(y,z)dz$.

II.1 Motivation

We will consider as a motivating example the edge-triangle model, which is a $2$-parameter exponential random graph model obtained from (10) by setting $H_{1}$ to an edge and $H_{2}$ to a triangle. Take $\beta=ru$ for $u=(u_{1},u_{2})\in{\mathbb{S}}$. Suppose that $u_{1}>0$ and $u_{2}<0$, i.e., the $2$-dimensional parameter vector $u$ is pointing towards the $4$th quadrant. Then ${\mathcal{K}}^{u}$ (which agrees with ${\mathcal{K}}^{\beta}$) consists of graphons $f$ that minimize $(u_{1}/u_{2})e+t$, where $e=t(H_{1},f)$ denotes the edge density and $t=t(H_{2},f)$ denotes the triangle density of $f$, respectively. This implies that $f$, the maximizers of $T^{u}$ (and hence of $T^{\beta}$) must lie on the Razborov curve, which is the lower boundary of the feasible region $R$ of edge-triangle homomorphism densities. As $u_{1}e+u_{2}t$ is a linear function, $f$ must minimize over the convex hull $P$ of $R$. Since $R$ and $P$ only intersect at the points corresponding to Turán graphons, $f$ must be a Turán graphon. This important fact about the structure of ${\mathcal{K}}^{\beta}$ was further used to derive the maximizing graphons of ${\mathcal{F}}^{\beta}$ in YRF . Using the boundedness of $I$, we found that ${\mathcal{F}}^{\beta}$ consists of graphons that can be made arbitrarily close to Turán graphons when the magnitude of $\beta$ is sufficiently large, and exactly which Turán graphon is favored by ${\mathcal{F}}^{\beta}$ depends on the direction of $\beta$. However, as nice as these results are, there is ambiguity concerning the optimal Turán graphon along the critical directions of $\beta$, which correspond to normal lines of the convex hull $P$. The subtlety might be due to the fact that Turán graphons, though close to our optimal graphons in cut distance, are not best approximations to the maximizing set ${\mathcal{F}}^{\beta}$. Generalizing from Turán graphons HN , we say that a graphon $f$ is random-free if $f=\chi_{A}$ for some symmetric measurable subset $A$ of $[0,1]^{2}$. When $\beta$ is finite, Chatterjee and Diaconis CD2013 showed that the maximizing graphons in ${\mathcal{F}}^{\beta}$ are almost nowhere random-free – that is the set $\{x\in[0,1]^{2}:f(x)\in\{0,1\}\}$ has zero measure. We are thus interested in finding coarse-grained graphons that are not random-free but close to ${\mathcal{F}}^{\beta}$ in cut distance, as they are sufficient to distinguish between candidates for the ground state. Resorting to perturbation analysis in the $|\beta|\to\infty$ regime, we will propose a method that keeps track of only the most significant characteristics of the maximizing graphons and demonstrate its effectiveness.

II.2 Assumptions

Let $R$ denote the range of $T$ in ${\mathbb{R}}^{d}$ and suppose that the boundary $\partial R$ is piecewise analytic. For $w\in\partial R$, let $F_{R}(w)$ denote the set of all feasible directions of $R$ at $w$, so that $v\in F_{R}(w)$ implies that there is $\epsilon>0$ with $w+\epsilon v\in R$. The (internal) tangent cone $C_{R}(w)$ of $R$ at $w$ is then given by the closure of $F_{R}(w)$. See Figure 1 for an illustration of this concept in the edge-triangle model. Let $\beta=ru$ for $u\in{\mathbb{S}}^{d-1}$, where $r$ is sufficiently large. Suppose that $T^{u}$ (and hence $T^{\beta}$ as $|\beta|\to\infty$) is maximized at a set of random-free graphons $f=\chi_{A}$, i.e., ${\mathcal{K}}^{u}$ (which agrees with ${\mathcal{K}}^{\beta}$) consists of random-free graphons only. Further suppose that $f\in{\mathcal{K}}^{u}$ has the property that

[TABLE]

where $\Delta_{H_{i}}f$ is defined as in (16) and $B:=[0,1]^{2}\texttt{\char 92\relax}A$. We remark that these nice properties that we assumed are enjoyed beyond the edge-triangle model. Denote by $K_{k}$ a complete graph on $k$ vertices. Since the vertices of the convex hull $P$ of $R$ for $K_{2}$ and $K_{n}$ (any $n$), not just $K_{2}$ and $K_{3}$ as in the edge-triangle case, are given by Turán graphons, our argument will run through without much modification in these cases. Utilizing geometry of their respective convex hull EN , similar analysis may be extended to more general models. We point out in particular that if $f$ is a Turán graphon, then (17) holds for any $H_{i}$ due to symmetry.

II.3 Perturbation of random-free graphons

Consider a maximizing graphon $h\in{\mathcal{F}}^{\beta}$. When the magnitude of $\beta$ is sufficiently large, since $T^{\beta}$ diverges while $I$ stays bounded, ${\mathcal{F}}^{\beta}$ and ${\mathcal{K}}^{\beta}$ can be made arbitrarily close in the cut metric. From our assumptions in Section II.2, there thus exists a random-free graphon $f=\chi_{A}$ that is close to $h$ under the cut distance. We will further construct a non-random-free graphon $X_{h}$ that is close to $h$, simple enough, and yet still retains important information of $h$. Most importantly, we will show that for any finite simple graph $H_{i}$, $t(H_{i},X_{h})$ approximates $t(H_{i},h)$ at least as well as $t(H_{i},f)$ asymptotically. The following proposition is specific for random-free graphons and will be useful for our investigation.

Proposition II.1**.**

Let $||\cdot||_{1}$ denote the $L^{1}$-norm. For a random-free graphon $f$, $d_{\Box}(f,h)=O(||f-h||_{1})$.

Proof.

By (6), it is clear that $d_{\Box}(f,h)\leq||f-h||_{1}$. For the other direction, suppose $f=\chi_{A}$ for some symmetric measurable subset $A\subset[0,1]^{2}$. If $A$ is a rectangle, then $||f-h||_{1}\leq 3d_{\Box}(f,h)$. The general conclusion follows once we recognize that any subset $A$ of $[0,1]^{2}$ may be approximated within $\epsilon>0$ by a finite union of disjoint rectangles. ∎

Let us first expand $t(H_{i},h)$ around $f$. Denote a perturbation of $f$ by $g=h-f$. Corresponding to the regime $|\beta|\to\infty$, as explained earlier, $f$ may be chosen so that $|g|\in{\mathcal{W}}$. We have

[TABLE]

Under our assumption (17), we compute the first variation:

[TABLE]

Define

[TABLE]

Using that $f=\chi_{A}$ and $g=h-f$, (II.3) reduces to

[TABLE]

The remainder terms (if they exist) give the homomorphism density after all possible ways of removing at least two edges from $H_{i}$, and so are bounded above by either $\int_{[0,1]^{3}}|g(x,y)||g(x,z)|dxdydz$ or $\int_{[0,1]^{4}}|g(x,y)||g(u,v)|dxdydudv$. Estimating the latter is easy.

[TABLE]

which is of negligible asymptotic order when compared to the first variation (II.3). The former corresponds to the $2$-star density of a graphon $|g|$ with edge density $||g||_{1}$. By AK , for small $||g||_{1}$, the $2$-star density is bounded below by $||g||_{1}^{2}$ and above by $||g||_{1}/2+O(||g||_{1}^{2})$, and the upper bound is achieved when $|g|$ is an anticlique of the form

[TABLE]

where $c=1-\sqrt{1-||g||_{1}}$. The lower bound for this remainder term is of desirable asymptotic order, but the upper bound is of the same order as the first variation. Recall that the maximizing graphon $h$ for a finite $\beta$ is uniformly bounded away from [math] and $1$ CD2013 . The graphon $|g|=|h-f|$ is thus likely quite different from an anticlique. This implies that the $2$-star density of $|h-f|$ is of higher order than $||h-f||_{1}$ and does not achieve the upper bound. The phenomenon was confirmed for example by simulations for the edge-triangle model KRRS2 .

Define the averaged perturbation by $X_{h}=a\chi_{A}+b\chi_{B}$, where $a$ and $b$ are given in (20). $X_{h}$ may be viewed as a flattened out version of $h$. Since

[TABLE]

$X_{h}$ is close to $f$ (and hence $h$) under the $L^{1}$ distance, and by Proposition II.1, also under the cut distance. Denote by $g^{\prime}=X_{h}-f$. We perform the same expansion for $t(H_{i},X_{h})$ around $f$ as in the last paragraph:

[TABLE]

Following similar reasoning as in (II.3) and using the definition of $X_{h}$,

[TABLE]

This says that $t(H_{i},h)$ and $t(H_{i},X_{h})$ agree except for the remainder terms. As for $t(H_{i},h)$, the remainder terms for $t(H_{i},X_{h})$ (if they exist) are bounded above by either $\int_{[0,1]^{3}}|g^{\prime}(x,y)||g^{\prime}(x,z)|dxdydz$ or $\int_{[0,1]^{4}}|g^{\prime}(x,y)||g^{\prime}(u,v)|dxdydudv$. Since $||g^{\prime}||_{1}=||g||_{1}=(1-a)|A|+b|B|$ by (24), the latter is of asymptotic order $||g||_{1}^{2}$; while the former gives

[TABLE]

and so is also of asymptotic order $||g||_{1}^{2}$. We conclude that $t(H_{i},X_{h})$ gives at least as good an asymptotic approximation for $t(H_{i},h)$ as $t(H_{i},f)$, and a better one when the error term related to the $2$-star density of $|h-f|$ may be dropped.

II.4 Maximizing graphons

Let $\beta=ru$ for $u\in{\mathbb{R}}^{d}$ and $\|u\|=1$. Under our assumptions, ${\mathcal{K}}^{u}$ (which agrees with ${\mathcal{K}}^{\beta}$) is a set of random-free graphons. Take $h\in{\mathcal{F}}^{\beta}$. Then for $\epsilon>0$ and $r$ large (corresponding to $|\beta|\to\infty$), $d_{\Box}(f,h)<\epsilon$ for some $f\in{\mathcal{K}}^{u}$, i.e., $f=\chi_{A}$ for some symmetric measurable subset $A$ of $[0,1]^{2}$. By (18) and following analysis,

[TABLE]

Similarly, by (25) and following analysis,

[TABLE]

In both equations above, the entropy $I$ is bounded in contrast with the energy contribution $T^{\beta}$. Except that $h$ is close to $f$ (and hence $X_{h}$) in cut distance, we do not have enough information regarding the structure of $h$. Rather than maximizing $\psi^{\beta}(h)$ over all possible graphons $h$ directly, we will maximize $\psi^{\beta}(X_{h})$ over $2$-parameter families $0\leq a,b\leq 1$. From the heuristics in Section II.3, this is an effective method to approximate the optimal graphon. More than that, in certain situations (for example the edge-triangle model to be discussed in detail in Section III), we will see that keeping $2$ parameters is not only effective but also sufficient. Notice that

[TABLE]

We rewrite $X_{h}$:

[TABLE]

Maximizing each first variation term, we have

[TABLE]

Under this choice of $a$ and $b$ and provided we can ignore the remainder terms, $\psi^{\beta}(X_{h})$ is strictly bigger than $rT^{u}(f)$, which is the random-free graphon approximation corresponding to $a=1$ and $b=0$. Let us verify that $\psi^{\beta}(X_{h})$ is indeed strictly bigger than $rT^{u}(f)$ by rigorously managing the error term in (31). As shown earlier in Section II.3, the remainder terms are of higher order:

[TABLE]

For any $\epsilon>0$, there exists large enough $r$ so that $||X_{h}-f||_{1}$ is sufficiently small, making $0<1-a<\epsilon$ and $0<b<\epsilon$. Applying these $\epsilon$ bounds in (33) gives

[TABLE]

Maximizing each first variation term as previously, this yields $\psi^{\beta}(X_{h})>rT^{u}(f)$. We conclude that, as expected, the addition of one more parameter improves the random-free graphon estimation.

III The edge-triangle model

Denote by $K_{k}$ a complete graph on $k$ vertices. The edge-triangle model is a $2$-parameter exponential random graph model obtained by taking $H_{1}$ an edge ($K_{2}$) and $H_{2}$ a triangle ($K_{3}$) in (10). Consider the set $R=\{(t(K_{2},f),t(K_{3},f)),f\in{\mathcal{W}}\}$ of all realizable values of the edge ($e$) and triangle ($t$) homomorphism densities as the graphon $f$ varies over the entire graphon space ${\mathcal{W}}$. See Figure 1. The upper boundary curve of $R$ is given by the equation $t=e^{3/2}$, and can be derived using the Kruskal-Katona theorem (see Section 16.3 of Lov ). The lower boundary curve is trickier. The trivial lower bound of $t=g_{1}(e):=0$, corresponding to the horizontal segment, is attainable at any $0\leq e\leq 1/2$ by graphons describing the possibly asymptotic edge density of subgraphs of complete bipartite graphs (Turán graphon with $2$ classes). For $e\geq 1/2$, the optimal bound was obtained by Razborov Razborov , who established, using the flag algebra calculus, that for $(k-1)/k\leq e\leq k/(k+1)$ with $k\geq 2$,

[TABLE]

All the curve segments $g_{k}(e)$ describing the nontrivial part of the lower boundary of $R$ are strictly concave. For $k=1,2,\ldots$, we set $v_{k}:=(e_{k},t_{k})=(t(K_{2},f^{K_{k+1}}),t(K_{3},f^{K_{k+1}}))$, where explicitly,

[TABLE]

is the Turán graphon with $k+1$ classes. Thus

[TABLE]

For $k=1,2,\ldots$, let $L_{k}$ be the line segment joining vertices $v_{k}$ and $v_{k+1}$ of neighboring Turán graphons. These infinitely many line segments form the convex hull $P$ of $R$, and the length of $L_{k}$ decreases monotonically to zero as $k$ gets large.

The normal vectors to $L_{k}$,

[TABLE]

are the critical directions of the edge-triangle model. Let $\beta=ro_{k}$ and take $r\rightarrow\infty$. While the vectors $o_{k}$ (38) are not normalized as in our derivation (see Section II), this can be easily adjusted by adapting $r$. Then Turán graphons with $k+1$ and $k+2$ classes both belong to ${\mathcal{K}}^{o_{k}}$ (which agrees with ${\mathcal{K}}^{\beta}$), and we concluded in YRF that a typical graph sampled from the model may behave like either a Turán graphon with $k+1$ classes or a Turán graphon with $k+2$ classes, with no clear preference. Though already quite informative, as explained earlier in Section II.1, this result remains somewhat unsatisfactory because it does not indicate whether both such graphons are actually realizable in the limit and in what manner. Notice that underneath our investigation, there is an ordered double asymptotic framework, in the sense that the network size $n$ goes to infinity first followed by the divergence of the parameters $\beta$. In hope of resolving this rather subtle ambiguity within the edge-triangle model, the “other” order was also examined in YRF , where we first let the magnitude of $\beta$ increase to infinity so as to isolate a simpler sub-model and then study its limiting properties as $n$ grows. Both ordered asymptotics imply a nearly identical convergence in probability in the cut metric along the noncritical directions. Under the “original” as well as the “reversed” order, there exist (possibly different) subsequences of the form $\{n_{i},\beta_{1,i},\beta_{2,i}\}$, with $n_{i}\rightarrow\infty$, $\beta_{1,i}\rightarrow\infty$ and $\beta_{2,i}\rightarrow-\infty$ for $i=1,2,\ldots$, where the edge-triangle model converges to some Turán graphon specified by the direction of the parameters $(\beta_{1},\beta_{2})$. Additionally, under the “reversed” asymptotics, a detailed categorization of the limiting behavior of the edge-triangle model was obtained: When $\beta$ diverges along the critical direction $o_{k}$, a typical sampled graph more likely resembles a Turán graphon with $k+2$ classes than with $k+1$ classes. Now that we are equipped with refined perturbation analysis, we would like to sharpen our results under the “original” asymptotics and inquire whether the same type of discontinuity in natural parametrization exists.

Let us make a further remark before carrying out the detailed calculations. In the physics literature, people are often interested in cases where the parameter $\beta$ depends on $n$ (some averages need to be satisfied for every $n$). In these models, in place of the normalization constant (free energy density), the relative entropy plays a central role. Analogous (but more complicated) maximization problems and concentration of measure results have been established, which lead to classifications of ensemble equivalence between the microcanonical ensemble and the canonical ensemble. The perturbative methods explored in the current paper are expected to apply in these general parameter situations. In some cases, the perturbation would still be around random-free graphons, and our argument will run through without much adaptation HMRS . In some other cases however, the perturbation would be around graphons admitting more intricate structures, and serious future work is neededGHR ; PN .

III.1 Perturbation analysis

Let $\beta=ro_{k}$. We will compare $\psi^{\beta}(X)$ and $\psi^{\beta}(Y)$, where $X$ is the flattened out graphon close to $f^{K_{k+1}}$ (Turán graphon with $k+1$ classes) and $Y$ is the flattened out graphon close to $f^{K_{k+2}}$ (Turán graphon with $k+2$ classes). Both $X$ and $Y$ are constructed with the optimal perturbation values (32). We set $u=o_{k}$ in our calculations below. For $i=1,2$, we compute $a_{i}^{X}$, $b_{i}^{X}$, $a_{i}^{Y}$, and $b_{i}^{Y}$ for $f^{K_{k+1}}$ and $f^{K_{k+2}}$ in (17), where $H_{1}=K_{2}$ and $H_{2}=K_{3}$. As pointed out earlier in Section II, $a_{1}^{X}=b_{1}^{X}=a_{1}^{Y}=b_{1}^{Y}=1$. For notational convenience, from now on we denote $f^{K_{k+1}}$ by $T_{k}$ and $\chi_{[0,1]^{2}}-f^{K_{k+1}}$ by $D_{k}$. Then $T_{k}$ and $D_{k}$ are indicator functions associated with sets of measure $k/(k+1)$ and $1/(k+1)$, respectively. Reconfirming our assumption in (17),

[TABLE]

which gives

[TABLE]

[TABLE]

This yields

[TABLE]

[TABLE]

We check the remainder terms for $\psi^{\beta}(X)$ after first order perturbation. Similar analysis will work for $\psi^{\beta}(Y)$ and we skip the details. Consider the (internal) tangent cone $C_{R}(v_{k})$, consisting of graphons of the type $X_{h}=aT_{k}+bD_{k}$. See Figure 2. Then

[TABLE]

For $(x,y)\in[0,1]^{2}$, we have

[TABLE]

which implies that

[TABLE]

[TABLE]

This yields

[TABLE]

Using (32) and (41) for $a$ and $b$, (45) gives $\psi^{\beta}(X)$. After first order perturbation as in (31), the remainder terms are bounded by

[TABLE]

The following lemma is useful for our asymptotic derivation.

Lemma III.1**.**

Let $I$ be defined as in (12). As $x\rightarrow\infty$,

[TABLE]

and as $x\rightarrow-\infty$,

[TABLE]

Proof.

We recognize that after simplification, the left hand side of (47) becomes $\log(1+e^{-x})/2$, and the left hand side of (48) becomes $x/2+\log(1+e^{-x})/2$. The rest is immediate. ∎

Applying Lemma III.1 to the first order perturbation terms, we see that conforming to our heuristic analysis in Sections II.3 and II.4, the remainder terms are indeed of negligible order. For $r$ large enough, we have

[TABLE]

[TABLE]

Since $\sum_{i=1}^{2}b_{i}^{X}u_{i}=-1/(3k+5)$ has the smallest absolute value among all equations in (41), $\psi^{\beta}(X)>\psi^{\beta}(Y)$ both in terms of first order perturbation (49) and the exact value.

Denote by $\psi^{\beta}_{\text{opt}}(X_{h})$ and $\psi^{\beta}_{\text{opt}}(Y_{h})$ the exact optimal value of $\psi^{\beta}$ within the (internal) tangent cone $C_{R}(v_{k})$ and $C_{R}(v_{k+1})$, respectively. Let us compare $\psi^{\beta}_{\text{opt}}(X_{h})$ with $\psi^{\beta}_{\text{opt}}(Y_{h})$. From our heuristic argument, for any $\epsilon>0$ and $r$ large enough, (34) is satisfied:

[TABLE]

Following similar reasoning, we also find a bound in the other direction:

[TABLE]

Maximizing each first variation term and applying Lemma III.1 as previously, this says that $\psi^{\beta}(X_{h})$ is asymptotically bounded below by

[TABLE]

and above by

[TABLE]

Similar analysis works for $\psi^{\beta}(Y_{h})$ and we skip the details. Since $\epsilon>0$ can be taken arbitrarily small and $\sum_{i=1}^{2}b_{i}^{X}u_{i}=-1/(3k+5)$ has the smallest absolute value among all equations in (41), $\psi^{\beta}_{\text{opt}}(X_{h})>\psi^{\beta}_{\text{opt}}(Y_{h})$. As discussed earlier in Section I, by Theorems 3.1 and 3.2 in Chatterjee and Diaconis CD2013 , we conclude that when $\beta$ diverges along the critical direction $o_{k}$, a typical sampled graph more likely resembles a Turán graphon with $k+1$ classes than with $k+2$ classes. See Table 1 and Figure 3.

Theorem III.2**.**

Consider the edge-triangle exponential random graph model, obtained by setting in (10) $H_{1}$ an edge and $H_{2}$ a triangle. For $k\geq 1$, let $\beta=ro_{k}$, where $o_{k}$ is the critical direction (38) and $r$ is sufficiently large. Then in the large $n$ limit, a typical graph drawn from the model behaves like a Turán graphon with $k+1$ classes,

[TABLE]

III.2 Geometric interpretation

We proceed further and examine the effect of infinitesimal perturbation on the associated edge and triangle densities. From (42) and (44),

[TABLE]

[TABLE]

Let us present another perspective on this calculation incorporating assumption (17), which may be employed to derive infinitesimal variations for more complicated homomorphism densities. The idea appeared in our heuristic analysis before and we make it explicit here. For any $H_{i}$ so that (17) is satisfied,

[TABLE]

[TABLE]

Using $a_{1}=b_{1}=1$, $a_{2}=3(k-1)/(k+1)$ and $b_{2}=3k/(k+1)$ (40), we recover the partial derivatives calculated above. In particular, we recognize that $-(a_{1},a_{2})$ points along the left tangent line and $(b_{1},b_{2})$ points along the right tangent line, while the critical direction $u=o_{k}$ is the normal vector to the line segment $L_{k}$ that connects neighboring vertices $v_{k}$ and $v_{k+1}$. This offers a geometric justification of the ($\pm$) signs in (41).

Recall that the lower boundary of attainable edge-triangle densities is a piecewise algebraic curve with infinitely many concave pieces $g_{k}$ (35), and the connection point of $g_{k}$ and $g_{k+1}$ is $v_{k}=(e_{k},t_{k})$. See Figure 1. We compute

[TABLE]

which implies that the left and right hand side derivatives at $e_{k}=k/(k+1)$ are respectively given by

[TABLE]

The partial derivative vectors in (55) thus delineate the boundary of the (internal) tangent cone $C_{R}(v_{k})$. In other words, $C_{R}(v_{k})$ spans all possible infinitesimal variations at the Turán graphon with $k+1$ classes. See Figure 2. Since the optimizing graphon associated with large enough $\beta$ must lie within the tangent cone of some Turán graphon, it may be represented by a linear combination of Erdős-Rényi and Turán graphons. Even though the graphon representation may not be unique, optimizing over all possible combinations provides insight into the structure of the maximizing set. Keeping track of the Erdős-Rényi and Turán characteristics in the edge-triangle model is thus not only an effective but also sufficient method to estimate the normalization constant, and gives evidence of discontinuity of the natural parametrization along the critical directions $o_{k}$ in the limit as $n$ and then $r$ tend to infinity. This demonstrates the occurrence of discontinuous phase transitions in the edge-triangle model.

Acknowledgements

The authors are very grateful to the anonymous referee for the invaluable suggestions that greatly improved the quality of this paper. Mei Yin thanks Sukhada Fadnavis for helpful conversations. Rajinder Mavi was supported by a postdoctoral fellowship from the Michigan State University Institute for Theoretical and Mathematical Physics. Mei Yin’s research was partially supported by NSF grant DMS-1308333.