arXiv:1709.02674v2 [math.CO] 14 Nov 2017

Uniform generation of random graphs with power-law
degree sequences
Pu Gao∗
School of Mathematics
Monash University
jane.gao@monash.edu

Nicholas Wormald†
School of Mathematics
Monash University
nick.wormald@monash.edu

Abstract
We give a linear-time algorithm that approximately uniformly generates a random
simple graph with a power-law degree sequence whose exponent is at least 2.8811.
While sampling graphs with power-law degree sequence of exponent at least 3 is fairly
easy, and many samplers work efficiently in this case, the problem becomes dramatically
more difficult when the exponent drops below 3; ours is the first provably practicable
sampler for this case. We also show that with an appropriate rejection scheme, our
algorithm can be tuned into an exact uniform sampler. The running time of the exact
sampler is O(n2.107 ) with high probability, and O(n4.081 ) in expectation.

1

Introduction

We consider the following problem. Given a sequence of nonnegative integers d with even
sum, how can we generate uniformly at random a graph with degree sequence d? Motivation
to do this can come from testing algorithms, to test null hypotheses in statistics (see Blitzstein
and Diaconis [2] for example), or to simulate networks. If the degree sequence comes from
a ‘real-world’ network then it will often follow a power law (defined precisely below), and
often the ‘exponent’ of this power law is between 2 and 3. We give a linear-time algorithm
for generating such graphs approximately uniformly at random in the sense that the total
variation distance between the output a distribution from the uniform one is o(1) as the
number of vertices goes to infinity. There has been no such approximate sampler before with
good guaranteed performance, i.e. small error in distribution and low time complexity, for
power-law degree sequences with γ below 3. We will also give an exact uniform sampler, for
γ less than but sufficiently close to 3, which has a reasonably low time complexity: n2.107
with high probability and n4.081 in expectation.
Research on uniform generation of graphs with prescribed degrees has a long history. Tinhofer [20] studied the case of bounded regular degrees and described an algorithm without
bounding how far the output is from the uniform distribution. Soon afterwards, a uniform
∗
†

Research supported by NSERC and ARC DP160100835.
Research supported by ARC DP160100835.

1

sampler using a simple rejection scheme (described in Section 2) arose as an immediate
consequence of several enumeration
√ methods [4, 5]. This algorithm works efficiently only
when the maximum degree is O( log n), where n is the number of vertices. A major breakthrough that significantly relaxed the constraint on the maximum degree was by McKay and
Wormald [17]. Their algorithm efficiently generates random graphs whose degree sequence
satisfies ∆ = O(m1/4 ), where ∆ denotes the maximum degree and m is the number of edges
in the graph. The expected running time of this algorithm is O(∆4 n2 ). Very recently, the
authors of the current paper√developed a new algorithm which successfully generates random
d-regular graphs for d = o( n), with expected running time O(nd3 ). This new algorithm is
related to the McKay-Wormald algorithm but contains significant new features which make
it possible to cope with larger degrees without essentially increasing the runtime.
While the uniform generation of graphs whose degrees are regular, or close to regular,
has already been found challenging, generating graphs whose degrees are very far from being
regular seems much more difficult. However, a desire to generate such graphs materialised in
recent years. Because of the important roles of the internet and social networks in modern
society, much attention has been paid to graphs with real-world network properties. One
of the most prominent traits of many real-world networks is that their degree distribution
follows the so-called power law, usually with parameter γ between 2 and 3 (i.e. the number of
vertices with degree i is roughly proportional to i−γ ). Graphs with such degree distributions
are sparse but have vertices with very large degrees. There is more than one definition of
a power law degree sequence examined in the literature, and we consider the common one
resulting when the degrees are independently distributed with a power law with parameter
γ. For this, the maximum degree is roughly n1/(γ−1) .
Generating random graphs with power-law degrees for γ > 3 is quite easy using the
above-mentioned simple rejection scheme or the McKay-Wormald algorithm. However, a
radical change occurs when the power-law parameter γ drops from above 3 to below. A
notable difference is that the second moment of the degree distribution changes from finite
to infinite. As a result, the running time of the simple rejection scheme changes from linear
to super-polynomial time. The running time required for the McKay-Wormald algorithm
also becomes super-polynomial when
γ < 3. Another major difference is the appearance of
√
vertices with degree well above m. To date, no algorithm has been proved to uniformly
generate such graphs in time that is useful in practice, with any reasonable performance
guarantee on the distribution.
Some known approximation algorithms have polynomial running times without being
provably practical. Jerrum and Sinclair [14] gave an approximate sampler using the Markov
Chain Monte Carlo method for generating graphs with so-called P-stable degree sequences.
Sufficient conditions for P-stability and non-P-stability are given by Jerrum, McKay and
Sinclair [13], but none of these conditions are met by power-law degree sequences with
2 < γ < 3. On the other hand, the authors found enumeration
formulae for power-law
√
degree sequences [8] that imply P-stability when γ > 1 + 3, and one can show using a
related but simpler argument that P-stability holds for all γ > 2. However the degree of the
polynomial bound for the running time of this algorithm is too high to make it suitable for any
practical use with a provable approximation bound. Recently Greenhill [11] used a different
Markov chain, an extension of the one used for the regular case in [7], to approximately
generated non-regular graphs. The mixing time of that chain is bounded by ∆14 m10 , and
2

the algorithm only works for degree sequences that are not too far from regular. In particular,
it is not applicable to power-law degree sequences of the type we consider that have γ below
3. In [2] it is explained how to use sequential importance sampling to estimate probabilities
in the uniform random graph with given degrees, using statistics from non-uniform samplers.
However, there are no useful theoretical bounds on the number of samples required to achieve
desirable accuracy of such estimates. As stated in [2, Section 8], “Obtaining adequate bounds
on the variance for importance sampling algorithms is an open research problem in most cases
of interest.” Using a uniform sampler completely nullifies this problem.
In this paper, we present a linear-time approximate sampler PLD* for graphs with powerlaw degree sequences, whose output has a distribution that differs from uniform by o(1) in
total variation distance. This approximate sampler, like others such as in [19, 16, 3, 22],
has weaker properties than the ones from MCMC-based algorithms. Its approximation error
depends on n and cannot be improved by running the algorithm longer. However, it has
the big advantage of a linear running time. None of the previously proposed asymptotic
approximate samplers [19, 16, 3, 22] is effective in the non-regular case when ∆ is above
m1/3 , whereas the type of power-law degree sequences we consider can handle ∆ well above
m1/2 . So PLD* greatly extends the family of degree sequences for which we can generate
random graphs asymptotically uniformly and fast.
Using an appropriate rejection scheme, we tune PLD* into an exact uniform sampler PLD
which, for γ > 2.8811, has a running time that is O(n4.081 ) in expectation, and O(n2.107 )
asymptotically almost surely (a.a.s.). Most of the running time of PLD is spent computing
certain probabilities in the rejection scheme. The reason that the expected running time
is much greater than our bound on the likely running time is that our algorithm has a
component that computes a certain graph function. This particular component is unlikely
to be needed on any given run, but it is so time-costly that, even so, it contributes immensely
to the expected running time.
We next formally define power law sequences. There have been different versions in the
literature. The definition by Chung and Lu [6] requires the number of vertices with degree
i to be√O(i−γ n), uniformly for all i. Hence the maximum degree is O(n1/γ ), which is well
below n. However, it has been pointed out [18] that this definition leads to misleading
results, as in many real networks there are vertices with significantly higher degrees. A more
realistic model [18, 21] is to consider degrees composed of independent power-law variables.
In this model the maximum degree is of order n1/(γ−1) , which matches the statistics in the
classical preferential attachment model by Barabási and Albert [1] for modelling complex
networks. In this paper, we will use a more relaxed definition [8] which includes the latter
version of power-law sequences. Let X be a random variable
power law distribution
P∞ with
−γ
−1
−γ
with parameter γ > 1, i.e. P(X = i) = ci , where c = i=1 i , then for γ > 1 we have
P(X ≥ x) = O(i1−γ ). As vertices with degree 0 can be ignored in a generation algorithm,
we assume that the minimum degree is at least 1.
Definition 1. A sequence d = (d1 , . . . , dn ) is power-law distribution-bounded (plib for short)
with parameter γ > 1 if the minimum component in d is at least 1 and there is a constant
K > 0 independent of n such that the number of components that are at least i is at most
Kni1−γ for all i ≥ 1.
In plib sequences, the maximum element can reach as high as n1/(γ−1) , making both
3

enumeration and generation of graphs with such degree sequences more challenging than the
Chung-Lu version mentioned above. It is easy to see that if d is composed of independent
power-law variables with exponent γ then d is a plib sequence with parameter γ 0 for any
γ 0 < γ.
Our main result for the approximate sampler PLD*, defined in Sections 2 and 3.4, is as
follows.
Theorem 2. Assume d =P(d1 , . . . , dn ) is a plib sequence with parameter γ > 21/10 +
√
61/10 ≈ 2.881024968 and 1≤i≤n di is even. Then PLD* runs in time O(n) in expectation
and generates a graph with degree sequence d whose distribution differs from uniform by o(1)
in total variation distance.
A general description of our uniform sampler, which we call PLD, will be given in Section 2, and its formal definition in Sections 3 and 6. Our result for the exact uniform sampler
is the following.
Theorem
3. Assume d = P
(d1 , . . . , dn ) is a plib sequence with parameter γ > 21/10 +
√
61/10 ≈ 2.881024968 and 1≤i≤n di is even. Then PLD uniformly generates a random
simple graph with degree sequence d, with running time O(n4.081 ) in expectation, and O(n2.107 )
with high probability.
We remark that the same statements are immediately implied for the Chung-Lu version
of power law sequences (which we called power-law density-bounded in [8]), however for that
version one can easily modify our methods to reduce the bound on γ somewhat. This
√ is
because the maximum degree in the Chung-Lu version is much smaller (well below n),
and moreover, all multiple edges are likely to have bounded multiplicity. Hence, the major
difficult issues we face in this paper do not appear in the Chung-Lu version of power law
sequences.
We remark that our argument easily gives slightly improved time complexity bounds for
any particular γ in the range given in Theorem 3; to simplify the presentation we have taken
a uniform bound on a complicated function of γ.
We describe the general framework of both PLD and PLD* in Section 2. The formal
definition of PLD is given throughout Sections 3 and 6, together with key lemmas used for
bounding the running time. We focus on the design of the algorithm including its definition,
and the specification of various parameters involved in the algorithm definition, leaving many
of the straightforward proofs to the Appendix, particularly those similar to the material in [9].
We prove Theorems 3 and 2 in Sections 4 and 5.

2

General description

As in [17] and [9], we will use the pairing model. Here each vertex i is represented
as a bin
P
containing exactly di points. Let Φ be the set of perfect matchings of these 1≤i≤n di points.
Each element in Φ is called a pairing, and each set of two matched points is called a pair.
Given a pairing P , let G(P ) denote the graph obtained from P by contracting each bin into
a vertex and represent each pair in P as an edge. Typically, G(P ) is a multigraph, since it
can have edges in parallel (multiple edges) joining the same vertices. An edge that is not
multiple we call a single edge. The pairs of P are often referred to as edges, as in G(P ). A
loop is an edge with both end points contained in the same vertex. The graph is simple if it
4

has no loops or multiple edges.
The simple rejection scheme repeatedly generates a random P ∈ Φ until G(P ) is simple.
Let Gd denote the set of simple graphs with degree sequence
d. An easy calculation shows
Q
that every simple graph in Gd corresponds to exactly ni=1 di ! pairings in Φ. Hence, the
output of the rejection scheme is uniformly distributed over Gd . The problem is that, when
the degrees become very non-regular, it takes too many iterations for the rejection scheme to
find a simple G(P ). With power-law degree sequences as in Theorems 2 and 3, the rejection
scheme is doomed to run in super-polynomial time.
The McKay-Wormald algorithm [17] does not reject P if G(P ) is a multigraph. Instead,
it uses switching operations (e.g. see Figures 1 and 2) to switch away the multiple edges.
These switchings slightly distort the distribution of the pairings away from uniform, and
so some rejection scheme is used to correct the distribution. The algorithm in [9] included
major rearrangements in set-up and analysis in order to incorporate new features that reduce
the probability of having a rejection, and thereby extended the family of degree sequences
manageable. Both [17] and [9] only deal with multiplicities at most three. However, it is
easy to see that for a plib degree sequence with parameter γ < 3, there are multiedges
of multiplicity as high as some power of n. Previous exact generation algorithms have
never reached the density at which unbounded multiplicity occurs in the pairing model,
and new considerations for the design of PLD* and PLD stem from this. One may think
that switchings for low multiplicities such as in Figures 1 and 2 work for any multiplicity.
However, these switchings would cause difficulties for high multiplicities because (a) the error
arising in the analysis is too hard to control, and (b) incorporating unbounded multiplicities
into the counting scheme would result in super-polynomial running time.
One new feature in the present paper is that we treat vertices of large degree differently.
We will specify a set H of such vertices, which we call heavy. Other vertices are called light.
An edge in P is heavy if its both end vertices are heavy. In particular, a loop is heavy if it
is at a heavy vertex. Non-heavy edges are called light. The algorithm PLD* contains two
stages. In the first stage, it uses some switchings, defined in Section 3.2, to turn every heavy
multiple edge ij into either a non-edge, or a single edge, and to switch away all heavy loops.
It can be shown that with a non-vanishing probability, the final pairing is in the following
set


L(P ) ≤ BL ; D(P ) ≤ BD ; T (P ) ≤ BT ;
A0 = P ⊆ Φ : G[H] (P ) is simple,
,
G(P ) contains no other types of multiple edges
(1)
where G[H] (P ) denotes the subgraph of G(P ) induced by H, L(P ), D(P ) and T (P ) denote
the number of light simple loop, light double edges, and light triple edges in P , and BL , BD
and BT are prescribed parameters specified in (19) in Section 6. PLD* restarts until the final
pairing P is in A0 , and this is the first stage of PLD*. We will show that the distribution of
the output of the first stage of PLD* is almost uniform over A0 . Algorithm PLD is PLD*
accompanied with a rejection scheme so that the output of Stage 1 is uniformly distributed
over A0 . After Stage 1, PLD* enters the second stage, in which it repeatedly removes light
loops, and then light triple edges using random “valid” switchings in Figures 1 and 2, which
were also used in [9]. Valid switchings mean that the performance of the switching does not
destroy more than one multiple edge (or loop), or create new multiple edges (or loops). To
remove the double edges, we will use parameters ρIII (i), specified in Section 6.2. In each
5

u0

u1

u4

u2

u1

u0

u4

u2

u3

u3

Figure 1: Switching for a light loop

u4

v4

u4

v4

u3

v3

u3

v3

u1

v1

u1

v1

u2

v2

u2

v2

Figure 2: Switching for a light triple edge

u2

u1

u3

u4

u2

u1

u3

u4

u2

u1

u3

u2

u1

u3

v2

v1

v3

v4

v2

v1

v3

v4

v2

v1

v3

v2

v1

v3

Figure 3: Type III

Figure 4: Type I

iteration, if the current pairing P contains i double edges, then with probability ρIII (i), it
performs a random valid switching as in Figure 3, and with probability 1−ρIII (i), it performs
a random valid switching as in Figure 4. Note that the switching in Figure 4 is the typical
type (called type I) used to switch away a double edge. The switching in Figure 3 is a new
type (called type III) to remedy the distribution distortion caused by merely performing
type I switchings.
We will prove:
(a) all parameters ρIII (i) can be computed in time o(n);
(b) the number of iterations in the second stage is o(n) in expectation and with high
probability;
(c) starting from a uniformly random P ∈ A0 , the second stage of PLD* outputs an
asymptotically uniform random graph with target degree sequence d.
Even though the type III switching is new, and some other new considerations will be
brought in to cope with various complications, the main ideas involved in the second stage
of PLD* are similar to those in [9]. Indeed, we will tune PLD* into PLD with introduction
of more types of switchings, and with various sorts of rejections. The output of PLD is a
uniformly random graph with the target degree sequence (see Lemma 25). We can show
that the probability that PLD ever performs a rejection or other types of switchings is o(1).
This allows us to couple PLD and PLD* and deduce that the output distributions of PLD*
differs from PLD, which is uniform, by o(1), and thus derive our main theorem. The full
details of PLD and PLD* for the second stage will be presented in Section 6.
Our major innovative ideas lie in the first stage of PLD and PLD*. We will later specify
a set of pairings Φ0 ⊆ Φ (in Section 3.1) such that with a non-vanishing probability a
random pairing in Φ is in Φ0 . Roughly speaking, Φ0 excludes pairings in which some vertex
is incident with too many heavy multiple edges or heavy loops. The algorithm starts by
repeatedly generating a random pairing P ∈ Φ until P ∈ Φ0 . Then, PLD* enters two
phases (both in Stage one) sequentially where heavy multiple edges and then heavy loops
are eliminated using switchings.
6

Since our analysis is tailored for the power-law degree sequence, we start by stating some
properties P
of power-law degree sequences. In the rest of this paper, we will assume that d is
such that i di is even and Definition 1 is satisfied for some fixed K > 0 and γ > 2.5, that
is
X
Nj ≤ Kni1−γ , for all 1 ≤ i ≤ ∆
(2)
j≥i

where Nj = |t : dt = j|. A stronger lower bound on γ will be imposed later. We will
also assume without loss of generality in the proof that γ < 3, since a plib sequence with
parameter γ ≥ 3 is also a plib sequence with parameter γ 0 for every γ 0 < γ.
Recall that ∆ denotes the maximum degree. Without loss of generality we may assume
that ∆ = d1 ≥ · · · ≥ dn ≥ 1. For every 1 ≤ i ≤ n, the number of vertices with degree at
least di is at least i, so (2) implies that i ≤ Kn(di )1−γ . This gives
di ≤ (Kn/i)1/(γ−1) for all 1 ≤ i ≤ n.
(3)
P
Qk−1
For a positive integer k, define Mk = 1≤i≤n [di ]k where [x]k denotes i=0 (x − i). It follows
immediately from (3) using γ > 2 that
∆ ≤ (Kn)1/(γ−1)
M1 = Θ(n)
(4)
k
k/(γ−1)
Mk = O(∆ ) = O(n
) for every fixed integer k ≥ 2.
Let 0 < δ < 1/2 be a constant to be specified later. Define
h = n1−δ(γ−1) ,
H = {i : 1 ≤ i ≤ h}, L = [n] \ H.

