# Convergence criteria for sampling random graphs with specified degree
sequences

Upasana Dutta<EMAIL_ADDRESS>Department of Computer Science,
University of Colorado, Boulder, CO, USA  Aaron Clauset
<EMAIL_ADDRESS>Department of Computer Science, University of
Colorado, Boulder, CO, USA  BioFrontiers Institute, University of Colorado,
Boulder, CO, USA  Santa Fe Institute, Santa Fe, NM, USA

###### Abstract

The configuration model is a standard tool for generating random graphs with a
specified degree sequence, and is often used as a null model to evaluate how
much of an observed network’s structure is explained by its degrees alone.
Except for networks with both self-loops and multi-edges, we lack a direct
sampling algorithm for the configuration model, e.g., for simple graphs. A
Markov chain Monte Carlo (MCMC) algorithm, based on a degree-preserving
double-edge swap, provides an asymptotic solution to sample from the
configuration model without bias. However, accurately detecting convergence of
this Markov chain on its stationary distribution remains an unsolved problem.
Here, we provide a concrete solution to detect convergence and sample from the
configuration model without bias. We first develop an algorithm for estimating
a sufficient gap between sampled MCMC states for them to be effectively
independent. Applying this algorithm to a corpus of 509 empirical networks, we
derive a set of computationally efficient heuristics, based on scaling laws,
for choosing this sampling gap automatically. We then construct a convergence
detection method that applies a Kolmogorov–Smirnov test to sequences of
network assortativity values derived from the Markov chain’s sampled states.
Comparing this test to three generic Markov chain convergence diagnostics, we
find that our method is both more accurate and more efficient at detecting
convergence.

In the analysis and modeling of networks, random graph models are widely used
as both a substrate for numerical experiments and as a null model or reference
distribution to evaluate whether some network statistic is typical or unusual.
In network science, the configuration model is a special kind of random graph
model, which conditions on a specified degree sequence, and allows researchers
to assess and explore the structural or dynamical consequences of a network’s
degree structure Molloy and Reed (1998, 1995); Milo  _et al._ (2004); Bender
and Canfield (1978); Łuczak (1989); Chatterjee  _et al._ (2011); Ring  _et
al._ (2020). As a result, the configuration model is among the most widely
used random graph models across network science, and it provides the basis for
many theoretical results, e.g., the broad implications of heavy-tailed degree
distributions Newman (2002); Newman  _et al._ (2001); Leicht  _et al._ (2006);
Amaral  _et al._ (2000); Volz (2008); Miller and Kiss (2014); Chatterjee  _et
al._ (2009).

The configuration model for networks with self-loops and multi-edges is the
most well known Newman (2010). In fact, however, there are eight different
configuration models, depending on whether the random graph to be generated is
vertex-labeled or stub-labeled—that is, whether or not it matters which “stub”
on a vertex $i$ an edge from vertex $j$ lands on—and whether it is allowed to
have self-loops and/or multi-edges (Fig. 1). The distinction between these
different flavours of the configuration model are not academic: Fosdick et al.
Fosdick  _et al._ (2018) showed that the corresponding distributions of
network statistics can differ so much that an incorrect choice of graph space
can lead to spurious or even opposite conclusions about the significance of
some empirical networks’ observed characteristics.

Figure 1: A small network in each of eight distinct configuration model graph
spaces, corresponding to all combinations of allowing (c-d, g-h) or not
allowing (a-b, e-f) multi-edges, and allowing (e-f, g-h) or not allowing (a-b,
c-d) self-loops, in either a vertex-labeled space (1st and 3rd columns) or a
stub-labeled space (2nd and 4th columns).

Drawing a graph from a configuration model is straightforward if the graph
space is stub-labeled and allows both self-loops and multi-edges. In this
case, a simple stub matching algorithm Newman (2010); Blitzstein and Diaconis
(2011), which chooses a uniform random matching of all the edge stubs of the
network, suffices and can be run in $O(m)$ time, where $m$ is the number of
edges. For the other graph spaces, however, it is common to generate a network
using the stub matching algorithm first, and then simply removing the self-
loops or collapsing the multi-edges in the generated network. In asymptopia,
such self-loops and multi-edges are a vanishing fraction of all edges when the
graph is sparse, as is usually desired. However, these modifications change
the resulting network in highly non-random ways, because high-degree nodes are
more likely to participate in self-loops and multi-edges than are low-degree
nodes. Hence, using the stub matching algorithm in this common way ultimately
violates the underlying assumptions of the configuration model, and produces
non-uniform draws from the target graph space. Or, to put it more bluntly: the
random graph matching algorithm is only applicable to graph spaces where
multi-edges and self-loops are allowed. An alternative solution to drawing
from the other seven graph space, e.g., that of simple graphs, is to sample
them using a Markov chain Monte Carlo (MCMC) algorithm based on double-edge
swaps. Although a number of such algorithms have been defined, the Fosdick et
al. MCMC is known to asymptotically converge on the target uniform
distribution for each of the eight graph spaces (with rare exceptions in graph
spaces that allow loops but not multi-edges Nishimura (2018)). However,
practical guidance on the time required for convergence with finite-sized
networks remains unknown. This is the question we investigate and answer here.

Past work does provide theoretical upper bounds on the convergence or mixing
times for some versions of the double-edge swap MCMC, for very specific types
of degree distributions. For example, the mixing time $\tau(\epsilon)$ for
sampling $d$-regular undirected graphs with fixed degree sequences is
polynomially bounded as $\tau(\epsilon)\leq
d^{15}n^{8}(dn\log(dn)+\log(\epsilon^{-1})))$ Cooper  _et al._ (2007), where
$n$ is the number of vertices and $d$ is the degree. A similar polynomial
bound exists for directed graphs, as $\tau(\epsilon)\leq
50d^{25}n9(dn\log(dn)+\log(\epsilon^{-1}))$ Greenhill (2011). For half regular
bipartite graphs with fixed degree sequence, the mixing time is upper bounded
by a polynomial in $n$ Erdös  _et al._ (2013). The mixing time
$\tau(\epsilon)$ for undirected graphs with a minimum degree of at least one
and a maximum degree satisfying $3\leq d_{\max}\leq 1/3\sqrt{M}$ is
polynomially bounded in $n$ for any degree sequence as
$\tau(\epsilon)\leq\frac{1}{10}d_{\max}^{14}M^{9}(M\log(M)+\log(\epsilon^{-1}))$
Greenhill (2014), where $M$ is the sum of the degrees. None of these results
cover the case of generic degree sequences and hence have limited utility for
working with realistic degree structure.

Because convergence for arbitrary MCMCs cannot be calculated theoretically
Tierney (1994), convergence is instead detected in an online fashion by
examining the local statistics of the sequence of states the MCMC visits. Such
convergence tests have a common structure: (1) a statistic calculated from
each state the MCMC visits, (2) a choice of how often to compute this
statistic, referred to here as the sampling gap $\eta$, and (3) a statistical
test to assess when this calculated statistic has converged to its steady
state or asymptotic distribution. Many such tests have been developed for
general MCMCs Garren and Smith (1993); Dixit and Roy (2017); Johnson (1996);
Heidelberger and Welch (1983); Liu  _et al._ (1992); Mykland  _et al._ (1995);
Ritter and Tanner (1992); Roberts (1992); Yu (1994); Yu and Mykland (1994);
Zellner and Min (1995), and several are included in popular Python and R
packages Gelman  _et al._ (1995); Raftery and Lewis (1995); Geweke (1991).
However, none of these methods is designed for the configuration model and, as
we will show, they do not perform well when applied to it. Refs. Brooks and
Gelman (1998); Roy (2020); Brooks and Roberts (1998) present a detailed
comparative analysis of several MCMC convergence detection techniques, and
highlight their limitations by analyzing their theoretical biases and
practical implementations. Relatedly, Cowles and Carlin Cowles and Carlin
(1996) study 13 different general convergence diagnostics and find that every
method can fail to detect the type of convergence they were designed to
identify.

Here, we develop an efficient and accurate convergence method specifically for
the double-edge swap MCMC for sampling from the configuration model. This
method requires only an input degree sequence and a choice of graph space, and
no other free parameters. First, we develop an algorithm for estimating a
sampling gap between MCMC states so that the sampled states are effectively
independent. We then apply this algorithm to a corpus of 509 real-world and
semi-synthetic networks spanning all eight graph spaces, and distill the
experimental results into a simple set of decision rules, based on scaling
laws, for selecting the sampling gap automatically, along with a set of
conditions under which the scaling laws can be applied. We then specify a test
for detecting MCMC convergence based on the distribution of sampled network
assortativity statistics, and show through a series of experiments that this
test is effective at detecting convergence. Further, we compare this method to
several generic MCMC convergence detection methods and show that our method is
both more accurate and more efficient when applied to real-world networks.
Although our method can be applied to any of the eight graph spaces, it is
most useful for graph spaces other than stub-labeled loopy multigraphs, where
the well-known random stub-matching algorithm provides a computationally
efficient direct solution. Finally, to make the methods described here more
accessible to the community, we provide our own implementations in a Python
package, which can be found here.

##  I Materials and Methods

The Markov chain Monte Carlo algorithms described by Fosdick et al. guarantee
that the resulting stationary distribution is uniform over graphs with the
specified degree sequence and graph space, except for very few loopy graphs
that do not allow multi-edges Nishimura (2018). A graph space is chosen by
specifying whether self-loops are allowed or not, whether multi-edges are
allowed or not, and whether the graph is vertex- or stub-labeled (Fig. 1). In
a vertex-labeled graph, each vertex has a distinct label, whereas in a stub-
labeled graph, each stub has a distinct label and hence each vertex can be
identified by the unique set of stubs attached to it. When a graph is vertex-
labeled, the permutations of the stubs attached to a particular vertex are
indistinguishable, and hence when multi-edges are present, “crossing” two
“parallel” edges does not create another distinct graph. For example, a co-
authorship network can be viewed as a vertex-labelled multigraph because it
would be nonsensical to match author $A$’s first collaboration with author $B$
to author $B$’s second collaboration with author $A$. On the other hand, when
a graph is stub-labeled, crossed edges and parallel edges have different
meanings because the different orderings of the stubs attached to each vertex
imply distinct patterns of interaction in the underlying system. For example,
in a friendship network where each node represents a cohort of students at a
school, and stubs represent students in that cohort, every possible matching
of stubs represents a distinct set of friendships. Hence, such a network would
be stub-labeled. See Ref. Fosdick  _et al._ (2018) for additional examples of
vertex-labeled and stub-labeled networks.

In the MCMC we study here for sampling from the configuration model, each step
of the chain is a degree-preserving double edge swap (Fig. 2a), which
transitions from one graph $G_{t}$ to another graph $G_{t+1}$. At each step in
the chain, we choose two edges $\\{(x,y),(w,z)\\}$ uniformly at random and
replace them with either $\\{(x,z),(w,y)\\}$ or $\\{(x,w),(y,z)\\}$ with equal
probability. Hence, a double-edge swap slightly rewires a graph, while keeping
the degrees of the graph preserved.

