Efficient Network Generation Under General Preferential Attachment

Preferential attachment (PA) models of network structure are widely used due to their explanatory power and conceptual simplicity. PA models are able to account for the scale-free degree distributions observed in many real-world large networks through the remarkably simple mechanism of sequentially introducing nodes that attach preferentially to high-degree nodes. The ability to efficiently generate instances from PA models is a key asset in understanding both the models themselves and the real networks that they represent. Surprisingly, little attention has been paid to the problem of efficient instance generation. In this paper, we show that the complexity of generating network instances from a PA model depends on the preference function of the model, provide efficient data structures that work under any preference function, and present empirical results from an implementation based on these data structures. We demonstrate that, by indexing growing networks with a simple augmented heap, we can implement a network generator which scales many orders of magnitude beyond existing capabilities ($10^6$ -- $10^8$ nodes). We show the utility of an efficient and general PA network generator by investigating the consequences of varying the preference functions of an existing model. We also provide"quicknet", a freely-available open-source implementation of the methods described in this work.


INTRODUCTION
There is a clear need for scalable network generators, as the ability to efficiently generate instances from models of network structure is central to understanding both the models and the real networks that they represent. Ideally, researchers of communication and social networks should be able to generate networks on the same scale as the real networks they study, and many interesting networks, such as the World Wide Web and Facebook, have millions to billions of nodes. Furthermore, network generation is the primary tool both for empirically validating the theoretical behavior of models of network structure and for investigating behaviors that are not captured by theoretical results. The generation of very large networks is of particular importance for these tasks because theoretically derived behavior is often asymptotic.
However, the generation of large networks is difficult because of its high complexity. In the case of preferential attachment (PA), arguably the most widely used generative model of networks, a non-local distribution over node degrees must be both sampled from and updated at each time-step. If we naively index this distribution, we will need to update every node at every time-step, which implies that generating a network will have complexity of at least O(|V | 2 ).
PA models are of particular interest because they account for the scale-free distribution of degree observed in many large networks [15]. For instance, scale-free degree distributions have been observed in the World Wide Web [5,6,7], the Internet [9,12,20], and telephone call graphs [1,2], bibliographic networks [11] and social networks [18].
Preferential attachment models generate networks by sequentially introducing nodes that prefer to attach to nodes with high degree. While many extensions to this model class exist, all members share the same basic form: At each timestep, sample a node from the network with probability proportional to its degree; introduce a new node to the network; and add an edge from the new node to the sampled node. This behavior has important implications for implementation. First, PA is inherently sequential, because the next action taken depends on the state of the network, and the state of the network changes at each time-step. This implies that the algorithm is not easily parallelized. Second, network nodes must be indexed such that they can be efficiently sampled by degree, and, because we are introducing a new node at each time-step, the index must also support efficient insertion. Third, the relevant distribution over nodes is non-local, in that the introduction of a new node and edge affects the probability of every node in the network through the normalization factor.
Much of the work in modeling network structure has focused on the asymptotic regime. A model is defined, and a limiting degree distribution (as |V | approaches infinity) is obtained analytically. Less effort has been focused on generating finite networks. In the following sections, we provide a robust framework for generating networks via PA. This framework easily scales to millions of nodes on commod-Price(n,λ): existing node ← sample in degree(G, λ) new node = add node(G) add edge(G, new node, existing node) end for return G Krapivsky(n, p, λ, µ): new node = add node(G) add edge(G, new node, existing node) else existing node tail ← sample in degree(G, λ) existing node head ← sample out degree(G, µ) add edge(G, existing node tail, existing node head) end if end while return G Figure 1: Generating a network with n nodes under Price and Krapivsky's models. G0 is some small seed network. λ and µ are scalers which give the fitness of nodes for incoming and outgoing edges, respectively.
ity hardware. We also provide "quicknet", a freely-available open-source C implementation of the framework 1 .
The remainder of the paper is structured as follows. In Section 2, we present an analysis of the complexity of generating networks from PA models. Section 3 describes candidate methods for efficiently implementing preferential attachment generators and presents results from a simulator which implements them. Section 4 describes several applications of a PA network generator which scales to many millions of nodes. We describe related work in Section 5 and present conclusions and future work in Sections 6 and 7, respectively.