(5)
(6)

Immediately from (3),

dh ≤ K 1/(γ−1) nδ ,

(7)

δ
and so every vertex in L, the set of
Plight vertices, has degree O(n ).
Finally, for k ≥ 1 define Hk = v∈H [dv ]k and Lk = Mk − Hk .

3

First stage of PLD and PLD*: heavy multiple edge
reductions

Recall that PLD (PLD*) first generates a uniformly random pairing P0 ∈ Φ0 , where Φ0 ⊆ Φ
is the set of pairings satisfying some initial conditions. Our next step is to specify Φ0 .
Suppose firstly that M2 < M1 . In this case, we define Φ0 to be the set of simple pairings
in Φ. Janson [12] showed that under this condition, there is a probability bounded away
from 0, uniformly for all n, that P0 ∈ Φ0 . So if P0 is not simple, PLD (PLD*) rejects it,
and otherwise, the remaining phases are defined to pass P0 directly to the output. In this
case, PLD (PLD*) is equivalent to the simple rejection scheme. Henceforth in the paper, we
assume
M2 ≥ M1 .
(8)

7

3.1

Specifying Φ0

Given i < j, let Wi,j (P ) denote the number of points in vertex i that belong to heavy loops
or heavy multiple edges in P with one end in i and the other end not in j. Let Wi (P )
denote the number of pairs that are in a heavy non-loop multiple edge, with one end in i.
For any
p i ≤ j, let mi,j = mi,j (P ) denote the number of pairs between i and j in P . Define
η = M22 H1 /M13 and

Φ0 =









P ∈Φ:

mi,j Imi,j ≥2 Wi,j ≤ ηdi
∀1 ≤ i, j ≤ h with i 6= j;
mi,i Wi ≤ ηdi
∀1 ≤ i ≤ h;
P
2
2
m
I
≤
4M
/M
;
2
1
P1≤i<j≤h i,j mi,j ≥2
m
≤
4M
/M
.
i,i
2
1
1≤i≤h






(9)





The following lemma ensures that the probability of an initial rejection is at most 1/4 +
o(1).
Lemma 4. Assume that η = o(1), i.e. M22 H1 /M13 = o(1). If P is a random pairing in Φ
then the probability that P ∈ Φ0 is at least 1/4 + o(1).
Next we define the switchings used in Stage 1.

3.2

Heavy switchings

Given a pairing P , let m(P, i, j) denote the number of pairs between i and j in P , and set
m(P, i) = m(P, i, i), the number of loops at i.
We first define the switchings we are going to use to switch away heavy multiple edges
between a fixed pair of heavy vertices (i, j), i < j. The switching that we are going to define
will change the multiplicity of ij from some m ≥ 1 to 0 (i.e. no pair between i and j). We
will only use this switching directly if m(P, i, j) ≥ 2, and we will also use an inverse version
with a certain probability to change the multiplicity from 0 to 1. The inverse switching is
needed for the present problem because i and j may be quite likely to have an edge between
them.
Definition 5. (Heavy m-way switching at (i, j) on P .)
This is defined for an ordered pair (i, j) of vertices with m(P, i, j) = m ≥ 1. Label
the endpoints of the m pairs between i and j as 2g − 1 and 2g, 1 ≤ g ≤ m, where points
1, 3, . . . , 2m − 1 are contained in vertex i. Pick m distinct light pairs x1 , . . . , xm . Label
the endpoints of xg as 2m + 2g − 1 and 2m + 2g. The switching operation replaces pairs
{2g − 1, 2g} and {2m + 2g − 1, 2m + 2g} by {2g − 1, 2m + 2g − 1} and {2g, 2m + 2g}. This
is a valid switching if no new heavy multiple edges or heavy loops are created. See Figure 5.
If S is a heavy multiple edge switching creating P 0 from P , then the inverse of S is an
inverse heavy multiple edge switching creating P from P 0 .
Definition 6. (Heavy m-way loop switching at i.)
This is defined if m(P, i) = m ≥ 1. Label the endpoints of the m loops at i as 2g − 1
and 2g, 1 ≤ g ≤ m. Pick m distinct light pairs in P , labelling the endpoints of the g-th pair
2m+2g −1 and 2m+2g. The switching replaces pairs {2g −1, 2g} and {2m+2g −1, 2m+2g}

8

2m + 2g − 1

2m + 2g − 1

2m + 2g

2m + 2g

i

2g − 1

i

2g − 1

2g − 1

j

2g

2g

2g − 1

2g

j
2g

2m + 2g − 1

Figure 5:
edge

Switching for heavy multiple

Figure 6:
loop

2m + 2g

2m + 2g − 1

2m + 2g

Switching for heavy multiple

by {2g − 1, 2m + 2g − 1} and {2g, 2m + 2g}. This switching is valid if no new heavy multiple
edges or heavy loops are simultaneously created. See Figure 6.
Note from the above definitions that heavy multiple edge and heavy loop switchings are
permitted to create or destroy light multiple edges or light loops. This was not the case
for the switchings introduced in [8] for enumeration purposes, and is another feature of our
arguments not present in [17] and [9]. A heavy multiple edge switching on the pair (i, j) still
does not affect other heavy multiple edges.

3.3

PLD: Phases 1 and 2

Now we are ready to define the first two phases of PLD. We define PLD* in Section 3.4.
3.3.1

Phase 1: heavy multiple edge reduction

Given a pairing P , we define an h × h symmetric matrix M(P ) as follows. For each
1 ≤ i ≤ j ≤ h, if ij is not a multiple edge or loop in G(P ), then the ij and ji entries
in M(P ) are set as ♣; otherwise, these entries are set as the multiplicity of the edge ij.
Therefore, in this matrix, all the diagonal entries take values in {♣} ∪ N≥1 and all offdiagonal entries take values in {♣} ∪ N≥2 where N≥k denotes the set of integers that are at
least k. We call M(P ) the signature of P . Given such a matrix M, we let C(M) denote
the set of pairings that have M as their signature. Phase 1 of PLD uses the heavy m-way
switchings of Definition 5 to switch away heavy multiple edges. Each switching step converts
a pairing uniformly distributed in C(M) for some M to a pairing uniformly distributed in
C(M0 ) where M0 is obtained from M by setting a symmetric pair of its off-diagonal entries
to ♣. By the end of Phase 1, if no rejection occurs, the resulting pairing has a signature
with all off-diagonal entries being ♣, and hence has no heavy multiple edges.
Recall Φ0 in (9). When Phase 1 starts, P0 is uniformly distributed in Φ0 . Let M0 =
M(P0 ).
Recall that H denotes the set of heavy vertices. Let J be a list of all pairs (i, j) with
i < j and both i, j ∈ H; i.e. J is an enumeration of all pairs of heavy vertices. This defines
an ordering of elements in J : (i0 , j 0 )  (i, j) if (i0 , j 0 ) appears before (i, j) in the list J . For
9

each k ≥ 1, let Mk denote the matrix obtained from M0 by setting the ij and ji entries of
M0 to ♣, for each ij among the first k entries in the list J .
Let k ≥ 1 and suppose (i, j) be the k-th element of J . For any non-negative integer m,
let C(Mk , i, j, m) denote the set of pairings with ij having multiplicity m, and all other edges
between heavy vertices satisfying the signature Mk (ignoring its ij entry). Given m ≥ 1
and a pairing P ∈ C(Mk , i, j, m), let fi,j (P ) denote the number of heavy m-way switchings
at (i, j) where m is the multiplicity of ij in P . Given m ≥ 1 and P 0 ∈ C(Mk , i, j, 0), let
bi,j (P 0 , m) denote the number of inverse heavy m-way switchings applicable to P 0 ; these
convert P 0 to some P ∈ C(Mk , i, j, m).
Note that Wi,j (P ) is determined by the signature of P ; hence we may write Wi,j (M(P ))
for Wi,j (P ). For any m ≥ 1, define
mfi,j (Mk , m) = m!M1m ;

(10)

mfi,j (Mk , m) = m!(M1 − H1 − 2m)m ;

(11)

mbi,j (Mk , m) = [di − Wi,j ]m [dj − Wj,i ]m ;

mbi,j (Mk , m)

(12)
2

= [di − Wi,j ]m [dj − Wj,i ]m − mh [di − Wi,j ]m−1 [dj − Wj,i ]m−1 ,

(13)

where Wi,j = Wi,j (Mk ) and Wj,i = Wj,i (Mk ).
By simple counting arguments (see Lemma 28) we have that for every m ≥ 1 and for
every P ∈ C(Mk , i, j, m) and P 0 ∈ C(Mk , i, j, 0),
mfi,j (Mk , m) ≤ fi,j (P ) ≤ mfi,j (Mk , m), mbi,j (Mk , m) ≤ bi,j (P 0 , m) ≤ mbi,j (Mk , m).

Phase 1 of PLD is defined inductively as follows. For each 1 ≤ k ≤ |J |, let Pk denote
the pairing obtained after the k-th step of Phase 1. To define this step, if m(Pk−1 , i, j) ≤ 1,
then Pk = Pk−1 ; otherwise, let m = m(Pk−1 , i, j) and do the following sub-steps:

(i) Choose a random heavy m-way switching at (i, j) on P . Let P 0 be the pairing created
by S;
(ii) Perform an f-rejection with probability 1 − fi,j (P )/mfi,j (Mk , m); then perform a brejection with probability 1 − mbi,j (Mk , m)/bi,j (P 0 , m);
(iii) If no f or b-rejection took place, set Pk∗ = P 0 ;
(iv) Choose a random inverse heavy 1-way switching S 0 at (i, j) on P 0 ; let P 00 denote the
pairing created by S 0 ;
(v) With probability 1/(1 + mbi,j (Mk , 1)/mfi,j (Mk , 1)), set Pk = Pk∗ = P 0 ; with the remaining probability, perform an f-rejection (w.r.t. S 0 ) with probability 1−bi,j (P 0 )/mbi,j (Mk , 1)
and perform a b-rejection with probability 1 − mfi,j (Mk , 1)/fi,j (P 00 ). If no f- or brejection occurred, set Pk = P 00 .
Lemma 7. The probability of an f -rejection or b-rejection during Phase 1 is o(1), if δ >
(3 − γ)/(γ − 2), δ > 1/(2γ − 3), δ < 1/2 and 1 − δ > 1/(γ − 1).
3.3.2

Phase 2: heavy loop reduction

The heavy loop reduction is similar to the heavy multiple edge reduction but simpler. Let
I be an ordering of elements in H. At step k of Phase 2, for k ≥ 1, the algorithm switches
away all loops at the k-th vertex in I using the heavy loop switching, if there are any. In
10

each step, sub-steps (i) and (ii) are the same as in heavy multiple edge reduction, except that
fij (P ) and bij (P 0 , m) are replaced by fi (P ) and bi (P 0 , m), and mfij (Mk , m) and mbij (Mk , m)
are replaced by mfi (Mk , m) and mbi (Mk , m), defined as below:
(14)
mfi (Mk , m) = 2m m!M1m ;
f
m
m
(15)
mi (Mk , m) = 2 m!(M1 − H1 − 2m) ;
b
mi (Mk , m) = [di − Wi ]2m ;
(16)
mbi (Mk , m) = [di − Wi ]2m − mh2 [di − Wi ]2m−2 .
(17)

Moreover, Pk is set to P 0 as long as no f or b-rejection happend. There are no sub-steps
(iii)–(v).
It is easy to show that for any P ∈ C(Mk , i, i, m), mfi (Mk , m) ≤ fi (P ) ≤ mfi (Mk , m);
and for any P ∈ C(Mk ), mbi (Mk , m) ≤ bi (P, m) ≤ mbi (Mk , m).
Lemma 8. If 2/(γ − 1) < 1 + δ(γ − 2) and δ(γ − 1) > 2/3, then the probability that a
rejection happens during Phase 2 is o(1).
3.3.3

Uniformity

Recall
Φ2 = {P ⊆ Φ0 : G[H] (P ) is simple}.

We have the following lemma confirming the uniformity of the output of the first two phases.
Lemma 9. Let P0 be a random pairing in Φ0 and P 0 be the output after the first two phases,
assuming that no rejection occurred. Then P 0 has the uniform distribution on Φ2 .

3.4

Stage 1 of PLD*

The first stage of PLD* is just PLD without rejection. I.e. Phase 1 of PLD* consists of
substeps (i), (iv) and (v) without f-rejection or b-rejection. Phase 2 of PLD* repeatedly
switches away heavy loops using random heavy loop switchings.

4

Proof of Theorem 3

Uniformity
By Lemma 9, PLD generates a uniformly random pairing P ∈ Φ2 after the first stage. In
Section 6 we show that the second stage of PLD outputs a uniformly random graphs with a
plib degree sequence with parameter γ; see Lemma 25.
Running time
The initial generation of P ∈ Φ takes O(n) time, as there are only O(n) points in the
pairing model.
If (γ, δ) satisfies all conditions in Lemmas 4, 7 and 8 then PLD only
√ restarts O(1) times
in expectation in Stage 1. In Section 6.4 we show that if γ > 21/10 + 61/10 ≈ 2.881024968
then there exists δ satisfying
n
1
1
2 − 2/(γ − 1) − (2γ − 3)/(γ − 1)2 o
< δ < min 1 −
,
,
(18)
2γ − 3
γ−1
3−γ
11

and that with any such choice of δ, the probability of any rejection occurring in the second
stage is o(1). It only remains to bound the time complexity of PLD assuming no rejection
occurs.
We first bound the running time in the first two phases. Each step of the algorithm
involves computing fi,j (P ) and bi,j (P, m) in Phase 1, and fi (P ) and bi (P, m) in Phase 2.
In the next lemma we bound the computation time for such functions. The proof is in the
Appendix.
Lemma 10. (a) The running time for computing fi,j (P ) and fi (P ) in Phases 1 and 2,
given P , is O(H1 + |H|m(P, i, j)) and O(H1 + |H|m(P, i, i)) respectively.
(b) The running time for computing bi,j (P, m) and bi (P, m) in Phases 1 and 2, given P ,
is O(∆ + m).
Next, we bound the time complexity of Phases 1 and 2.
Lemma 11. The expected running time for the first two phases is o(n).
Proof. By the definition of Φ0 in (9),
X
E(mi,j Imi,j ≥2 ) = O(M22 /M12 )
i<j≤h

for every P ∈ Φ0 . Let P0 ∈ Φ0 be the pairing obtained before entering Phase 1. Recall that
M(P0 ) is the signature of P0 and let mi,j denote the ij entry of M(P0 ) for (i, j) ∈ J . Note
that the heavy multiple edge switchings are applied only for (i, j) ∈ J such that mi,j ≥ 2. By
Lemma 10, for each such (i, j), the running time is O(mi,j ) for sub-step (i), O(H1 + |H|mi,j )
for (ii), O(1) for (iii) and (iv), and O(H1 + |H|) for (v). Hence, the total time required for
switching away each heavy multiple edge is O(H1 + |H|mi,j ) and hence, the time complexity
for Phase 1 is
X
X
O((H1 + |H|mi,j )Imi,j ≥2 ) = O(H1 )
E(mi,j Imi,j ≥2 ) = O(H1 M22 /M12 ),
1≤i<j≤h

1≤i<j≤h

as |H| = O(H1 ). The above is bounded by o(n) for any γ in the range of Theorem 3 and
any choice of δ satisfying (18). Similarly, the time complexity for Phase 2 is o(n).
The analysis for the running time of Phases 3–5 is similar to that of [9]. In Section 6.2 we
show that parameters ρτ (i) such as ρIII (i) can be computed inductively in at most BD = o(n)
steps. Each switching step involves counting the number of certain structures (involving a
bounded number of edges) in a pairing (e.g. the number of 6-cycles). Using direct brute-force
counting methods yields the time complexity claimed in Theorem 3. The details are given
in the Appendix under the heading “Running time for Phases 3–5”.

5

Proof of Theorem 2

The approximate sampler PLD* is simply PLD without rejections, except for rejections when
checking P ∈ Φ0 and P ∈ A0 . As a result, there is no need to compute the number of certain
structures, as done in PLD.
Proof of Theorem 2. We have already discussed above that it takes PLD, and thus PLD* as
well, O(n) time to find P ∈ A0 before entering the last three phases. The computation of
ρIII (i) takes at most BD = o(n) units of time. It is easy then to see then that each iteration
12

of PLD* in these three phases takes only O(1) time, as it only involves switching several
pairs, and it is shown in Section 6 that the total number of iterations in these three phases
is O(BL + BD + BT ) = o(n). Hence, PLD* is a linear time algorithm in expectation. It only
remains to show that the output of PLD* is close to uniform.
Clearly PLD and PLD* can be coupled so that they produce the same output if no
rejection occurs, and no types other than τ = I, III are chosen in PLD in the last three
phases. By Lemma 25, any rejection occurs in Phase 3–5 with probability o(1). We also prove
(see (40) and Lemma 15) that the probability of ever performing a switching of type other
than I and III is o(1) (the probability is O(ξBD ) which is guaranteed o(1) by Lemma 17).
Thus, the probability that PLD and PLD* have different output is o(1). Let A be an
arbitrary subset of Gd . The probability that PLD outputs a graph in A whereas PLD* does
not is o(1), as that happens only if they have different outputs. Hence,
PP LD (A) − PP LD∗ (A) = o(1) for all A ⊆ Gd .
That is, the distribution of the output of PLD* differs by o(1) in total variation distance
from the distribution of the output of PLD, which is uniform.

6

Second stage: phases for light switchings

The second stage of PLD* has been described already in Section 2. In this section, we tune
it into PLD with a rejection
scheme, as well
P as introducing more types of switchings.
P
Recall that Mk = 1≤i≤n [di ]k , Hk = i∈H [di ]k and Lk = Mk − Hk . Define
4L2
2L3 M3
4L2 M2
BL =
BT =
.
(19)
; BD =
2
M1
M1
M13
and recall from (1) that