Figure 2: (a) The two degree-preserving double-edge swaps on a pair of edges
$\\{(x,y),(w,z)\\}$ results in either
$\\{(x,y),(w,z)\\}\rightarrow\\{(x,z),(w,y)\\}$ or
$\\{(x,y),(w,z)\\}\rightarrow\\{(x,w),(y,z)\\}$ as shown. (b) If the vertices
x, y, w, z are not distinct, double edge swap
$\\{(x,y),(w,x)\\}\rightarrow\\{(x,x),(w,y)\\}$ can introduce a self-loop, (c)
Similarly, if a third edge already exists among the vertices $x,y,w,z$, a
double edge swap $\\{(x,y),(w,z)\\}\rightarrow\\{(x,w),(y,z)\\}$ can introduce
a multi-edge.

In order to sample correctly from the target distribution, the space of graphs
we choose imposes restrictions on which double-edge swaps are permitted. If
the change proposed by a particular swap moves the chain $G_{t}$, $G_{t+1}$,
$G_{t+2}$, $\dots$ outside the specified graph space, e.g., by introducing a
self-loop (Fig. 2b) or a multi-edge (Fig. 2c) forbidden by the given graph
space, or if certain other technical conditions are satisfied (see Ref.
Fosdick  _et al._ (2018)), then the proposed change is rejected and the graph
is resampled. Because the Markov chain traverses a sequence of graphs where
graph $G_{t}$ differs from $G_{t+1}$ by exactly one double edge swap, graphs
close to each other in the sequence are highly serially correlated.
Additionally, if graphs are resampled very often, the serial correlation is
even stronger. This correlation naturally decays between more distantly
sampled points in the chain, but the distance required to achieve a certain
level of correlation must grow with the size of the network, because each
double-edge swap changes only a constant number of edges at a time.

The Fosdick et al. MCMC algorithm is guaranteed to sample graphs with a given
degree sequence uniformly from any of the eight desired graphs spaces (with
very few exceptions in the loopy graph space Nishimura (2018)), after it has
converged to its stationary distribution. Fosdick et al. provide a proof that
this convergence is achieved in the asymptotic limit of $t\to\infty$. However,
in practice, there exists some finite time $t_{*}$ at which the MCMC has
effectively reached this asymptotic state. Our goal is to detect the earliest
time at which this is true, but no earlier.

Here, we develop a complete solution to the practical task of sampling from
the configuration model. Our solution divides this problem into three parts.

  1. 1.

We select a network-level summary statistic that quantifies a sufficiently
non-trivial aspect of a network’s structure, so that we may transform a
sequence of graphs $G_{t}$, $G_{t+1}$, $G_{t+2}$, …from the MCMC into a
standard scalar time series $x_{t}$, $x_{t+1}$, $x_{t+2}$,….

  2. 2.

We develop an algorithm for choosing a “sampling gap” $\eta_{0}$ for the given
network such that the values $x_{t}$ and $x_{t+\eta_{0}}$ from the MCMC are
effectively independent.

  3. 3.

We define a statistical test to compare two consecutive sets of sampled values
$\\{x_{t}\\}_{1}$ and $\\{x_{t}\\}_{2}$ within a sliding window of the MCMC
and decide if they are plausibly drawn from the same underlying (stationary)
distribution Pr($x_{t}$).

Throughout, we make extensive use of a set of real-world and semi-synthetic
networks which span social, biological, technological domains Aaron Clauset
and Sainz (2016). We use these networks both to evaluate the methods we
develop, and to develop a set of computationally lightweight heuristics for
automatically parameterizing our solution. These real-world networks include
103 simple graphs, 51 loopy graphs, 32 multigraphs and 110 loopy multigraphs,
and range in size from n = 16 to 30,269 nodes with a variety of edge densities
and degree distributions. Real networks with self-loops but no multi-edges and
networks with multi-edges but no self-loops are relatively rare compared to
those with either both or neither of them. To obtain a sufficient number of
such networks for our numerical experiments on these graph spaces, we add or
delete self-loops from other networks in our empirical corpus to obtain semi-
synthetic networks (see Sections II.2 and II.4).

##  II Results

To characterize the progression of the MCMC through a graph space, we select
the degree assortativity $r$ of a network as the cognizant network summary
statistic $x_{t}=f(G_{t})$ . The degree assortativity quantifies the tendency
of nodes with similar degrees to be connected, and ranges over the interval
$[-1,1]$. Mathematically, $r$ is calculated as the normalized covariance of
the degrees across all the edges of the network, given by,

$r=\displaystyle\frac{\sum_{xy}(A_{xy}-k_{x}k_{y}/2m)k_{x}k_{y}}{\sum_{xy}(k_{x}\delta(x,y)-k_{x}k_{y}/2m)k_{x}k_{y}},$
(1)

where $A_{xy}$ is the adjacency matrix entry for nodes $x$ and $y$, $k_{x}$ is
the degree of node $x$, $m$ is the number of edges in the network, and
$\delta(x,y)=1$ if $x=y$ and 0 otherwise.

The degree assortativity takes the value $r=1$ if the graph is composed of
only cliques, because in that case, the degree of every node is the same as
that of its neighbors. The degree assortativity takes the value $r=-1$ if the
graph is composed of only equal-sized stars, i.e., trees with exactly one
internal node and $\ell$ leaves $(\ell\geq 2)$, because in that case, the
highest degree nodes of the network connect only to nodes with the lowest
degree. While it is common to expect that $r=0$ is a random graph, the precise
form of the null distribution, and hence its expected value, depends on the
chosen graph space Fosdick  _et al._ (2018).

There are, of course, many alternative network-level summary statistics that
could be used instead of degree assortativity, e.g., the clustering
coefficient $C$, the mean geodesic path length $\langle\ell\rangle$, the mean
betweenness centrality, etc. However, most such statistics are expensive to
calculate repeatedly over the course of the MCMC, which would limit the
scalability of the sampling algorithm. In contrast, the degree assortativity
admits a computationally efficient update equation such that it can be
calculated quickly after every double-edge swap of the Markov chain, which
allows us to sample longer chains and larger networks.

Suppose that a double-edge swap
$\\{(x,y),(w,z)\\}\rightarrow\\{(x,w),(y,z)\\}$ is performed (Fig. 2a). Using
the definition of degree assortativity, its change from this swap can be
written as

$\Delta
r=\displaystyle\frac{(k_{x}k_{w}+k_{y}k_{z}-k_{x}k_{y}-k_{w}k_{z})\times
4m}{(\sum_{i}k_{i}\times\sum_{i}k_{i}^{3})-(\sum_{i}k_{i}^{2})^{2}},$ (2)

where $k_{i}$ is the degree of the $i^{th}$ node and $m$ is the number of
edges in the network.

Given a particular network, the denominator in eq. (2) is a constant, and
hence can be calculated once at a cost of $O(m)$, when the MCMC is first
initialised, and stored for reference later. The numerator only requires the
degrees $k_{x}$, $k_{y}$, $k_{w}$, $k_{z}$ of the four vertices involved in
the swap, and the number of edges $m$. Hence it takes only $O(1)$ constant
time to update $r$, given the assortativity of the current graph $G_{t}$ and
the degrees of the nodes chosen for the double-edge swap. Calculating the
initial assortativity $r_{0}$ is more expensive, but must be done only once.
Its value can be calculated as,

$r_{0}=\displaystyle\frac{S_{1}S_{\ell}-S_{2}^{2}}{S_{1}S_{3}-S_{2}^{2}},$ (3)

where $S_{1}=\sum_{i}k_{x},S_{2}=\sum_{i}k_{x}^{2},S_{3}=\sum_{i}k_{x}^{3}$
and $S_{\ell}=\sum_{xy}A_{xy}k_{x}k_{y}=2\sum_{(x,y)\in E}k_{x}k_{y}$. The
expression in eq. (1) contains only $O(n+m)$ terms and hence is substantially
more efficient than eq. (3), which contains $O(n^{2})$ terms Newman (2010).
Only in case of dense networks, where $m=\Theta(n^{2})$, are the two
calculations equally inefficient.

To illustrate the evolution of the degree assortativity $r$ over the course of
the Markov chain, we apply the Fosdick et al. MCMC to a modest sized stub-
labeled simple graph with $n=16,062$ nodes and $m=25,593$ edges(Fig. 3). In
the early part of the Markov chain, the quick movement away from the initial
value $r_{0}$ is highly apparent, as the double-edge swaps initially randomize
the empirical correlations in each network’s structure. As each chain
progresses, the assortativity converges towards a value, and then displays
non-trivial fluctuations around it. The key problem that convergence detection
seeks to solve is deciding when these excursions are sufficiently random that
we may declare the Markov chain to have reached its stationary distribution.
In making this decision, we note that it is better to detect convergence too
late rather than too early: a late decision merely wastes time in the form of
extra steps in the Markov chain, while an early decision results in sampling
from the wrong distribution of graphs.

Figure 3: Network degree assortativity $r$ as a function of the number of
double-edge swaps, $s$ performed in the Markov chain, for a stub-labeled
simple graph with $n=16,062$ nodes and $m=25,593$ edges

By construction, the states the Markov chain visits are serially correlated,
as each double-edge swap changes at most four adjacencies. The magnitude of
this serial correlation must therefore increase with the size of the network
$m$, because it takes more steps in the Markov chain to randomize a larger
number of edges. As the task of detecting convergence is one of deciding when
the fluctuations in the structure of the states visited by the Markov chain
are indicative of a stationary distribution, large and long-lasting serial
correlations pose a significant problem by creating the appearance of non-
random structure in the chain. A solution to this problem is to only sample
states that are sufficiently well separated in the chain that they are
statistically independent. To identify the spacing between states that yields
such a sample, we develop an algorithm based on the autocorrelation function
and statistical sampling theory, which can be applied to any network for
obtaining an appropriate sampling gap $\eta_{0}$ between states. Given such a
choice of $\eta_{0}$, we expect the sample
$X_{\eta_{0}}=x_{t},x_{t+\eta_{0}},x_{t+2\eta_{0}},...,x_{t+T\eta_{0}}$ to
behave as a set of $T$ independent draws from the MCMC’s stationary
distribution when the MCMC has converged.

The autocorrelation function of a time series measures the pairwise
correlation of values $x_{t}$ and $x_{t+\Delta}$ as a function of the gap or
lag $\Delta$ that separates them (see Appendix A). Our algorithm (Algorithm 1)
uses a test for independence based on the autocorrelation function to
determine whether a sample $X_{\eta}$ from the Markov chain can be considered
to be composed of draws that are independent and identically distributed
(iid).

