A Deterministic Partition Function Approximation for Exponential Random Graph Models

1 Introduction

Social networks are becoming central components to many aspects of marketing, recruiting, web search, and education programs [Kempe et al., 2003; Cetintas et al., 2011; Mancilla-Caceres et al., 2012]. Careful use of social network analysis in those areas is the key to future advances. For that reason, many researchers and practitioners model their relevant social networks and learn them from data [Carrington et al., 2005]. Many of those social networks are large, and precise modeling of them is difficult. A family of simple models called Exponential Random Graph Models (ERGM) [Robins et al., 2007] is commonly used by researchers and practitioners for this purpose.

An ERGM defines a distribution over all graphs of n nodes. Coefficients and subgraph statistics, such as number of edges, triangles, and k-stars, are then used to specify ERGM distributions [Robins et al., 2007]. The model captures the correlation of network sub-structures and enables various inferences on complex networks. For example, we can tell whether tran-

Jaesik Choi

Ulsan National Institute of Science and Technology, Ulsan, Korea jaesik@unist.ac.kr

Eyal Amir

University of Illinois at Urbana-Champaign Urbana, IL, USA eyal@illinois.edu

sitivity is prominent in a network by fitting an ERGM with related subgraphs features, such as triangles.

Learning ERGMs from data is done by Maximum Likelihood Estimation (MLE). Unfortunately, such learning is hard even for networks of modest sizes (e.g. 40 nodes) because calculating normalizing constants (partition functions) precisely for such models is intractable. For this reason most current techniques involve sampling using Markov Chain Monte Carlo (MCMC) [Handcock et al., 2003; van Duijn et al., 2009]. This often results in intractable computation or highly imprecise results for these modest-sizeor-larger networks [Bhamidi et al., 2008; Snijders, 2002; Handcock, 2003; Hunter et al., 2012].

Recently, Chatterjee and Diaconis [Chatterjee and Diaconis, 2011] derived an analytic approximation to the loglikelihood function of ERGM by applying a new large deviation principle result on Erd¨os-R´enyi models. Chatterjee and Diaconis's theoretical analysis proves a strong convergence result on their approximation, but it only applies to ERGMs with uniformly non-negative/non-positive coefficients for subgraph features.

In this paper, we introduces a new deterministic approximation solution which generalizes Chatterjee and Diaconis's approximation for estimating ERGM partition functions. The two approximations are asymptotically equivalent for large n. However, instead of complex large deviation analysis, we derive the new approximation using much simpler analysis based on weaker convergence results. Moreover, our analysis lifts the unrealistic condition on uniformly non-negative/nonpositive subgraph feature coefficients, enabling application of the approximation to the whole subgraph feature coefficient space. Finally, we provide a constructive and pragmatic approach to the problem, which enables us to evaluate the accuracy of the approximation empirically against common stochastic approximation methods in realistic settings.

Specifically, we present a quadratic time (or linear time wrt the number of random variables) deterministic approximation to the log partition function of ERGMs. Asymptotic properties of the subgraph statistics space enable this new approximation. The approximation works as follows: given (coefficient) parameters θ of an ERGM, find that edge-count u (between 0 and (n2 )) that maximizes ˜ γ(θ, u) = θT ρ(u)+C(n, u) (See (8) for definition), where ρ(u) is a vector of subgraph statistics approximated for graphs with u edges and function C(n, u) approximates the logarithm of the number of graphs with subgraph statistics close to ρ(u). Once the maximizing u is found, we estimate the log partition function lnZ(θ) by ˜ γ(θ, u). The approximation works because this ρ(u) captures the subgraph statistics of a large (asymptotically) mass of graphs of n nodes. So, in a sense, many graphs would look similar from a subgraph statistics perspective. We also reveal that θ = Θ(n2 ) is a necessary condition for θ to be asymptotically relevant in an ERGM for network of order n.

Our results show that the new method performs well experimentally comparing to existing sampling methods [Gelman and Meng, 1998; Handcock et al., 2003] on synthetic data and real-world social networks. Our results also show that the new algorithm yields reliable approximation for many models when the size of the network is larger than 30.