L(P ) ≤ BL ; D(P ) ≤ BD ; T (P ) ≤ BT ;
A0 = P ∈ Φ2 :
.
G(P ) contains no other types of multiple edges
Roughly speaking, the only possible multiple edges in a pairing in A0 are loops, double edges
and triple edges, and there are not too many of them. Lemma 9 guarantees that the resulting
pairing P after the first two phases is uniformly distributed in Φ2 , and the following lemma
guarantees that PLD restarts O(1) times in expectation until entering the second stage, i.e.
the last three phases. The proof of the lemma is presented in Appendix.
Lemma 12. Assume L4 = o(M12 ) and L4 M4 = o(M14 ). With probability at least 1/2, a
random P ∈ Φ2 is in A0 .
The second stage consists of Phases 3, 4 and 5 during which light multiple edges (loops,
triple edges and double edges) are eliminated in turn. These phases are similar to those in [9]
for regular graphs, except that here we use more types of switchings to boost several types
of structures due to the non-regular degree sequence. In each phase, the number of a certain
type of light multiple edges will be reduced until there are no multiple edges of that type.
There will be f and b-rejections designed so that after each phase, the output is uniform
conditioned on the event that the number of multiple edges of a certain type is zero. In each
1
phase, the set of pairings is partitioned into classes ∪ii=0
Si , where Si is the set of pairings
with exactly i multiple edges of the type that is treated in the phase; i1 ∈ {BL , BD , BT } is a
pre-determined integer whose validity is guaranteed by the definition of A0 . Let Ai be the
13

output after each phase for i = 3, 4, 5. In Phase 3 (and 4), the number of loops (triple edges)
will be reduced by one in each step. Hence, A4 will be the set of pairings in A0 containing
no loops or triple edges. The analysis for Phases 3 and 4 is quite simple, almost the same
as that in [17] and [9]. We postpone this until Sections 6.3. In Phase 5 the number of light
double edges is reduced. For this, we will use the boosting technique introduced in [9], but
we use more booster switchings to cope with non-regular degrees. We briefly describe the
framework of Phase 5 here, pointing out the major differences compared with [9], with the
detailed design, definitions and analysis given in the rest of this section. Phase 5 consists
of a sequence of switching steps, where each step typically eliminates one double edge using
a type I switching operation. This is the same operation used in [17]. However, in some
occasional switching steps, some other types of switchings are performed. The algorithm first
determines certain parameters ρτ (i), 0 ≤ i ≤ BD . In each switching step, PLD chooses a
switching type τ , e.g. τ = I, from a given set of switching types, with probability ρτ (i), where
i is the number of double edges in the current pairing. After choosing type τ , PLD performs
a random switching of type τ . Rejections can happen in each switching step with small
probabilities, and these probabilities are carefully designed to maintain the property that
the expected number of times a pairing is reached, throughout the course of the algorithm,
depends only on the number of double edges the pairing has. (This feature was introduced
in [9] and permits the use of switchings that occasionally increase the number of double
edges, and is a marked departure from the method in [17].) As a consequence, the output
after Phase 5 has uniform distribution. It follows that the output of PLD is uniform over
simple graphs with degree sequence d.
The parameters ρτ (i) are designed to guarantee the uniform distribution of the output
of PLD, and meanwhile to have a small probability of any rejection during Phase 5. They
are determined as the solution of a system of equations for which we give an efficient (o(n)
time) computation scheme. The analysis is a little different from [9] due to the different
types of switchings being used. Switchings in [9] are classified into different classes (class
A and class B), whereas the switchings in the present paper can all be classified as class A
under that paradigm. However, the analysis of b-rejections is more complicated than in [9],
and as a convenience for the argument we need to introduce “pre-states” of a pairing, and
“pre-b-rejections”. These concepts will be discussed in detail in Section 6.1.1.
As in [9], under various conditions on Mk , Hk and Lk for several small values of k, we
bound the probability of any rejection during Phases 3–5 by o(1) (see Lemmas 14, 17, 18, 19,
20, 22, and 24).
√ These conditions set constraints on (γ, δ). Eventually we will show that if
γ > 21/10 + 61/10 ≈ 2.881024968 then there exists δ > 0 satisfying all conditions required
in these lemmas; see Lemma 25.
As mentioned before, we postpone the discussions of Phases 3 and 4 till Section 6.3 and
we start our discussions of Phase 5 by assuming that the algorithm enters Phase 5 with a
uniform random pairing P ∈ A4 , where A4 is the set of pairings in A0 containing no loops
or triple edges. By the definition of A0 , the number of double edges in P is at most i1 = BD ,
and they are all light. In Phase 5, Si denotes the set of pairings in A4 containing exactly i
double edges. Then, S0 is the set of pairings in Φ that produce simple graphs.

14

6.1

Phase 5: double edge reduction

We will use five different types of switchings in the phase for the double edge reduction. We
introduce one of them now, and the rest after some discussion. The most commonly used
type is I, which was shown in less detail in Figure 4. We formally define it as follows.
Type I switching. Take a double edge z; note that one end of z must be a light vertex. Label
the end points of the two pairs in z by {1, 2} and {3, 4} so that points 1 and 3 are in the same
vertex and this vertex is light. Choose another two pairs x and y and label their ends by
{5, 6} and {7, 8} respectively. Label the vertices that containing points 1 ≤ i ≤ 8 as shown
in Figure 7 (same as Figure 4 but with points and vertices labelled as in the definition). If
(a) ui , vi , 1 ≤ i ≤ 3 are six distinct vertices; (b) both x and y are single edges; (c) there are
no edges between u1 and ui for i ∈ {2, 3} and no edges between v1 and vi for i ∈ {2, 3}, then
replace the four pairs {2i − 1, 2i}, 1 ≤ i ≤ 4 by {1, 5}, {3, 7}, {2, 6} and {4, 8} as shown in
Figure 7.
5

1

u2

6

v2

3

7

u1

2

v1

4

5

u3

8

1

u2

v3

6

v2

3

7

u1

2

v1

4

u3

8

v3

Figure 7: type I
The inverse of the above operation is called an inverse type I switching. Hence, to perform
an inverse switching, we will choose a light 2-star and label the points in the two pairs by
{5, 1} and {3, 7} so that 1 and 3 are in the same light vertex; choose another 2-star (not
necessarily light) and label the points in the two pairs by {6, 2} and {4, 8} so that 2 and 4
are in the same vertex. Label all involved vertices as shown in Figure 7. If (a) all the six
vertices are distinct; (b) both 2-stars chosen are simple; (c) there is no edge between ui and
vi for all 1 ≤ i ≤ 3, then replace all four pairs {1, 5}, {3, 7}, {2, 6} and {4, 8} by {2i − 1, 2i},
1 ≤ i ≤ 4.
Note that from the definition of the type I switching, there can be multiple switchings of
type I that switch a pairing P to another pairing P 0 , caused by different ways of labelling
the points in the switched pairs.
When a type I switching is performed, the number of double edges will decrease by one.
Given a pairing P , let F1 (P ) denote the set of mappings from the numbers 5, 1, 3, 7 to the
points of P such that {5, 1} and {3, 7} are pairs of P and 1 and 3 are in the same vertex,
which is light. We call ({5, 1}, {3, 7}) a light ordered 2-star. Define F2 (P ) similarly from
the numbers 6, 2, 4, 8 without the restriction to light vertices. Let F(P ) = F1 (P ) × F2 (P ).
The elements of F(P ) are called doublets. As |F1 (P )| = L2 and |F2 (P )| = M2 , we have
|F(P )| = M2 L2 , which is independent of P . If a type I switching switches a pairing P 0 to
a new pairing P , then it creates a new element in F(P ) (but not all elements in F(P ) are
created this way).
Note the labels on the points in each type I switching induce, in an obvious way, a unique
doublet in the resulting pairing. Now, for a given pairing P , define Z0 (P ) to be the set of
doubles in F(P ) that can be created in this way. Then, the number of ways that P can
15

be created via a type I switching is |Z0 (P )|. Since Z0 (P ) ⊆ F(P ), we immediately have
|Z0 (P )| ≤ M2 L2 . For a power-law distribution-bounded degree sequence with γ < 3, |Z0 (P )|
varies a lot among pairings P ∈ Si . As we discussed before, a big variation of this number
will cause a big probability of a b-rejection. It turns out that five main structures cause
the undesired variation of this number: Hi , 1 ≤ i ≤ 5 in Figure 8. In order to reduce the
probability of a b-rejection, we boost pairings P such that G(P ) does not contain sufficient
number of copies of H0 , or in other words, G(P ) contains many copies of Hi for some
1 ≤ i ≤ 5. These force the use of switchings of type III, IV, V, VI and VII respectively (we
leave out type II as there is a type of switching called type II used for the d-regular case which
we do not use in this paper; but it is possible that using this type, together with strengthening
other parts can lead to a weaker constraint on γ). These switchings will be defined in
Section 6.1.3 and are shown in Figures 3 and 9 – 12. Let Λ = {I, III, IV, V, V I, V II} be
the set of switching types to be used. All types of switchings, except for type I, will switch
a pairing to another pairing containing the same number of double edges. In other words, a
type τ (τ 6= I) switching does not change the number of double edges in a pairing.
The basic method in this part of the algorithm and analysis is similar to that in [9] but
there are new ingredients, which include the use of pre-states. We briefly describe how and
why they are used.
Given a pairing P , let Z(P ) be the subset of F(P ) that includes just those doublets
induced (in the way described above) by applying a switching of any type to some other
pairing to create P . Then, each switching (of any type) converting some P 0 to P induces a
doublet z ∈ Z(P ). However, except for type I, other type of switchings create extra pairs
along with z. For instance, a type III switching creates one extra pair whereas a type IV
switching creates four extra pairs (see Figures 3 and 9). These extra pairs modify (slightly)
the expected number of times that a given z ∈ Z(P ) is created. In order to account for this,
we will associate with each pairing P a set of pre-states. These pre-states are basically a copy
of P plus a designated z ∈ Z(P ), together with several pairs, depending on how P is reached
via a switching. Then, going from a pre-state to P , we may perform a pre-b-rejection. This
rejection trims off the fluctuation caused by the extra pairs created together with z, ensuring
that every z is created with the same expected number of times.
For every pairing P ∈ Si , we want to make sure that the expected number of times that
P is created in the phase is independent of P ; we do so by equalising the expected number
of times that each z ∈ Z(P ) is created, after the pre-b-rejection; denote this quantity by
qi (i.e. this quantity is independent of P ). Then the expected number of times that P is
created, before the b-rejection, is qi |Z(P )|. Then, the probability of a b-rejection depends
only on the variation of |Z(P )| over all P ∈ Si . The reason to use type III,IV,V,VI and VII
aside from type I is to ensure that the variation of |Z(P )| is sufficiently small in a class Si .
Let Z ∗ (P ) = F(P ) \ Z(P ). We will prove that |Z ∗ (P )| is sufficiently small compared with
|Z(P )|.
Next, we partition Z(P ) to ∪5i=0 Z(P ) according to the different restrictions put on to
z ∈ Z(P ) when different types of switchings are applied. Let Zi (P ) (1 ≤ i ≤ 5) denote the
set of z ∈ F(P ), as shown in Figure 7, satisfying
(a) all these four pairs {5,1},{3,7},{6,2},{4,8} are contained only in single edges; v1 is a
light vertex; and all six vertices ui and vi , 1 ≤ i ≤ 3, are distinct;
16

(b) for each 1 ≤ i ≤ 5, z corresponds to a copy of Hi in Figure 8 in G(P ).

H0

H1

H2

H3

H4

H5

Figure 8: A pair of 2-paths
Let Zi (P ) = |Zi (P )| and Z ∗ (P ) = |Z ∗ (P )|. Let Sτ+ (P ) denote the number of type τ
switchings that converts a pairing to P . Then, SI+ (P ) = Z0 (P ). Clearly,
5
X
Zi (P ) + Z ∗ (P ) = |F(P )| = M2 L2 .
(20)
i=0

P
We will show later that the variation of 5i=0 Zi (P ) caused by Z ∗ (P ) is sufficiently small for
the range of γ we consider.
Note that each type τ ∈ Λ corresponds to a unique integer t so that the created element
z ∈ Z(P ), formed after a type τ switching creating a pairing P , is in Zt (P ). We define
t = t(τ ) to be this unique integer.
For instance, t(I) = 0, and by looking at Figures 3, 9 and 8, we have t(III) = 1 and
t(IV ) = 2.
Before defining all types of switchings, we first give a general description of the whole
approach.
6.1.1

A brief outline

Note that a pairing can be reached more than once in Phase 5. Let σ(P ) denote the expected
number of times that pairing P is reached. Our goal is to design the algorithm in such a
way that σ(P ) = σ(i) for every P ∈ Si and for every 0 ≤ i ≤ i1 ; this ensures that the output
of Phase 5, if no rejection happens, is uniform over S0 .
Choosing a switching type.
As in [9], in Phase 5, the algorithm starts with a uniform distribution on A4 . In each
step, the algorithm starts from a pairing P ∈ Si obtained from the previous step, and chooses
a switching type according to a pre-determined “type distribution”. This distribution only
depends on i, i.e. the class where P is in, and is independent of P . We will use ρτ (i)
to
P denote the probability that type τ switching is chosen for a pairing in Si . We require
τ ρτ (i) ≤ 1 for every i. If this inequality is strict, then with the remain probability, the
algorithm terminates. Such a rejection is called a t-rejection.
An f-rejection.
If a type τ is chosen, then the algorithm chooses a random type τ switching S that can be
performed on P ∈ Si . An f-rejection may happen in order to equalise the expected number
of times a particular type τ switching is performed, for every P ∈ Si . Let fτ (P ) denote the
number of type τ switchings that can be performed on P , then, we will perform an f-rejection

17

with probability
fτ (P )
,
(21)
mτ (i)
where mτ (i) will be specified later as an upper bound for maxP ∈Si fτ (P ). Therefore, every
switching S is chosen and is not f-rejected with probability ρτ (i)/mτ (i). Hence, the expected
number of times that S is performed during Phase 5 is
ρτ (i)
,
σ(i)
mτ (i)
assuming that each pairing in Si is reached with σ(i) expected number of times.
1−

Pre-states and pre-b-rejection.
We first give a formal definition of pre-states of a pairing P 0 ∈ Sj . Given pairing P 0 ∈ Sj ,
a switching type τ and an arbitrary z ∈ Zt(τ ) (P 0 ), a pre-state of P 0 , corresponding to τ and
z is a copy of P 0 in which z is designated together with a set of pairs X such that an inverse
bτ (P, z) denote the set of pre-states of
type τ switching can be performed on X and z. Let B
bτ (P, z)|. Define
P corresponding to switching type τ and z ∈ Zt(τ ) (P ) and let bbτ (P, z) = |B
bbτ (P ) = minz∈Z bbτ (P, z). We will define later m
b τ (i) as a lower bound for minP ∈Si bbτ (P ).
t(τ )
In particular, we will show that for every 0 ≤ i ≤ i1 − 1, bbI (P, z) = 1 for any P ∈ Si and for
any z ∈ Zt(I) (P ) = Z0 (P ) (see (22)), which immediately implies that
m
b I (i) = 1, for all 0 ≤ i ≤ i1 − 1.
(22)
Now, let P 0 be the pairing that S converts P into and let z be the new pair of 2-stars
created by S. Then, z ∈ Zt(τ ) (P 0 ). Assume P 0 ∈ Sj . Note that S designates a pre-state of
P 0 . We will perform a pre-b-rejection before entering P 0 , to ensure that every z ∈ Zt(τ ) (P 0 )
are created with equal expected number of times, for every P 0 ∈ Sj .
A pre-b-rejection will be performed in the following way. Note that our algorithm has
been designed so that given τ , P ∈ Si and z ∈ Zt(τ ) (P ), all pre-states of P corresponding to
τ and z will be reached with equal expected number of times, which is
ρτ (i)
,
yi := σ(i)
mτ (i)
which is independent of P , P 0 and z. Then, a pre-b-rejection will be performed with probability
m
b τ (j)
1−
.
bbτ (P 0 , z)
Now, for any z ∈ z ∈ Zt(τ ) (P 0 ), the expected number of times that z is created is
X
m
b (j)
ρτ (i)
yi τ
= yi m
b τ (j) = σ(i)
m
b (j),
bbτ (P 0 , z)
mτ (i) τ
0
b
s∈B (P ,z)

(23)

τ

which is independent of z. This ensures that every z ∈ Zt(τ ) (P 0 ) are created with equal
expected number of times, as desired. In fact, we will equalise this quantity for every
switching type τ as well. This will be done by choosing ρτ (i) properly. We will later let q(j)
denote this equalised quantity, which depends only on j, i.e. the class P 0 is in.
A b-rejection
18

If S is not f-rejected or pre-b-rejected, then a b-rejection will be performed with probability
m(j)
,
1−
b(P 0 )
where b(P 0 ) = |Z(P 0 )| and m(j) will be a specified lower bound for minP 0 ∈Sj |Z(P 0 )|. If no
b-rejection occurred, then P is switched to P 0 and that completes a switching step. The
expected number of times that P 0 is reached via a switching is
q(j)m(j),
(24)
as the expected number of times a z ∈ Z(P 0 ) is created is equalised to qj , as we mentioned.
This confirms that the expected number of times that P 0 is reached is dependent only on
the state in which P 0 lies.
The Markov chain.
This defines a Markov chain run on A4 . Each state in the Markov chain is a pairing in A4
together with {O, R}; O is the state representing an output and R is the state representing
a rejection. A switching step from P to P 0 is a transition from P to P 0 in the chain. If a
rejection (of any kind) happens in a step, then the chain enters state R and the algorithm
terminates without an output. If P ∈ S0 , then we interpret a type I switching on P as
outputting P . So if the chain reaches a pairing in S0 , then it enters O in the next step with
probability ρI (0).
III, IV, V, V I, V II
I
Si+1

I
Si

Si−1

Recall the following quantities that we are going to specify later:
mτ (i)
an upper bound for max fτ (P );
P ∈Si

m
b τ (i)
m(i)

a lower bound for min bbτ (P ) = min min bbτ (P, z);
P ∈Si

P ∈Si z∈Zi(τ )

a lower bound for min b(P ).
P ∈Si

