NISTIR 8066 

Improving the Computational 

Efficiency of the Blitzstein-Diaconis 

Algorithm for Generating Random 

Graphs of Prescribed Degree 

Elizabeth Moseman 

This publication is available free of charge from: 

http://dx.doi.org/10.6028/NIST.IR.8066
NISTIR 8066 

Improving the Computational 

Efficiency of the Blitzstein-Diaconis 

Algorithm for Generating Random 

Graphs of Prescribed Degree 

Elizabeth Moseman 

Applied and Computational Mathematics Division 

Information Technology Laboratory 

This publication is available free of charge from: 

http://dx.doi.org/10.6028/NIST.IR.8066 

June 2015 

U.S. Department of Commerce 

Penny Pritzker, Secretary 

National Institute of Standards and Technology 

Willie May, Under Secretary of Commerce for Standards and Technology and Director 
Abstract
 

When generating a random graph, if more structure is desired than is given in the popular 

Erdos–Renyi model, one method is to generate a degree sequence first then create a graph  ̋

with this degree sequence. Blitzstein and Diaconis[1] (among others) developed a sequential 

algorithm to create a random graph from a degree sequence. This algorithm is assured to 

always terminate in a graph with the desired degree sequence; unfortunately, it is slow. This 

work focuses on the subroutine of the previous algorithm which determines the candidate 

edges, improving the runtime of the overall algorithm from O(mn2) to O(mn). 

Key words: random graphs; degree sequence; graph algorithms 

ii
 
Contents 

1 The Blitzstein–Diaconis Algorithm 1
 

2 Creating the Candidate List Efficiently 2
 

3 Impact 5
 

iii
 
For modeling purposes, one frequently wishes to generate a random graph. Many models for 

random graph generation, including the frequently used Erdos–Renyi model, will generate a graph  ̋

on a specified number of vertices, but the degrees of the vertices in the resulting graph will be 

concentrated among very few values. This makes such models of limited use in modeling real 

world networks, which often exhibit degree distributions that are not concentrated in the same 

manner. In an effort to generate random graphs which more closely resemble real world graphs, 

effort has been made on creating random graphs with a specified degree sequence. 

The simplest (and most computationally efficient) algorithm is that of Bayati, Kim and Saberi[2]. 

In this algorithm, each edge is added sequentially from among the allowable edges. In the analysis 

of the algorithm, it is shown that this algorithm has running time O(mdmax) where m is the number of edges and dmax is the maximum degree of any vertex. One of the main drawbacks of this 

algorithm is that it will not always generate a graph; sometimes, the order of the edge selection 

results in no edges being allowable before all the necessary edges have been generated. The probability of this occurring is shown to be asymptotically small when the maximum degree is bounded 

as a function of the total number of edges. However, from the standpoint of actually using such 

an algorithm, asymptotics are unsatisfying. In addition, this result only holds for graphs where 

dmax = O(m1/4−τ ) with τ a positive constant, and does not hold for the degree sequences that 

arise in many real world networks. 

There are three algorithms that generate random graphs of prescribed degree that ensure successful termination. Mihail and Vishnoi[3] created an algorithm that transforms the problem to 

the well-studied problem of finding a maximum matching. This is done at the expense of making 

the problem bigger: the matching problem is performed on graph with O(n2) vertices, where the 

original graph had only n. 

The most widely used algorithm is that of Gkantsidis, Mihail, and Zegura[4], which is based on 

a Monte Carlo Markov chain. In it, a graph is created which realizes the degree sequence, and then 

edge swaps are performed in order to randomize the graph. The main drawback of this method is 

that there are no results on the mixing time (the number of swaps to perform) in general. 

This technical report focuses on the algorithm of Blitzstein and Diaconis[1]. This algorithm 

also adds edges to the graph sequentially, at each step only a few edges are allowed. More details 

of the algorithm will be given in the following section. The contribution here is a more efficient 

method of calculating the candidates in this algorithm, improving the runtime from O(n2m) to 

O(nm). 

1 The Blitzstein–Diaconis Algorithm 

Beginning with a degree sequence, it is graphical if there is some graph which realizes the degree 

sequence. As edges are added to a graph sequentially, each vertex has a desired degree (the degree 

given by the degree sequence), a current degree (how many edges already contain that vertex) and a 

residual degree (the difference between the desired degree and the current degree). The main idea 

of the Blitzstein–Diaconis algorithm is that if edges are added to the vertex of smallest residual 

degree, and the remaining residual degree sequence is graphical, then the graph can always be 

completed with the desired degree sequence. This is codified as Algorithm 1. 

In this algorithm and later sections, we use the notation 8i,jd to indicate the degree sequence 