COMPLEXITY
In this section, we describe the two existing PA models as examples. We then provide a more formal definition of a PA model, which is followed by an analysis of the complexity of generating networks from PA models. Figure 1 describes Price's algorithm. Briefly, at each timestep, a node is sampled from the network with probability proportional to its in-degree; a new node is introduced to the network; and a directed edge is added from the new node to the sampled node. Notice that a node is added at each time-step, so that the generation of a network with |V | nodes takes |V | steps. At each step, the algorithm of Price's model is followed with probability p, and a "preferential edge step" is taken with 1 https://github.com/hackscience/quicknet probability 1 − p. During a preferential edge step two nodes, no and ni, are sampled from the network by out-and indegree, respectively, and an edge is added from no to ni. Note that a node is no longer added at every step; rather, a node is added at a given step with probability p. This implies that the number of iterations required to generate a network with |V | nodes is a random variable with expected value |V |/p. |V |/p is Θ(|V |) ∀p, so asymptotically this is no different than Price's model. More generally, the number of iterations required to generate a network with |V | nodes via a PA model is Θ(|V |).

Definitions
Let Gt = (Vt, Et) be the network that results from t iterations of a PA simulation. Vt is the set vertices (or nodes) within the network and Et is the set of edges between elements of Vt. Let T (Gt) be the worst-case time complexity of generating Gt; that is, the worst-case time complexity of a preferential attachment simulation of t iterations.
Recall that the number of iterations required to generate a network with |V | nodes via PA is Θ(|V |). Accordingly, we will omit t and frame our discussion of complexity T (G) in terms of |V |.
Let A = {a1, a2, ..., a |A| } be a set of attributes that can be defined on a network node. Let Xv = {xva 1 , xva 2 , ..., xva |A| } ∈ R |A| be a setting of A for node v ∈ V , and let λva i ∈ R be the fitness of node v for attribute ai. Let be a set of functions, where fa i ∈ f maps xva i ∈ R and λva i ∈ R to a preference mass µva i ∈ R + . The "preference mass" µva i is a non-negative real value that is proportional to the probability of selecting v by ai under the PA model. We will refer to the elements of f as the "preference functions" of the PA model.
A PA model has one or more preference functions. Price's model, for example, has a single linear preference function. Krapivsky's model has two: one for in-degree and another for out-degree. A "linear preferential attachment model" only admits linear preference functions of the form g(x, λ) = c1x + λ, a "quadratic preferential attachment model" only admits quadratic preference functions of the form g(x, λ) = c2x 2 + c1x + λ, and so on.

Derivation of Generation Complexity
We obtain a trivial lower bound on T (G) by noting that, in order to generate G, we must at the very least output |V | nodes, so T (G) = Ω(|V |).
A discussion of the upper bound follows. Recall that the salient problem in generating networks from a PA model is indexing the network's nodes in such a way that sampling, insertion, and incrementation can be accomplished efficiently. Tonelli et al. [19] provide a clever method for accomplishing all three tasks in constant time, provided that the preference function is linear and the fitness is both uniform across all nodes and constant. Given constant insertion and sampling, the generation of a network with |V | nodes takes O(|V |) time. Considering that the lower bound is Ω(|V |), we have the asymptotically tight bound of T (G) = Θ(|V |).
However, this method does not extend to nonlinear preferential attachment. We can improve performance by shifting to data structures which provide O(log|V |) insertion, sam-  Figure 2: General algorithm to sample from the augmented tree structure. η is the mass observed thus far and u is a sample from the standard uniform distribution.
pling, and incrementation, giving an overall complexity of We accomplish this with a set of augmented tree structures. Each tree supports a preference function of the model by indexing the preference mass assigned to each node in the network by that preference function. Each item in the tree indexes a node in the network. The tree items are annotated with the preference mass of the network node under the preference function, and the subtree mass, which is the total preference mass of the subtree that has the item as root; see Figure 3. Note that we refer to "items" in the tree rather than the more typical "nodes"; this is to avoid confusion between elements of the tree and elements of the network. We can sample from such a structure by recursively comparing the properly normalized subtree mass of a given item and its children to a uniform random draw; see Figures 2 and 3. Note that, at each iteration of a standard PA simulation, we must sample a node, update that node's mass, and insert a new node. In what follows we show that each of these steps can be accomplished in asymptotically logarithmic time.