In the next section, we will deduce a system involving ρτ (i) and σ(i) using the above quantities, and use the solution of that system to specify ρτ (i), the type distribution function to
be used in the algorithm.
6.1.2

The system with ρτ (i) and σ(i)

Note that Pt is the pairing obtained after the t-th step of Phase 5. Initially, each pairing in
A4 can be P0 with probability 1/|A4 |.
Note that while performing switchings, a pairing can be visited more than once. Recall
that σ(P ) denotes the expected number of times that a given pairing P is visited. We wish
to ensure that σ(P ) = σ(i) for every P ∈ Si .
Given a pairing P ∈ Si , we want to equalise the expected number of times that any pair
of 2-stars in Z(P ) is created, after the pre-b-rejection. Recall that q(i) denotes this quantity
(see the discussion below (23)). We first consider switchings of type I. For such a switching
19

S converting P 0 to P , P 0 must be in Si+1 . By (23) (with i in (23) being i + 1 and j in (23)
being i and τ = I) and (22),
ρI (i + 1)
ρI (i + 1)
m
b I (i) = σ(i + 1)
.
(25)
q(i) = σ(i + 1)
mI (i + 1)
mI (i + 1)
For any switching type τ 6= I, a switching of type τ will switch a pairing to another pairing
in the same class. So, for each τ ∈ Λ \ {I}, let S be a switching of type τ that converts a
pairing P 0 ∈ Si to P . By (23),
ρτ (i)
m
b (i).
q(i) = σ(i)
mτ (i) τ
Combining this with (25), we must have
q(i)mτ (i)
σ(i + 1) ρI (i + 1) mτ (i)
ρτ (i) =
=
,
σ(i)m
b τ (i)
σ(i) mI (i + 1) m
b τ (i)
for each τ ∈ Λ \ {I}.
Finally, by (24) and (25),
ρI (i + 1)
+ 1/|A4 |,
σ(i) = m(i)q(i) + 1/|A4 | = σ(i + 1)
mI (i + 1)
where 1/|A4 | is the probability that P0 = P . The boundary condition is q(i1 ) = 0.
Letting xi = σ(i)|A4 |, we have
xi1 = 1;
(26)
m(i)
+ 1 ∀0 ≤ i ≤ i1 − 1;
(27)
xi = xi+1 ρI (i + 1)
mI (i + 1)
xi+1
mτ (i)
ρτ (i) = ρI (i + 1)
∀τ ∈ Λ \ {I}, 1 ≤ i ≤ i1 − 1;
(28)
xi mI (i + 1)m
b τ (i)
X
ρτ (i) ≤ 1.
(29)
ρτ (i) ≥ 0; ρτ (i1 ) = 0 (τ 6= I);
τ

It follows immediately that
1
xi+1
mI (i + 1)
≤
.
(30)
xi
ρI (i + 1) m(i)
We will show thatPthere is a solution ρ∗τ (i) and x∗i to the above system with ρ∗I (i) close
to 1 for every i and τ ∈Λ ρ∗τ (i) = 1 − o(i1 ). Then we will set the type probability ρτ (i) in
the algorithm to be ρ∗τ (i) for every τ and i. A lemma in [9] shows that with this setting, the
expected number of times that a pairing in Si is reached is x∗i /|A4 |, independent of P ∈ Si .
Therefore, for each P ∈ S0 , the probability that P is outputted is equal to (x∗0 /|A4 |)ρI (0),
and thus
the output of PLD is uniformly distributed in S0 .
(31)

Before proving the existence of such a solution, we need to define all types of switchings and
b τ (i) that have appeared in the system.
then specify quantities mτ (i), m(i) and m

6.1.3

Switchings

In the definition of a switching (of each type), we will choose a set of pairs and label their end
points in a certain way; then we replace this set of pairs by another set, without changing
20

the given degree sequence. The type I switching has been formally defined in detail. The
definition of other types is clear by the illustrations in Figures 3, and 9 – 12. All pairs
involved in the switchings and the inverse switchings are required to be contained only in
single edges.
A type I switching creates only a pair of 2-stars in Z0 (P ) for some P and no other extra
pairs. Hence, given any P ∈ Si (i ≤ i1 − 1) and any z ∈ Z0 (P ), there is only one pre-state
of P corresponding to type I and z. In other words, for every 0 ≤ i ≤ i1 − 1, bbI (P, z) = 1 for
any P ∈ Si and for any z ∈ Zt(I) (P ) = Z0 (P ). This confirms (22), i.e.
m
b I (i) = 1, for all 0 ≤ i ≤ i1 − 1.

u2

u1

u3

u2

u1

u3

v2

v1

v3

v2

v1

v3

Figure 9: type IV

u2

u1

u3

u2

u1

u3

v2

v1

v3

v2

v1

v3

Figure 10: type V

u2

u1

u3

u2

u1

u3

v2

v1

v3

v2

v1

v3

Figure 11: type VI

6.1.4

Specifying mτ (i), m(i) and m
b τ (i)

Define
mI (i) = 4iM12 ; mIII (i) = M3 L3 ; mIV (i) = 2M23 L2 ;
mV (i) = 2M22 M3 L3 ; mV I (i) = M25 L2 ; mV II (i) = M24 M3 L3 .

21

(32)
(33)

u2

u1

u3

u2

u1

u3

v2

v1

v3

v2

v1

v3

Figure 12: type VII
It is easy to verify that for every P ∈ Si and τ ∈ Λ, fτ (P ) ≤ mτ (i) for each 0 ≤ i ≤ i1 . The
factor 2 in mIV (i) and mV (i) accounts for the two cases that (a) u2 is adjacent to v2 ; (b) u3
is adjacent to v3 .
Define
2
2
U1 = (Kn)1/(γ−1) (Kn)(γ−2)/(γ−1) = (Kn)(2γ−3)/(γ−1) ;
(34)
γK k/(γ−1) 1−δ(γ−k−1)
Uk =
n
for k ≥ 2.
(35)
k+1−γ
The Uk (k ≥ 1) are upper bounds for the number of various types of structures in P ∈ A0 .
For instance, U1 is an upper bound for the number of 2-paths with one end at a given vertex
(see the Appendix for details). These Uk will be used to define parameters such as m(i) and
to bound the probabilities of f-rejections, pre-b-rejections and b-rejections in Phase 5.
b τ (i). Define
Next we specify m(i) and m
m(i) = M2 L2 − 8i(dh M2 + d1 L2 ) − (2id21 d2h + 4iU12 + 8M2 U1 + L4 ).
(36)
and

m
b I (i) = 1; m
b III (i) = M1 − 2U1 − 4i; m
b IV (i) = (M1 − 2U1 − 4i)3 ;
m
b V (i) = (M1 − 2U1 − 4i)4 ; m
b V I (i) = (M1 − 2U1 − 4i)6 ;
m
b V II (i) = (M1 − 2U1 − 4i)7 .

(37)
(38)
(39)

Simple counting arguments yield the following lemma, whose proof is presented in the
Appendix.
Lemma 13. (a) For every 0 ≤ i ≤ i1 , P ∈ Si , b(P ) ≥ m(i).
(b) For every τ ∈ Λ, 0 ≤ i ≤ i1 , P ∈ Si and z ∈ Zt(τ ) , bbτ (P, z) ≥ m
b τ (i).

6.2

Specifying ρτ (i)

Let ξ = ξn > 0 be a function of n to be specified later. Consider the system (26)–(29). We will
showP
that there is a unique solution (x∗i , ρ∗τ (i)) to the system such that ρ∗I (i) + ρ∗III (i) = 1 − ξ
and τ 6=I,III ρ∗τ (i) < ξ for all i.
We first specify a recursive algorithm to compute ρ∗I (i) and ρ∗III (i) and x∗i satisfying (26)–
(28) so that
ρ∗I (i) + ρ∗III (i) + ξ = 1, for all i.
Base case: ρ∗I (i1 ) = 1 − ξ and ρ∗III (i1 ) = 0. From this we can compute x∗i1 .
Inductive step: assume that we have computed ρ∗I (i + 1), ρ∗III (i + 1) and x∗i+1 . For i: we

22

can compute ρ∗I (i), ρ∗III (i) and x∗i by solving the following system.
m(i)
+ 1;
x∗i = x∗i+1 ρ∗I (i + 1)
mI (i + 1)
x∗i+1
mIII (i)
∗
∗
ρIII (i) = ρI (i + 1) ∗
;
xi mI (i + 1)m
b III (i)
ρ∗I (i) + ρ∗III (i) = 1 − ξ.

Putting

C1 = x∗i+1 ρ∗I (i + 1)

mIII (i)
m(i)
+ 1; C2 = ρ∗I (i + 1)x∗i+1
,
mI (i + 1)
mI (i + 1)m
b III (i)

we have
x∗i = C1 ;
ρ∗I (i) = 1 − ξ −

C2
;
C1

C2
.
C1
Hence, all ρ∗I (i), ρ∗III (i) and x∗i can be uniquely determined, once we fix ξ. We only have to
verify that the solution satisfies (29) by choosing ξ appropriately. Define
32M22
ξ=
.
(40)
M13
Lemma 14. Assume δ + 1/(γ − 1) < 1 and ξ = o(1). Then,
X
mτ (i)
max
≤ ξ.
0≤i≤i1
m
b τ (i)m(i)
ρIII (i) =

τ ∈Λ\{I,III}

Proof. Recall (32) – (39) and it is easy to see that m(i) ≥ M2 L2 /2 and each m
b τ (i) is at
least a half of the first term in their expressions. Hence, for every 0 ≤ i ≤ i1 ,


X
2
mτ (i)
4M23 L2 4M22 M3 L3 2M25 L2 2M24 M3 L3
≤
+
+
+
.
m
b τ (i)m(i)
M2 L2
M13
M14
M16
M17
τ ∈Λ\{I,III}

Using the inequality that L3 ≤ L2 dh = O(L2 nδ ), and the fact that M2 = O(n2/(γ−1) ) by (4),
and our assumptions that δ + 1/(γ − 1) < 1 and ξ = o(1), the last three terms in the above
formula are dominated by the first term. Hence,
X
mτ (i)
32 M23 L2
≤
= ξ.
m
b τ (i)m(i)
M2 L2 M13
τ ∈Λ\{I,III}

This holds uniformly for all 0 ≤ i ≤ i1 and this completes the proof of this lemma.
Then we have the following lemma.
Lemma 15. Assume ξ = o(1) and M3 L3 /M2 L2 M1 = o(1). System (26)– (29) has unique
solution (x∗i , ρ∗τ (i)) with ρ∗I (i) + ρ∗III (i) = 1 − ξ for every i. Moreover,




M3 L3
M3 L3
∗
∗
ρI (i) = 1 − O
; ρIII (i) = O
;
M2 L2 M1
M2 L2 M1


mτ (i)
∗
ρτ (i) = O
for all τ ∈ Λ \ {I, III}.
m(i)m
b τ (i)
23

Proof.
We have proven the uniqueness. We only need to verify that the solution given
above satisfies the inequalities in (29). It is obvious that both xi > 0 for each 0 ≤ i ≤ i1 .
This implies that ρIII (i) > 0. Next, we bound ρIII (i) and ρτ (i) for τ ∈ Λ \ {I, III}. Since
C1 ≥ x∗i+1 ρ∗I (i + 1)m(i)/mI (i + 1) we have
mIII (i)
C2
≤
.
ρIII (i) =
C1
m
b III (i)m(i)
By (32), (36) and (37),


M3 L3
ρIII (i) = O
= o(1).
M2 L2 M1
This verifies that 0 < ρIII (i) < 1 and 0 < ρI (i) < 1 for every i.
Moreover, by (28) and (30), for each τ ∈ Λ \ {I, III},
x∗ ρ∗ (i + 1)
mτ (i)
mτ (i)
ρτ (i) ≤ i+1 I ∗
≤
.
xi
mI (i + 1)m
b τ (i)
m
b τ (i)m(i)
Therefore,
X
X
ρτ (i) = 1 − ξ +
ρτ (i) ≤ 1 − ξ + ξ = 1,
τ ∈Λ

τ ∈Λ\{I,III}

by (40) and so (29) is satisfied. Hence, (x∗i , ρ∗τ (i)) is the unique solution of (26)– (29) satisfying
ρI (i) + ρIII (i) = 1 − ξ for every 0 ≤ i ≤ i1 .
In the rest of the paper, we fix ξ as in (40) and we set the type distribution ρτ (i) to be
∗
ρτ (i), the unique solution specified in Lemma 15.
6.2.1

Probability of a t-rejection

Lemma 16. The expected number of iterations during the phase of double edge reduction
is at most 10BD . Moreover, with probability 1 − o(1), the phase of double edge reduction
terminates within 8BD iterations.
Proof.
Let t ≥ 8BD be an integer. Assume the phase does not stop before step t. Let
x denote the number of steps among the first t steps that the number of double edges does
not decrease. Since the number of double edges in P0 is at most i1 , if Phase 5 lasts at least
t steps then i1 + x − (t − x) ≥ 0, since the number of double edges in each of the x steps
will increase by at most one whereas in the other t − x steps it decreases by exactly one. It
follows then that x ≥ (t − i1 )/2. The probability that the number of double edges (if there
is any) does not decrease in any step is bounded by 1 − min1≤i≤i1 {ρI (i)} < 1/10. Hence, the
probability that x ≥ (t − i1 )/2 is at most


t
(1/10)(t−i1 )/2 ≤ 2t (1/10)7t/16 < 0.8t ,
(t − i1 )/2
as (t − i1 )/2 ≤ (t − t/8)/2 = 7t/16. Putting t = 8i1 , the above probability is clearly o(1).
Hence, with probability 1 − o(1), this phase lasts for at most 8i1 steps. The expected number
of steps in this phase is at most
X
8i1 +
0.8t ≤ 10i1 .
t≥8i1

Recall that BD from (19) and recall that i1 is set to BD in Phase 5.
24

Lemma 17. Assume BD ξ = o(1). Then with probability 1 − o(1), no t-rejection happens in
Phase 5.
Proof.
By the definition of ρτ (i), the probability of a t-rejection in each iteration is at
most ξ. By Lemma 16, a.a.s. the number of iterations in this phase is at most 8i1 . Hence,
the probability of a t-rejection during this phase is at most 8i1 ξ + o(1) = o(1).
6.2.2

Probability of an f-rejection, a pre-b-rejection and a b-rejection.

By bounding the difference between fτ (P ) from mτ (i) we can bound the probability of an
f-rejection; similarly, by bounding the difference between bb(P, z) and m
b τ (i), and bounding
the difference between b(P ) and m(i), we can bound the probability of a pre-b-rejection and
a b-rejection in Phase 5. The arguments are similar to that in [9]. Thus we just state the
assumptions needed to have a small probability of a rejection and leave the proofs in the
Appendix.
Lemma 18. Assume that BD ξ = o(1), BD (M3 d2h + L3 d21 ) + L6 + M3 U3 + M2 U1 U2 = o(M13 )
2
and BD n(2γ−3)/(γ−1) −1 = o(1). The probability of an f-rejection during Phase 5 is o(1).
Lemma 19. Assume U1 M2 L2 /M12 + (M2 L2 /M12 )2 = o(M1 ). Then the probability of a preb-rejection during Phase 5 is o(1).
Lemma 20. Assume (dh M2 + d1 L2 + d21 d2h + U12 )M2 L2 /M12 + M2 U1 + L4 = o(M12 ). The
probability of a b-rejection during Phase 5 is o(1).

6.3

Phases 3 and 4: reduction of light loops and light triple edges

The switchings used in these two phases are the same as in the corresponding phases in [9],
with almost the same analysis. Thus, we briefly describe the switchings, define the switching
steps, and specify key parameters used in these two phases. Then we give the lemmas
which give bounds on f (P ) and b(P ), and bounds on the probabilities of f-rejections and
b-rejections. The proofs are almost the same as in [9] and thus will be presented in the
Appendix.
When Phase 3 starts, the input is a pairing uniformly distributed in A0 : the set of pairings
P such that G[H] (P ) is simple. Set i1 = BL in this phase. Then by (19), i1 = O(L2 /M1 ). In
Phase 3, Si is the set of pairings in A0 containing exactly i loops.
We use the switching shown in Figure 1 to reduce the number of loops.
Define
m` (i) = 2iM12 ;
m` (i) = L2 M1 − 2dh M1 (2i + 4BD + 6BT + idh /2) − L2 (2i + 4BD + 6BT + 6d1 + 2U1 ).
Lemma 21. For each P ∈ Si ,
m` (i) − O(iM1 (i + BD + BT + d1 + U1 )) ≤ f (P ) ≤ m` (i);