obtained from d by subtracting 1 from di and dj but leaving the other entries unchanged. The 

1
 
Algorithm 1: Blitzstein–Diaconis algorithm for generating a random graph with a given 

degree sequence 

Input: a graphical degree sequence d = (d1, . . . , dn)

Output: the edge list E for a graph with degree sequence d

1 E ← ∅

2 while d = 0 do 

3 Choose i with di a minimal nonzero entry. 

4 while di > 0 do 

5 

6 

Compute a candidate list J = {j = i | {i, j} ∈/ E and 8i,jd is graphical}. 

Pick j ∈ J with probability proportional to its degree in d. 

7 

8 

E ← E ∪ {{i, j}}

d ← 8i,jd. 

9 end 

10 end 

largest bottleneck in this algorithm comes from calculating the candidate list, J, before each edge 

is added. Since this is also where the improvements come in, some time should be spent on this 

step. 

To test for graphicality, Blitzstein and Diaconis use the Erdos–Gallai conditions for graphical­  ̋

ity. 

Definition 1.1. For a monotone decreasing sequence d = ) (d1, . . . , dn , the kth Erdos–Gallai  ̋

condition is 

kn k

k

k(k − 1) + min(k, di) − di ≥ 0. (1) 

i=k+1 i=1

This definition differs from the standard in the ordering of the terms, but this usage will make 

later definitions easier. In general, an integer sequence is graphical if it satisfies all of the Erdos–  ̋

Gallai conditions for k = 1, . . . n. 

In their algorithm, Blitzstein and Diaconis utilize a Θ(n) subroutine to test all the Erdos–Gallai  ̋

conditions for each candidate, one at a time. Since there are n potential candidates, the running 

time for generating the candidate list is O(n2) and since a candidate list is generated for every edge 

this makes the total running time O(n2m). Anecdotally, Blitzstein and Diaconis note that it took 

13 seconds to generate graphs for a specified 33 node degree sequence on a 1.33 GHz machine 

when the algorithm was coded in R. Matlab1 code for the same procedure took 0.05 seconds on 

the same 33 node degree sequence. 

2 Creating the Candidate List Efficiently 

To improve the running time of the candidate calculations, we need to observe that the values 

obtained from calculations of the Erdos–Gallai conditions for  ̋ d will be very similar to the values 

1Certain commercial equipment, instruments, or materials are identified in this paper to foster understanding. Such 

identification does not imply recommendation or endorsement by the National Institute of Standards and Technology, 

nor does it imply that the materials or equipment identified are necessarily the best available for the purpose. 

2
 
of the same condition for 8i,jd. We will see later that the values differ by at most two, and it is 

easy to tell by exactly how much they differ. But first, we need to address the issue of order in the 

degree sequence. 

In testing the Erdos–Gallai conditions efficiently, we assume that the sequence fed into the  ̋

tester is ordered from greatest to least. This is important because there is a Θ(n log n) lower bound 

on a sorting algorithms when nothing is known of the sequence. However, the sorted order of d

and the sorted order of 8i,jd may differ. Assume that d is in sorted order. Notice that trailing zeros 

do not effect whether d is graphical or not, so we assume that n is the number of nonzero entries 

and thus dn is a minimal nonzero entry. 

Definition 2.1. Let d be a sorted, graphical degree sequence with no trailing zeros. We say that 

j ≤ n − 1 is a candidate if 8j,nd is graphical. 

We first address the case where 8j,nd is not in sorted order. Notice that since there are no entries 

less than dn, decrementing to dn − 1 makes no difference to the sorting, so the only possibility is 

dj − 1 < dj+1. Since dj ≥ dj+1, this forces dj = dj+1. If, however, j = n − 1, then even 

dj = dj+1 = dn will not cause 8j,nd to be unsorted. 

Lemma 2.2. If dj = dj+1 with j < n − 1, then j is a candidate if and only if j + 1 is a candidate. 

Proof. This is a necessary consequence of the sorted order of 8j,nd. In particular, in sorted order 

8j,nd and 8j+1,nd are the same sequence so if the common sequence is graphical then both are 

candidates and otherwise both are not. 

The following theorem addresses the candidacy of the remaining j values. 

Theorem 2.3. Let d = (d1, . . . , dn) be a sorted, graphical degree sequence with dn > 0. Fix k

with 1 ≤ k < n and 1 ≤ j < n. Define d

 ̃= 8j,nd, 

kn k

k

ek = k(k − 1) + min(k, di) − di and (2) 

i=k+1 i=1

kn k

k

e ̃k = k(k − 1) + min(k, d

 ̃

i) − d

 ̃

i. (3) 

i=k+1 i=1

Suppose, without losing any generality via Lemma 2.2, dj > dj+1 or j = n − 1 so that d

 ̃ is also 

sorted. Then 

⎧

⎪⎪ 1 j ≤ k and dn > k ⎪⎪⎨

0 j ≤ k and dn ≤ k

e ̃k − ek = 0 j > k and dn > k (4) 

−1 j > k and dj > k ≥ dn

⎪⎪⎪⎪⎩ −2 j > k and k ≥ dj

Proof. We observe that e ̃k and ek are very similar. We consider two cases. 

Case j ≤ k: Then 

e ̃k − ek = min(k, dn − 1) − min(k, dn) + 1 (5) 

3
 
 

where the 1 is due to the change in dj to dj−1. If dn > k, we have min(k, dn−1) = min(k, dn) = k

so that e ̃k − ek = 1. Otherwise, min(k, dn − 1) = dn − 1 and min(k, dn) = dn so that e ̃k − ek = 0. 

Case j > k: We see that 

e ̃k − ek = min(k, dn − 1) − min(k, dn) + min(k, dj − 1) − min(k, dj ). (6) 

When dn > k, the sorted nature of d ensures dj ≥ dn > k so that min(k, dn − 1) = min(k, dn) =

min(k, dj − 1) = min(k, dj ) = k and e ̃k − ek = 0. When dj > k ≥ dn, we have min(k, dn − 1) =

dn − 1, min(k, dn) = dn, and min(k, dj − 1) = min(k, dj ) = k so that e ̃k − ek = −1. Finally, 

when k ≥ dj ≥ dn, we have min(k, dj − 1) = dj − 1 and min(k, dj ) = dj so that e ̃k − ek = −2. 

Algorithm 2: Calculating the max candidate.
 

Input: a sorted, graphical degree sequence d = (d1, . . . , dn), with dn > 0

Output: a max candidate j' with 8j,nd graphical for any j ≤ j'

1 Calculate ek for k = 1, . . . n

2 for k = 1 to n do 

3 if ek ≥ 2 then 

4 jk ← n

5 else if ek = 1 then 

6 jk ← max(k, min(j | dj ≤ k))

7 else if ek = 0 then 

8 if dn ≤ k then 

9 jk ← k

10 else 

11 jk ← max(k, min(j | dj ≤ k))

12 end 

13 end 

14 end 

15 j' ← mink(jk)

16 while dj' = dj'+1 do 

j' ← j 17 ' − 1

18 end 

Each condition will have a maximum index j so that 8j,nd passes that condition, and the 

minimum of these indices will be the maximum candidate j'

. Then the candidates are J = {j ≤

j' | {j, n} ∈/ E and dj = dj'+1}. This process is shown in Algorithm 2. This is then a Θ(n)

process for calculating the candidate list. If we use the theorem to adjust the ek once a candidate 

has been selected, the constant may be improved but we cannot calculate ek in better than linear 

time because all of them may need to be updated. 

4
 
Figure 1: Running time for trials with the unaltered (green) and altered (blue) algorithms. Average 

of 10 trials at each n. As a reference, a single trial with the unaltered algorithm with n = 10000

ran in 5300 s, or approximately 88 minutes. 

3 Impact 

To show the practicality of the algorithm, degree sequences were chosen to follow a power law, 

so that the probability that a given vertex has degree k is proportional to k−β for some constant β. 

We chose to work with β = 2.1 since many graphs found in practice seem to have this exponent. 

Elementary calculations yield the results that the expected degree on these graphs is β−1 and the 2−β

expected maximum degree is at least n1/(β−1) and thus the Bayati, Kim and Saberi algorithm [2] 

will converge for β > 5, but not in the range that we are interested in. Figure 1 shows the average of 

10 running times for random degree sequences generated according to a power law with β = 2.1.

2 

2Code was generated in Matlab and run on a 2.4GHz AMD Athlon 64 processor. 

5
 
References 

[1] J. Blitzstein and P.	 Diaconis (2010). “A sequential importance sampling algorithm for 

generating random graphs with prescribed degrees”. Internet Mathematics 6 (4), 489. 

http://dx.doi.org/10.1080/15427951.2010.557277. 

[2] M.	 Bayati, J. H. Kim, and A. Saberi (2010). “A sequential algorithm for generating random graphs”. Algorithmica 58 (4), 860–910. 

http://dx.doi.org/10.1007/s00453-009-9340-1. 

[3] M. Mihail and N. Vishnoi (2002). “On generating graphs with prescribed vertex degrees for 

complex network modelling”. In ARACNE. 

[4] C. Gkantsidis, M. Mihail, and E. Zegura (2003). “The markov chain simulation method for 

generating connected power law random graphs”. In ALENEX. 

6