The rest of the paper is organized as follows: Section 2 reviews ERGM, Section 3 describes the components of the approximation and key theoretical results, Section 4 describes our experimental evaluation, Section 5 reviews related work and Section 6 concludes the paper.

2 Background

An ERGM defines the following distribution over order-n graphs g ∈ G:

where φ(g) is the feature vector for graph g ∈ G; the parameter θ is a real vector; partition function Z(θ) is a normalizing constant.

The feature vector φ(g) may include any network and nodal attributes of g , and the edge statistics is almost always included [Robins et al., 2007]. In this work, we focus on undirected graphs and subgraph statistics features. Specifically, for a set of subgraph structures of interest {L1, . . . ,Lr}, the feature vector of undirected graph g can be defined with subgraph densities as below:

Here t(g,Li) counts the number of subgraphs in g that are isomorphic to Li; Kn is the order-n complete graph, therefore t(Kn,Li) = ( n vi )t(Kvi ,Li) is a constant for any Li of order vi. Notice that the simplest subgraph K2, or edge, is almost always included as a feature in ERGMs. Its role in the model is similar to that of the intercept term in most linear regression models [Robins et al., 2007; Hunter, 2007]. For the rest of the paper, we assume K2 is always included in the feature subgraphs.

Example: Figure 1 illustrates a simple example network of order 7. It has seven edges, one triangle, eleven 2-stars, five 3stars and no rectangle. The third column shows the subgraph densities of the network. For example, the 7-node labeled graph can have at most (73) × 1 = 35 triangles, therefore the triangle density is 1/35 ≃ 0.029.

Given a network g, the MLE of parameter vector θ is:

In this paper, we are interested in approximation of the log partition function lnZ(θ).

3 Approximating Log Partition Functions

In this section, we derive a deterministic approximation to the log partition function lnZ(θ) and analyze the behavior of θ for networks of different sizes. We first introduce the counting function for graphs with the same feature vector, which leads to a set of edge-number induced lower bounds of lnZ(θ) and its approximation. Figure 2 gives an overview of the new approach.

3.1 Graph counting in the feature space

We introduce the key concept of graph counting function for the feature space of ERGM. Let φ : G → H be the function that maps a graph in g ∈ G to the subgraph density space H of the graphs. For h ∈ H, we define counting function #(h) = |Gh| where Gh = {g ∈ G|φ(g) = h}|, i.e. the number of graphs in G having h as subgraph densities. We re-write the partition function (1) as a compact form using counting function 1 :

Notice that when θ=0, each term in (4) simply counts the graphs with given subgraph configuration, and the normalizing constant becomes the total number of graphs |G|. Later we will show how the graph counting interpretation helps in computing lnZ(θ).

Let L1,L2, . . . ,Lr be simple graphs of interest and vi be the number of nodes for Li. The following lemma provides an upper bound to |H|. Under the assumption ∀i, n ≫ vi and n ≫ r, the lemma establishes reasonable error bounds for several arguments in the rest of the paper:

1 Note that all isomorphic graphs have the same subgraph densities, but the reverse it not true. Two non-isomorphic graphs may also have the same subgraph densities.

Figure 2: The algorithm has two approximations: a© γ(θ, u) is an edge-count-u-induced lower bound for lnZ(θ), Lemma 3 shows the error is bounded in O(lnn); b©We propose ˜ γ(θ, u) as an approximation to the unknown γ(θ, u), following Lemma 4 and Lemma 5.

Given some set S, and any function f : S → R, a well known computation trick of computing ln ∑x∈S expf(x) is to use maxx∈S f(x) as an approximation if |S| is small. Specifically, we have the following bounds:

Lemma 2. Let f be a function on S and x∗ = argmaxx∈Sf(x), it holds that:

Direct application of Lemma 2 to lnZ(θ) yields a sloppy approximation because the huge size of G. Thanks to Lemma 1, the following approximation to (4) has a much tighter error bound:

Now, we discuss how to estimate the first term of (5).

3.2 Edge-Count Induced Lower Bounds

In this section, we first derive an alternative representation to the approximation in (5), then develop an estimator for the approximation.

Let Gu ⊂ G be the set of graphs with u edges, Hu ⊂ H be the set of subgraph statistics induced by Gu, and #u(h) be the restricted counting function which only counts graphs in Gu, i,e. #u(h) = |{g ∈ Gu|φ(g) = h}|. For any θ and u, we have the following lower bound to (5):