m` (i) ≤ b(P ) ≤ L2 M2 .

Let Pt denote the pairing obtained after step t of Phase 3. At step t ≥ 1, choose uniformly
at random a switching S applicable on Pt−1 . Assume that Pt−1 ∈ Si and P is the pairing
obtained if S is performed. Perform an f-rejection with probability 1 − f (Pt−1 )/m` (i) and a
b-rejection with probability 1 − m` (i − 1)/b(P ). If no rejection happened, then set Pt = P .
Repeat until Pt ∈ S0 .
25

As proved in [17, 9], by the end of Phase 3, the output is uniformly distributed in
A3 = {P ∈ A0 : L(P ) = 0}.
Lemma 22. Assume dh M1 (BD + BT + BL dh ) + L2 (d1 + U1 ) = o(M12 ). The probability of an
f-rejection or a b-rejection during Phase 3 is o(1).
Phase 4 starts with a pairing P0 uniformly distributed in A3 . We will use the switching
shown in Figure 2 to reduce the triple edges.
In this phase, i1 is set to BT defined in (19) and hence A3 = ∪0≤i≤i1 Si where Si is the
set of pairings in A3 containing exactly i triple edges.
Recalling Uk from (34) and (35), define
mt (i) = 12iM13
mt (i) = M3 L3 − 3M3 (4BD d2h + 6id2h ) − 3L3 (4BD d21 + 6id21 ) − L6
−16M3 U3 − 3M2 U1 U2 .
(41)
Lemma 23. For each P ∈ Si ,
mt (i) − O(iM12 (i + BD + d1 + U1 )) ≤ f (P ) ≤ mt (i);

mt (i) ≤ b(P ) ≤ M3 L3 .

The algorithm in Phase 4 is similar to that in Phase 3. At step t, the algorithm chooses
uniformly at random a switching S applicable on Pt−1 . Assume that Pt−1 ∈ Si and P is the
pairing obtained if S is performed. Perform an f-rejection with probability 1 − f (Pt−1 )/mt (i)
and a b-rejection with probability 1 − mt (i − 1)/b(P ). If no rejection happened, then set
Pt = P . Repeat until Pt ∈ S0 . The output then is uniformly distributed in A4 = {P ∈ A0 :
L(P ) = T (P ) = 0}.
Lemma 24. Assume (d1 +U1 )L3 M3 = o(M14 ) and M3 d2h (BD +BT )+L3 (BD d21 +BT d21 )+L6 +
M3 U3 + M2 U1 U2 = o(M13 ). The overall probability that either an f-rejection or a b-rejection
occurs during Phase 4 is o(1).

6.4

Fixing δ

We show that there exists δ such that the various conditions on γ and δ together with the
assumptions of Lemmas 4, 29,√ 30, 8, 12, 14, 17, 18, 19, 20, 22, and 24 are simultaneously
satisfied when γ > 21/10 + 61/10. After removing redundant constraints (see the full
details in Appendix), these conditions finally reduce to
f1 (γ) < δ < f2 (γ), 5/2 < γ < 3
(42)
where

f1 (γ) = max

1
4/(γ − 1) − 2 3 − γ
,
,
γ−2
γ − 2 2γ − 3


,

1
2
2 − 2/(γ − 1) − (2γ − 3)/(γ − 1)2 2 − 3/(γ − 1) o
f2 (γ) = min
,1 −
,
,
,
2
γ−1 7−γ
3−γ
4−γ
It is easy to see that (42) is feasible whenever
√
21/10 + 61/10 < γ < 3,
√
where 21/10 + 61/10 ≈ 2.881024968.
n1

26

For γ in that range, we may choose any
1
2 − 3/(γ − 1)
<δ<
.
2γ − 3
4−γ

6.5

Uniformity

√
Lemma 25. Assume that d is a plib sequence with parameter γ > 21/10 + 61/10 ≈
2.881024968 with even sum of the components. Then, the output of PLD is uniform over the
set of graphs with degree sequence d. Moreover, the probability of any rejection occurring in
Stage 2 (Phases 3–5) is o(1).
Proof. The uniformity follows by Lemma 9 and (31). By choosing δ satisfying (18), the
probability of any rejection occurring in Stage 2 is o(1) by Lemmas 4, 29, 30, 8, 12, 14, 17,
18, 19, 20, 22, and 24.

References
[1] A.L. Barabási and R. Albert. Emergence of scaling in random networks. science 286.5439
(1999): 509–512.
[2] J. Blitzstein and P. Diaconis. A sequential importance sampling algorithm for generating
random graphs with prescribed degrees. Internet mathematics 6.4 (2011): 489–522.
[3] M. Bayati, J. H. Kim, and A. Saberi. A sequential algorithm for generating random
graphs. Algorithmica, 58(4):860–910, 2010.
[4] E. A. Bender and E. R. Canfield. The asymptotic number of labeled graphs with given
degree sequences. J. Combinatorial Theory Ser. A, 24(3):296–307, 1978.
[5] B. Bollobás. A probabilistic proof of an asymptotic formula for the number of labelled
regular graphs. European J. Combin., 1(4):311–316, 1980.
[6] F. Chung and L. Lu. The volume of the giant component of a random graph with given
expected degrees. SIAM J. Discrete Math., 20(2):395–411, 2006.
[7] C. Cooper, M. Dyer, and C. Greenhill. Sampling regular graphs and a peer-to-peer
network. Combin. Probab. Comput., 16(4):557–593, 2007.
[8] P. Gao and N. Wormald. Enumeration of graphs with a heavy-tailed degree sequence.
Advances in Mathematics, 287:412–450, 2016.
[9] P. Gao and N. Wormald. Uniform generation of random regular graphs. Proc. FOCS’15:
1218–1230, IEEE.
[10] P. Gao and N. Wormald.
1511.01175.

Uniform generation of random regular graphs.

arXiv:

[11] C. Greenhill. The switch markov chain for sampling irregular graphs. Proc. SODA’15:
1564–1572, SIAM.
27

[12] S. Janson. The probability that a random multigraph is simple. Combin. Probab.
Comput., 18(1-2):205–225, 2009.
[13] M. Jerrum, B. McKay and A. Sinclair. When is a graphical sequence stable? Random
graphs, Vol. 2 (Poznań, 1989), 101–115.
[14] M. Jerrum and A. Sinclair. Fast uniform generation of regular graphs. Theoret. Comput.
Sci., 73(1):91–100, 1990.
[15] J. H. Kim and V. H. Vu. Sandwiching random graphs: universality between random
graph models. Adv. Math., 188(2):444–469, 2004.
[16] J. H. Kim and V. H. Vu. Generating random regular graphs. Combinatorica, 26(6):683–
708, 2006.
[17] B. D. McKay and N. C. Wormald. Uniform generation of random regular graphs of
moderate degree. J. Algorithms, 11(1):52–67, 1990.
[18] M. E. Newman. The structure and function of complex networks. SIAM review 45.2
(2003): 167–256.
[19] A. Steger and N. C. Wormald. Generating random regular graphs quickly. Combin.
Probab. Comput., 8(4):377–396, 1999. Random graphs and combinatorial structures
(Oberwolfach, 1997).
[20] G. Tinhofer. On the generation of random graphs with given properties and known
distribution. Appl. Comput. Sci., Ber. Prakt. Inf 13 (1979): 265–297.
[21] R. Van Der Hofstad. Random graphs and complex networks. Cambridge Series in
Statistical and probabilistic Mathematics 43 (2016).
[22] J. Y. Zhao. Expand and Contract: Sampling graphs with given degrees and other
combinatorial families. arXiv preprint arXiv:1308.6627 (2013).

28

Appendix
Bounds on Hk and Lk
By (3),
h
X
H1 =
O((n/i)1/(γ−1) ) = O(n1/(γ−1) h1−1/(γ−1) ) = O(n1−δγ+2δ ).

(43)

i=1

For each k ≥ 2, (4) gives

Hk ≤ Mk = O(nk/(γ−1) ).

(44)

Since γ > 5/2, (43) implies H1 = o(n) and so L1 = M1 − H1 = Θ(n). By (3) and (5),
n
X
Lk ≤
((Kn/i)1/(γ−1) )k = O(n1−δ(γ−1−k) ) for each fixed k ≥ 2.
(45)
i=n1−δ(γ−1)

Several parts of our analysis will require the following upper bound.
Lemma 26. Assume
1 − δ > 1/(γ − 1).

(46)

Then

X
i<j≤h


min

(di dj )k (di dj )
,
n
n2


k+1



2−γ

log n + n−2δ(γ−2)
k=0
O n
3−γ
= O (n
log n)
k=1



2k/(γ−1)−1
k+2−γ
O n
+n
log n k ≥ 2.

(47)

Proof. We wish to obtain an upper bound for
X

min (di dj )k /n, (di dj )k+1 /n2 ,
i<j≤h

for fixed integer k ≥ 0. This function takes the first term if di dj ≥ n and the second term
otherwise. We use f ≤a g to denote f = O(g) for simplicity. By (3), we know that di dj = O(n)
if ij > n3−γ . Thus, we obtain an upper bound for the minimum function by taking the first
term if ij ≤ n3−γ and taking the second term otherwise. Hence,


h
h
X
X
X
(di dj )k X X
(di dj )k+1
(di dj )k (di dj )k+1
,
+
≤
min
a
n
n2
n
n2
3−γ
i=1 3−γ
i=1
i<j≤h
j≤min{h,n

≤a n2k/(γ−1)−1

/i}

h
X

X

/i<j≤h

(ij)−k/(γ−1)

i=1 j≤min{h,n3−γ /i}

+n2(k+1)/(γ−1)−2
3−γ

n

h
X

X

i=1

n3−γ /i<j≤h

(ij)−(k+1)/(γ−1) .

(48)

By (5) and (46), we have n
< h. Hence, using (5) and (48), we obtain the desired
bound as in Lemma 26.
We bound Wi,j and Wi in the following lemma, which has been used to specify Φ0 .

29

Lemma 27. Let P be a pairing chosen uniformly at random from Φ. Then

 2 

X 
M2 H1
Wi,j Wj,i
,
+
= O
E mi,j Imi,j ≥2
3
d
d
M
i
j
1
i<j≤h

 2 
X 
Wi
M2 H1
,
E mi,i
= O
3
d
M
i
1
i≤h
X

E mi,j Imi,j ≥2

i<j≤h

X
i≤h



M22
= (1/2 + o(1)) 2 ,
M1

E (mi,i ) ≤ (1/2 + o(1))

M2
,
M1

where mi,j denotes the the number of pairs between i and j and mi,i denotes the number of
loops at i in P and Wi,j = Wi,j (P ), Wi = Wi (P ).
Proof. We will prove the first equation; the other equations follow with a similar argument.
Note that


X
X X  mi,j mi,w Imi,j ≥2,mi,w ≥2 
Wi,j
E mi,j Imi,j ≥2
≤
E
d
di
i
1≤i<j≤h
1≤i<j≤h 1≤w≤h


X
mi,j mi,i Imi,j ≥2,mi,i ≥1
E
+
di
1≤i<j≤h
Note also that E(mi,j Imi,j ≥2 ) ≤ E(mi,j ) ∼ di dj /M1 . Conditional on any possible value of mi,j ,
which is o(M1 ) always, the expectation of mi,w Imi,w ≥2 is at most the conditional expectation
of [mi,w ]2 which is O([di ]2 [dw ]2 /M12 ) and the conditional expectation of mi,i Imi,i ≥1 is at most
O(d2i /M1 ). Hence,


 2 
X
X X  d2 d2 dj 
X
Wi,j
di dj
i w
E mi,j Imi,j ≥2
=
O
+
O
3
di
M1
M12
1≤i<j≤h
1≤i<j≤h w≤h
1≤i<j≤h
= O(H2 H1 M2 /M13 + H2 H1 /M12 ) = O(M22 H1 /M13 ), (49)
where the last equation above holds by the fact H2 ≤ M2 and by our assumption (8).
By symmetry we have


X
Wj,i
E mi,j Imi,j ≥2
= O(M22 H1 /M13 ),
d
j
1≤i<j≤h
which together with (49) yields the first equation in the lemma.
Proof of Lemma 4. It follows as a direct corollary of Lemma 27.
Lemma 28. For every k ≥ 1, every m ≥ 1, and for every P ∈ C(Mk , i, j, m),
mfi,j (Mk , m) ≤ fi,j (P ) ≤ mfi,j (Mk , m);
for every P ∈ C(Mk , i, j, 0),
mbi,j (Mk , m) ≤ bi,j (P, m) ≤ mbi,j (Mk , m).

where (i, j) is the k-th element in J .
Proof. The upper bound for fi,j (P ) is obvious: there are m! ways to ordering the m
pairs between i and j and there are at most M1m ways to pick another m pairs and label
30

their end points. For the lower bound, note that the total number of pairs incident with
a heavy vertex is at most H1 and thus, the number of choices to pick the k-th pair (k ≤
m) that has not been picked by the first k − 1 pairs, and label its end points is at least
M1 − H1 − 2(k − 1) ≥ M1 − H1 − 2m. This proves the lower bound for fi,j (P ). Now, consider
bi,j (P 0 , m). Recall that Wi,j = Wi,j (P ) is the number of pairs that are in a multiple edge
and with one end in i and the other wind in some w 6= j in P . To perform an inverse
switching on P , pick m ordered points in i and their mates and m ordered points in j and
their mates. Switch these 2m pairs as indicated in Figure 5. One forbidden case is that
one of these 2m pair is contained in a heavy multiple edge. Hence, an upper bound for
bi,j (P, m) is [di − Wi,j ]m [dj − Wj,i ]m . The only other forbidden case is that the k-th pair
{2k − 1, 2m + 2k − 1} incident with i and the k-th pair {2k, 2m + 2k} incident with j are
such that both 2m + 2k − 1 and 2m + 2k are contained in heavy vertices for some 1 ≤ k ≤ m.
The number of ways to choose such two pairs is at most mh2 [di − Wij ]m−1 [dj − Wj,i ]m−1 ,
as m is the number of choices for k; [di − Wij ]m−1 [dj − Wj,i ]m−1 is an upper bound for the
number of choices for the other 2(m − 1) pairs and h2 is an upper bound for the number of
choices for {2k − 1, 2m + 2k} and {2k, 2m + 2k − 1}; both are contained in a single edge, as
there are only h heavy vertices. Note that we do not need to consider the case that these
pairs are contained in a multiple edge as this case has already been excluded in the upper
bound [di − Wi,j ]m [dj − Wj,i ]m . The lower bound follows.
Note that Lemma 7 follows from the following two lemma.
Lemma 29. The probability of an f -rejection during Phase 1 is o(1), if δ > (3 − γ)/(γ − 2),
δ > 1/(2γ − 3), δ < 1/2 and 1 − δ > 1/(γ − 1).
Lemma 30. The probability of a b-rejection during Phase 1 is o(1), if δ > (3 − γ)/(γ − 2),
δ > 1/(2γ − 3), δ < 1/2 and 1 − δ > 1/(γ − 1).
Proof for Lemma 29. Let (i, j) be the k-th element in J and P be a random pairing in
C(Mk−1 ). Let Yi,j be the random variable denoting the multiplicity of {i, j} in P . Note that
Yi,j is the ij entry of Mk−1 and thus the ij entry of M0 = M(P0 ) by the definition of Mk−1 .
An f-rejection can happen in sub-step (ii) or (v).
For (ii), the probability of an f-rejection is
!
X
fi,j (P )
Im(P,i,j)=m .
P(m(P, i, j) = m) 1 − f
m
(M
,
m)
k
i,j
m≥2
By (10) and (11), if m(P, i, j) = m then mfi,j (Mk , m)−fi,j (P ) ≤ mfi,j (Mk , m)−mfi,j (Mk , m) =
O(m!mM1m−1 n1−δγ+2δ ). Hence, the above probability is at most
X

P(m(P, i, j) = m) · O mn−δγ+2δ = O(n−δγ+2δ )E(Yi,j IYi,j ≥2 ).
m≥2

We will bound E(Yi,j IYi,j ≥2 ) by

O min



di dj d2i d2j
, 2
n
n


.

Then the above probability is bounded by

O(n−δγ+2δ ) min di dj /n, d2i d2j /n2 .
For (v), the probability that Pk is not set as Pk∗ is at most one. The probability of an
31

f-rejection at sub-step (v) is
P(Yi,j ≥ 2)(1 − bi,j (P

0

)/mbi,j (Mk , 1))


= min

di dj d2i d2j
, 2
n
n



· O(h2 /di dj ),

by (12) and (13).
Hence, the overall probability of an f-rejection is at most
X

O(n−δγ+2δ + h2 /di dj ) min di dj /n, d2i d2j /n2
i<j≤h

= O(n−δγ+2δ+3−γ log n) + O(h2 n−2δ(γ−2) ) = O(n3−γ−δ(γ−2) log n + n2−2δ(2γ−3) ),
by (5) and (47) and using δ < 1/2.
To ensure that the overall probability of an f-rejection is o(1), we need δ > (3−γ)/(γ −2),
δ > 1/(2γ − 3) besides the conditions δ < 1/2 and 1 − δ > 1/(γ − 1).
Proof for Lemma 30. A b-rejection can happen in sub-step (ii) or (v) in a switching step.
For (ii), the probability that a b-rejection happens is
!
X
mbi,j (Mk , m)
.
P(m(P, i, j) = m) 1 −
bi,j (P 0 , m)
m≥2
By (12), (13) and (9),


mbi,j (Mk , m)
mh2
=O
.
1−
bi,j (P 0 , m)
di dj
Hence the above probability is at most a constant times
X
mh2
h2
P(Yi,j = m)
=
E(Yi,j IYi,j ≥2 ).
di dj
di dj
m≥2
The overall probability of a b-rejection in sub-step (ii) is at most a constant times

!
X h2
X h2
di dj d2i d2j
E(Yi,j IYi,j ≥2 ) ≤ O
min
, 2
d
d
n
n
i dj
i dj
i<j≤h
i<j≤h

≤ O h2 n−2δ(γ−2)
by (47)

≤ O n2−2δ(2γ−3)
by (5).
For (v), the probability that a b-rejection happens is at most
!
mfi,j (Mk , 1)
P(m(P, i, j) ≥ 2) 1 −
.
fi,j (P 00 )

By (10) and (11),


mfi,j (Mk , 1)
H1
1−
=O
= O(n−δγ+2δ ),
00
fi,j (P )
M1
by (43). Hence, the above probability is at most a constant times



di dj d2i d2j
−δγ+2δ
−δγ+2δ
P(Yi,j ≥ 2) · n
≤O n
min
, 2
.
n
n

32

Now the overall probability of a b-rejection for (v) is at most a constant multiple of


X

di dj d2i d2j
−δγ+2δ
, 2
n
min
≤ O n−δγ+2δ+3−γ log n ,
n
n
i<j≤h
which is o(1) for γ and δ that satisfy the conditions of the lemma.
Proof for Lemma 8. The argument is analogous to, but simpler, than those in Lemmas 29
and 30, so we only give a simple sketch. The probability of an f -rejection in Phase 2 is
X
X
E(Yi,i ) · O(n−δγ+2δ ) =
O([di ]2 /n) · O(n−δγ+2δ ),
i≤h

i≤h

P
where Yi,i is the number of loops at i. Noting that i≤h [di ]2 ≤ M2 , the above probability is
O(n2/(γ−1)−δγ+2δ ), which is o(1) if 2/(γ − 1) < 1 + δ(γ − 2).
The probability of a b-rejection in Phase 2 is at most a constant multiple of
 3
X
h2
h
h2 X
2
O(di /n) 2 = O
,
E(Yi,i ) 2 =
di
di
n
i≤h
i≤h
which is o(1) by (5) if δ(γ − 1) > 2/3.
Proof for Lemma 9. Recall that G(P ) denotes the multigraph obtained from P . For any
vertex set S ⊆ [n], let G[S] (P ) denote the subgraph of G(P ) induced by S. We first confirm
that the output of Phase 1, if no rejection occurs, is uniformly distributed in Φ1 , defined by
Φ1 = {P ∈ Φ0 : G[H] (P ) contains no multiple edges ij with i 6= j}.
(50)

+
Given a pairing P ∈ C(Mk , i, j, 0), let Si,j,m
(P ) denote the set of inverse heavy m-way
+
0
switchings that convert P to some pairing P in C(Mk , i, j, m). Then |Si,j,m
(P )| = bi,j (P, m).
−
Conversely, let Si,j (P ) denote the set of heavy switchings that convert P ∈ C(Mk , i, j, m) to
−
some P 0 in C(Mk , i, j, m); then |Si,j
(P )| = fi,j (P ).
For each k ≥ 0, let (i, j) be the k-th element in J . We prove by induction on k that for
every k ≥ 0, the pairing obtained after the k-th step of Phase 1 is uniformly distributed in
C(Mk ).
The base case k = 0 is trivially true as P0 is uniformly distributed in C(M0 ), conditional
on the value of M0 . Assume k ≥ 1 and that the statement is true for k − 1. Note that
C(Mk−1 ) denote the set of pairings P such that the multiplicity of (i, j) in P is 0 or 1 for all
(i0 , j 0 )  (i, j).
Consider the k-th step in Phase 1. Let Pk∗ be as in (iii). Then, for any pairing P ∈
C(Mk , i, j, 0),
X
X
mbi,j (Mk , m)
1
P(Pk∗ = P | m(Pk−1 , i, j) ≥ 2) =
mfi,j (Mk , m) bi,j (P, m)
m≥2 0 +
P ∈Si,j,m (P )

=

X mbi,j (Mk , m)
mfi,j (Mk , m)
m≥2

Note that the final expression above does not depend on P .

33

.

(51)

Now, for any pairing P ∈ C(Mk , i, j, 0):
P(Pk∗ = P | m(Pk−1 , i, j) ≥ 2)
P(Pk = P ) = P(Pk−1 = P ) +
1 + mbi,j (Mk , 1)/mfi,j (Mk , 1)
X mbi,j (Mk , m)
1
= P(Pk−1 = P ) +
;
1 + mbi,j (Mk , 1)/mfi,j (Mk , 1) m≥2 mfi,j (Mk , m)

(52)

for any pairing P ∈ C(Mk , i, j, 1):
P(Pk = P ) = P(Pk−1 = P )
f
X
mbi,j (Mk , 1)/mfi,j (Mk , 1)
bi,j (P 0 ) mi,j (Mk , 1)
0
∗
+
P(Pk = P | m(Pk−1 , i, j) ≥ 2)
1 + mbi,j (Mk , 1)/mfi,j (Mk , 1) mbi,j (Mk , 1) fi,j (P )
P 0 ∈S − (i,j)(P )
= P(Pk−1 = P ) +
= P(Pk−1 = P ) +

mfi,j (Mk , 1)
1
mfi,j (Mk , m) 1 + mbi,j (Mk , 1)/mfi,j (Mk , 1) mbi,j (Mk , 1) fi,j (P )
P 0 ∈S − (i,j)(P ) m≥2
X

X mbi,j (Mk , m)

mbi,j (Mk , 1)/mfi,j (Mk , 1)

X mbi,j (Mk , m)

1

1 + mbi,j (Mk , 1)/mfi,j (Mk , 1) m≥2 mfi,j (Mk , m)

,

since |S − (i, j)(P )| = fi,j (P ).
By (52) and (53) and noting that C(Mk ) = C(Mk , i, j, 0) ∪ C(Mk , i, j, 1), for any P ∈
C(Mk ),
X mbi,j (Mk , m)
1
P(Pk = P ) = P(Pk−1 = P ) +
1 + mi,j (Mk , 1)/mi,j (Mk , 0) m≥2 mfi,j (Mk , m)
X mbi,j (Mk , m)
1
1
=
,
+
|C(Mk−1 )| 1 + mi,j (Mk , 1)/mi,j (Mk , 0) m≥2 mfi,j (Mk , m)
which is independent of P . Hence, Pk is uniformly distributed in C(Mk ).
With an analogous but simpler analysis as in Section 3.3.1, we can show that the output
of Phase 2, if no rejection happens, is uniformly distributed on Φ2 . This completes the proof.
Proof of Lemma 12. Given a pairing P , let P[H] be the set of pairs in P whose end vertices
are both in H. Conditional on (P0 )[H] be a given simple pairing P 0 , P0 can be generated by
pairing up all remaining points in H to points in L and then take a uniform pairing over all
remaining points in L. Let Φ(P 0 ) be this conditional probability space. It suffice to prove
that for every P 0 , with probability at least 1/2, a random pairing in Φ(P 0 ) contains at most
BL loops, BD double edges, mk,0 multiple edges with multiplicity k, for 3 ≤ k ≤ k0 and no
other types of multiple edges.
Let d0i = di − di (P 0 ) where di (P 0 ) is the degree of i in P 0 . In other words, d0i is theP
number
of remaining points in vertex i after removing all points used in P 0 . Let H10 = i∈H d0i .
Clearly, d0i ≤ di and so H10 ≤ H1 = O(n1−δγ+2δ ) = o(n) by (43) and by our assumption that
γ > 5/2.
Recall that the generation of Φ(P 0 ) has two stages: in the first stage, each remaining
point in H is paired randomly to a point in L; in the second stage, a uniform pairing is taken
over all remaining L01 := L1 − H10 points in L.
34

(53)

Let X denote the number of double loops in Φ(P 0 ): two loops at the same vertex. Since
P 0 is simple, the only possible way a double loop can be created is that for some vertex
i ∈ L, there are four points pj , 1 ≤ j ≤ 4, in i such that p1 is paired to p2 and p3 is paired
to p4 . The probability for this to happen is
1
−2
∼ L−2
1 = O(n ),
0
(L1 − 1)(L01 − 3)
as L1 = M1 − o(n) = Θ(n) and H10 ≤ H1 = o(n). The total number of ways to choose
{pj : 1 ≤ j ≤ 4} is at most
X di 
3
.
4
i∈L
Hence,



L4
EX = O
= O(L4 /M12 ) = o(1).
L21
by (45) and the assumption that δ(5 − γ) < 1.
Let Y denote the number of multiple edges with multiplicity at least 4. Such multiple
edges can be created either in the first stage or in the second stage; one of the end vertices
must be light.
EY = O(M4 L4 /L41 ) = O(n4/(γ−1)+1−δ(γ−5)−4 ) = o(1).
by (44) and (45).
Let L = L(P ) denote the number of loops in P for P ∈ Φ(P 0 ). Then,

X di
L2
2
EL =
∼
.
0
L
−
1
2M
1
1
i∈L
Hence, with probability at most 1/8 + o(1) will L exceeds BL .
Let Yk denote the number of multiple edges with multiplicity k. It only remains to verify
that with probability at least 1/4 + o(1), Y2 ≤ BD and Y3 ≤ BT .
  !
P
P
di dj
k!
i∈[n]
j∈L k
k
EYk ≤ (1 + o(1))
k
L1
Lk Mk
.
k!M1k
Hence, with probability 1 = o(1), there are no multiple edges with multiplicity greater than
3 by our assumptions in the lemma. Moreover, with probability at most 1/8 + o(1), the
number of double edges is more than BD , and with probability at most 1/12 + o(1), the
number of triple edges is more than BT . By taking the union bound, the probability that P
does not satisfy conditions in (b) is less than 1/2 + o(1).
Proof of Lemma 13. We first prove part (b). It is easy to see that, for any fixed z ∈ Zt(τ ) , the
only condition for the choices of the additional pairs is that (a) each pair is contained in a
single edge; (b) for each pair, their end vertices are not allowed to be adjacent to two certain
vertices. There are at most 4i choices to violate (a) and at most 2U1 choices to violate (b)
since the number of 2-paths starting from a given vertex is bounded by U1 by (54). Therefore,
the number of choices for each pair is at least M1 − 4i − 2U1 . The lower bounds in (b) follows
by the illustrations of the inverse switchings of each type shown in Figures 3–12.
≤ (1 + o(1))

35

For (a), it is sufficient to bound Z ∗ (P ) uniformly by 8i(dh M2 + d1 L2 ) + (2id21 d2h + 4iU12 +
8M2 U1 + L4 ) by (20). Recall that Z ∗ (P ) is the size of Z ∗ (P ) = F(P ) \ Z(P ). Elements in
Z ∗ (P ) must fall into at least one of the following cases.
(a1) one of the two stars contains a pair that is contained in a double edge;
(a2) uj and vj are adjacent with a double edge, for some j = 1, 2, 3;
(a3) the two stars share a vertex.
For (a1), there are 4i choices to choose a point that is contained in a double edge. Label
this point 1 or 3 or 2 or 4. If it is labelled 1 (or 3), there are at most dh ways to choose the
point 3 (or 1) and there are at most M2 ways to choose the other 2-star. If it is labelled 2
(or 4), there are at most d1 ways to choose the point 4 (or 2) and then there are at most
L2 ways to choose the other light 2-star. Hence, the number of choices for (a1) is at most
4i · 2 · dh · M2 + 4i · 2 · d1 · L2 = 8i(dh M2 + d1 L2 ).
For (a2), we discuss two cases: j = 1 and j = 2, 3. If j = 1, there are at most 2i ways to
choose vertices u1 and v1 . There are at most d2h ways to choose points and label them by 1
and 3 in u1 and at most d21 ways to choose points and label them by 2 and 3 in v1 . So the
number of choices in this case is at most 2id21 d2h . If j = 2, then there are 2i ways to choose
vertices u2 and v2 ; then at most U12 ways to fix a 2-path starting at u2 and a 2-path starting
at v2 . Therefore, the number of choices for j = 2 is at most 2iU12 . Similarly, this holds for
j = 3. Hence, the total number of choices for (a2) is at most 2id21 d2h + 4iU12 .
For (a3), we discuss three cases: (i) v1 ∈ {u2 , u3 }, the number of such choices is at most
2M2 U1 ; (ii) v1 = u1 , then u1 must be light and so the number of such choices is at most L4 ;
(iii) vj ∈ {u1 , u2 , u3 } for some j = 2, 3, then the number of such choices is at most 6L2 U1 .
Hence, the total number of choices for (a3) is at most 2M2 U1 + L4 + 6L2 U1 ≤ 8M2 U1 + L4
as L2 ≤ M2 .
This implies the lower bound for b(P ) as desired.
Before bounding the probabilities of various types of f- and b-rejections, we bound the
number of certain structures in P , uniformly for all P ∈ A0 . These bounds turn out to be
Uk defined in (34) and (35).
2-paths at v
For an integer i ≥ 1, an i-path in a pairing P is a sequence of i pairs: {p1 , p2 }, . . . , {p2i−1 , p2i }
such that p2j and p2j+1 are in the same vertex for each 1 ≤ j ≤ i − 1.
Let Ai (P, v) denote the number of i-paths starting at v in P and let
Ai (d) = max max Ai (P, v).
P ∈Φ v∈[n]

We will bound A2 = A2 (d).
For any u ∈ [n], the degree of u is at most d1 = ∆. The number of 2-paths from u is at
most the sum of degrees of vertices incident with u, which is at most
d1
∆
X
X
(γ−2)/(γ−1)
(Kn/i)1/(γ−1) = (Kn)1/(γ−1)
i−1/(γ−1) ≤ (Kn)1/(γ−1) d1
.
i=1

i=1

1/(γ−1)

Using d1 ≤ (Kn)

) by (4), we have
2

2

A2 ≤ (Kn)1/(γ−1) (Kn)(γ−2)/(γ−1) = (Kn)(2γ−3)/(γ−1) = U1 .
36

(54)

Light k-blooms at v
A light k-bloom at vertex v consists of k + 1 pairs {p2i−1 , p2i }, 1 ≤ i ≤ k + 1 such that
p1 ∈ v; p2 , p2i−1 for 2 ≤ i ≤ k + 1 are in the same vertex that is light. See Figure 13 for
an example where k = 2. Note that a light k-bloom at v is also a light (k + 1)-star at u
where u is the vertex containing points p2 and p2i−1 , 2 ≤ i ≤ k + 1. Let Bk (P, v) denote
p4
p1

v

p2
p3
p5

p6

Figure 13: 2-bloom
the number of light k-blooms at v in P ; and let Bk (d) = maxP ∈Φ maxv∈[n] Bk (P, v). We will
bound Bk = Bk (d) for k ≥ 2. Similar to the calculations for A2 ,
h+∆
h+∆
X
X
X
(Kn/i)k/(γ−1) ≤ (Kn)k/(γ−1)
i−k/(γ−1)
dki ≤
Bk ≤
i=h+1

i=h+1

i≥h+1

γ
≤
(Kn)k/(γ−1) h(γ−k−1)/(γ−1) .
k+1−γ
By (5), it follows that
γK k/(γ−1) 1−δ(γ−k−1)
n
= Uk for k ≥ 2.
(55)
Bk ≤
k+1−γ
Proof of Lemma 18. For a given pairing P , we first bound the probability that a random
switching performed on P is f-rejected. Assume P ∈ Si ; the probability that a type τ
switching is performed on P and is f-rejected, is


fτ (P )
.
ρτ (i) 1 −
mτ (i)
Summing over all τ ∈ Λ we obtain that assuming P arises in the algorithm, the probability
that a random switching performed on P is f -rejected is


X
fτ (P )
ρτ (i) 1 −
.
(56)
mI (i)
τ
Now we bound the probability that an f-rejection happens during Phase 5. Recall that Pt is
the pairing obtained after step t of Phase 3. Then, by (56), the probability of an f-rejection
during Phase 3 is

!
X X X
X
fτ (P )
P(Pt = P )
ρτ (i) 1 −
.
m
(i)
τ
τ
t≥0 1≤i≤i1 P ∈Si
P
Note that t≥0 P(Pt = P ) is the expected number of times that P is reached, which is σ(i)

37

for every P ∈ Si . Hence, the above probability is
!

X X
X
fτ (P )
.
σ(i)
ρτ (i) 1 −
m
(i)
τ
τ
1≤i≤i P ∈S
1

(57)

i

Note that for each τ ,




X X
X
fτ (P )
fτ (P )
σ(i)ρτ (i) 1 −
=
σ(i)ρτ (i)|Si | · E 1 −
,
mτ (i)
mτ (i)
1≤i≤i P ∈S
1≤i≤i
1

1

i

where the above expectation is taken on a random pairing P ∈ Si . It will be sufficient to
show that the product of this expectation and ρτ (i) is o(1/i1 ) for every 1 ≤ i ≤ i1 and for
every τ , since then

!


X X
X
X X
fτ (P )
fτ (P )
σ(i)
ρτ (i) 1 −
=
σ(i)ρτ (i)|Si | · E 1 −
m
(i)
mτ (i)
τ
τ
τ 1≤i≤i1
1≤i≤i1 P ∈Si
X X
= o(1/i1 )
σ(i)|Si |
1≤i≤i1 P ∈Si

P

P

and the lemma follows since 1≤i≤i1 P ∈Si σ(i)|Si | is the expected number of steps that the
phase takes, which is O(i1 ) by Lemma 16.
For τ ∈ Λ \ {I, III} we have ρτ ≤ ξ and so


fτ (P )
= O(ρτ (i)) = o(1/i1 )
(58)
ρτ (i)E 1 −
mτ (i)
by our assumption BD ξ = o(1).
For τ = III, by Lemma 15,


M3 L3
ρIII (i) = O
.
M2 L2 M1
The basic operation of type III switching is as in Figure 3, repeated in Figure 14 for convenience.

u2

u1

u3

u4

u2

u1

u3

u4

v2

v1

v3

v4

v2

v1

v3

v4

Figure 14: Type III again
To perform a type III switching on P , we randomly and independently choose a 3-star
and a light 3-star. Label the points in the stars and the related vertices as in the figure (here
u1 is light). Perform an f-rejection if a type III switching cannot be performed with the given
choice. Now E(1 − fτ (P )/mτ (i)) is the probability that an f-rejection is performed.
An f-rejection happens if
(a) one of the 3-stars contains a double edge;
(b) the two 3-stars share a vertex;
38

(c) there is an edge between uj and vj for some j = 1, 2, 3, 4.
The number of choices for (a) is at most 3M3 · 4i1 d2h + 3L3 · 4i1 d21 .
For (b), there are four cases: (b1), u1 = v1 . Since u1 is light, the number of such choices
is at most L6 ; (b2), v1 ∈ {u2 , u3 , u4 }. Recall the light 2-bloom defined above Figure 13. The
number of choices for (b2) is at most 3M3 B2 ≤ 3M3 U2 ; (b3), vj = u1 for some j = 2, 3, 4.
The number of such choices is at most 3M2 B3 ≤ 3M2 U3 ; (b4), vj ∈ {u2 , u3 , u4 } for some
j = 2, 3, 4, the number of such choices is at most 9M3 B2 ≤ 9M3 U2 . Hence, the total number
of choices for (b) is at most L6 + 12M3 U2 + 3M2 U3 .
For (c), there are two cases: (c1), u1 is adjacent to v1 . The number of such choices is
at most M3 B3 ≤ M3 U3 ; (c2), uj is adjacent to vj for some j = 2, 3, 4. The number of such
choices is at most 3M2 A2 B2 ≤ 3M2 U1 U2 . So the total number of choices for (c) is bounded
by M3 U3 + 3M2 U1 U2 .
Hence, the probability of an f-rejection in each step is bounded by
M3 i1 d2h + L3 i1 d21 + L6 + M3 U3 + M2 U1 U2
,
M3 L3
by noting that M3 U2 = O(M3 U3 ) and M2 U3 = O(M3 U3 ). Since ρIII (i) = O(M3 L3 /M2 L2 M1 )
and i1 = O(M2 L2 /M12 ), (58) is satisfied for τ = III as long as
M2 L2
+ L6 + M3 U3 + M2 U1 U2 = o(M13 ),
(M3 d2h + L3 d21 )
2
M1
which is guaranteed by the hypotheses of the lemma.
Finally, for τ = I, we use the bound ρI (i) ≤ 1. An f-rejection occurs if u2 and v2 is chosen
such that they are not distinct from u1 and v1 , or there is an existing edge between u1 and
u2 or an existing edge between v1 and v2 ; same for u3 and v3 . The total number of such
2
invalid choices for u2 and v2 is O(∆ + A2 ), where A2 is bounded by O(n(2γ−3)/(γ−1) ) in (54).
2
Hence, the probability of an f-rejection in each step is bounded by O(n(2γ−3)/(γ−1) − 1),
which satisfies (58) by the assumption of the lemma.
Proof of Lemma 19. Note that a pre-b-rejection only happens when type τ 6= I switchings are
performed. At a given step t, we estimate the probability that a given switching S = (P 0 , P, z)
is performed and is pre-b-rejected at step t. Here S converts P 0 to P and z is the ordered
pair of 2-stars created by S. Note that there can be several switchings converting P 0 to P ,
corresponding to different z, due to a different labelling of points in the involved pairs that
are switched. Let τ = τ (P 0 , P ) denote the type of S. Assume that P ∈ Si and so P 0 ∈ Si .
The probability that S is performed, conditional on Pt−1 = P 0 , is
ρτ (i)
.
mτ (i)
Conditional on that, the probability that S is pre-b-rejected is 1 − m
b τ (i)/bb(P, z).
+
0
Let Sτ (P ) denote the set of pairings P that can be switched to P via a type τ switching.

39

The probability that a pre-b-rejection occurs during Phase 5 is at most
!
XX X X
X
X
m
b
ρ
(i)
(i)
τ
1− τ
P(Pt−1 = P 0 )
bb(P, z)
m
(i)
τ
t≥1 τ 6=I 1≤i≤i1 P ∈Si P 0 ∈Si ∩Sτ+ (P ) z
!
X ρτ (i)
X XX
X
m
b τ (i)
σi
=
1−
,
bb(P, z)
m
(i)
τ
+
1≤i≤i1 P ∈Si τ 6=I P 0 ∈Si ∩Sτ (P ) z
P
as t≥1 P(Pt−1 = P 0 ) = σ(P 0 ) which is σj for every P 0 ∈ Sj . The z in the summation are
determined by (P 0 , P ), and there can be several of them, due to different labelling of points
in the pairs that are switched. For each τ 6= I, define Uτ to be the same as m
b τ (i) but without
b
terms −2U1 − 4i. Note that Uτ is a trivial upper bound for b(P, z). Now, by Lemma 13, for
every τ 6= I,


U1 + i
m
b τ (i)
=O
.
1−
bb(P, z)
M1
P
Noting that given τ and z ∈ Zt(τ ) (P ), P 0 ∈Si ∩Sτ+ (P ) 1 = bb(P, z) ≤ Uτ , and noting also that
|Zt(τ ) (P )| ≤ M2 L2 for any τ , the above probability of a pre-b-rejection becomes
!
!
X XX
X X ρτ (i)
ρτ (i) U1 + i
U1 + i1 X
O
M2 L2 Uτ σi
= O M2 L2
Uτ
σi
.
mτ (i) M1
M1 τ 6=I
mτ (i)
1≤i≤i P ∈S τ 6=I
1≤i≤i P ∈S
1

1

i

i

Note that mτ (i) does not depend on
 τ 6= I; so we may write mτ for mτ (i) and
Pi for every
rewrite the above probability as O
τ 6=I h(τ ) , where

since by Lemma 16,

(U1 + i1 )i1 Uτ ∗
ρ , ρ∗τ = max ρτ (i),
1≤i≤i1
M1
mτ τ
1≤i≤i1
P ∈S σi ≤ 10i1 and so
X X i
X X
σi ρτ (i) ≤ ρ∗τ
σi = O(i1 ρ∗τ ).

h(τ ) = M2 L2
P
P

1≤i≤i1 P ∈Si

1≤i≤i1 P ∈Si

b τ (i) = Ω(Uτ ) for each i, by Lemma 15, ρ∗τ = O (mτ /M2 L2 Uτ ).
Since m(i) = Ω(M2 L2 ) and m
Hence,




(U1 + i1 )i1
mτ
(U1 + i1 )i1 Uτ
O
=O
, for each τ 6= I.
h(τ ) = M2 L2
M1
mτ
M2 L2 Uτ
M1
Noting that i1 = O(M2 L2 /M12 ), now h(τ ) = o(1) for every τ 6= I by the hypotheses of the
lemma and so the probability of a pre-b-rejection in Phase 5 is o(1).
Proof of Lemma 20. Now, we bound the probability of a b-rejection during Phase 5. At a
given step t, we estimate the probability that a given switching S = (P 0 , P ) is performed
and is b-rejected at step t. Let τ = τ (P 0 , P ) denote the type of S. Assume that P 0 ∈ Sj and
P ∈ Si . Hence, j = i + 1 if τ = I and j = i otherwise.
The probability that S is performed, conditional on Pt−1 = P 0 , is
ρτ (j)
.
mτ (j)
Conditional on that, the probability that S is b-rejected is 1 − m(i)/b(P ). Let Sτ (P 0 , P )
denote the set of type τ switchings that converts P 0 to P . So the probability that a b-rejection
40

occurs during Phase 3 is at most
XX X X X



m(i)
ρτ (j)
1−
P(Pt−1 = P )
m
(j)
b(P )
τ
0
+
t≥1 τ 1≤i≤i1 1≤j≤i1 P ∈Si P 0 ∈Sj ∩Sτ (P ) S∈Sτ (P ,P )


X XX X
X
X
ρτ (j)
m(i)
=
σj
1−
,
mτ (j)
b(P )
0
+
0
1≤i≤i P ∈S τ 1≤j≤i
1

i

1

X

X

0

P ∈Sj ∩Sτ (P ) S∈Sτ (P ,P )

P(Pt−1 = P 0 ) = σ(P 0 ) which is σj for every P 0 ∈ Sj .
Note that after fixing i, for every τ and j such that Sj ∩ Sτ+ (P ) is non-empty for P ∈ Si ,
ρτ (j) were chosen so that
σj ρτ (j)
= q(i).
mτ (j)
We also have
X X
X
X
b(P ) =
1.
as

P

t≥1

τ

So the above expression is
X X
1≤i≤i1 P ∈Si

1≤j≤i1 P 0 ∈Sj ∩Sτ+ (P ) S∈Sτ (P 0 ,P )

q(i)(b(P ) − m(i)) =

X
1≤i≤i1

|Si |q(i)E(b(P ) − m(i)),

where the expectation is taken over a random P ∈ Si .
The above is
X
X
σ(i + 1)
E(b(P ) − m(i))
|Si |
|Si |q(i)E(b(P ) − m(i)) ≤
mI (i + 1)
1≤i≤i1
1≤i≤i1
!
X
Eb(P ) − m(i)
= O
|Si |σ(i)
m(i)
1≤i≤i1
P
by (25) and (30). By Lemma 16, 1≤i≤i1 |Si |σ(i) = O(i1 ) and so the above is o(1) if we can
show that for every 0 ≤ i ≤ i1 ,
Eb(P ) − m(i)
= o(1/i1 ).
(59)
m(i)
By (20), (36) and Lemma 13 we have
b(P ) = M2 L2 − O(i1 (dh M2 + d1 L2 + d21 d2h + U12 ) + M2 U1 + L4 ).

We also have i1 = O(M2 L2 /n2 ). Hence,


i1 (dh M2 + d1 L2 + d21 d2h + U12 ) + M2 U1 + L4
Eb(P ) − m(i)
· i1 = O
.
m(i)
M12
By the hypotheses of the lemma, the probability of a b-rejection during Phase 5 is o(1).
Proof of Lemma 21. The proof is a slight modification of [10, Lemma 15]. For the upper
bound of f (P ), the number of ways to choose a loop and label its end points in P is 2i.
Then, the number of ways to choose another two pairs and label their end points is at most
M12 . This shows that f (P ) ≤ m` (i) for any P ∈ Si . For the lower bound, note that each pair
of y and z cannot be contained in a loop or a multiple edge, the total number of such choices
is bounded by O(iM1 (i + BD + BT )). Moreover, u1 and u2 (same for u3 and u4 ) cannot be
the same vertex as u0 or be adjacent to u0 . The total number of forbidden choices like this
41