1:$G_{0}$ (initial graph), $T$ (the size of the list $S_{T}$), $D$ (the number
of parallel MCMC runs), and $\alpha_{c}$ (the significance level).

2:The sampling gap $\eta_{0}$ for graph $G_{0}$.

3:Let m be the number of edges in $G_{0}$.

4:Run the MCMC for $1000m$ swaps (burn-in).

5:Let $x_{1}$ be the first degree assortativity value after the burn-in.

6:$\eta=0$.

7:$f_{\eta}=1$.

8:$Q$ = $\Big{(}100-\frac{(100-C)}{2}\Big{)}\%$ quantile of N(0, 1).

9:while $f_{\eta}>\alpha_{c}$ do

10: $\eta=\eta+1$

11: $f=0$

12: for $k\leftarrow 1$ to $D$ do

13: Construct $S_{T}$ such that $s_{i}\\!=x_{\eta\times(i-1)+1}$, $1\leq i\leq
T$.

14: Compute the autocorrelation function $R_{h}$ for $S_{T}$, $0<h<T$.

15: $f_{k}$ = fraction of autocorrelations beyond
$\pm\displaystyle\frac{Q}{\sqrt{n}}$ in $R_{h}$.

16: $f=f+f_{k}$

17: $f_{\eta}=f/D$

18:$\eta_{0}=\eta$

Algorithm 1 Choosing the sampling gap $\eta_{0}$

To construct this algorithm, we exploit the fact that for a sufficiently large
choice of gap $\eta$, the sampled values in $X_{\eta}$ will behave like an iid
draw from a stationary distribution. Our goal is to find the smallest value
$\eta_{0}$ for which this behavior occurs. When it does, because the sampled
values have finite variance, the autocorrelations of $X_{\eta_{0}}$for any lag
$\Delta>0$ will follow a Normal distribution N($0,1/T$) for large sample size
$T$ Brockwell  _et al._ (1991). Hence, to assess if $X_{\eta}$, $\eta>0$,
comprises of iid draws from a stationary distribution (the null hypothesis),
we can use a standard autocorrelation hypothesis test that assesses the
normality of the autocorrelation values of $X_{\eta}$. This hypothesis test
requires choosing a significance level $\alpha_{c}$ and deriving the
corresponding thresholds outside of which no more than a fraction $\alpha_{c}$
of the autocorrelation values should fall, if the sampled values in $X_{\eta}$
are effectively independent. In our experiments, we choose $\alpha_{c}$=0.05,
and measure the fraction $f_{\eta}$ of autocorrelations of the sample
$X_{\eta}$ that fall outside the corresponding range of $\pm 1.96/\sqrt{T}$
(see Appendix A). Finally, to select the minimum gap $\eta_{0}$, we scan
upward through the values of $\eta>0$ until we fail to reject the null
hypothesis in the autocorrelation test.

To initialise this algorithm, the MCMC is first run for $1000m$ swaps. This
“burn-in period” ensures that in expectation, every edge has been swapped 2000
times, which we take as a reasonable degree of randomization before samples
are taken for any experiment. After this burn-in period, the algorithm creates
a sequence of lists of the form $S_{T}=[s_{1},s_{2},s_{3},\dots,s_{T}]$ where
$s_{i}=x_{\eta(i-1)+1}$ for $1\leq i\leq T$, that is, each $S_{T}$ is a
candidate set of sampled values from the Markov chain, for a particular choice
of gap $\eta$. As $\eta$ increases, the serial correlations in the Markov
chain, and hence also the fraction $f_{\eta}$ of statistically significant
autocorrelations will decay toward zero. Since the Markov chain is stochastic,
we calculate $f_{\eta}$ by averaging over $100$ Markov chains. The first value
of $\eta$ at which $f_{\eta}\leq\alpha_{c}$ is returned as the effective
choice of sampling gap $\eta_{0}$ by which to construct a sample with the
target properties, where $\alpha_{c}$ is the chosen critical threshold.

Unspecified, however, is the number of sampled values $T$ to be used in the
test. Intuitively, the smaller the value of $T$, the faster the sampling, but
the less accurate our estimate of $\eta_{0}$ may be. We conduct a set of
numerical experiments to identify the number of sampled values $T$ that
appears in practice to be sufficient to obtain a stable estimate of
$\eta_{0}$. Specifically, we construct a set of lists $S_{T}$, for
$T=\\{20,30,50,100,150,200,250\\}$ as described above. For example, when
$T=20$ and $\eta=1$, $S_{T}$ is a list of 20 consecutive degree assortativity
values from the Markov chain after the burn-in, and when $T=100$ and
$\eta=2000$, $S_{T}$ contains 100 values sampled every 2000 steps from the
chain after the burn-in, so that the last sample is 200,000 double-edge swaps
after the first sample.

Fig. 4 shows the fraction of significant autocorrelation $f_{\eta}$ as a
function of the sampling gap $\eta$, illustrating how $f_{\eta}$ decays
towards zero with increasing $\eta$. With a fixed sample size $T$, the value
of $\eta$ for which the fraction of significant autocorrelations becomes less
than or equal to 5% is recorded as the sampling gap $\eta_{0(T)}$ for the
network. We observe that the sampling gap $\eta_{0(20)}$ obtained for sample
$S_{20}$ is less than both $\eta_{0(30)}$ and $\eta_{0(50)}$, obtained for
samples $S_{30}$ and $S_{50}$, respectively. Hence, the sampling gap
$\eta_{0(20)}$ is large enough by our criteria, for a sample containing $20$
values, but not large enough when the sample contains $50$ values. To resolve
such a discrepancy, we take the largest length of sub-samples so that adding
more values to it does not induce a significant change in the $\eta_{0}$ value
for this network. Fig. 4 shows that while the sample size is less than $50$,
increasing it also increases the corresponding value $\eta_{0}$. However, Fig.
4 inset shows that this pattern disappears once the sample size exceeds $50$,
meaning the value of $\eta_{0}$ stabilizes for sufficiently large sub-samples.
Hence, to obtain a reliably stable estimate of $\eta_{0}$, in our sampling gap
experiments, we fully automate the selection of $\eta_{0}$ by running the
algorithm for $T=200$ and $T=250$ and taking the mean of the two corresponding
estimates $\eta_{0(200)}$ and $\eta_{0(250)}$ returned by the algorithm.

The autocorrelation test described above also makes an assumption of
stationarity of the test sample $X_{\eta}$ being tested, and may not perform
as expected if this behavior is uncommon in practice. Using our corpus of 509
empirical and semi-synthetic networks, we evaluated this assumption Dickey and
Fuller (1979) on all samples taken for the experiments after the “burn-in”,
and found that for 92% of the samples tested across all values of sampling gap
$\eta$, the stationarity assumption is met.

Figure 4: Fraction of significant autocorrelations $f_{\eta}$ as a function of
sampling gap $\eta$ shown for different list sizes for a network with $n=1589$
nodes and $m=2742$ edges. The inset is a zoomed-in view of the outer plot
between sampling gaps $850$ and $1000$, with the same axes as the outer plot.

Although for any particular network, one could run the above sampling gap
estimation algorithm to choose an appropriate value for $\eta_{0}$, the
procedure itself is computationally expensive. To avoid this cost, we now
develop a set of heuristics and decision criteria by which to quickly and
automatically choose $\eta_{0}$, given only the fixed degree sequence and the
choice of the graph space to sample from.

###  II.1 Simple graphs

For simple networks, the Markov chain’s transition probabilities are the same
for both stub-labeled and vertex-labeled spaces (see Appendix B.1). Hence, the
sampling gaps that our algorithm estimates for particular real-world simple
networks will be the same, in either the vertex-labeled and the stub-labeled
spaces.

We begin by running the sampling gap algorithm on each of the 103 simple
networks in our corpus. Of these, the vast majority (93.2%) exhibit a nearly
linear pattern in how the estimated sampling gap $\eta_{0}$ varies with the
size of the network m (Fig. 5a,e). This pattern is very close to the simple
form of $\eta_{0}=m/3+300$. Although there are a few networks for whom this
heuristic would overestimate the sampling gap, it is generally better to
choose a sampling gap $\eta_{0}$ too large—sampling a pair of networks in the
Markov chain that are even more uncorrelated—than to choose it too small and
risk non-trivial serial correlations. However, for a small number of networks
(7), the estimated sampling gap $\eta_{0}$ is far larger than this simple
heuristic predicts. This behavior occurs mainly for graph spaces that prohibit
multi-edges. In that case, the network’s density plays a crucial role in
governing how often the change proposed by a double-edge swap is rejected, and
hence a graph is resampled by the MCMC.

The reason the MCMC produces these unusually large $\eta_{0}$ estimates stems
from the Markov chain’s boundary enforcement criterion: if the Markov chain
proposes a double-edge swap in a simple graph that would create a self-loop or
a multi-edge, thereby exiting the specified graph space, the algorithm rejects
this change and resamples its current state. Each such resampling extends the
time needed for the serial correlation to decay.

Let $\rho$ is the iid probability of an edge existing between a uniformly
random pair of nodes, i.e., the network density. (We note that for simple
graphs $\rho=\langle k\rangle/n-1$ and for loopy graphs $\rho=\langle
k\rangle/n$, where $\langle k\rangle$ is the mean degree). Suppose two edges
$(x,y)$ and $(w,z)$ are chosen uniformly at random for a double-edge swap.
Each of the two kinds of swaps (Fig. 2a) is then chosen with equal
probability. In both the cases, the probability of rejecting the swap
$\\{(x,y),(w,z)\\}\rightarrow\\{(x,z),(w,y)\\}$ or the swap
$\\{(x,y),(w,z)\\}\rightarrow\\{(x,w),(y,z)\\}$ due to the creation of a
multi-edge after the swap is governed by the likelihood that an edge already
exists where the swap proposes to place one:

Pr(rejection) $\displaystyle=\textrm{Pr\big{(}multi-
edge}\big{\rvert}\\{(x,y),(w,z)\\}\rightarrow\\{(x,z),(w,y)\\}\big{)}$
$\displaystyle\hskip 12.80365pt+\textrm{Pr\big{(}multi-
edge}\big{\rvert}\\{(x,y),(w,z)\\}\rightarrow\\{(x,w),(y,z)\\}\big{)}$
$\displaystyle=\displaystyle\frac{1}{2}\textrm{Pr\Big{(}at least one of
}(x,z)\textrm{ and }(w,y)\textrm{ exists}\Big{)}$ $\displaystyle\hskip
12.80365pt+\displaystyle\frac{1}{2}\textrm{ Pr\Big{(}at least one of
}(x,w)\textrm{ and }(y,z)\textrm{ exists}\Big{)}$
$\displaystyle=\displaystyle\frac{1}{2}\Big{(}1-\textrm{Pr\big{(}neither
}(x,z)\textrm{ nor }(w,y)\textrm{ exists}\big{)}\Big{)}$ $\displaystyle\hskip
12.80365pt+\displaystyle\frac{1}{2}\Big{(}1-\textrm{Pr\big{(}neither
}(x,w)\textrm{ nor }(y,z)\textrm{ exists}\big{)}\Big{)}$
$\displaystyle=\displaystyle\frac{1}{2}\Big{(}1-(1-\rho)^{2}\Big{)}+\displaystyle\frac{1}{2}\Big{(}1-(1-\rho)^{2}\Big{)}$
$\displaystyle=2\rho-\rho^{2},$

