Unified Representation and Lifted Sampling for Generative Models of Social Networks

1 Introduction

Many complex systems are modeled using networks-where nodes represent entities and edges represent some type of relation between the entities. An open problem in network science is how to create generative models that accurately capture the characteristics of real world, sparse networks, in order to better understand the properties of the systems, and to facilitate statistical analysis.

Recent research has explored various models to accurately characterize network structure (e.g., [Chung and Lu, 2002; Hoff, 2008; Leskovec et al., 2010; Benson et al., 2014]). A number of these generative network models (GNMs) use a common procedure to generate the edges in the network using a matrix of probabilities. The approaches differ in the number and structure of parameters used to specify the edge probabilities. While several GNMs can capture important characteristics of real world networks such as power law degree distributions and community structure, it has been more difficult to develop methods to capture longer-range dependencies that impact global characteristics. To address these limitations, new hierarchical GNMs have been proposed with more complex dependencies between the edge probabilities, e.g., mixed Kronecker product graph model (mKPGM) [Moreno et al., 2013], Block two-level Erd˝os-R´enyi Model (BTER) [Seshadhri et al., 2012], bipartite stochastic block model (biSBM) [Larremore et al., 2014], and hierarchical graph models [Peixoto, 2014].

Moreover, while GNMs specify a generative process, it can still be computationally intensive to naively sample large networks. There have been efficient sampling methods proposed for some GNMs [Leskovec et al., 2010; Yun and Vishwanathan, 2012], but because models are specified procedurally it is difficult to guarantee correctness, and consequently the efficient sampling methods can generate improbable network structures [Moreno et al., 2014]. Thus, it is still an open question how to develop accurate and efficient sampling methods. Another challenge is to identify commonalities across GNMs. Again, since models are specified procedurally, it is difficult to discern the impact of algorithmic differences-both w.r.t. the number and structure of parameters-except by comparing the structure of the networks that are generated from the methods.

In this paper, we show that hierarchical GNMs can be abstracted as Bayesian networks (BNs) with random variables (RVs) representing the existence of edges in the network, and hierarchies of latent RVs representing relations among groups of edges. With this representation we can identify common properties of unrelated models, such as mKPGM and BTER, and use the knowledge to guide the design of new GNMs. The BN representation also facilitates the development of a general sampling method that is both efficient and accurate. The transformation of a GNM to BN representation makes it easy to identify parametric symmetries allowing us to sample groups of edge RVs rather than sequentially-using insights from lifted inference [Poole, 2003; Jha et al., 2010; Van den Broeck and Darwiche, 2013; Mittal et al., 2015]. Moreover, we can maximize computational efficiency and guarantee correctness by using insights from context specific independence [Boutilier et al., 1996; Poole and Zhang, 2003] to dynamically sparsify the set of RVs to be sampled based on the context of parent RV values.

Previous work had provided important insights about the relationship among some GNMs (e.g. [Jacobs and Clauset, 2014]) and BNs have been used to model GNMs (e.g.

[Schein et al., 2016; Liang et al., 2016; Ho et al., 2011; Neiswanger et al., 2014]), However, to our knowledge, this work is the first attempt to model a family of existing GNMs in a principled way that is general enough to facilitate the creation of new classes of GNMs. We do that by taking advantage of BNs to model GNMs but without enforcing restrictions on the topology of a BN except for a hierarchical conditional dependence which we describe in this paper.

In summary, the key contributions of this work are:

1. Demonstration that hierarchical GNMs can be represented in a universal, succinct way that supports the development of new models.
2. Development of a general efficient, and provably correct, sampling algorithm that applies to all GNMs in the family of hierarchical GNMs. Our algorithm is based on two important properties:
3. (a) Deterministic context specific dependence-which sparsifies the RV sampling space.
4. (b) Parametric symmetries-which, due to GNM parameterizations, can be exploited with lifted sampling.
5. Design, and evaluation, of a new GNM to illustrate the utility of the BN representation; in particular, of the related sampling algorithm.

2 Background: Generative NetworkModels

First, we review some details about generative network models. Let G=(V, E) be a network with set of vertices V and edges E ⊂ V × V. We define a generative network model (GNM) M with parameters Θ as follows.

Definition 1. Generative network model (GNM)