is bounded by iM1 (d1 + U1 ). This verifies the lower bound for f (P ).
Next we prove the bounds for b(P ). To perform an inverse switching on P , we pick a
light 2-star and label the points in the pairs by {3,1} and {2,5}, with 1 and 2 being in the
same light vertex u0 , together with another pair x whose end points are labeled by 4 and 6.
There are L2 ways to pick the light 2-star and M1 ways to pick and label x. This yields the
upper bound for b(P ). We need to exclude choices so that
(a) the 2-star contains a loop, or a double edge or a triple edge; or there is a loop at u0 ;
(b) x is a loop or is contained in a multiple edge;
(c) x and the 2-star share a vertex;
(d) u1 and u2 are adjacent or u3 and u4 are adjacent.
We first show that the number of choices for (a) is at most 2(2idM1 +4BD dM1 +6BT dM1 ).
We first count the choices so that {1,3} is a loop. There are 2i ways to choose a loop and
label their end points by 1 and 3. Then, there are at most dh ways to choose 2 and 5 as u0 ,
the vertex containing point 1 has degree at most dh . Therefore, the number of ways to choose
the 2-star and x so that {1,3} is a loop is at most 2idh M1 . Similarly, the number of such
choices so that {1, 3} is contained in a double edge or a triple edge is at most 4BD dh M1 and
6BT dh M1 respectively as there are at most BD double edges and BT triple edges. The extra
factor 2 accounts for the case that {2,5} is a loop or is contained in a double edge or triple
edge. The number of choices so that there is a loop at u0 is at most id2h M1 as there are i ways
to choose u0 . Hence, the number of choices for (a) is at most 2dh M1 (2i+4BD +6BT +idh /2).
The number of choices for (b) is at most L2 (2i + 4BD + 6BT ), for (c) is at most 6L2 d1
and for (d) is at most 2L2 U1 . Thus
b(P ) ≥ L2 M1 − 2dh M1 (2i + 4BD + 6BT + idh /2) − L2 (2i + 4BD + 6BT + 6d1 + 2U1 ) = m(i).
Proof of Lemma 22. By Lemma 21, the probability of an f-rejection in each iteration is