where each swap can produce a multi-edge with either of the two swapped edges.
We refer to the simplified expression of $2\rho-\rho^{2}=\omega$ as the
density factor of a network, which gives a simple density-based estimate of
Pr(rejection). A more precise estimate would exploit the moment structure of
the degree distribution, but such a formula is not necessary for our purposes.

Hence, in a space that does not allow multi-edges, higher network density
substantially increases the probability that the MCMC will resample states.
And, resampling states necessarily increases the estimated sampling gap
because resampling preserves the serial correlations in the chain that we seek
to avoid. In our corpus of 103 simple networks, the vast majority (86.4%) have
a density factor $\omega\leq 0.25$, six (5.8%) have one $0.25<\omega\leq 0.5$,
five (4.9%) have one $0.5<\omega\leq 0.75$ and the remaining four (3.9%) have
$\omega>0.75$. These frequencies reflect the fact that real-world networks
with very high density typically occur only rarely Broido and Clauset (2019),
unless the network is especially small.

Figure 5: Space specific estimated sampling gap of networks that passed all
the criteria of our proposed heuristics for (a) stub-labelled simple graphs,
(b) stub-labeled loopy graphs, (c) stub-labeled loopy multigraphs, (d) stub-
labeled multigraphs, (e) vertex-labeled simple graphs, (f) vertex-labeled
loopy graphs, (g) vertex-labeled loopy multigraphs, (h) vertex-labeled
multigraphs.

Across our simple network corpus, we find that networks with a density factor
$\omega\leq 0.25$ yield sampling gap estimates that scale linearly with the
number of edges m (Fig. 5 a,e), while networks with a larger density factor
are more likely to yield anomalously high sampling gap estimates. We use this
observation to divide networks into two categories, based on their calculated
density factors $\omega$. If a network’s density factor is $\omega>0.25$, the
sampling gap should be estimated directly via the algorithm described above
(Algorithm 1); otherwise, the sampling gap follows a simple pattern and may be
chosen automatically via a scaling law that fits the estimated gap values
(Fig. 5a,e). For any particular network, choosing too small a value of
$\eta_{0}$ would cause the Markov chain to mix more poorly, and so we manually
fit a scaling law of the form $\eta_{0}=m/3+300$ to the empirical data, which
provides an upper bound on nearly all of the estimated gap values. For the
seven (7.8%) networks for which this function is not an upper bound, the error
is modest: 1.9% and 1.1% in the stub and the vertex-labeled spaces,
respectively. We note that if the graph is small, e.g., $m<1000$, and has a
relatively high density factor $\omega>0.25$, then the scaling law still
yields a sufficiently high value of $\eta_{0}$, i.e., one comparable to what
the sampling gap estimation algorithm typically produces.

###  II.2 Loopy graphs

Real-world “loopy” networks, those with self-loops but no multi-edges, are
uncommon. We construct a reasonable empirical corpus of 154 loopy graphs by
obtaining 51 such networks from public repositories Aaron Clauset and Sainz
(2016) and generating an additional 103 loopy graphs by adding self-loops to
our previous corpus of simple networks. To convert a simple network into a
loopy network, we first measured the fraction of nodes with a self-loop in
both our corpus of 110 real-world loopy multigraphs and the 51 real-world
loopy graphs to obtain an empirical distribution of loopiness. For each simple
graph in our corpus, we then chose a fraction iid from this loopiness
distribution and added a self-loop to each node with the corresponding
probability. We note that this loopiness distribution is strongly bimodal: 44%
have at most 5% of nodes with self-loops, while 39.1% have at least 95% of
nodes with self-loops.

Although the MCMC algorithm proposed by Fosdick et al. correctly samples
graphs with fixed degree sequence post convergence, there are certain
conditions Nishimura (2018) that the degree sequence must satisfy for this to
hold true in a loopy graph space (i.e graph space that allows loops but no
multi-edges). These conditions ensure that the double-edge swap algorithm
reaches every valid loopy graph with the given degree sequence. Nearly all
empirical loopy degree sequences satisfy these conditions, and in our corpus
of loopy graphs, all the 154 networks’ degree sequences do so. However, for
the rare empirical degree sequences that do not, a different MCMC sampling
algorithm exists, which augments double-edge swaps with triangle-loop swaps
for connecting the space of loopy graphs Nishimura (2018). The solution
proposed in this paper only serves those loopy graphs that satisfy the
necessary conditions for the Fosdick et al. MCMC. Designing a convergence
detection method for the MCMC sampler that uniformly samples loopy graphs with
any degree sequence Nishimura (2018) is an interesting direction of future
work.

For our loopy network corpus, we repeat our analysis by first running the
sampling gap algorithm on each network in the corpus. We note that because
multi-edges are not permitted in the loopy space, the density factor
$\omega=2\rho-\rho^{2}$ has the same relevance here as it does with simple
graphs. In our loopy corpus, 90.9% of networks have a density factor
$\omega\leq 0.25$, and the same scaling law of $\eta_{0}=m/3+300$ is an upper
bound in both stub- and vertex-labeled spaces (see Appendix B.2) on nearly all
of the observed sampling gap values for these networks (Fig. 5b,f). Of the 154
loopy networks, 10 (6.5%) and 18 (11.7%) networks in the stub and vertex-
labeled spaces, respectively, have low density factors but the scaling law is
not an upper bound for them. In these cases, the error from using the scaling
law to choose $\eta_{0}$ remains modest, being 2.3% and 1.1%, respectively.

###  II.3 Loopy multigraphs

For stub-labeled loopy multigraphs, it is common to directly construct
networks with a fixed degree sequence by choosing a uniformly random matching
on the set of edge “stubs” given by that sequence Newman (2010); Blitzstein
and Diaconis (2011). Furthermore, this algorithm is computationally cheap,
taking $O(m)$ time, compared to an MCMC approach. However, stub matching
cannot be used to correctly sample from the vertex-labeled loopy multigraph
space. There, we must instead use the Fosdick et al. MCMC, with a convergence
criterion. For completeness, we analyze the MCMC’s behavior in both the stub-
labeled and vertex-labeled loopy multigraph settings. In the former setting,
our results may be useful in cases where correlated draws are desired.

In the loopy multigraph space, the Markov chain’s transition probabilities in
the stub-labeled and vertex-labeled spaces are different. In the stub-labeled
case, the Markov chain never resamples a state, while in the vertex-labeled
space, some transitions are forbidden and hence the Markov chain resamples a
state under certain conditions (see Appendix B.3). Due to this resampling, the
sampling gap required in the vertex-labeled space should be greater than the
gap in the stub-labeled space.

We investigate these differences by running the sampling gap algorithm on our
corpus of 110 real-world loppy multigraphs (Fig. 5c,g). We find that, indeed,
for a given network size, the estimated sampling gap tends to be higher in the
vertex-labeled case than in the stub-labeled case. However, for 42 of 110
(38%) loopy multigraphs in the vertex-labeled case, the estimated sampling gap
is much higher than the simple linear pattern observed in the majority of
cases. This peculiar behaviour of the MCMC for these networks can be
understood by considering how a network’s degree distribution influences the
MCMC’s dynamics.

Suppose that two edges of the network $(x,y)$ and $(w,z)$ share exactly one
common endpoint, implying that only 3 of the 4 vertices are distinct. If a
double edge swap proposes the change
$\\{(x,y),(w,z)\\}\rightarrow\\{(x,w),(y,z)\\}$, then half the time the degree
assortativity will remain unchanged after the swap, i.e., $\Delta r=0$.
Because these edges are chosen uniformly, the higher the degree of a node, the
greater the probability of that node appearing twice in the set
$\\{x,y,w,z\\}$. The presence of such high-degree nodes in the network
increases the likelihood of MCMC steps with $\Delta r=0$. The more of these
$\Delta r=0$ steps that occur, the greater the serial correlation in the
Markov chain, and the greater the gap required to obtain independent samples.
In the vertex-labeled space, in addition to choosing edges having one common
endpoint, the Fosdick et al. MCMC resamples states for other reasons Fosdick
_et al._ (2018), which also creates a $\Delta r=0$ step, and therefore further
increases the serial correlations in the Markov chain.

Chung-Lu random graphs are a vertex-labeled model of simple graphs that places
a constraint on the value of the maximum degree. This constraint serves to
limit the likelihood of producing a network with multi-edges, and has a form
$\max_{i}k_{i}^{2}\leq 2m$ Chung and Lu (2002). In our setting, this
constraint allows for a simple heuristic for deciding how skewed a degree
distribution can be before state resamplings are likely to generate too-strong
serial correlations in the vertex-labeled space. In our sampling gap
experiments, we find that every loopy multigraph in our corpus that exceeds
this maximum degree threshold does indeed require an unusually large sampling
gap (even larger than $m$) in the vertex-labeled space. In contrast, networks
in our corpus that fall below $\frac{1}{3}$ of this threshold (a safe margin)
produce estimated sampling gaps that exhibit the same familiar linear
relationship with $m$. In this case, a scaling law of $\eta_{0}=m/2.7+100$
(constructed in the same way as for previous graph spaces) produces a
conservative upper bound on $\eta_{0}$ (Fig. 5g). And, if a network exceeds
this maximum degree safe margin, the sampling gap of the network in the
vertex-labeled space should instead be chosen using our sampling gap algorithm
(Algorithm 1).

In the stub-labeled case, the corresponding scaling law obtained from all the
110 loopy multigraphs is $\eta_{0}=m/3+300$ (Fig. 5c)111The panels Fig. 5c and
Fig. 5g show the estimated sampling gaps for 110 and 68 networks in the stub-
labeled and vertex-labeled spaces, respectively. In the latter, we exclude the
42 networks whose maximum degree exceeds the threshold given in the text.. For
the nine (8.1%) loopy multigraphs for which the scaling law is not an upper
bound in the stub-labeled space, the resulting error is 1.9%. And, for the
nine (13.2%) loopy multigraphs that satisfy our maximum degree constraint and
yet have an estimated sampling gap higher than our proposed scaling law for
the vertex-labeled space, the error is 14.1%.

###  II.4 Multigraphs

Real world networks with multi-edges but no self-loops are also uncommon. We
construct an empirical corpus of 142 such networks by obtaining 32 of them
from public repositories and generating 110 multigraph networks by removing
the self-loops of our loopy multigraphs. We selected this approach for
generating semi-synthetic networks because the set of possible graphs with and
without various self-loops is exponentially smaller, and thus easier to sample
from, than the set of graphs with or without various multi-edges.