A GNM is a statistical model M with parameters Θ that define, either explicitly or via a construction process, a size |V|×|V| matrix P of probabilities. Each cell [i, j] ∈ P corresponds to the binary random variable Aij and the value Pij represents the Bernoulli probability that the [un]directed edge ei,j exists in the network (i.e., if Aij =1 then ei,j ∈ E and P (Aij =1) = Pij). Thus, P models the structure of the network through the set of binary random variables Aij ∀i, j ∈ {1, . . . , |V|}.

The size of Θ differs for each GNM and can vary in the range [1, |V|2 ]. This restriction in the number of parameters is intended to make sampling of random graphs efficient. However, most GNMs have a small number of parameters (usually |Θ|≪|V|) to avoid overfitting. Thus, multiple cells in P (and RVs A) could have the same probability values. We use < U,T > to represent this, where U is the set of unique probabilities that appear in the matrix P, i.e. U=unique(P)={π1, π2, . . . , πu , . . . πκ}, and T is the set of the list of positions Tuwhere each of the πu appears in P.

2.1 Basic Edge-Oriented GNMs

There is a range of basic GNMs such as Erd˝os-R´enyi Model (ER) [Erdos and Renyi, 1960], Stochastic Block Model (SBM) [Holland et al., 1983], Chung-Lu Model (CL) [Chung and Lu, 2002], etc. that build P with some mathematical operation over a set of parameter(s). For space reasons we will define only CL here since is relevant for our later discussion.

Chung-Lu Model (CL) [Chung and Lu, 2002]. In a CL model, Pij=wiwj/ ∑ wk for a sequence of expected nodedegrees w=(w1,. . .,w|V|). There are then, for the case of undirected networks, at most |V|(|V|-1)/2 probabilities πi in P (and associated U), in which case T would consist of lists of a single element each.

2.2 Hierarchical GNMs

Hierarchical GNMs are a super-classs of GNMs that use a hierarchical process to introduce dependencies among sets of edge RVs. While a basic GNM defines the matrix P directly, a hierarchical GNM uses intermediate latent variables at a certain level in the hierarchy to model the variable interactions of the next level in the hierarchy, modeling the entries in P indirectly. Thus, there are two levels of randomness associated with sampling in hierarchical GNMs: one for the intermediate (hierarchical) RVs that impact the generation of P, and one for edge sampling. Notice that basic GNMs can be regarded as hierarchical GNMs with a single layer. Most hierarchical GNMs are defined iteratively or procedurally, but not all iterative or procedural GNMs are hierarchical. We denote the levels in the hierarchy by λ ∈ {0, 1, . . . , φ - 1}, where φ is the number of levels in the hierarchy, λ=0 represents the root of the hierarchy, and λ =φ - 1 the edge layer. B[λ] is the adjacency matrix associated with the blocks at layer λ.

Block two-level Erd˝os-R´enyi Model (BTER) [Seshadhri et al., 2012]. In a preprocessing step, BTER creates groups of nodes based on the sequence of expected degrees. Then, sampling is done in two steps: (1) An ER model is used to link nodes within groups, where the probability is proportional to the smallest node-degree in the group; (2) A CL model is used to create link between groups, with probability of edges is proportional to the excess degree from (1). Step 1 ensures that edges within-blocks have higher probability than betweenblocks. Step 2 ensures a power law degree-distribution and since it uses a CL model, U & T have the same worst case scenario as CL.

Mixed Kronecker Product Graph Model (mKPGM) [Moreno et al., 2013]. mKPGM is a GNM that overcomes known limitations (see e.g. [Seshadhri et al., 2013]) of stochastic Kronecker graphs [Leskovec et al., 2010]. Given a b×b parameter matrix Θ, K and ℓ, it samples a network as follows.

1. Compute P [0] , a bℓ × bℓ matrix, equal to ℓ-1 Kronecker products of Θ with itself.
2. Construct G[0] from P [0] by sampling each edge RV independently ∼ Bernoulli(P [0] ij ).

3. For l=1 . . .K-ℓ: Set P [l] =G[l-1] ⊗Θ and sample G[l] . Note, |V[0] |= bℓ and |V[K-ℓ] |= bK. Since G[0] . . .G[K-ℓ-1] represent auxiliary graphs, where each edge influences a block of possible edges in the next iteration of the hierarchical sampling process (and they are not present in the final output network); we refer to them as block adjacency matrices: B[0] . . .B[K-ℓ-1] . To simplify notation we will refer to the final sampled network G[K-ℓ] =(V[K-ℓ] , E[K-ℓ] ) as G=(V,E). In mKPGMs, the size of U is b2 for each level l>0 of the hierarchy and T consists of b2 lists of edges.