i + BD + BT + d1 + U1
O
,
M1
and the probability of a b-rejection in each iteration is


dh M1 (i1 + BD + BT + i1 dh ) + L2 (i1 + BD + BT + d1 + U1 )
O
L2 M1
Since the number of loops decreases by one in each iteration, Phase 3 lasts for at most
i1 ≤ BL = O(L2 /M1 ) steps. Now the probability of an f-rejection or a b-rejection in Phase
3 is bounded by


dh M1 (BL + BD + BT + BL dh ) + L2 (i0 + BD + BT + d1 + U1 )
O
M12


dh M1 (BD + BT + BL dh ) + L2 (d1 + U1 )
=O
,
M12
by noting that L2 ≤ dh M1 . Hence, this probability is o(1) by our hypotheses.
Proof of Lemma 23. For f (P ): there are i ways to choose a triple edge and then 12 ways to
label their end points and then at most M13 ways to choose the other 3 pairs. This gives the
upper bound. The lower bound follows by an analogous argument as in Lemma 21.
42

For b(P ): in order to perform an inverse switchings on P , we choose an ordered pair of
3-stars, the first being light, and label their end points as in Figure 2 (here u1 is required
to be light). The number of such choices is L3 M3 . This gives the upper bound for b(P ).
Compare this inverse switching with the type III switching in Phase 5. It is easy to see that
b(P ) equals the number of type III switchings applicable on P . Hence, the same analysis
in Lemma 18 for the lower bound on the number of type III switchings on a pairing P (the
arguments below (58)) works here for the lower bound on b(P ), except that for case (a):
one of the 3-stars contains a multiple edge, we need to consider also contribution from triple
edges instead of just double edge. This gives 3M3 (4BD d2h + 6id2h ) + 3L3 (4BD d21 + 6id21 ). The
discussions for cases (b,c) are exactly the same. Thus, we obtain that
b(P ) ≥ M3 L3 − 3M3 (4BD d2h + 6id2h ) − 3L3 (4BD d21 + 6id21 ) − L6 − 12M3 U2 − 3M2 U3
−M3 U3 − 3M2 U1 U2 ≥ mt (i),

by noting that M3 U3 ≥ M2 U3 and M3 U3 ≥ M3 U2 .
Proof of Lemma 23. The number of iterations in Phase 4 is at most i1 = O(L3 M3 /M13 ). In
iteration t, assuming Pt−1 ∈ Si , the probability of an f-rejection is


i1 + BD + d1 + U1
f (Pt−1 )
=O
1−
mt(i)
M1
by Lemma 23.
Similarly, the probability of a b-rejection in iteration t is at most


M3 (BD d2h + id2h ) + L3 (BD d21 + id21 ) + L6 + M3 U3 + M2 U1 U2
mt (i − 1)
=O
1−
.
maxP ∈Si−1 b(P )
M3 L3
Hence, the probability of an f-rejection or a b-rejection during Phase 4 is


(BT + BD + d1 + U1 )L3 M3 M3 d2h (BD + BT ) + L3 (BD d21 + BT d21 ) + L6 + M3 U3 + M2 U1 U2
+
O
M4
M13
which is o(1) by our lemma hypothesis by noting that (BT +BD )L3 M3 = O((BT +BD )M3 M1 d2h ).
Proof for Lemma 10. We prove bounds for the computation time of fi,j (P ) and bi,j (P, m)
only; the arguments for fi (P ) and bi (P, m) are similar.
Given P , let m denote the number of pairs between i and j. Assume that m ≥ 1 and we
want to compute fi,j (P ) exactly. We say that an ordered pair (p1 , p2 ) in P is nice if either
(i) both p1 and p2 are contained in light vertices, or (ii) p1 is contained in a light vertex, and
p2 is contained in a heavy vertex w 6= j such that there is no pair between j and w, or (iii)
p2 is contained in a light vertex and p1 is contained in a heavy vertex w 6= i such that there
is no pair between i and w. To perform a heavy multiple edge switching, we need to order
the m pairs between i and j and switch them away together with another m nice ordered
pairs, as shown in Figure 5. There are m! ways to order the m pairs between i and j. Then,
fi,j (P ) equals to m! multiplied by the number of ways to sequentially choose m distinct nice
ordered pairs, (pk1 , pk2 ), . . . , (pk2m−1 , pk2m ) satisfying
(a) no pk2g−1 and pk2h−1 are in the same heavy vertex for some 1 ≤ g < h ≤ m;
(b) no pk2g and pk2h are in the same heavy vertex for some 1 ≤ g < h ≤ m.

We describe a computing scheme for this number. Let I be the set of heavy vertices that
are not adjacent to i and J be the set of heavy vertices that are not adjacent to j in P ;
43

note that I ∩ J is not necessarily empty. Let 0 ≤ ` ≤ m. We first compute the number
of choices such that there are exactly ` nice ordered pairs of cases (ii) or (iii), denoted by
(pk1 , pk2 ), . . . , (pk2`−1 , pk2` ), and exactly m − ` nicely ordered pairs of case (i). Let Z denote
the total number of nice ordered pairs of case (i) in P . For each vertex w ∈ I, it can contain
at most one point from pk2g−1 , for 1 ≤ g ≤ `, and for each vertex v ∈ J, it can contain at
most one point from pk2g , for 1 ≤ g ≤ `. For each w ∈ I, we can compute Xw , the number
of nice ordered pairs (p1 , p2 ) with p1 ∈ w. For each v ∈ J, we can compute Yv , the number
of nice ordered pairs (p1 , p2 ) with p2 ∈ v. Now, the generating function for the choices of the
nice ordered pairs satisfying conditions (a,b) is
Y
Y
(1 + Xw z) (1 + Yv z).
w∈I

v∈J

Thus, the number of such choices with exactly ` pairs of cases (ii) or (iii) is
 