Figure 6: Decision tree of the space specific sampling gaps given by our
heuristic, depending on the satisfiability of the constraints.

As with loopy multigraphs, the network’s maximum degree and whether the space
is stub-labeled or vertex-labeled govern the choice of sampling gap $\eta_{0}$
(see Appendix B.4), just as it does in the loopy multigraph space. Repeating
our analysis on this corpus of 142 multigraphs, we find that the estimated
sampling gaps of all 142 networks in stub-labeled space and of 78 (54.9%)
networks in the vertex-labeled space exhibit a familiar linear pattern in the
number of edges $m$ (Fig. 5d,h)222As in Fig. 5g, the estimated sampling gaps
shown in for networks whose maximum degree exceeds the threshold given in the
text are excluded from Fig. 5h. The scaling laws of $\eta_{0}=m/3+300$ in the
stub-labeled space, and $\eta_{0}=m/2.7+100$ in the vertex-labeled space for
networks whose maximum degree does not exceed the safe margin of
$\max_{i}k_{i}^{2}\leq\frac{1}{3}2m$, provide an upper bound on the estimated
sampling gaps for nearly all multigraphs in the corpus. For the 26 (18.3%)
loopy multigraphs where this is not true in the stub-labeled space, the
resulting error is 3.9%; and, for the seven (8.9%) loopy multigraphs that
satisfy our maximum degree constraint and yet the scaling law is not an upper
bound in the vertex-labeled space, the error is 11.2%.

###  II.5 A complete decision tree for choosing the sampling gap

Although an appropriate sampling gap $\eta_{0}$ can always be estimated using
Algorithm 1, our numerical experiments show that in many cases, a sufficient
gap may be chosen more efficiently using a simple scaling law in the number of
edges in the network (Fig. 5). However, there are some conditions on applying
these scaling laws in practice, depending on the particular graph space and
certain structural properties of the network, and in some cases, it may still
be necessary to choose $\eta_{0}$ via Algorithm 1.

The decision tree in Fig. 6 organizes the insights, conditions, and scaling
laws obtained from our numerical experiments into a single heuristic for
choosing the sampling gap, depending on the network’s particular properties
and specified graph space. In cases where a network does not fall into a
category in which a scaling law can be used to choose $\eta_{0}$, the serial
correlations in the Markov chain will tend to far exceed those predicted by
any of our scaling laws, and we recommend instead directly estimating a
network-specific sampling gap via Algorithm 1. Again, we note that the
heuristics developed in this paper are derived from the double-edge swap MCMC
Fosdick  _et al._ (2018) and do not apply to an extremely rare class of loopy
graphs for which the necessary conditions of the graph space connectivity
Nishimura (2018) are not satisfied. In those rare cases, one would have to use
the MCMC proposed by Nishimura et. al. Nishimura (2018).

We note that the appearance of scaling laws in this setting, and the
consistency of their form across different graph spaces is intriguing. This
pattern may reflect a currently unknown but common underlying topology both
within and across these different graph spaces. Investigating the origins of
this common pattern may yield deeper theoretical insights into or guarantees
for MCMC convergence, and is an interesting direction for future work.

From a more practical perspective, this common pattern provides a simple and
highly effective heuristic that works well across a wide variety of graph
types.

###  II.6 Convergence criterion

Given the choice of sampling gap $\eta_{0}$ (Fig. 6), we now specify a
convergence test for the Fosdick et al. MCMC. The test we propose compares the
states of two consecutive sequences of states from the Markov chain, and
assesses whether they are plausibly drawn from a common underlying
distribution, and hence are evidence that the Markov chain has converged on
its equilibrium.

The test samples states from the Markov chain to populate a window of two
lists, each composed of 30 degree assortativity values, denoted as
$L_{1}=\\{r_{t},r_{t+\eta_{0}},r_{t+2\eta_{0}},\dots r_{t+29\eta_{0}}\\}$ and
$L_{2}=\\{r_{t+30\eta_{0}},r_{t+31\eta_{0}},r_{t+32\eta_{0}},\dots
r_{t+59\eta_{0}}\\}$. We then apply the non-parametric two-sample
Kolmogorov–Smirnov (KS) test Massey Jr (1951) to the lists $L_{1}$ and $L_{2}$
to determine whether they are plausibly drawn from the same distribution.

This choice of 30 observations in a list is a trade off: it is large enough to
provide a reasonable degree of statistical power to the KS test, but small
enough to be computationally expedient. In some applications, a larger list
may be desired, but we caution against a shorter list, which may increase the
risk of incorrect early convergence due to low statistical power. Similarly,
there are a variety of other tests that could be used to assess the
statistical similarity of $L_{1}$ and $L_{2}$. For this task, the KS test has
three attractive properties: (i) it is a test of distributional similarity and
hence is more sensitive than comparing summary statistics like the mean or
variance, (ii) it is computationally light-weight, and (iii) it is non-
parametric, making it appropriate for Markov chains with unknown
distributional characteristics.

Figure 7: Network degree assortativity $r$ as a function of the number of
double-edge swaps, $s$ performed in the Markov chain, for (a) a stub-labeled
simple graph with $n=2\,642$ nodes and $m=3\,303$ edges, and (b) a vertex-
labeled loopy multigraph graph with $n=2\,837$ nodes and $m=11,407$ edges. The
vertical red line indicates the point at which our diagnostic detects
convergence.

If the KS test rejects the null hypothesis that the two samples are drawn from
the same underlying distribution, we discard the contents of the first list,
“slide” the window forward by a list-length (30 values), sample new values
from the Markov chain in order to populate the rest of the window, and finally
test again. We repeat this process of slide, repopulate, and test, until the
KS test fails to reject the null hypothesis, which we interpret as an
indicator that the MCMC has converged to its stationary distribution. Every
$\eta_{0}$ steps beyond this point provides an iid draw from this
configuration model.

To illustrate this procedure, Fig. 7 shows the convergence detected by our
this test for a simple network in the stub-labeled space with $n=2\,642$ nodes
and $m=3\,303$ edges, and for a loopy multigraph in the vertex-labeled space
with $n=2\,837$ nodes and $m=11,407$ edges.

###  II.7 Validating the diagnostic

We now turn to the question of whether our convergence test can correctly and
reliably detect when the MCMC has reached its stationary distribution. If the
procedure works as desired, then before convergence, the distribution of
assortativity values will be non-stationary, while after convergence is
detected, the distribution should appear stationary, meaning that running the
Markov chain longer should not alter the sampled assortativity distribution
(Fig. 8).

To assess this behavior, we tabulate the assortativity distributions of $500$
states sampled after $m/8$, after $m/4$, and after $m/2$ steps, and then again
at convergence and after $1000m$ steps (well beyond convergence). The latter
distribution represents the target distribution of the MCMC, and hence
provides a clear comparison against distributions sampled earlier in the
Markov chain. As a first test, we apply this assessment to a simple stub-
labeled network with $n=1589$ nodes and $m=2742$ edges.

Fig. 8 shows these five sampled distributions, along with the assortativity
coefficient for the original empirical network, which is the initial condition
of the Markov chain. In the pre-convergence phase, as the MCMC walk lengthens
from $m/8$, to $m/4$, to $m/2$ steps, the assortativity distribution for the
sampled networks moves progressively further away from the empirical value of
the initial network. This non-stationarity is indicative of a steady
decorrelation of the Markov chain relative to its initial state, as
progressively more double-edge swaps are applied to randomize the network’s
structure. Once convergence is detected, this non-stationary behavior is no
longer present, and instead the difference between the assortativity
distribution at convergence, and well beyond it, is negligible. The
stationarity of the two latter distributions indicates that the test of
convergence correctly detected the Markov chain’s convergence on the target,
uniform distribution. Finally, repeating this same experiment on networks
drawn from our corpuses for other graph spaces produces similar results.

Figure 8: Null distributions obtained after increasing numbers of double edge
swaps are applied to a simple network with $n=1589$ nodes and $m=2742$ edges
in the stub-labeled space, validating the correctness of our convergence
detection method. We find similar results for the other graph spaces.

###  II.8 Comparing convergence tests

Detecting convergence in a Markov chain Monte Carlo algorithm is a common and
non-trivial problem, especially for scalar time series. Although the
underlying states in the double-edge swap Markov chain are networks, our
approach for detecting convergence converts this graph sequence into a
sequence of scalar values. Hence, it is reasonable to wonder how existing
scalar-valued techniques perform in this setting.

To explore this question, we compare the method described above with three
commonly used scalar time series convergence tests: (i) the Geweke diagnostic
Geweke (1991), (ii) the Gelman-Rubin diagnostic Gelman  _et al._ (1995) and
(iii) the Raftery-Lewis diagnostic Raftery and Lewis (1995). (Each of these
methods is available via libraries in Python Fonnesbeck  and R Plummer  _et
al._ ).

A crucial consideration when selecting a convergence test is whether the
Markov chain in question satisfies the test’s underlying requirements and
assumptions Cowles and Carlin (1996). Each of these three alternative tests
make assumptions that may not hold for the configuration model. For instance,
the Geweke diagnostic assumes that the Geweke statistic derived from the
sampled states will be distributed as a standard normal variable in the
asymptotic limit (see Appendix C.1). The Gelman-Rubin diagnostic assumes that
the MCMC’s stationary distribution is normally distributed, and requires
multiple Markov chains to be initialized at highly dispersed initial states in
the sample space (see Appendix C.2). And finally, the Raftery-Lewis diagnostic
depends on the quantile of interest, and in some cases, the estimated
convergence rate may fall far below the rate required of the full chain Brooks
and Roberts (1999) (see Appendix C.3).

The Geweke and the Raftery-Lewis method assess one list at a time for
detecting convergence, while our method and the Gelman-Rubin method compares
two consecutive lists from the Markov chain to conduct its test for
convergence (see Appendix C). However, none of the other three diagnostics
specify their sampling list sizes, and that choice of parameter is left to the
user. In contrast, our method described above has a fixed list size of
$29\eta_{0}+1$ steps, which yields exactly 30 sampled values, each $\eta_{0}$
double-edge swaps apart.

We evaluate the performance of the four convergence tests along two
dimensions: accuracy and efficiency. First, we say that a test is more
accurate if the assortativity distribution at the time that the test detects
convergence is closer to the corresponding distribution at $1000m$ steps into
the Markov chain (a point sufficiently deep in the chain that we assume it
represents the converged distribution). Second, we say that a test is more
efficient if it uses fewer steps to detect convergence. An ideal convergence
test will perform well on both accuracy and efficiency measures. An efficient
but inaccurate test would use relatively little data, but would tend to
declare convergence too early, producing a distribution of assortativity
values that differs substantially from the target distribution. In contrast,
an inefficient but accurate test would run a very long Markov chain and
declare convergence long after the stationary distribution had been reached.