IMPLEMENTATION
The tree structure that we described in the previous section can be implemented in a number of different ways that have a generation time of O(|V |log|V |). They differ in their computational time for finite |V |. In this section, we empirically evaluate a set of realizations of the annotated tree structure. Specifically, we investigate a simple binary maxheap where priority is defined by node mass, and a set of binary treaps with various sort and priority keys.
Note that, in the discussion of the heap-based and treapbased implementations of the tree structure, we will often refer to a "sort invariant" and a "heap invariant". The sort invariant states that, for any three nodes Y ← X → Z where Y and Z are the left and right children of parent X, respectively, and a "sort key" k that is associated with each item, Y.k ≤ X.k ≤ Z.k. The heap invariant states  that for any three nodes Y ← X → Z (defined in the same fashion) and some "priority key" p associated with each item, X.p ≥ Y.p and X.p ≥ Z.p.
We first describe the binary maximum heap. We annotate each item in the heap with a node mass, which is defined by the preference function, and a subtree mass, which is initialized to the node mass. When inserting a new item i, i's node mass is added to all traversed items, so that the subtree mass is remains accurate upon insertion. Sampling is accomplished via the algorithm of Figure 2. Node mass may only increase, so we implement an augmented version of increase-key which maintains subtree mass under exchanges; see Figure 4 for a diagram of the exchange operation. The increase-key operation supports the Increment operation, which is described below. We set priorities to be equivalent to node masses so that the most probable nodes can be accessed more quickly.
The PA process is supported by the binary maximum heap via the operations Sample, Increment and Insert. As previously mentioned, sampling is performed via the algorithm of Figure 2. Increment increases an item's mass and then performs heap exchanges to account for any violation of the heap invariant; it is described in Figure 5. Insert adds an item to the index, appropriately updating the subtree masses of any parent items; see Figure 6.
Note that, if we were to annotate each item with a node's probability mass rather than preference mass, insertion would be a linear time operation. When a new node is introduced, the probability of every existing node decreases because the normalization factor increases. Thus, upon insertion, every item's probability mass would need to be updated. There

Insert(heap, item):
heap.add(item) node mass ← item.node mass while item has a parent do item ← parent(item) item.subtree mass += node mass end while Figure 6: The Insert operation of the heap-based tree structure. Note that each item has O(log|V |) ancestors, so Insert is a O(log|V |) operation.
are |V | items, so insertion becomes a Θ(|V |) operation in this situation. Conversely, the preference mass of each node is unaffected by the introduction of a new node. Insertion in this scenario is a O(log|V |) operation; see Figure 6.
We use Price's model as an illustrative example. Recall that, in Price's model, a new node is introduced at each time-step, and an edge from the new node to an existing node is added preferentially. We first identify an existing node via Sample. We then create a new node and add an edge from the new node to the existing node. Increment is called on the existing node to reflect the change in preference mass due to the new incoming edge. Finally, the new node is added to the index via Insert. Sample, Increment and Insert are O(log|V |) operations, which implies that a single iteration is O(log|V |), and that a simulation with |V | iterations is O(|V |log|V |). Generating a network with |V | nodes takes Θ(|V |) iterations, so T (G) = O(|V |log|V |).
The treap-based implementations were designed to make more efficient use of space. The heap was implemented via a dynamic array, which provides amortized constant-time insertion at the cost of some wasted space. We sought to avoid this wastage by instead using binary treaps, which are extensions of binary trees that maintain heap invariants over random priorities to ensure balance in expectation [3]. We find that the treap-based implementation consistently underperformed the heap-based implementation, so we devote less space to its description. Figure 7 shows the empirical run time of each of these structures as a function of generated network size. All networks were generated from Krapivsky's model. We find that the binary heap consistently took significantly less time than the treap-based methods to generate networks of several different sizes.

APPLICATIONS
We validate our generation model by generating sets of networks from the Krapivsky model and comparing the marginal degree distributions inferred from the generated networks with the asymptotic value predicted by the model. We then use the generator to explore some interesting questions. Specifically, we analyze the effect of changing the fitnesses of the Krapivsky model from a constant value to a random variable with various distributions. We also analyze the robustness of Krapivsky's model to superlinear preference functions.

Validating the Network Generator
We validate our framework by comparing the inferred exponents of the marginal distributions of generated networks with the known (theoretical, asymptotic) values for the exponents. We generated 10 networks with 10 7 nodes each. Figure 8 shows a plot of the base-10 logarithm of both degree an complementary cumulative distribution. The exponents of the marginal degree distributions were inferred via linear regression. We find, as expected, that they both exhibit power-law behavior (evident in the linearity) and that the inferred exponents of the distributions are in relatively good agreement with asymptotic theoretical values. Note that, while networks with 10 7 nodes are very large, they are still finite; we believe that this accounts for the small discrepancy between the inferred exponents and the theoretic values.