3 Representing GNMs as Bayesian Networks

In this section, we introduce a common BN representation for GNMs with hierarchical dependencies among the edgeprobabilities, and develop an associated, efficient sampling method. We also identify common properties of the GNMs based on their BN form and exploit these properties for efficient and provably correct sampling.

The BN representation is comprised of a hierarchical structure of edges and blocks. We denote the levels in the hierarchy by λ ∈ {0, 1, . . . , φ - 1}, where φ is the number of levels in the hierarchy, λ=0 represents the root of the hierarchy, and λ=φ-1 the edge layer. A block j at level λ has an associated RV Z [λ] j which represents its state: sampled (Z [λ] j =1) or not (Z [λ] j =0). Since blocks could have different sizes we keep a one dimensional subindex for each block RV (Zi) rather than use the rectangular structure (with two subindices Zij) of some GNMs (e.g., mKPGM in Fig. 2). Likewise, P could be regarded as a one-dimensional structure for this same purpose (notice that P(B[λ] j =1) = P [λ] j is the probability of sampling block j as defined by the GNM). This illustrates that our algorithm is easily applicable to hierarchical models with complex topology of the interacting layers and not just to tree-like BNs. The function pa (Z [λ] j ):=Z [λ-1] i for λ>0 returns the parent

RVs for Z [λ] j . We denote as Z[λ] the set {Z [λ] j }. It is important to notice that the BN is comprised of Z[0] , . . . ,Z[φ-1] and that the parameters of the BN are the associated probabilities. These probabilities are a function of, or derived from, P [λ] . However, this relation is not straightforward for all GNMs.

3.1 Transformation of GNMs to BNs

As discussed in the previous section, hierarchical GNMs are a superclass of GNMs that include edge-based GNMs as a special case, as well as more complex models that use a hierarchy of latent variables to represent dependencies among edge RVs. However, since most hierarchical GNMs are defined iteratively, the hierarchical structure is not immediately evident and a transformation is necessary. In particular for procedural models, transforming a GNM to BN form consists of two steps, as outlined in Figure 1: (1) Restructure the output of each step in the iterative sampling of a GNM as a set of RVs. This allows the process to be reorganized as levels in a hierarchy of RVs. (2) Use the hierarchy of RVs to build a BN using parameters for the BN derived from the GNM. Some care should be exercised for this transformation as the relations among RVs are not trivial.

General Transformation from Hierarchical GNM to BN A hierarchical model M, with parameter matrix Θ, can be represented as a BN N with parameters Θ′ obtained from Θ: MΘr⇀NΘ′ as follows:

1. Represent as Z [0] j the RV that models the block edge B[0] j
2. Define the probability of Z [0] j : P(Z [0] j =1) = P [0] [j]

(b) Specify pa (Z [λ]

1. For λ = 1, . . . , φ-1 (a) Represent as Z [λ] k the RV that models the block B[λ] k

k

) = Z [λ-1]

k′

to be the corresponding

P (Z [λ]

) = 0) = 0

1. Add the RVs Z[0] , . . . , Z[φ-1] (from 1 & 3.a) to the BN N and add their associated CPDs (from 2 & 3.b) to Θ′ .