As a common basis for evaluating the four tests, we select $50$ loopy
multigraphs uniformly at random from our empirical corpus. Because the three
alternative tests do not specify a sampling list size, we evaluate them under
two choices, which allows us to also assess whether a test’s performance is
sensitive to that choice. In the first, we fix their sampling list size at an
arbitrary $m/2$ steps, while our method uses its fixed value of
$29\eta_{0}+1$. In the second, we fix all sampling list sizes to be that of
our method. In both experiments, we apply the four tests to each of the
networks in our test corpus. To measure a test’s accuracy, we perform a
standard two-sample KS-test (with a significance level of $\alpha=0.05$)
between the assortativity distribution at the time of convergence detected by
a test and the distribution obtained after $1000m$ swaps, and record the
corresponding KS distance as our measure of accuracy. To measure a test’s
efficiency, we record the number of steps run in the Markov chain before the
test detects convergence.

Fig. 9a,b show the measured KS distance (accuracy) vs. the mean number of
double-edge swaps $\bar{s}$ used to detect convergence (efficiency). In the
first experiment (Fig. 9a), our approach exhibits a very high accuracy, in the
form of a low rate of early convergence ($2\%$ false alarms). The Geweke
diagnostic also performs well ($6\%$ rate), and both our method’s rates and
the Geweke diagnostic’s rates are statistically indistinguishable from the
$5\%$ rate expected from our choice of $\alpha$ (test of proportions, $p=0.41$
and $p=0.82$, respectively). In contrast, the Gelman-Rubin diagnostic exhibits
a modest and significant error rate ($18\%$, $p=0.04$), while the Raftery-
Lewis diagnostic performs very poorly ($94\%$, $p<10^{-5}$). We find that our
method’s accuracy is statistically indistinguishable from that of the Geweke
diagnostic ($p=0.30$), but significantly better than the Gelman-Rubin
diagnostic ($p<0.01$) and the Raftery-Lewis diagnostic ($p<10^{-5}$).

Figure 9: KS distance (accuracy) as a function of the mean number of double-
edge swaps (efficiency) applied on 50 stub-labeled loopy multigraphs, for the
four diagnostics. (a) Convergence tests every $m/2$ swaps for the three
alternative tests, and (b) convergence tests every $29\eta_{0}+1$ swaps for
all tests.

In the second experiment, where all tests use the same sampling list size
(Fig. 9b), our approach again exhibits a very high accuracy ($0\%$ false
alarms). The Geweke and Gelman-Rubin diagnostics also exhibit high accuracies
($2\%$ and $8\%$ false alarms, respectively), and all three methods’ rates are
indistinguishable from the expected $5\%$ false alarm rate based on the chosen
$\alpha$ (test of proportions, $p=0.10$, $p=0.41$, and $p=0.54$,
respectively). The Raftery-Lewis diagnostic again performs poorly ($98\%$,
$p<10^{-5}$). We find that in this setting our method’s accuracy is
statistically indistinguishable from that of the Geweke diagnostic ($p=0.31$),
but significantly better than the Gelman-Rubin diagnostic ($p=0.04$) and the
Raftery-Lewis diagnostic ($p<10^{-5}$). Comparing error rates across
experiments, we find that each alternative convergence test has a similar
error rate in both experiments. This consistency suggests that the particular
choice of sampling list size makes little difference in the resulting accuracy
of these three methods. In terms of accuracy, the Geweke diagnostic performs
similarly to our method, and both methods are highly accurate. However, the
Geweke diagnostic can be substantially less efficient (Fig. 9b), letting the
Markov chain run for many more steps than is necessary, often producing late
alarms, compared with our method.

Across both experiments, our assortativity-based test exhibits high accuracy
and efficiency, while other techniques tend to perform either slightly or
dramatically worse on one or both dimensions. Of the three, the Geweke
diagnostic performs most similarly, exhibiting equivalent accuracy, but can
also be substantially less efficient (Fig. 9b). An interesting direction of
future work would be to understand why, and whether a more efficient test than
ours can be constructed. Finally, we find that repeating these experiments on
the other graph spaces yields similar results.

##  III Discussion

The configuration model is among the most widely used models of random graphs
in the analysis and modeling of networks. For instance, it is both a standard
way to generate synthetic networks with a specified degree structure, e.g.
with a heavy tail or from an empirical network, and a null model by which to
evaluate whether some observed network statistic is typical or unusual, given
the degree structure. The most common method for sampling configuration model
networks—the popular stub-matching algorithm—correctly generates only stub-
labeled loopy multigraphs, and heuristics for converting these networks into
other types, e.g., simple graphs, can introduce structural artifacts that
contaminate their use Fosdick  _et al._ (2018). A Markov chain Monte Carlo
algorithm propsed by Fosdick et al. Fosdick  _et al._ (2018) based on a
degree-preserving double-edge swap can correctly sample from the configuration
model after its convergence, in any of the either graph spaces, with very rare
exceptions in the loopy graph space Nishimura (2018). However, the proof of
correctness for the Fosdick et al. MCMC does not provide criteria for
identifying when the Markov chain has converged on the target distribution.

The method we developed here provides a complete solution to the convergence
detection problem of the Fosdick et al. MCMC, for all eight graph spaces,
defined by whether the target graph is stub-labeled or vertex-labeled, and
whether it allows self-loops or not, and multi-edges or not. This solution
provides an efficient and accurate means to detect its convergence on the
target, uniform distribution over the specified graphs.

Our approach transforms the sequence of graphs sampled by the Markov chain
into a scalar-valued sequence of degree assortativity values for the graphs
(Fig. 3), and then uses a hypothesis test to detect when this time series
appears to have converged (Fig. 7). There are, of course, many alternative
network summary statistics that one could use instead of degree assortativity,
such as the clustering coefficient or average path length. For this task, the
global assortativity measure $r$ is more computationally efficient to
calculate on a potentially long sequence of networks than many alternatives,
and it provides a good measure of the overall “shape” of a state in the Markov
chain.

The strong serial correlations in the Markov chain, however, creates problems
both for drawing a graph from the configuration model, after convergence, and
for using hypothesis tests to detect convergence, which typically assume that
samples are iid. To mitigate this issue, we develop a novel algorithm for
estimating a sampling gap $\eta_{0}$ (Algorithm 1) by which to obtain
effectively uncorrelated draws from the MCMC. This algorithm is based on a
standard autocorrelation test to determine how far apart two sampled states
must be in order to be statistically independent. Here again, other choices
are certainly possible for how one might quantify the degree of
autocorrelation between two sampled states. We leave an exploration of these
possibilities, which may yield more efficient gap estimation algorithms, for
future work. Applying this gap estimation algorithm to a large corpus of $509$
real world and semi-synthetic networks, we identified a set of simple decision
rules based on the empirical scaling behavior of estimated gaps $\eta_{0}$
with the number of edges $m$. These rules allow a researcher to automatically
select an appropriate sampling gap for a given network usually without having
to run the gap estimation algorithm (Algorithm 1).

Our analysis of the estimated sampling gaps revealed several interesting
patterns, along with useful insights as to when the double-edge swap MCMC
tends to misbehave. For instance, in multigraphs, when the maximum degree is
particularly large, the MCMC tends to resample states at higher rates in the
vertex-labeled space, producing more serial correlations and hence delayed
convergence. A similar issue occurs, but for a different reason, when the
overall network density is especially high in networks where multi-edges are
not allowed. Additionally, across all eight graph spaces, we find that the
estimated sampling gaps $\eta_{0}$ of nearly all networks scale very nearly
linearly with the number of edges in the network $m$ (Fig. 5). We have no
explanation for this appearance of scaling laws here, or their common form
across different graphs spaces. However, these findings suggest that the
overall topologies of the different graph spaces may be fundamentally similar
in some unknown, but likely important, way. Further investigation of these
underlying similarities may lead to more reliable convergence criteria.

In cases where the conditions of the heuristics (Fig. 6) are not satisfied,
the sampling gap algorithm (Algorithm 1) can be used to estimate a bespoke
sampling gap $\eta_{0}$ for the particular network, albeit at a non-trivial
computational cost.

Given a set of effectively independent samples from the Markov chain, many
different statistical convergence tests could be applied. The non-parametric
two-sample Kolmogorov–Smirnov test we employ here has the advantage of being a
test of distributional differences, which provides greater sensitivity than
tests that merely compare summary statistics like the mean value. Moreover,
when applied to a set of values sampled according to the estimated gap, this
test can accurately detect convergence (Fig. 8). Compared to three common
tests of convergence for MCMCs, our approach is both more accurate and more
efficient—converging neither too early nor too late—than three commonly used
MCMC convergence tests (Fig. 9).

The convergence detection method we introduce thus resolves the open question
of how to accurately sample from the configuration model. Our study may also
provide insights on developing convergence detection methods for other Markov
chain sampling algorithms for networks Nishimura (2018), and we note that our
approach can likely be adapted to other variations of the configuration model,
e.g., on graphs of fixed core-value sequence Van Koevering  _et al._ (2021).

The methods developed here are only applicable to the classic configuration
model on dyadic networks, i.e., networks where edges are defined as pairs of
nodes. As such, different methods may be needed for correctly sampling from
either the hypergraph configuration model Chodrow (2020), in which edges are
polyadic, or the configuration models for simplicial complexes Young  _et al._
(2017), or for networks where edge weights cannot be interpreted as multi-
edges, although it can be adapted to work on directed networks. Finally, the
degree-preserving double-edge swap approach employed by the Fosdick et al.
MCMC cannot be applied to the loopy networks that do not satisfy the necessary
conditions for the loopy graph space to be connected Nishimura (2018).
However, such loopy networks occur extremely rarely.

In closing, we note that the lack of accurate and efficient convergence tests
for drawing graphs from the configuration model has encouraged researchers to
often use the fast random stub-matching algorithm first and then convert that
loopy-multigraphs into whatever graph type is desired, e.g., a simple network.
Fosdick et al. showed that this approach is risky when one desires a precise
comparison between an empirical network and the distribution defined by the
configuration model, as the conversion process can introduce problematic
statistical artifacts. The methods described here provide a complete solution
to sampling from the degree-preserving configuration model using the Fosdick
et al. MCMC algorithm, and we hope they will enable researchers to more
accurately and more broadly compare the rich variety of empirical patterns
found in networks with those of random graphs with a specified degree
structure.

###### Acknowledgements.

The authors thank Dan Larremore, Joel Nishimura, Vanja Dukic and Brian
Zaharatos for helpful conversations, and acknowledge the BioFrontiers
Computing Core at the University of Colorado Boulder for providing High
Performance Computing resources (NIH 1S10OD012300) supported by BioFrontiers
IT. The authors thank Joel Nishimura for providing the code for checking if
the loopy graph space is connected or not. This work was funded in part by Air
Force Office of Scientific Research Award FA9550-19-1-0329 (UD, AC).