Pareto Fitness
We use our network generator to investigate the effects of altering the Krapivsky model. Specifically, we generated networks from a variant where the fitnesses assigned to each node were sampled from a Pareto distribution, rather than assigning the same constant value to each node. Results can be seen in Figure 8. The distribution of in-degree fitness is λd λ m d λ+1 and has expected value λdm λ−1 . The parameter dm is set to (λ − 1) so that the expected value of the distribution simplifies to λ. The same form was used for the out-degree fitness. Note that this variant still exhibits scale-free behavior, that the inferred exponents are in better agreement with the predicted values than the exponents inferred from the simulation of the unaltered model, and that the variance of the inferred exponents is higher.

Normal Fitness
We also simulated a variant of the Krapivsky model where fitnesses were sampled from a truncated normal distribution. Results can be seen in Figure 8. In-degree fitnesses were sampled from N(λ, (λ/4) 2 ) and out-degree fitnesses from N(µ, (µ/4) 2 ). The variances were chosen such that the probability of sampling a negative fitness is very small (less than 10 −4 ); the distributions were truncated so that any negative samples were replaced with zero. Note that scale-free behavior is still observed and that the inferred exponents of the marginal distributions of in and out-degree are in very close agreement with the simulation of the original model.

Robustness to Superlinear Preference Functions
Super-linear preference functions increase the strength of the "rich-get-richer" effect. This can lead to situations where one node quickly overtakes all others and is thus a component of most of the edges in the network. In the extreme case, a star will form; all edges will be connected to the outlier node. We investigate the robustness of Krapivsky's model to super-linear preference functions by plotting the ratio dmax |E| , were dmax is the maximum degree, as a function of the preference function exponent α; see Figure 9. dmax |E| will approach 1 as the network approaches a star formation.
There is an interesting side effect to the transition from scale-free to star-structured networks. As the network becomes more star-like, the probability of selecting the most probable node tends to increase. The most probable node always sits at the top of the heap, so it can be accessed in constant time. So, the closer a network's structure is to a star formation, the larger the probability that an iteration of a PA algorithm will be constant time. For a star structured network in the limit, every iteration will be constant time and the generation of a network with |V | nodes will be Θ(|V |). This behavior is apparent for finite |V |; we have observed that the runtime of the generator tends to decrease as α increases.

RELATED WORK
This work is concerned with the problem of efficiently generating networks from PA models. Some examples of PA models include the models of Price (directed networks with scale-free in-degrees) [16], Barabasi and Albert (undirected networks with scale-free degrees) [4], Krapivsky et al. (directed networks with non-independent in and outdegrees which exhibit marginally scale-free behavior) [14], and Capocci et al. (like Krapivsky's model, but with reciprocation) [8].
There has been some prior work in efficiently generating networks from PA models. Tonelli et al. [19] provide a method for computing an iteration of the linear Yule- gives the proportion of edges that involve the maximum degree node. dmax |E| = 1 indicates a star formation. Note that there is a phase transition from scale-free to star-structured networks between α = 1.0 and α = 1.2. 100 networks with 10 6 nodes each were generated for each value of α. Error bars indicate the 95% confidence interval.
Simon cumulative advantage process in constant time. This method can naturally be extended to network generation through linear PA. However, the extension to nonlinear PA (not shown due to space constraints) is very inefficient in both time and space. Ren and Li [17] describe the simulation of a particular linear PA model, RX, but do not address the general problem of simulating networks from models with general preference functions. Hruz et al. [13] and D'Angelo and Ferreti [10] provide methods for parallelizing the simulation of linear PA, but do not treat the nonlinear case. To the best of our knowledge, our work is the first to address the efficient generation from PA models under possibly nonlinear preference functions.

CONCLUSION
We provide an efficient framework for simulating preferential attachment under general preference functions which scales to millions of nodes. We validate this framework empirically and show applications in the generation and comparison of large networks.

FUTURE WORK
We have shown that, for nonlinear preferential attachment, the complexity of generating a network with |V | nodes is both Ω(|V |) and O(|V |log|V |). Future work could provide asymptotically tighter bounds.

ACKNOWLEDGMENTS
This work was supported by MURI ARO grant 66220-9902 and NSF grant CNS-1065133.