Notice that the equality holds when K2 is a feature subgraph and u = argmaxuγ(θ, u), because in this case Hu ∩Hu′ = ∅ if u′ ̸= u, therefore {Hu} is a partition of H. Specifically, we have the following lemma:

Lemma 3. Given K2 is included in subgraph features, the following equation holds:

Lemma 3 shows maxu{γ(θ, u)} can be treated as a good approximation of lnZ(θ) with bounded error. However, the computation of γ(θ, u) is still non-trivial. For the rest of this section, we develop an approximation of γ(θ, u) by exploiting the asymptotic property of #u(h) in Gu.

Define h′ and h∗ as the optimum of γ(θ, u) and maximizer of #u(h) respectively:

The following bounds of γ(θ, u) hold for all θ and u:

Notice that the gap between the upper and lower bounds in (6) is θT (h′ - h∗ ). When θ = 0, we have h′ = h∗ ; As θ deviates from 0, h′ will also stray from h∗ and settle the conflict between the increasing linear term θT h and a very quickly diminishing ln#u(h) as we will soon discuss in the next section. We argue that h∗ can be used to approximate h′ in terms of estimating the log partition function lnZ(θ) if the model of interest is not dominated by a few graphs 2 . Compared to h′ and ln#u(h′ ), h∗ and ln#u(h∗ ) are much easier to estimate, therefore lead to a feasible approximation to γ(θ, u).

Estimating h∗ :

h∗ maximizes the counting function #u(h) on Gu, therefore for any randomly picked graph g ∈ Gu, h∗ is the most likely value of φ(g). If we define a uniform distribution on Gu, then h∗ is the mode of φ(Gu).

The process of generating graphs randomly from Gu is known as as Erd¨os-R´enyi (ER) random graphs model G(n,M) [Erd¨os and R´enyi, 1960]. Here n is the number of nodes in the graph and M = u is the number of edges. An alternative (and popular) definition of ER model is G(n, p) [Gilbert, 1959], in which an order-n graph is constructed by picking each edge independently with probability p. The distribution of subgraph statistics for ER models has been actively studied in probabilistic graph theory[Nowicki, 1989; D¨oring and Eichelsbacher, 2009].

Nowicki [Nowicki, 1989] proved that φ(g) is asymptotically normally distributed for g ∈ G(n, p). The following lemma extends that result to G(n,M) over Gu using Chebyshev's inequality:

Lemma 4. Let si be the edge count of Li, define function ρi(u) = (u/(n2 ))si . Given any edge density µ, write the edge count u = (n2 )µ as a function of n. Then for any real vector a = (a1, a2, . . . , ar)T and random graph g ∈ G(n,M = u), the following holds as n → ∞:

where ρ(u) = (ρ1(u), . . . , ρr(u))T and c is some constant.

Notice here ρi(u) is the expected density of Li in G(n, p = u/(n2 )). Lemma 4 suggests that the subgraph densities for most graphs of Gu are close to ρ(u)! In a sense, graphs in Gu form a cluster in terms of the subgraph statistics. Figure 3 illustrates the phenomenon using order 12 unlabeled graphs [Brouwer, accessed 2012 Sep 30]. This result suggests that ρ(u) is a good estimation of h∗ for large n.

2 In the cases where a few graphs are dominating the model, h′ will sway away from h∗ . Large entries in θ will lead to this scenario as we will show case in Section 4.3.

Estimating ln#u(h∗ ):

Eq (7) also hints using |Gu| to approximate #u(h∗ ) as φ(g) concentrates. With the help of Lemma 1, it turns out to be a very good estimation:

Lemma 5. Given an edge count u, it holds that ln#u(h∗ )=(n2 )H(u/(n2 ))-O(ln n) where H(x)=-xln x-(1-x)ln(1-x).

Estimating Lower Bound γ(θ, u):

Given estimations of h∗ ≃ ρ(u) and ln#u(h∗ ) ≃ (n2 )H (u/(n2 )), we immediately have the following approximation to γ(θ, u):