## Author contributions

AC conceived the research, UD and AC designed the analyses, and wrote the
manuscript. U.D. conducted the analyses. Competing interests: None.

## References

  * Molloy and Reed (1998) M. Molloy and B. Reed, Combinatorics Probability and Computing 7, 295 (1998).
  * Molloy and Reed (1995) M. Molloy and B. Reed, Random Structures & Algorithms 6, 161 (1995).
  * Milo  _et al._ (2004) R. Milo, N. Kashtan, S. Itzkovitz, M. E. Newman, and U. Alon, arXiv preprint cond-mat/0312028 (2004).
  * Bender and Canfield (1978) E. A. Bender and E. R. Canfield, Journal of Combinatorial Theory 24, 296 (1978).
  * Łuczak (1989) T. Łuczak, in  _Proceedings of the Symposium on Random Graphs, Poznań_ (1989) pp. 165–182.
  * Chatterjee  _et al._ (2011) S. Chatterjee, P. Diaconis, and A. Sly, Annals of Applied Probability 21, 1400 (2011).
  * Ring  _et al._ (2020) J. H. Ring, J.-G. Young, and L. Hébert-Dufresne, in  _International Conference on Network Science_ (2020) pp. 33–47.
  * Newman (2002) M. E. Newman, Physical Review E 66, 016128 (2002).
  * Newman  _et al._ (2001) M. E. Newman, S. H. Strogatz, and D. J. Watts, Physical Review E 64, 026118 (2001).
  * Leicht  _et al._ (2006) E. A. Leicht, P. Holme, and M. E. Newman, Physical Review E 73, 026120 (2006).
  * Amaral  _et al._ (2000) L. A. N. Amaral, A. Scala, M. Barthelemy, and H. E. Stanley, Proceedings of the National Academy of Sciences 97, 11149 (2000).
  * Volz (2008) E. Volz, Journal of Mathematical Biology 56, 293 (2008).
  * Miller and Kiss (2014) J. C. Miller and I. Z. Kiss, Mathematical Modelling of Natural Phenomena 9, 4 (2014).
  * Chatterjee  _et al._ (2009) S. Chatterjee, R. Durrett, _et al._ , The Annals of Probability 37, 2332 (2009).
  * Newman (2010) M. Newman,  _Networks: An Introduction_ (Oxford university press, 2010).
  * Fosdick  _et al._ (2018) B. K. Fosdick, D. B. Larremore, J. Nishimura, and J. Ugander, SIAM Review 60, 315 (2018).
  * Blitzstein and Diaconis (2011) J. Blitzstein and P. Diaconis, Internet Mathematics 6, 489 (2011).
  * Nishimura (2018) J. Nishimura, Journal of Complex Networks 6, 927 (2018).
  * Cooper  _et al._ (2007) C. Cooper, M. Dyer, and C. Greenhill, Combinatorics, Probability and Computing 16, 557 (2007).
  * Greenhill (2011) C. Greenhill, The Electronic Journal of Combinatorics 18, P234 (2011).
  * Erdös  _et al._ (2013) P. L. Erdös, I. Miklós, and L. Soukup, The Electronic Journal of Combinatorics 20, P16 (2013).
  * Greenhill (2014) C. Greenhill, in  _Proceedings of the twenty-sixth annual ACM-SIAM symposium on discrete algorithms_ (2014) pp. 1564–1572.
  * Tierney (1994) L. Tierney, The Annals of Statistics 22, 1701 (1994).
  * Garren and Smith (1993) S. Garren and R. L. Smith,  _Convergence Diagnostics for Markov Chain Samplers_ , Tech. Rep. (University of North Carolina, 1993) technical Report.
  * Dixit and Roy (2017) A. Dixit and V. Roy, Journal of Statistical Computation and Simulation 87, 2622 (2017).
  * Johnson (1996) V. E. Johnson, Journal of the American Statistical Association 91, 154 (1996).
  * Heidelberger and Welch (1983) P. Heidelberger and P. D. Welch, Operations Research 31, 1109 (1983).
  * Liu  _et al._ (1992) C. Liu, J. Liu, and D. B. Rubin, in  _Proceedings of the American Statistical Association, Statistical Computing Section_ (1992) pp. 74–78.
  * Mykland  _et al._ (1995) P. Mykland, L. Tierney, and B. Yu, Journal of the American Statistical Association 90, 233 (1995).
  * Ritter and Tanner (1992) C. Ritter and M. A. Tanner, Journal of the American Statistical Association 87, 861 (1992).
  * Roberts (1992) G. O. Roberts, Bayesian Statistics 4, 777 (1992).
  * Yu (1994) B. Yu, The Annals of Probability 22, 94 (1994).
  * Yu and Mykland (1994) B. Yu and P. Mykland,  _Looking at Markov samplers through CUMSUM path plots: a simple diagnostic idea_ , Tech. Rep. (Department of Statistics, University of California at Berkeley, 1994) technical Report 413.
  * Zellner and Min (1995) A. Zellner and C.-K. Min, Journal of the American Statistical Association 90, 921 (1995).
  * Gelman  _et al._ (1995) A. Gelman, J. B. Carlin, H. S. Stern, D. B. Dunson, A. Vehtari, and D. B. Rubin,  _Bayesian data analysis_ (CRC press, 1995).
  * Raftery and Lewis (1995) A. E. Raftery and S. M. Lewis, Practical Markov Chain Monte Carlo 7, 763 (1995).
  * Geweke (1991) J. Geweke,  _Evaluating the accuracy of sampling-based approaches to the calculation of posterior moments_ , Vol. 196 (Federal Reserve Bank of Minneapolis, 1991).
  * Brooks and Gelman (1998) S. Brooks and A. Gelman, Computing Science and Statistics , 30 (1998).
  * Roy (2020) V. Roy, Annual Review of Statistics and Its Application 7, 387 (2020).
  * Brooks and Roberts (1998) S. P. Brooks and G. O. Roberts, Statistics and Computing 8, 319 (1998).
  * Cowles and Carlin (1996) M. K. Cowles and B. P. Carlin, Journal of the American Statistical Association 91, 883 (1996).
  * Aaron Clauset and Sainz (2016) E. T. Aaron Clauset and M. Sainz, “The Colorado Index of Complex Networks,”  (2016).
  * Brockwell  _et al._ (1991) P. J. Brockwell, R. A. Davis, and S. E. Fienberg,  _Time series: Theory and Methods_ (Springer Science & Business Media, 1991).
  * Dickey and Fuller (1979) D. A. Dickey and W. A. Fuller, Journal of the American Statistical Association 74, 427 (1979).
  * Broido and Clauset (2019) A. D. Broido and A. Clauset, Nature Communications 10, 1 (2019).
  * Chung and Lu (2002) F. Chung and L. Lu, Annals of Combinatorics 6, 125 (2002).
  * Massey Jr (1951) F. J. Massey Jr, Journal of the American Statistical Association 46, 68 (1951).
  * (48) C. J. Fonnesbeck, “PyMC 2.3.8 documentation Modelchecking and diagnostics,” .
  * (49) M. Plummer _et al._ , “Package coda,” .
  * Brooks and Roberts (1999) S. P. Brooks and G. O. Roberts, Biometrika 86, 710 (1999).
  * Van Koevering  _et al._ (2021) K. Van Koevering, A. R. Benson, and J. Kleinberg, arXiv preprint (2021).
  * Chodrow (2020) P. S. Chodrow, Journal of Complex Networks 8, cnaa018 (2020).
  * Young  _et al._ (2017) J.-G. Young, G. Petri, F. Vaccarino, and A. Patania, Physical Review E 96, 032312 (2017).
  * Box  _et al._ (1976) G. E. Box, G. M. Jenkins, and G. C. Reinsel,  _Time series analysis: forecasting and control_ , Vol. 734 (John Wiley & Sons, 1976).

##  Appendix A The Autocorrelation function

The autocorrelation function of a sequence measures the degree to which its
values are serially correlated, so that the greater the serial correlation,
the larger the value of the autocorrelation function Box  _et al._ (1976).
Mathematically, the autocorrelation $R_{h}$ quantifies the average correlation
between a pair of values $x_{t}$ and $x_{t+h}$ in the sequence, separated by a
lag of $h$, and ranges from $-1\leq R_{h}\leq 1$. It is defined as
$R_{h}=\displaystyle\frac{C_{h}}{C_{0}}$, where
$C_{h}=\displaystyle\frac{1}{T}\displaystyle\sum_{t=1}^{T-h}{(x_{t}-\mkern
1.5mu\overline{\mkern-1.5mux\mkern-1.5mu}\mkern 1.5mu)}{(x_{t+h}-\mkern
1.5mu\overline{\mkern-1.5mux\mkern-1.5mu}\mkern 1.5mu)},$ and
$C_{0}=\displaystyle\frac{\displaystyle\sum_{t=1}^{T}(x_{t}-\mkern
1.5mu\overline{\mkern-1.5mux\mkern-1.5mu}\mkern 1.5mu)^{2}}{T}.$

Hence, the autocorrelation $R_{h}$ is the covariance of the sample and itself
at a lag of $h$, normalized by the variance of the sample. Here $T$ is the
maximum separation over which we consider the autocorrelation of the pairs of
states.

It can be shown (see Brockwell and Davis Brockwell  _et al._ (1991) p.222)
that if a sample of $n$ points comprises of only independent and identically
distributed (iid) noise with finite variance, the sample autocorrelations
$R_{h}$ of this sample, for $h>0$, are approximately iid N($0,1/n$) for large
$n$. Hence, $C\%$ of the sample autocorrelations will lie between the bounds
$\pm Q/\sqrt{n}$, where $Q$ is the $(100-\frac{(100-C)}{2})\%$ quantile of the
standard normal distribution.

##  Appendix B Dependence of sampling gap on the MCMC transition probabilities

For reference, we have presented the stub-labeled and the vertex-labeled
double-edge swap MCMC algorithm in Algorithm 2 and 3 adapted from Fosdick et.
al Fosdick  _et al._ (2018). In this appendix, we provide a discussion of how
the MCMC’s transition probabilities in different graph spaces govern the
respective sampling gaps we obtain using Algorithm 1.

1:initial graph $G_{0}$, graph space (simple, multigraph, or loopy multigraph)

2:sequence of graphs $G_{i}$

3:for $i<$ number of graphs to sample do

4: choose two edges at random

5: randomly choose one of the two possible swaps

6: if edge swap would leave graph space then

7: resample current graph: $G_{i}\leftarrow G_{i-1}$

8: else

9: swap the chosen edges, producing $G_{i}$

Algorithm 2 stub-labeled MCMC

1:initial graph $G_{0}$, graph space (simple graph, multigraph, or loopy
multigraph)