Y
Y
m
[Z]m−` [z ` ]
(1 + Xw z) (1 + Yv z),
`
w∈I
v∈J

m
where ` is the number of ways to choose the ` out of m pairs of nice ordered pairs to be of
cases (ii) or (iii). For the m − ` pairs of Q
case (i), there are
Q [Z]m−` ways to sequentially take
`
m − ` distinct pairs. There are then [z ] w∈I (1 + Xw z) v∈J (1 + Yv z) ways to sequentially
choose ` nice ordered pairs satisfying conditions (a,b). So,
X m
Y
Y
fi,j (P ) =
[Z]m−` · [z ` ]
(1 + Xw z) (1 + Yv z).
`
w∈I
v∈J
0≤`≤m
Thus, computing all of these numbers take O(H1 )-time. Similarly, computing Z takes O(H1 )time as Z = M1 − Z 0 , where Z 0 counts the number of pairs incident with at least one
heavy
vertex, and soQcomputing Z 0 takes only O(H1 )-time. Next, we show how to compute
Q
[z ` ] w∈I (1 + Xw z) v∈J (1 + Yv z) for each 0 ≤ ` ≤ m.
Denote I by w1 , w2 , . . . wg and J by v1 , v2 , . . . , vh . Let a = (a0 , . . . , am ) be a vector,
initialised by a0 = 1 and ak = 0 for all 1 ≤ k ≤ m. Note that a is the coefficient sequence of
P0 (z) = 1 (we only track the first m + 1 coefficients).
Q Let P1 (z) =Q(1 + Xw1 z) and P2 (z) = P1 (z)(1 + Xw2 z) and repeat. Thus, Pg+h (z) =
v∈J (1 + Yv z). We will update a so that after iteration k, a is the coefficient
w∈I (1 + Xw z)
sequence of Pk (z). By the inductive definition of Pk (z), it is easy to see that, for each
0 ≤ ` ≤ m,
 `
[z ]Pk−1 (z) + Xwk [z `−1 ]Pk−1 (z) for 1 ≤ k ≤ g
`
[z ]Pk (z) =
(60)
[z ` ]Pk−1 (z) + Yvk−g [z `−1 ]Pk−1 (z) for g + 1 ≤ k ≤ g + h

Note that (60) gives the updating rule for a in each iteration.
Since each
Q
Q iteration runs
`
in O(m)-time, the total number of steps for computing [z ] w∈I (1 + Xw z) v∈J (1 + Yv z) for
every 0 ≤ ` ≤ m is O(m(|I| + |J|)) = O(m|H|). Thus, the total number of steps required to
compute fi,j (P ) is O(H1 + m|H| + m) = O(H1 + m|H|).
Next, we bound the computation time for bi,j (P, m), where P is a pairing with no pairs
between i and j and m ≥ 1 is an integer. Recall that Wi,j (P ) denotes the number of points in
i that belong to heavy loops or heavy multiple edges with one end in i and the other end not
in j in P . Now for 1 ≤ k ≤ m, let Bk denote the set of choices of pairs {2g − 1, 2m + 2g − 1},
{2g, 2m + 2g} for 1 ≤ g ≤ m, as in Definition 5 for the inverse switching, such that (a) point
44

2g − 1 is not counted by Wi,j (P ) and point 2g is not counted by Wj,i (P ), and (b) switching
{2k − 1, 2m + 2k − 1} and {2k, 2m + 2k} to {2k − 1, 2k}, {2m + 2k − 1, 2m + 2k} would
create new heavy loops or heavy multiple edges. Then,
bi,j (P, m) = [di − Wi,j ]m [dj − Wj,i ]m − | ∪1≤k≤m Bk |.

It is easy to see that Wi,j and Wj,i can be computed in O(∆)-time. By the inclusionexclusion principle, we can express | ∪1≤k≤m Bk | by
m
X
X
X
|Bk |−
|Bk1 ∩Bk2 |+
|Bk1 ∩Bk2 ∩Bk3 |−· · ·+(−1)m−1 |B1 ∩B2 ∩· · ·∩Bk |.
k=1

1≤k1 <k2 ≤m

1≤k1 <k2 <k3 ≤m

If a set of 2m pairs are in Bk , then we must have that the vertex containing point 2m+2g −1,
denoted by w1 , and the vertex containing point 2m + 2g, denoted by w2 , are both heavy. Let
Y1 and Y2 denote the number of points in i and j respectively that are paired with a point in
a heavy vertex via a single edge in P . Then, |Bk | = Y1 Y2 [di − Wi,j − 1]m−1 [dj − Wj,i − 1]m−1
for every 1 ≤ k ≤ m. Hence,
m
X
|Bk | = mY1 Y2 [di − Wi,j − 1]m−1 [dj − Wj,i − 1]m−1 .
k=1

Clearly, both Y1 and Y2 can be computed in O(∆)-time.
Now for any 2 ≤ ` ≤ m, consider any |Bk1 ∩ Bk2 ∩ · · · ∩ Bk` |. It is easy to see that given
k1 , . . . , k` ,
|Bk1 ∩ Bk2 ∩ · · · ∩ Bk` | = [Y1 ]` [Y2 ]` [di − Wi,j − `]m−` [dj − Wj,i − `]m−` .

Thus,

 
m
|Bk1 ∩ Bk2 ∩ · · · ∩ Bk` | =
[Y1 ]` [Y2 ]` [di − Wi,j − `]m−` [dj − Wj,i − `]m−` .
`
1≤k1 <···<k` ≤m
P
Since Y1 and Y2 have been computed, computing each 1≤k1 <···<k` ≤m |Bk1 ∩ Bk2 ∩ · · · ∩ Bk` |
takes only O(1)-time. Hence, the computation of bi,j (P, m) takes O(∆ + m) steps.
Running time in Phases 3–5
We first discuss the time complexity of Phase 5, which significantly dominates that of
the other two phases.
Again, each step of the algorithm involves information on fτ (P ), bbτ (P, z) and b(P ). However, it is easy to see that it is never necessary to compute fτ (P ), due to the way the mτ (i)
is defined. For instance, mI (i) = 4iM12 which is the total number of ways to pick a double
edge and label its end points, and pick two pairs, with repetition allowed, and label their end
points. If such a choice does not yield a valid switching then an f-rejection is performed. It is
easy to see that the f-rejection is performed with the correct probability, as defined in (21).
Thus, it suffices to compute bbτ (P, z) and b(P ) only.
Consider the computation of b(P ). By (20), b(P ) = M2 L2 − Z ∗ (P ). Note that Z ∗ (P )
is the number of pairs of 2-paths involving at most 5 vertices, or pairs of 2-paths looking
like structures in Figure 8 but involving double edges. It has been shown in the proof of
Lemma 13 that Z ∗ (P ) ≤ 8i(dh M2 +d1 L2 )+(2id21 d2h +4iU12 +8M2 U1 +L4 ). By the assumption
in Lemma 20, we have Z ∗ (P ) = o(M12 ). Hence, computing Z ∗ (P ) takes no more than M12
units of time by brute-force search. So the total computation time for b(P ) is at most M12 .
Next consider the computation of bbτ (P, z) for τ ∈ {III, IV, V, V I, V II}. Each such
X

45

function counts the number of ways to choose a structure corresponding to one of those
in Figure 8, plus several other pairs with restrictions that between certain pairs of vertices
there should be no edge. An important quantity in these bounds will be an upper bound Ij
on the number of j-stars which have a given vertex v as a leaf. These are similar to light
(j − 1)-blooms at v, but without the lightness condition. For j = 2, it also corresponds to a
2-path starting at v. Thus,
2
I2 ≤ U1 = O(n(2γ−3)/(γ−1) ).
(61)

A j-star is P
a set of j pairs all of which include a point in some given vertex u. We clearly
j−1
have Ij ≤ ∆
where the degrees of the graph are d1 ≥ d2 ≥ · · · ; here di represents
i=1 di
the possible degree of vertex u. Hence for any j ≥ 3,
∆
∆
X
X
X
j−1
i−(j−1)/(γ−1) .
(62)
Ij ≤
di ≤
(Kn/i)(j−1)/(γ−1) ≤ (Kn)(j−1)/(γ−1)
i=1

i=1

i≥1

We first focus on τ = III, as a simple illustration of the method. Here, we have bbIII (P, z)
countings the choices for pairs representing H1 in Figure 8 and an extra pair that is not too
close to the H1 . Using brute force, we can find all copies of H1 in time O(L2 I3 ) by breadthfirst search starting at light vertices v (recall Lk below (7)) and investigating all 3-stars
with v as a leaf. For each copy of H1 , we can count the valid choices of the extra pair
u4 v4 in Figure 3 by inclusion-exclusion: make a list of each possible forbidden event (vertex
coincidences, and edges present that are not allowed) and then, for each subset S of the
forbidden events, count that choices for which all events in S hold. For S = ∅, the number
is 2M1 , without any computation. For all other S, the counting is done by breadth-first
search. For the S involving a vertex coincidence, there is a trivial bound of O(∆); the pair
u4 v4 can be chosen in this many ways starting at the specified vertex. Once the subgraph
with the vertex coincidence is constructed, any other events specified in S can be checked
to see if this particular subgraph satisfies them all. Similarly, it is easy to see that all with
an edge coincidence can be bounded above by O(I2 ). Hence O(L2 I3 I2 ) is the complexity for
this case.
Now consider the other types τ . Each of these involves choosing a 2-path with a light
vertex in the middle, as well as two other pairs u1 u2 and u1 u3 , together with a set of up
to seven edges other edges with various forbidden events. We use an inclusion-exclusion
method as for type III.
We treat type VII in detail, as the most complicated case. We can check every occurrence
of H1 as in the type III case, in time O(L2 I3 ), and check each for the existence of the extra
pairs u2 v2 and u3 v3 in constant time. Then take any subset S of all the single events that
must be avoided, and check for each copy of H1 how many choices of the remaining seven
pairs satisfy all events in S.
Note that for S to be satisfied, there must exist a set of edges inducing a graph of a
given isomorphism type, and containing H at a certain location. Pick such an isomorphism
type, W say. The first thing to notice is that if we delete the vertices in H1 , W falls into
components that can be treated separately: for each component we can find the number of
ways to place that component, plus the specified edges in W that join in to H, and multiply
together the counts for the various components. From this it is easy to see that the worst
case is when there is only one such component. (If that is not clear now, it is at the end
46

of this proof.) Consider the graph W0 consisting of W − H1 together with one edge of W
from H1 to W − H1 . (If there is no such edge, then the number of possibilities for W can
be computed in the following manner independently of H1 , and it will be easy to see by the
time we reach the end of this proof that this case is insignificant.) Take a spanning tree T
of W0 . We just need to find all possible attachments of T to the copy of H, as each can be
investigated in constant time for whether the required extra edges are present. Since there
was only one component in W , the seven edges involve at most eight vertices, and so T has
at most eight edges. One leaf, call it w, of T , is a vertex of H0 . We can decompose T , as an
abstract tree, into stars, as follows. If we pull off all leaves of T other than W , the new leaves
are centres of stars. We remove the edges of these stars from T and repeat with the tree
formed from the remaining edges. In this way, T is decomposed into stars that all contain
at least two edges each, except perhaps there might be one star containing a single edge
incident with w. Moreover, each star contains the edge that goes from its centre towards
w. These stars can be placed in the graph, starting with the one containing w, and then
successively building outwards. The number of possibilities for a single edge is at most ∆,
and for a k-star it is at most Ik . The time taken to do this is thus the maximum of R8 and
∆RQ
7 , where Rk is the maximum, over all partitions π1 + · · · + π` of k, with each πi ≥ 2,
of i Iπi . Using the bounds for Ij in (61) and (62) it is easy to see that the maximum is
achieved by n7/(γ−1) , the upper bound for I8 .
The argument is easy to adapt to the other types. For each type, we obtain the following
bounds on the running time in a single switching step. For types IV and VI, M2 L2 is an
upper bound for ui and vi , 1 ≤ i ≤ 3, and the rest of the arguments are the same as above,
adapted to the smaller number of extra edges involved. For convenience, we repeat the
bound obtained above for type III.
III: L2 I3 I2 ,

IV: M2 L2 I4 ,

V: L2 I3 I5 ,

VI: M2 L2 I7 ,

VII: L2 I3 I8 .
Let
denote max0≤i≤i1 ρτ (i). Then the expected total cost for computing bbτ (P, z) during
the algorithm, for a given τ , is bounded by the above quantity times BD ρ∗τ . With the
value of ρ∗τ in Lemma 15, it is easy to see after some computations that the type VII term,
BD M23 M3 L3 I3 I8 /M17 dominates. Recall that BD = 4L2 M2 /M12 from (19). For the range of
γ in Theorem 3 and for sufficiently small δ satisfying (18), this shows the total running time
is bounded above by
ρ∗τ

O(M24 M3 L2 L3 I3 I8 /M19 ) = O(n4.081 )
for the range of γ in Theorem 3 and for sufficiently small δ satisfying (18).
Moreover, since BD ρ∗τ = o(1) for τ ∈ {IV, V, V I, V II} and a.a.s. Phase 5 lasts only
O(BD ) iterations by Lemma 16, the probability that at least one switching of one of these
types is performed in one run of PLD is o(1). Therefore, a.a.s. the computation time of
bbτ (P, z) for τ ∈ {III, IV, V, V I, V II} is bounded by
BD L2 I3 I2 ρ∗III = O(M3 L2 L3 I2 I3 /M13 ) = O(n2.89 ).
This can be improved by the judicious use of data structures. Before running through all
possible H1 , we can compute and record the number of single edges incident with every
vertex, the number of 2-paths starting at any vertex, and the numbers of 2-paths and of
47

3-paths connecting any two vertices. The time taken for this pre-computation is at most the
time taken to generate all 3-paths, which is O(M2 ∆). Multiplying this by BD still gives easily
O(n2 ). Using this information, it is easy to compute, for a given copy of H1 , the number of
choices for the extra edge satisfying a given set of constraints corresponding to S, in constant
time. Then the expected complexity for this case is reduced to O(BD L2 I3 ρ∗III ) = O(n2.107 ).
It is easy to bound the expected computation time in Phase 3 (reduction of loops) using
brute-force searching, by O(n2.087 ). For Phase 4, we can consider for each pair of 2-stars
(L2 M2 possibilities), exactly how many ways to extend it to a pair of 3-stars (ordered of
course), using the same data structures as introduced just above for type III switchings and
a similar inclusion-exclusion scheme. The resulting complexity is O(BT L2 M2 ) = O(n2.107 ).
Hence, we obtain the time complexities claimed in Theorem 3.
We remark that slight improvements in the running time bounds can be gained by making
more use of data structures as we did for type III. In [17], very significant gains were possible
in the case of regular graphs by this means, but in the present case no such large further
gains are easily apparent. So, to keep the argument simple, we don’t pursue this here.
Constraints on (γ, δ)
Various constraints on Mk , Lk and Uk were placed to enforce a small probability of a
rejection during the algorithm. In this section, we list all these constraints and we deduce a
condition on γ under which we may find an appropriate δ so that all these constraints are
satisfied. That completes the definition of L and H and thus completes the definition of the
algorithm.
We always assumed 5/2 < γ < 3; we also assumed 1 − δ > 1/(γ − 1) in (46). Collecting
constraints from Lemma 4, Lemmas 29, 30 and 8, we have the following constraints for δ
and γ.
5/2 < γ < 3
(63)
4
− 2 < δ(γ − 2)
(64)
γ−1
3−γ
< δ
(65)
γ−2
1
1
< δ<
(66)
2γ − 3
2
1
< 1−δ
(67)
γ−1
2
< 1 + δ(γ − 2)
(68)
γ−1
2
δ(γ − 1) >
.
(69)
3
Lemma 12 requires L4 = o(M12 ) and L4 M4 = o(M14 ). Using (45) and (4), these are
implied by
δ(5 − γ) < 1
(70)
4
< 3.
(71)
δ(5 − γ) +
γ−1

48

Recall from (40) that
32M22
.
M13
Lemma 14 requires that δ + 1/(γ − 1) < 1 and ξ = o(1). The first condition is taken into
consideration in (67). The constraint ξ = o(1) requires 4/(γ − 1) < 3 and thus γ > 7/3
which is guaranteed by (63).
Recall from (19) that
4L2 M2
2L3 M3
4L2
BT =
,
BL =
; BD =
2
M1
M1
M13
Lemma 17 requires BD ξ = o(1) which is guaranteed by
6
δ(3 − γ) +
< 4.
(72)
γ−1
Lemma 18 requires BD ξ = o(1), BD (M3 d2h + L3 d21 ) + L6 + M3 U3 + M2 U1 U2 = o(M13 ) and
2
BD n(2γ−3)/(γ−1) −1 = o(1), which are implied by (72) and
5
< 4,
(73)
δ(5 − γ) +
γ−1
4
δ(7 − 2γ) +
< 3,
(74)
γ−1
δ(7 − γ) < 2,
(75)
3
δ(4 − γ) +
< 2,
(76)
γ−1
2
2γ − 3
δ(3 − γ) +
+
< 2.
(77)
γ − 1 (γ − 1)2
Lemma 19 requires U1 M2 L2 /M12 + (M2 L2 /M12 )2 = o(M1 ) which is implied by (77) and
4
< 3.
(78)
2δ(3 − γ) +
γ−1
Lemma 20 requires (dh M2 + d1 L2 + d21 d2h + U12 )M2 L2 /M12 + M2 U1 + L4 = o(M12 ) which is
implied by (70), (71),
2γ − 3
2
+
< 2 (already implied by (77))
γ − 1 (γ − 1)2
and
3
2δ(3 − γ) +
< 2
(79)
γ−1
2
2(2γ − 3)
δ(3 − γ) +
+
< 3.
(80)
γ−1
(γ − 1)2
Lemma 22 requires dh M1 (BD + BT + BL dh ) + L2 (d1 + U1 ) = o(M12 ) which is implied
by (70),
2
δ(4 − γ) +
< 2 (implied by (76))
γ−1
3
δ(5 − γ) +
< 3 (implied by (71))
γ−1
ξ=

49

and

1
< 1
(81)
γ−1
2γ − 3
δ(3 − γ) +
< 1.
(82)
(γ − 1)2
Lemma 24 requires (d1 + U1 )L3 M3 = o(M14 ) and M3 d2h (BD + BT ) + L3 (BD d21 + BT d21 ) +
L6 + M3 U3 + M2 U1 U2 = o(M13 ), which is implied by (73), (75), (77), (76), (74) and
4
δ(4 − γ) +
< 3
(83)
γ−1
2γ − 3
3
+
δ(4 − γ) +
< 3
(84)
γ − 1 (γ − 1)2
6
< 5
(85)
δ(6 − γ) +
γ−1
5
2δ(4 − γ) +
< 4.
(86)
γ−1
By (63), 1/2 < 1/(γ −1) < 1, (2γ −3)/(γ −1)2 > 1/(γ −1) and (2γ −3)/(γ −1)2 < 1, and
so (71) implies (70) and (73), (79) implies (78), (77) implies (80), (82) implies (81) and is
implied by (77), (76) implies (83). Using (67), we have that (71) implies (85). Next we show
that (76) implies (71), (72), (74), (79) and (86). Using γ > 2 we have that (76) implies (79)
and (71) implies (74). By (63), 1/(γ − 1) < 2/3 and so 3/(γ − 1) < 2. Using this together
with γ < 3 we have that (79) implies (72). Using (67) that δ+1/(γ −1) < 1 we have that (76)
implies (71), and (79) implies (86). Thus, we have proved that, using (63) and (67), (76)
implies (71), (72), (74), (79) and (86). Moreover, since 5/2 < γ < 3, it is easy to see that
(2γ −3)/(γ −1)2 < 1 and so (76) implies (84). Since γ < 3, we have 2/(γ −1) > 1 and so (64)
implies (68). Finally, it is easy to see that for all 5/2 ≤ γ < 3, 1/(2γ − 3) > 2/3(γ − 1) and
so (66) implies (69). After removing the redundant constraints (68) – (74) and (78) – (86),
system (63)–(86) reduces to the following
δ(3 − γ) +

which leads to (42).

5/2
4
−2
γ−1
3−γ
γ−2
1
2γ − 3
1
γ−1
3
δ(4 − γ) +
γ−1
δ(7 − γ)
2
2γ − 3
+
δ(3 − γ) +
γ − 1 (γ − 1)2
50

< γ<3
< δ(γ − 2)
< δ
< δ<

1
2

< 1−δ
< 2
< 2
< 2,