k = 1|pa (Z [λ]

k

Notice that when φ (number of levels) is 1 then the hierarchical GNM collapses to a traditional edge-based GNM. Thus, hierarchical GNMs include edge-based GNMs as a special case. Finally, Θ′ is fully defined by Θ because all P are derived from it. We now discuss two specific transformations.

Transforming mKPGM to BN representation

|V[0] |

Figure 2: Example mKPGM in BN form. The left subplot, descendants for 1 RV per level show the relations among the variables. The right subplot shows the BN using plate notation.

An mKPGM model M with parameter matrix Θ, K, and ℓ (Sec. 2.2) can be represented as a BN N with a tree structure and parameters Θ′ obtained from Θ: MΘ r⇀NΘ′ as follows.

1. Represent as Z [0] ij the RV that models the block edge B[0] ij ∀i, j ∈[1, bℓ ]
2. For λ = 1, . . . , φ, where φ =K - ℓ do:
3. Define the CPD of Z [0] ij : P(Z [0] ij =1) = P [0] [i, j] for all i, j, where P [0] =⊗ℓk=1Θ
4. (a) Represent as Z [λ] kl the RV that models the block edge B[λ] kl ∀k, l∈[1, bℓ+λ ]
5. (b) Specify pa (Z [λ] kl ) = Z [λ-1] ⌈k/b⌉,⌈l/b⌉ to be the parents of RV Z [λ] kl , where x=mod (k-1, b)+1 and y=mod (l1, b)+1, then, define the CPD of Z [λ] kl , ∀k, l∈[1, bℓ+λ ]: P (Z [λ] kl = 1|pa (Z [λ] kl ) = 1) = Θ[x, y] P (Z [λ] kl = 1|pa (Z [λ] kl ) = 0) = 0
6. Add the RVs Z[0] , Z[1] , ..., Z[λ=K-ℓ] (from steps 1 and 3.a) to the BN N and add their associated CPDs (from steps 2 and 3.b) to Θ′ .

An example BN representation of an mKPGM is shown in Fig. 2 for λ=0, 1, 2, 3. Here λ=0 corresponds to B[0] in the mKPGM sampling process. There is a total of (bℓ )2 =|V[0] |2 RVs each of them represented by a Z [0] ij (a double subindex is used to indicate the position of the RV in the block/edge matrix). Each of the RVs has b2 descendants at λ>0.

Transforming BTER to BN representation

A BTER model M with parameters Θ=({nd}d∈N, {cd}d∈N), where {nd}d∈N is a sequence of target degrees and {cd}d∈N is the target clustering coefficient per degree, can be represented as a BN N with a tree structure and parameters Θ′ obtained fromΘ: MΘr⇀NΘ′ as follows.

As discussed earlier, BTER creates groups of nodes with possibly the same degree and uses an ER model to generate blocks of edges linking those nodes together. These are called affinity-blocks. To produce a power law degree distribution, there are many affinity-blocks of edges linking nodes with small degree and fewer blocks for nodes with larger degree. Recall that each node has an expected degree; then, bulk nodes are nodes where the minimum degree of its affinity block is equal to its expected degree. Otherwise, if the minimum degree is smaller, they are called filler nodes. Contrary to mKPGMs, BTER has an heterogenous number of descendent edges. The nodes of degree one have no edges between them because they are used as filler nodes; i.e. they are used to link blocks of degrees > 1.

Let nb, rd, db be the number of edges in affinity block B[0] ij , the ratio of the number of filler nodes to the total nodes, and the minimum degree in the block, respectively. Let, Bin-1 be the inverse CDF of a binomial. Then, the mapping of BTER sampling to BN representation is as follows:

1. Represent as Z [0] ij the RVmodeling the affinity blocks B[0] ij
2. Obtain the parameters nb, rd, db of BTER as described in Algorithm 1 of [Kolda et al., 2014]
3. Define the CPD of Z [0] ij : P(Z [0] ij = 1) = P [0] [i, j] ∀ i, j where P [0] is probability of sampling edges from degree equal to the degree of B[0] ij
4. Let pa(·) be the affinity block containing ·. Define P (Z [φ] kl =1|pa(Z [φ] kl )=0)=0 ∀k, l ∈ [1, . . . , |V|]. Otherwise, for edges ek,l s.t. pa(Z [φ] kl ) = 1 and b its affinity block with clustering coefficient cdb
5. (a) Represent its associated RV as Z [φ] kl

• (c) Else: Given σtype = rd 2 if type = filler, and σtype = 1-rd 2 if type = bulk, define the CPD of Z [φ] kl : P(Z [φ] kl = 1|pa(Z [λ] kl ) = 1) = Bin-1 (∑i∈Vd σtypei [di - (c1/3 db · db)]) for all k, l and nodes i ∈ Vd of degree d
• Add the RVs Z[0] ,Z[1] to the BNN(steps 1, 5.a.i, 5.b) and their associated CPDs (steps 3, 4, 5.a.2, & 5.c) to Θ′ .

As we can see BTER and mKPGM are similar in the sense that both models constrain the sampling so very few edges are sampled with higher probabilities. This allows for skewed degree distribution. The difference in the depth of the hierarchies is now easily seen, mKPGM groups edges in blocks and super-blocks to abstract their relations as separate latent SBM models, while BTER uses shallower hierarchies for redistribution of probabilities according to a power-law.

3.2 GNM-BNCommon Properties

We identify 2 properties of GNMs transformed to BN form.

1. Group Symmetries-The first key insight about our BN representation is that symmetries appear in the BN from the parameterization of the GNM. Since there are fewer than |V|2 unique Θs, RVs can be grouped and sampled together. Moreover, since the models are designed to generate sparse networks, it is more efficient to sample the number of edges and then select their locations rather than sequentially sample each edge RV, as proved in [Moreno et al., 2014]. As in earlier lifted techniques, that were based on finding identical structures for which a common computation was performed only once ([Koller and Pfeffer, 1997; Pfeffer et al., 1999]), we propose an algorithm that identifies RVs with the same conditionals and group them together in a single representation. Then, sampling will be done only once by using a Binomial (this determines the number of edges to sample) and randomly selecting the edges from the locations with the same unique probability (as illustrated in Fig. 3). Because in GNMs the symmetries are parametric similarities among a set of independent RVs there is no need to assess the dependencies among individual RVs or among groups of RVs.
2. DCSD Sparsification-The second property of the BN representation for GMNs is deterministic context-specific dependence (DCSD). As seen in the earlier transformations, this property arises because some values of the variables are sampled with value 0. Since sampling at the next layer of the hierarchy is conditional on the sampling of the previous layer, the sampling space becomes more and more sparse as we move through the hierarchy, as shown in Figure 4. Because a

BN form of Fig. 1 b

Figure 3: Sampling of GNM-BN: Symmetries of RV values are exploited for efficiency and correctness.

Algorithm GNM-BN-Sampling

1: Input: N 2: Ouput: G = (V, E) {defined by its adjacency matrix A} 3: V = {1, . . . , |V|} 4: Obtain the set U[0] and the list T[0] using N 5: for λ = 0, . . . , φ-1 do 6: for i = 1, . . . , |U[λ] | do 7: Take πi = U[λ] i {unique prob in P[λ] } 8: Take τi = |T[λ] i | {# of positions in P[λ] where πi appears} 9: Sample Nbi ∼ Bin(τi, πi) 10: --→ pos =RandomNbi positions fromT[λ] i 11: if λ < φ - 1 then 12: Set Z[λ] j = 1, for all j ∈ --→ pos 13: Obtain the set U[λ+1] and the list T[λ+1] using Z[λ] 14: else 15: Set -→u =⌈ --→ pos |V| ⌉ and -→v=mod(--→ pos+|V|-1, |V|) + 1 16: Set Au,v =1 where u= -→uw, v= -→vw for w = 1, . . . , Nbi

Line 4 obtains the set of unique probabilities U[0] and the list T[0] (i.e. the corresponding positions where those probabilities appear). This information is used for sampling at the root layer λ = 0. The algorithm consists of using each πi ∈ U[λ] and the corresponding list of positions for randomly sampling blocks. Instead of sampling from a Bernoulli for each position of the block, we sample from a binomial which models the distribution of successes in |T[λ] i | Bernoulli trials with probability πi. Sampling from this distribution gives the number of sampled blocks (Nbi ) for layer λ whose positions are randomly picked and stored in the vector --→ pos. This vector is used to sample the blocks, i.e. to set their RVs Z[λ] j to one, where j ∈ --→ pos (lines 11-13). Sampling of edges (lines 15-16) follows the same idea except the positions --→ pos are transformed to two indices: the row and column, u ∈ -→u , and v ∈ -→v of the adjacency matrix.

Complexity: The average-case complexity of the algorithm is O(φ|E|), where φ is the depth of the hierarchy (typically a constant) and |E| = E[|E|]. This can be proved as follows. Lines 9-10 have average-case complexity O(E[|Z[λ-1] |]). This complexity is O(τiπi) (τi=# of positions in P [λ] where πi appears). The average complexity of the loop for (λ<φ-1)

value in any given U[λ] may be 0, we don't need to sample at all for the associated T . Thus this sparsification can be used for more efficient sampling. Formally:

Definition 2. Deterministic context-specific dependence (DCSD): Let N be a BN that generates a network G through a hierarchical sampling process. Let P (Z [λ] j ) be the probability that a block is sampled. Then N is "DCSD" iff it partitions all Zj, such that:

Remark 1. Z [λ] j is sampled iff its parent is sampled with value 1. However, if a superblock is not sampled, sampling of subblocks or edges is inhibited. Note that while P(Euv|pa(Euv)=0)=0, the marginal P(Euv)>0.

DCSD is related to context-specific independence (independence in a BN due to specific realizations of RVs) [Boutilier et al., 1996], but in this case the specific context precludes further sampling of subblocks and edges. We show below how these two properties can be used to develop an efficient and provably correct general sampling algorithm for the GNM-BNs.

3.3 GNM-BN Sampling

RV Reorganization: As we described before, a GNM can be represented as a set of unique probabilities πu ∈ U and their positions Tu. We can apply this organization to the RVs Z[0] , ...,Z[λ=φ-1] (and their probabilities) in N from the transformations Mr⇀N explained before. Thus, N can be represented by the set of ordered pairs {(U[λ] , T[λ] )}φ-1 λ=0 , where U[λ] ={π1, ..., πu , . . .}=unique(P(Z [λ] )) and T[λ] is a vector with the list of positions T[λ] u where each of the probabilities πu appear at level λ.

Algorithm: We propose an algorithm for sampling from models trasnformed to BN form that is efficient (as proved below). Our algorithm is easily applicable to hierarchical models with complex topology. The input of the algorithm is the BN N and its parameters. N has an associated set of unique probabilities U[λ] and list of positions T[λ] as we described above.

...

...

...

...

is approximately O(E[∑φ-1 λ=0 |Z[λ] |])=O(φ|E|) (φ =# layers) Since the average-case complexity for λ = φ - 1 is O(|E|).

Scalability: Because of the linear (in the number of edges) average case complexity of the algorithm, a BN-GNM is scalable for sufficiently sparse networks. Additionally, a BNGNM is amenable to be parallelized. After symmetries are identified, sampling could be separated in independent tasks due to the independence of the parameters.

Lemma 1. Given a valid model in BN representation NΘ′ with probability PN , and the algorithm GNM-BN-Sampling with probability PS , then ∀G[λ] and 0 < λ ≤ φ - 1:

Lemma 2. Given a valid model in BN representation NΘ′ with probability PN , and the algorithm GNM-BN-Sampling with probability PS , then ∀G[λ] and 0 < λ ≤ φ - 1:

Proof. By lemma 1 and sampling independence of Z[λ] then PS (Z[λ] |pa(Z[λ] )) = ∏∀k,l PS (Z [λ] kl = 1|pa(Z [λ] kl )) =

Theorem 1 (Correctness). Given a valid model in BN representation NΘ′ with probability PN , and the algorithm GNMBN-Sampling with probability PS , then ∀G PN (G) = PS (G).

Proof. PN (G)=∑G[λ-1]∈Γ[λ-1] PN (G[λ] |G[λ-1] )PN (G[λ-1] ) for the set Γ of all possible graphs at λ - 1. By Lemma 2 PN (G) = ∑G[λ-1]∈Γ[λ-1] ,Z[λ] PS (Z[λ] |pa(Z[λ] ))PN (G[λ-1] )

By the same procedure, it can be proven recursively that: PN (G[1] )=∑G[0]∈Γ′[0] ,Z[1] PS (Z[1] |pa(Z[1] ))PN (G[0] ) for the set Γ′ of all possible graphs at level 0. Since PS(pa(Z[0] )) = 1, by Lemma 1: PN (G[1] ) = ∑G[0]∈Γ′[0] ,Z[1] PS (Z[1] |pa(Z[1] ))PS (pa(Z[1] )). Then, ∀G PN (G)=PS (G)

Remark 2. Given a valid mKPGM model M(Θ) with probability PM , there is a BN N such that the algorithm GNM-BNSampling, with sampling probability PS , samples correctly from PM , i.e. ∀G PS (G) = PM (G).

Let PS be the probability of the GNM-BN sampling. Given the mapping MΘr⇀NΘ′ , as shown in the transformation of mKPGM to BN. Since the CPD are the same for the sampling and the BN representation, then, by Thm 1, ∀G PM (G)= PN (G)=PS (G)

Remark 3. Given a valid BTER model M(Θ) with probability PM , and the algorithm GNM-BN-Sampling with probability PS , then from Thm. 1 it also follows that ∀G PM (G) = PS (G).

4 Designing NewGNMs: a Case Study

To take advantage of the properties of the BN representation, we propose a new GNM, namely the BN-Graphlet model. We will show that the new GNM can model graph characteristics that are not possible with several existing GNMs. We consider undirected networks but BN-Graphlet could be easily extended for directed graphs as well.

In order to take advantage of our framework, and in particular our sampling algorithm, we design a GNM that is directly interpretable as a BN hierarchy with common parameters per layer and with a generative process that exploits sampling conditioned on the previous layer. The key idea in BN-Graphlet is to redistribute symmetries and sparsity based on several graph patterns, to convey (1) locations where edge probability mass is concentrated in the desired graphs, (2) primitives regarding interactions among nodes in the graph, and (3) hierarchical dependencies among primitives. To do so we consider (a) the degree distribution, (b) minimal graphlets (mingraphlets, i.e. triangles, wedges, and edges) that encode node interactions per degree, and (c) a graph decomposition due to the mingraphlets.

The reason to use mingraphlets is that larger graphlets can be constructed via algebraic transformations of the minimal graphlets. Since each of these characteristics are dependent on the previous one, each will be encoded in a level of the BN hierarchy.

Let ∆ ∼ fdeg(θ∆) be the sequence of degrees drawn from the degree distribution fdeg. BN-Graphlet creates groups of nδi nodes of degree δi for δi ∈ ∆. Then we can partition the adjacency matrix via the joint degree distribution, and use block bij =B[0] ij to refer to links among nodes with degrees δi, δj respectively. Let P [0] [i, j]=0 for blocks with no links; otherwise, P [0] [i, j]=1.

Next, the model partitions a single block in subgroups of nodes that link only triangles, only wedges, and the remainder, i.e., every block will be subdivided into nine subblocks B[1]

B[1]

·· . To simplify notation, we use

k=[1..3]l=[1..3] (i, j) to refer to appropriate cell in B[1]

associated with its parent i, j in B[0] . Then: P [1] [k, l] = ⎧⎨⎩ p△bij if k=1, l=1 i.e., both nodes link only to triangles p∧bij if k=2, l=2 i.e., both nodes link only to wedges 1 o.w.

Here p△bij and p∧bij are respectively: the probability that the two nodes incident to an edge in block B[0] ij participate only in triangles, and the probability that the two nodes incident to an edge in block B[0] ij participate only in wedges.

Lastly, the final adjacency matrix is constructed from submatrices of size nδi × nδj for each degree pair δi, δj. The nodes of each sub-matrix will be associated with roles based on their parent B[1] kl (i, j). For 3x3 submatrices of edges em,n, where m, n have the role only triangles, the off-diagonal edge probabilities in B[2] mn will be all 1; for blocks where the nodes have the role only wedges, there will be two edges with probability 1s and one with probability 0. For all other edges, the probability will be κ1/3 bij , where κbij is the clustering coefficient centered at a node in block B[0] ij .

It is straightforward to see that the BNN of BN-Graphlet is directly derived from the blocks B[λ] ij and their probabilities, contrary to mKPGM and BTER where transformations were needed. This is because we use our framework directly to build the GNM. The equivalence can be directly obtained with the following procedure:

Estimation: We need to estimate θ∆, the sequence κbij , and probabilities p△bij and p∧bij. θ∆ is learned as the MLE of the fitted degree distribution. κbij is the clustering coefficient per block of the data graph. The parameters p△bij and p∧bij are estimated as the MLE of the respective multinomial.

Figure 5: Variation of mKPGM and BN-Graphlet graph properties in synthetic networks. (LCC: largest connected component)

4.1 Experiments and Analysis

We performed two sets of experiments. First, we encoded BN-Graphlet using the new representation and show, using synthetic networks, that it can model different networks compared to mKPGM. Because mKPGM was shown in [Moreno et al., 2013] to be capable of modeling a wider range of network characteristics than Chung Lu and KPGM models, in our experiments we compare to mKPGM. We generated networks for a wide range of parameter values Θ and plotted the characteristics of the sampled graphs. For mKPGM, we used 22,060 different values of Θ (b = 2). The parameters were generated using every possible combination of Θ, such that θ11, θ12 ∈ {0.01 : δ : 1.00}, θ12 = θ21, θ22 ∈ {θ11 : δ : 1}, and 2.1 ≤ ∑ij θij ≤ 2.4 We utilized δ = 0.015 and θ12 ∈ {0.01 : δ : 1} to avoid repetition of the parameters due to isomorphisms of Θ. For each Θ setting, we generated 75 undirected networks with mKPGM (K=9, ℓ=6). For BN-Graphlet we considered θ∆ that lead to networks with a number of edges in the range [800, 2650] for a fair comparison with mKPGM models. For κbij we generated random values under the restriction that the global clustering coefficient is realistic. The values of p△bij and p∧bij were assigned using grid search.

We evaluate the characteristics of the generated networks using: (1) average cluster coefficient, (2) average geodesic distance, (3) number of non isolated nodes, and (4) size of the largest connected component. Figure 5 reports the results. BN-Graphlet produces the lowest geodesic distance (small world phenomena), the highest cluster coefficient (community structure), and largely reduces the number of isolated nodes (larger connected component). Thus BN-Graphlet can model networks that are not easily modeled with mKPGM.

Second, we show that BN-Graphlet can model real world networks better than the other models. We fitted the models to three real datasets: the CoRA citations network (comprises 11,881 AI papers with 31,482 citations between them), the National Longitudinal Study of Adolescent Health (AddHealth)

network (1155 nodes and 7884 edges), the astrophysics arXiv network (9,987 nodes and 25,973 edges), and Facebook wall postings (449,748 nodes and 1,016,621 edges).

We compared the BN-Graphlet with BTER and mKPGM. Figure 6 shows the 3-dimensional Kolmogorov-Smirnov (KS3D ) distance of the learned to the true network characteristics: hop-plot, degree, and clustering coefficient. In all datasets, BN-Graphlet obtains the lowest error, except in Facebook where its error is the second lowest; showing it can jointly model the hop-plot, degree, and clustering coefficient of real networks consistently better than BTER or mKPGM.

5 RelatedWork and Discussion

GNMs have been developed in various ways, using different motivations and frameworks. There has been some work to identify the common features of GNMs. For example, [Jacobs and Clauset, 2014] provides an overview of different families of GNMs and gives some important insights about relationships among them. [Rohe et al., 2011] considers the relation between the latent space model [Hoff et al., 2002] and stochastic block models [Holland et al., 1983]. However, these works do not consider complex GNMs with a hierarchical sampling process involving intermediate latent variables.

Previous work have used BNs to create GNMs [Schein et al., 2016; Liang et al., 2016; Ho et al., 2011; Neiswanger et al., 2014]. However, these works proposed a specific BN with the purpose of solving a specific network problem (e.g., node actions, link predictions, etc). Our contribution is to propose a general representation that is not dependent on any dataset, problem, or BN topology, but rather only on an assumption of conditional independence during sampling. Using this minimal building block, we have proposed a sampling method that is universal across GNMs that fulfill this assumption.

Also related to our work is lifted inference. Lifted inference algorithms identify and exploit abundant symmetries in graphical models, in order to avoid repeated computations and speed up probabilistic inference [Koller and Pfeffer, 1997; Poole, 2003; Jha et al., 2010; Van den Broeck and Darwiche, 2013; Mittal et al., 2015; Sen et al., 2009]. Lifted Importance Sampling (LIS) was proposed in relational learning for probabilistic inference [Gogate and Domingos, 2011], and extended by [Gogate et al., 2012]. We note that their task is different since the symmetries are exploited to improve the precision of LIS (i.e., reduce variance) while our work exploits parametric symmetries for efficient network sampling, which provides guaranteed performance (both time complexity and correctness) for GNMs. [Venugopal and Gogate, 2014] showed a way to deal with situations where symmetries are broken. As in their work, we show ways to find symmetries but in the realm of social network models.

In general, like these existing lifted-techniques, our lifted sampling identifies structural patterns in GNMs to sample a set of random variables as a group rather than treating each random variable separately. This leads to an efficient and correct sampling method. Unlike existing lifted-techniques, we use symmetries to create a meta-model that encodes existing social network models and allows for creation of new ones.

In our work, we take advantage of the symmetries in the intrinsic graphical model of the GNMs for efficient sampling and then exploit the hierarchical representation to construct a generative model that produces sampled graphs with characteristics that are not easily achieved with current GNM methods. In the experiments, we show the range of networks that can be modeled with this new generative model. To our knowledge, our work is one of the first applications of lifted techniques that exploits the symmetries existing in social network models.

6 Conclusion

In this paper, we showed how hierarchical GNMs that are specified procedurally, can be represented as generally as BNs. Our unifying BN view provides the following advantages: (1) universal representation, which can highlight the similarities and differences between GNMs, and (2) efficient, provably correct and universal sampling, due to the adherence of the algorithm to a directed BN structure and an efficient sampling algorithm. The representation also facilitates the creation of new GNMs. To illustrate the benefits of the BN representation and associated sampling method, we proposed the hierarchical BN-Graphlet model and showed it more accurately captures the structure of four real-world networks.

Acknowledgements

We thank the anonymous reviewers for their useful comments. This research is supported by NSF under contract numbers: IIS-1546488, IIS-1618690, CCF-0939370, and by "CONICYT + PAI/Concurso nacional de apoyo al retorno de investigadores/as desde el extranjero, convocatoria 2014 + folio 82140043."