2:sequence of graphs $G_{i}$

3:for $i<$ number of graphs to sample do

4: choose two distinct edges $(u,v)$ and $(x,y)$ uniformly at random

5: if Unif$(0,1)<0.5$ then

6: $u,v\leftarrow v,u$

7: if edge swap would leave the graph space then

8: resample current graph: $G_{i}\leftarrow G_{i-1}$

9: if $\exists$ 4 distinct vertices in $u,v,x,y$ then

10: $A\leftarrow w_{uv}w_{xy}$

11: $B\leftarrow(w_{ux}+1)(w_{vy}+1)$

12: else if $\exists$ 3 distinct vertices in $u,v,x,y$ then

13: if $u=v$ or $x=y$ then

14: $A\leftarrow 2w_{uv}w_{xy}$

15: $B\leftarrow(w_{ux}+1)(w_{vy}+1)$

16: else

17: $A\leftarrow w_{uv}w_{xy}$

18: $B\leftarrow 2(w_{ux}+1)(w_{vy}+1)$

19: else if $\exists$ 2 distinct vertices in $u,v,x,y$ then

20: if only one of $(u,v)$ or $(x,y)$ is a self-loop then

21: $G_{i}\leftarrow G_{i-1}$

22: continue

23: else if both $(u,v)$ and $(x,y)$ are self-loops then

24: $A\leftarrow 2w_{uu}w_{xx}$

25: $B\leftarrow\frac{1}{2}(w_{ux}+2)(w_{ux}+1)$

26: else

27: $A\leftarrow\frac{1}{2}w_{uv}(w_{uv}-1)$

28: $B\leftarrow 2(w_{uu}+1)(w_{vv}+1)$

29: else

30: $G_{i}\leftarrow G_{i-1}$

31: continue

32: $P\leftarrow\min\Big{(}1,\displaystyle\frac{B}{A}\Big{)}$

33: if Unif$(0,1)<P$ then

34: swap $(u,v),(x,y)\leadsto(u,x),(v,y)$ to produce $G_{i}$

35: else

36: $G_{i}\leftarrow G_{i-1}$

Algorithm 3 vertex-labeled MCMC

###  B.1 Simple graphs

In the stub-labeled space, if a double-edge swap would cause the Markov chain
to leave the graph space, e.g., by creating a self-loop or a multi-edge, the
swap is rejected and the current graph $G_{t-1}$ is resampled by the MCMC
(Algorithm 2). Otherwise, the swap is accepted, producing the new graph
$G_{t}$. In the vertex-labeled space, if a double-edge swap would cause the
Markov chain to leave the graph space, the current graph $G_{t-1}$ is also
resampled (Algorithm 3). Otherwise, the algorithm checks the number of
distinct vertices and self-loops in the chosen double edges. In simple graphs,
the number of distinct vertices among $u,v,x$ and $y$ would either be 3 or 4,
and the multiplicity $w_{ab}$ of any edge $(a,b)$ would be 1, because multi-
edges are forbidden in the simple graph space. Hence, in Algorithm 3, the
variable $B/A$ will be greater than or equal to 1, and the variable $P$ in the
algorithm will always be 1. As a result, if a double-edge swap does not leave
the graph space in the vertex-labeled graph space, it will always be accepted,
producing the new graph $G_{t}$. Hence, for simple graphs, the transition
probabilities are the same for both the stub-labeled and the vertex-labeled
space, and thus so are the induced null distributions too.

###  B.2 Loopy graphs

In the stub-labeled space, if a double-edge swap introduces a multi-edge in
the graph (which is forbidden), the current graph $G_{t-1}$ is resampled.
Otherwise, the swap is accepted, producing the new graph $G_{t}$. In the
vertex-labeled space, if a double-edge swap would cause the Markov chain to
leave the graph space, $G_{t-1}$ is resampled, otherwise, Algorithm 3 performs
the following check. If none of the edges chosen for the swap are self-loops,
then $P=1$ and the swap is accepted. But, if either of the chosen edges is a
self-loop, then in Algorithm 3, we set $A=2$ and $B\in{1,2,4}$. When $B=1$,
then $P=1/2$, and the swap will be rejected half the time. For all other
values of B, we have $P=1$, and hence the swap is accepted. In our corpus of
161 multigraphs and loopy multigraphs, 86 (53%) of these networks have less
than 5% of edges as self-loops, and 129 (80%) of them have less than 20% of
edges as self-loops, Because of the relative scarcity of self-loops, the
chance of a double-edge swap resulting in $P=1$, when the swap does not leave
the graph space, is high for most networks. For this reason, the Markov chains
for loopy graphs in the stub-labeled and vertex-labeled spaces are very
similar, but slightly different in some cases.

###  B.3 Loopy multigraphs

In the stub-labeled space, every double-edge swap is allowed, and hence never
rejected. Hence, the graph $G_{t-1}$ is never resampled. However, in the
vertex-labeled space, whether the proposed double-edge swap is accepted or
rejected depends wholly on the multiplicity of the edges and the number of
distinct vertices chosen for the double-edge swap. In Algorithm 3, there can
be several ways in which the value of the variable $P$ is far below 1 and
hence the graph $G_{t-1}$ will be resampled. For this reason, the MCMCs for
stub-labeled and vertex-labeled loopy multigraphs behave very differently,
with a substantial number of resamplings done in the vertex-labeled space, and
none in the stub-labeled space.

###  B.4 Multigraphs

In both the stub-labeled and the vertex-labeled spaces, a graph is resampled
only if a double-edge swap would introduce a self-loop in the network.
However, in the vertex-labeled space, resampling can occur for an additional
reason (Algorithm 3), depending on the multiplicity of the edges chosen for
the double edge swap. In this case, the Markov chain in the stub-labeled and
the vertex-labeled spaces have very different transition probabilities, which
leads to substantially different sampling gaps between the two spaces.

##  Appendix C Commonly known MCMC convergence tests

PyMC Fonnesbeck  is a Python package for Bayesian statistical modeling and
advanced MCMC algorithms, and it provides implementations of three commonly
used convergence tests for generic MCMCs the Geweke diagnostic Geweke (1991),
the Gelman-Rubin diagnostic Gelman  _et al._ (1995), and the Raftery-Lewis
diagnostic Raftery and Lewis (1995).

###  C.1 The Geweke diagnostic

The Geweke convergence diagnostic Geweke (1991) compares the mean and the
variance of samples from the beginning and the end of a single chain of a MCMC
walk. Geweke Geweke (1991) suggests that the beginning section be the first
10% of the chain and the ending section be the last 50%. The test uses the
Geweke statistic, defined as the difference between the means of the two
samples divided by the standard error. These statistics are estimated using
the spectral densities of the two samples evaluated at zero, which takes into
account the autocorrelations in the samples. The method uses the departure of
the statistic from the standard normal assumption as an indicator of
convergence failure. However, the test statistic is calculated under the
assumption that when the MCMC reaches its stationary distribution, the two
chain samples will be distributed according to a standard normal, in the
asymptotic limit. The test also assumes that the spectral density of the time
series has no discontinuities at frequency zero Cowles and Carlin (1996). Of
course, whether these assumptions are satisfied depends on the process that
produces the time series . The Geweke diagnostic is known to be very sensitive
to the chosen spectral window Cowles and Carlin (1996). Hence, the behavior of
the Geweke test will likely depend on the degree to which the properties of
the stationary distribution of the particular MCMC are known.

###  C.2 The Gelman-Rubin diagnostic

The Gelman-Rubin diagnostic Gelman  _et al._ (1995) compares the within-
sequence variability and between-sequence variability of multiple (at least
two) sequences obtained from MCMC walks with starting points sampled from an
overdispersed distribution. This test is based on the idea that when an MCMC
has not yet converged, the variance within each chain is much less than that
between the chains, because prior to convergence, the MCMC samples states non-
uniformly. If $\theta$ denotes the sequence of an MCMC chain, then first, an
estimate of the marginal posterior of Var($\theta$), is calculated as a
function of the within-sequence variance $(W)$ and between-sequence variance
$(B)$ of the multiple sequences. If the MCMC walks have not yet converged, $B$
would overestimate Var($\theta$) because the walks’ starting values were
chosen to be overdispersed, whereas $W$ would underestimate Var($\theta$)
because the MCMC walks have not yet saturated the states in the stationary
distribution. The Gelman-Rubin statistic is then given by the square root of
the ratio of Var($\theta$) and $W$. In the stationary distribution of the
MCMC, both Var($\theta$) and $W$ should approach the true variance of the MCMC
chain, and hence values of the Gelman-Rubin statistic close to 1 indicate
convergence. In general, practitioners often use a cutoff of 1.1 as an
indicator of convergence. This convergence test assumes that the stationary
distribution of the MCMC is normal. Cowles and Carlin Cowles and Carlin (1996)
suggest that both the assumption of a normal approximation of the target
distribution, and the requirement of multiple MCMC chains with highly
dispersed initial conditions, may not be reasonable in most practical
situations.

###  C.3 The Raftery-Lewis diagnostic

The Raftery-Lewis diagnostic Raftery and Lewis (1995) is based on two-state
Markov chain theory and standard sample size formulae for binomial variance.
In particular, it calculates the burn-in period of the Markov chain and the
total number of subsequent iterations required to accurately estimate $u$, a
$q$-quantile of the MCMC’s posterior distribution. The value of $q$, the
margin of error $r$, and the probability $s$ of obtaining the estimate in the
interval $(q-r,q+r)$, are all user-defined parameters, although the default
value of $s$ is 95%. The method also calculates a thinning interval $k$, which
is the number of iterations that should be skipped to produce a chain of
independent samples from the Markov chain (analogous to $\eta_{0}$ in our
method).

From the Markov chain $\\{\theta_{t}\\}$, the Raftery-Lewis method first
constructs a 0-1 binary chain $\\{Z_{t}\\}$, and then chooses the thinning
interval $k$ to be the smallest natural number for which a first-order Markov
chain model of the thinned out chain is statistically preferred over the
second-order Markov chain model. For convergence detection purposes, we are
interested only in the burn-in period’s value, which provides an estimate of
the number of steps needed before the MCMC reaches its stationary distribution
Brooks and Roberts (1998). A detailed analysis of the Raftery-Lewis method by
Brooks and Roberts Brooks and Roberts (1999) shows that the method depends
strongly on the quantile of interest $q$ and it does not provide information
about the chain as a whole Brooks and Roberts (1999). Furthemore, in certain
cases, the convergence rate estimated by this method is far below the
convergence rate of the full chain Brooks and Roberts (1999). As a general
rule, they find that the routine value of $q=0.025$ commonly used in the
literature should not be used, as it tends to underestimate the true
convergence time. Instead, they suggest that for practical purposes, the
diagnostic could be applied to several different $q$ values and then choose
the quantile that estimates the largest burning length.