Lemma 4 hints that as h′ deviates away from h∗ , ln#u(h′ ) will diminish rapidly as θT h′ increases linearly. When the gradient of the linear term is small, h∗ tends to be a good approximation of h′ . If the gradient is steep, h′ will lean towards the extreme entry in Hu that maximizes the linear term but leads to a minimal #u(h′ ), i.e. only one graph (or a few graphs) in Gu has feature vector h′ , but it dominates all other graphs.

Before the discussion of the approximate algorithm and its behavior as n →∞, we first investigate the behavior of fixed θ for networks of different order n.

3.3 Necessary Condition for Asymptotically Relevant Parameters

The analysis in Section 3.2 also reveals that ERGM with fixed θ may have very different behavior for networks of different order n. To see this, let u∗ = argmaxu γ(θ, u), we show that as n → ∞, any fixed θ becomes asymptotically irrelevant in maxu γ(θ, u), because u∗/(n2 ) will converge towards 1/2 and h′ (θ, u∗ ) will converge towards ρ(u∗ ):

Lemma 6. Let u∗ = argmaxu γ(θ, u) with θ fixed, h′ (θ, u∗ ) converges to ρ(u∗ ) asymptotically as n → ∞.

Lemma 6 implies that the effects of any fixed θ will diminish to a one-dimensional function ρ(u∗ ) as n increases. The shifting of model behavior for different n is closely related to the discussion of instability of ERGM sufficient statistics [Schweinberger, 2011], and more recently the result of ERGM's inconsistency under sampling [Shalizi and Rinaldo, 2013].

Lemma 6 further suggests that for large n, θ needs to be in Θ(n2 ) to be relevant in the model. This property leads to a functional representation of θ:

Here ξ is a "scale-free" meta parameter. This representation will give us a reasonable estimation on the scale of θ, which is necessary perform point-wise MLE for our approximation explained in the next section.

3.4 Approximate Algorithm

The estimation of edge-count induced lower bound immediately leads to an approximation of lnZ(θ): Edge Count Search (ECS) approximation:

Here ξ is the "scale-free" parameter defined in Eq (9). Algorithm 1 reports a straightforward implementation of (10), which simply searches through all the u to maximize ˜ γ(θ, u). Notice that the algorithm requires no extra parameters, which makes the ECS approximation very easy to apply compared to current MCMC sampling methods.

Assume the number of subgraph features r ≪ n, the time complexity of Algorithm 1 is in O(n2 ), which is linear in terms of the number of random variables (i.e. edges) of the model, and quadratic in terms of the order of the network.

A straightforward approximation to the log-likelihood function ℓ(g|θ) is to replace lnZ(θ) with ECS(ξ, n):

The decision of approximating h′ with h∗ in Section 3.2 leads to a simple algorithm. As n → ∞, ECS approximation (10) converges to another closely related approximation proposed by Chatterjee and Diaconis [Chatterjee and Diaconis, 2011], who show that for certain θ the approximation will converge to the true log partition functions. In this paper, we also show that lnZ(θ) will converge to ECS(θ, n) as n → ∞ in Lemma 6.

In next section, we will resort to experiments to verify the effectiveness of the proposed algorithm.

`BR (£103 )

Figure 4: Scatter plot of log-likelihood estimations for ECS and BR on networks of n = 160. Many BR estimations fail UB test (11). Otherwise, ECS and BR estimations are very close (top right).

Algorithm 1 Our new ECS Approximation to the log partition function lnZ(θ)

Input: model parameter θ and number of nodes n Output: estimation of lnZ(θ)

Initialize ECS ← -∞

for u ← 0 to n(n - 1)/2 do

ECS ← max{˜ γ(θ, u), ECS}

end for

4 Experimental Results

In this section, we first use two tasks on synthetic data set to evaluate the performance of ECS approximation: estimating log-likelihood functions and MLE estimation. Then we perform experiments on a real-world social network dataset to show case the quality and stability of our algorithm. We implement the commonly used triad model (edge, 2-star, triangle) for the experiments. To mitigate the degeneration problem in case study, we also used model (edge, altkstar(λ = 1.5)). Here, altkstar is the set of all possible k-stars subgraphs for a network with their weights being set by a function of λ and k. For a fixed n, this feature is mathematically equivalent to including all k-star subgraphs into ERGM, except some constraints on their parameters.

We compare the output of ECS with the state of the art MCMC sampling algorithm for ERGMs: Bridge Sampling [Gelman and Meng, 1998; Handcock et al., 2003; Hunter et al., 2012]. 3 The details of our experimental settings are not described here due to lack of space.

To alleviate the interference of the well known stability problem from sampling based methods on ERGMs [Handcock, 2003; Bhamidi et al., 2008], we employ the following upper bound to the log likelihood function as an indicator of bad approximation:

3 For large n > 20, finding the exact ERGM log partition function is intractable because the number of all graphs with n nodes are O(2n2 ).

The bound holds for any θ and g, because lnZ(θ) must be larger than the log potential of empty graph, which is 0, and of complete graph, which is ∑ri=1 θi. Notice that by design, ECS will never violate the bound. Because for any θ, we have γ(θ, 0) = ˜ γ(θ, 0) and γ(θ, (n2 )) = ˜ γ(θ, (n2 )). We apply this upper bound test (UB test) for all log-likelihood estimations.

4.1 Estimating log-likelihood functions

We sample synthetic networks from a wide range of parameters to evaluate their log-likelihoods. We first generate a 6×6×6 grid of ξ in [-5, 5]×[-5, 5]×[-5, 5], and drop the tuples in which all values have the same sign. We ended up with 162 different ξs. Then for each ξ, we sampled networks for different n ∈ {30, 40, 60, 80, 100, 120, 140, 160}. The total number of sampled graphs is 1,296. We estimate the log-likelihood for each sampled network using both Bridge Sampling and ECS.

Figure 4 reports the scatter plot of the results of both methods for n = 160. Points close to the dashed line suggest ECS and BR produce similar results; points far away from the dashed line suggests the estimation results are very different. For each estimation of BR, we also check whether it exceeds the UB test. If the estimation exceeds the log-likelihood upper bound, we mark the data point with a cross (×); otherwise we mark with a blue circle.

From 4 we can tell when BR estimation fails the UB test, the difference between ECS and BR results are almost negligible. However, there is a significant portion (about 30%) of BR estimation results turn out to be unrealistic, while ECS keeps producing results consistent to (11).

To further compare ECS estimations with the legit BR estimations, we report their relative differences for models that BR estimation pass the upper bound test: reldiff=|(ℓECS-ℓBridge)/ℓBridge|. Figure 5 reports the mean and variance of the relative difference for networks of which BR estimation passes the UB test. The plot shows both the mean and variance decrease as n increases. As n increases, the estimations become very close.

4.2 MLE estimation

In this section, we use ECS as a sub-routine to perform MLE estimations on network data. Because the number of subgraph features in triad model is r = 3, it is practical to perform grid search over a restricted sub-space of ξ. Using real social networks, we compare the performance of this simple ECS-MLE with MCMC-MLE [Snijders, 2002], which uses Bridge Sampling as a sub-routine [Handcock et al., 2003].

We first generated a 6×6×6 grid of ξ in in [-3, 3]×[-3, 3]×[-3, 3]. For each ξ, we sampled one network of n = 60. The purpose of this procedure is not to generate samples that represents the underlying ξ well, but to diversify the sampled graphs. Then we fit the triad model with the sampled network using both MCMC-MLE and ECS-MLE. For ECS-MLE, we performed grid search in a slightly enlarged parameter space with finer granularity.

To evaluate, we estimated the log-likelihood of the network on both fitted models using Bridge Sampling. As we observed in Figure 4, BR tends to generate unrealistic estimations. If a BR estimation fails the UB test, we take the upper bound as the log-likelihood estimation instead. Notice that the adjusted value is still valid for the purpose of comparing two MLE algorithms.

4.3 Case Study on Real Network Data

To study the stability and quality of ECS-MLE results, we fit two different models with four networks with more than 40 nodes: one kapferer2 from statnet; and the other three networks, prison, dolphins, and sanjuansur, from CMU CASOS.

We simulate graphs from the model and compare the node degree distribution to the original network [Hunter et al., 2008]. The quality of the fitted models is measured by the total variation distance 4 .

Both for triad models and alternating k-stars models, we set the initial value of grid search by using Maximum PseudoLikelihood Estimate (MPLE) as in the MCMC-MLE. For triad models, our ECS algorithm searched the grid of two parameters of (edge, 2-star, triangle) ranging the initial value plus offsets in [-5, 5]×[-5, 5]×[-5, 5], with granularity of 0.2. For alternating k-stars, the ECS searched the grid of three parameters of (edge, altkstar) ranging, initial value plus offsets in [-15, -15]×[-15, 15], with the same granularity.

Figure 6 presents the node degree distributions of the simulated network models. The first column shows the network plot of the four networks. The following two columns show the node degree distributions of the triad and the altkstar models fitted by the MCMC-MLE. The final two columns show the node degree distributions fitted by our ECS algorithm. The bold line is the degree distribution for the original network. The curve with error bar is the degree distribution for simulated networks. A good match between the two line suggests the simulated networks implies the quality of the fitted model.

4 ‖f-g‖=supA∈B (f(A)-g(A)), B is a class of Borel sets.

Table 1 shows the total variation distance of node degree distributions of the original and the fitted. We marked the minimum total variation among four approaches for each network. As you can see the ECS-MLE results show better performance compared with the MCMC-MLE results in overall.

5 RelatedWork

Modeling social network structures has been actively studied in machine learning community. Latent variable models, such as matrix factorization [Hoff, 2008], block modeling [Airoldi et al., 2008; Kemp et al., 2006; Ho et al., 2012] and others [Miller et al., 2009; Lloyd et al., 2012], represent the relational data with latent variables. Among those, Ho et al. [Ho et al., 2012] proposed triangular motifs as network representation, which is closely related to ERGM's subgraph features. In comparison, ERGM posts a simple model with intuitive feature specifications that fits for many network analysis tasks.

Computing normalizing constants for complex and highdimensional models, such as ERGMs, is intractable. MCMC simulations are arguably among the most effective methods. Gelman and Meng [Gelman and Meng, 1998] proposed the path sampling formulation to unify acceptance ratio method and thermodynamic integration from theoretical physics for estimating the (ratios of) normalizing constants. Annealed importance sampling (AIS) [Neal, 2001], popular in deep learning literature [Salakhutdinov and Murray, 2008], can also be viewed as one form of thermodynamic integration. Although effective in many applications, Bhamidi et al. [Bhamidi et al., 2008] shows the mixing time for any local Markov chain in low temperature regimes of ERGMs is exponentially slow, rendering these methods computationally intractable in many cases. In comparison, ECS approximation is deterministic, therefore avoids the sampling.

ECS approximation is a variational inference algorithm. In this category, there are many other techniques, such as pseudo-log-likelihood [Strauss and Ikeda, 1990], mean field approximation and Bethe approximation [Wainwright and Jordan, 2008]. In the context of ERGM, these methods have been reported to be inferior to sampling based methods [van Duijn et al., 2009], and are usually used to generated initial states for sampling based algorithms [Hunter et al., 2012]. ECS distinguishes from others by exploiting the asymptotic property in the feature space of the model. This macroscopic view goes beyond the conditional independence in local structures of the model, and may be more effective for complex high-dimensional models like ERGMs.

Many studies have empirically shown that the triad model tends to produce degenerated models when fitted by the MCMC-MLE [Hunter et al., 2008; 2012]. We also observe that the MCMC-MLE fitted triad model was completely degenerated in kapferer2 and sanjuansur. However, our ECSMLE could find an accurate model with no degeneration. To mitigate this problem, the alternating k-stars [Robins et al., 2007; Hunter et al., 2012] have been introduced. In our experiments, the MCMC-MLE fitted altkstar outperforms the MCMC-MLE triad. Although some of the MCMC-MLE result avoided degeneration, ECS-MLE results either with triad and altkstar show the better match to the degree distribution in original network.

6 Conclusions

In this paper, we propose a novel deterministic approximation to the log partition functions of ERGMs. Computing the partition functions (or the ratio of them) is essential in learning ERGMs. Our results show the new method is able to overcome some of the stability issues faced by sampling based methods without losing accuracy. The new algorithm does not depends on extra parameters, making it easy to implement and apply compared to sampling.

We also show that the proposed approximation can be used to build an effective MLE algorithm for ERGMs. In the future, we plan to address various types of MLE problems in EMRGs by using the proposed approximation principles.

Acknowledgments

Wen Pu and Eyal Amir were supported by NSF EAR grant 09-43627, IIS grant 09-17123, IIS grant 09-68552, and a DARPAR grant as part of the Machine Reading Program under AFRL prime contract no. FA8750-09-C-0181. Jaesik Choi and Yunseong Hwang were supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) grant funded by the Ministry of Science, ICT & Future Planning (MSIP) (NRF2014R1A1A1002662), the NRF grant funded by the MSIP (NRF-2014M2A8A2074096). The authors would like to thank David Hunt, Michael Schweinberger, Dan Roth, Gerald DeJong and anonymous reviewers for their constructive suggestions and feedbacks.

A Appendix

A.1 Proof of Lemma 1

Proof. Subgraph count for Li in any g is bounded by 0 ≤ t(g, Li) ≤ t(Kn, Li) ≤ ( n vi )vi!, therefore

A.2 Proof of Lemma 4

Before the proof of Lemma 4, we need some preparations. [Nowicki, 1989] proved that a vector of subgraph counts in G(n, p) are asymptotically normally distributed with a degenerated co-variance matrix with rank 1, as the order of the graph n → ∞. In other words, the subgraph counts are asymptotically linearly dependent on each other. Formally, let φ(g′ ) = {φ1(g′ ), φ2(g′ ), . . . , φr(g′ )} be the densities of subgraphs L1,L2, . . . ,Lr (i.e. φi(g′ ) = t(g′ ,Li) t(Kn,Li) ) for g′ ∈ G(n, p), the sizes (number of edges) of these subgraphs are s1, s2, . . . , sr, and u ∼ Bin((n2 ), p) is the edge count of g′ , we have the following theorem:

Theorem 1. [Nowicki, 1989] For g′ ∈ G(n, p), and real vector a = (a1, a2, . . . , ar)T , the following asymptotic property holds:

where ρ(u, p) = (ρ1(u, p), . . . , ρr(u, p)), and ρi(u, p) = sipsi-1 · u (n2) - (si - 1)psi .

In theorem 1, if we set p = u/(n2 ), then ρi(u, u/(n2 )) = ( u (n2) )si , which becomes the expected density of Li in G(n, p = u/(n2 )).

Next step is to extend the above property from G(n, p) to G(n, M).

Corollary 1. For g ∈ G(n,M = u), as n → ∞, it holds that

Proof. Following theorem 1, let ρi(u) = ρi(u, u/(n2 )), as n →∞, the following holds for g′ ∈ G(n, p = (u/(n2 ))):

Because ∑u p(u) = 1 and p(u) > 0, the claim holds.

Let c be some positive constant, apply Chebyshev's inequality to the linear combination aT φ(g), we get:

Now we start to prove Lemma 4.

Proof of Lemma 4. We first define function ε(u):

As we know

The equality holds if and only if ε(u) = 0. We can get the following property after applying it to corollary 1: as n → ∞

Therefore, as n → ∞

The last step used corollary 1.

We slacks (13) using (15), and rewrite the inner expectation term using (14):

Using (17), we can get

Therefore, apply corollary 1, as n → ∞, we get

A.3 Proof of Lemma 5

Proof. Let Gu be the set of graphs with edge count u, since h∗u is the maximizer, we have

Together with the trivial #(h∗u ) ≤ |Gu|, we can get:

Apply Stirling's approximation on ln |Gu|:

Therefore the claims hold.

A.4 Proof of Lemma 6

Proof. By definition h′ (θ, u∗ ) is vector of densities ranging in [0, 1], therefore the product |θT h′ (θ, u∗ )| is bounded by constant ∑i |θi|.

Let α = maxu H(u/(n2 )) ≈ H(1/2), Lemma 5 shows that

Because Eq (6) implies

as n → ∞ we have:

Let b be some real vector, assume there is some ϵ > 0 so that as n → ∞ we have:

In this case, Lemma 4 implies ln#u∗(h′ (θ, u∗ )) → 0. However, given that ln#u∗(h∗ (u∗ )) is in Ø(n2 ), it contradicts with the definition of h′ (θ, u). Therefore, h′ (θ, u) → ρ(u∗ ).