# State-Dependent Kernel Selection for Conditional Sampling of Graphs

James Scott
Department of Mathematics
Imperial College London
London SW7 2AZ
<EMAIL_ADDRESS>
&Axel Gandy
Department of Mathematics
Imperial College London
London SW7 2AZ
<EMAIL_ADDRESS>

###### Abstract

This paper introduces new efficient algorithms for two problems: sampling
conditional on vertex degrees in unweighted graphs, and sampling conditional
on vertex strengths in weighted graphs. The algorithms can sample conditional
on the presence or absence of an arbitrary number of edges. The resulting
conditional distributions provide the basis for exact tests. Existing samplers
based on MCMC or sequential importance sampling are generally not scalable;
their efficiency degrades in sparse graphs. MCMC methods usually require
explicit computation of a Markov basis to navigate the complex state space;
this is computationally intensive even for small graphs. We use state-
dependent kernel selection to develop new MCMC samplers. These do not require
a Markov basis, and are efficient both in sparse and dense graphs. The key
idea is to intelligently select a Markov kernel on the basis of the current
state of the chain. We apply our methods to testing hypotheses on a real
network and contingency table. The algorithms appear orders of magnitude more
efficient than existing methods in the test cases considered.

_K_ eywords Contingency Table, Degree Sequence, Exact Test, Markov Chain Monte
Carlo, Random Network

##  1 Introduction

Inference on graphs conditional on vertex level data arises in many diverse
disciplines including network science, psychometrics, community ecology and
categorical data analysis. Two difficult, yet important problems involve
sampling from the set of graphs (weighted graphs) conditional on prescribed
vertex degrees (strengths). In many settings, researchers additionally need to
condition on an arbitrary set of known edges/non-edges in the graph.

The resulting conditional distributions are used to perform hypothesis tests.
Approximate tests rely on distributional approximations to the null. These
approximations can be arbitrarily poor in sparse graphs, and so alternative
approximations, based on sampling, become particularly important.
Unfortunately, most proposed samplers are inefficient in sparse graphs.

Additionally, if the null distribution conditions on known edges/non-edges, it
is difficult to construct a connected Markov chain on the relevant state
space. Existing methods either specialise to particular patterns of fixed
edges, or in the general case, use techniques from computational algebra to
compute a Markov basis. These methods are computationally intensive and are
impractical for graphs with more than a few vertices.

We propose a new class of MCMC methods that use state-dependent mixing of
Markov kernels. This technique allows us to construct samplers that require
little tuning to the problem at hand, and do not require computation of a
Markov basis. The samplers are irreducible in the face of arbitrary patterns
of fixed edges/non-edges. They are efficient both in extremely sparse and
dense graphs, and appear orders of magnitude more efficient than existing
methods in the test cases we consider.

The first focus of this paper is on uniformly sampling unweighted graphs
conditional on prescribed vertex degrees. In the context of hypothesis
testing, the vertex degrees are often sufficient statistics for nuisance
parameters in the null distribution; in such applications, conditioning allows
researchers to perform exact tests.

Consider social network analysis. A social network equipped with a dichotomous
relation can be expressed as a simple digraph. Vertices represent actors, with
edges representing the applicability of the relation between actors.
Frequently researchers are interested in testing the presence of reciprocity
in the network; defined loosely as a preference for mutual dyads in the
corresponding digraph.

Holland and Leinhardt (1981) introduce an exponential family model under which
the UMPU test for reciprocity conditions on the observed degree sequences. In
this case, conditioning serves to remove unknown parameters from the null, and
the resulting distribution is then uniform on the reference set.

The complex interactions that result from conditioning render analytic
analysis of the null distribution difficult or impossible. Efforts have been
made to develop recursive formulas to enumerate all graphs in the reference
set (see Wasserman and Faust, 1994), however these are impractical for even
moderately sized graphs.

If we can sample graphs (nearly) uniformly, then we can approximate the null
distribution of an arbitrary test statistic. Thus, the literature has focused
on simulation, whose methods can broadly be divided into two camps; Markov
Chain Monte Carlo (MCMC) (Rao et al., 1996; Roberts, 2000; Milo et al., 2002;
McDonald et al., 2007; Verhelst, 2008) and sequential importance sampling
(SIS) (Snijders, 1991; Zhang and Chen, 2013; Chen et al., 2005; Bayati et al.,
2010).

A digraph can be represented by its adjacency matrix. Thus, the sampling
problem is equivalent to uniformly sampling zero-one tables with given margins
and a structurally zero diagonal. Sampling zero-one tables with given margins
is applied in community ecology to test for patterns in co-occurrence tables,
and in psychometrics to test the Rasch hypothesis (see Gustafsson, 1980).
Thus, there exists a substantial parallel literature along these lines.

Most MCMC algorithms proposed for sampling graphs are adaptations of methods
proposed for zero-one tables. Typically, they use a combination of ‘switch’
moves (Ryser, 1963) and additional moves to maintain irreducibility in the
face of structural zeros. Rao et al. (1996) and McDonald et al. (2007)
consider ‘compact alternating hexagon’ and ‘hexad’ updates respectively. Most
proposed methods suffer from poor mixing in unbalanced matrices, rendering
them impractical for moderate to large graphs. Additionally, they are not
extensible to arbitrary patterns of fixed edges and non-edges.

SIS builds the graph sequentially, at each iteration choosing a candidate edge
with probability proportional to the vertex degrees. Early methods for this
application include (Snijders, 1991; Chen et al., 2005). Most of these
samplers get stuck, and the probability of restarting approaches 1 as the
degree sequences grow. Bezáková et al. (2012) provide examples where such
algorithms are slow. More recent methods avoid the issue of restarting and
often come with better theoretical guarantees (Bayati et al., 2010; Blitzstein
and Diaconis, 2011; Zhang and Chen, 2013).

The second focus of this paper is on sampling integer-weighted graphs given
vertex strengths. This can be used to conduct network tomography in the case
of a star network topology. However, the motivating application is
approximating the null distribution for evaluating exact tests on two-way
contingency tables. This is a classical problem in statistics which is
important because standard asymptotics justifying approximate tests (notably
Pearson’s $\chi^{2}$ test of independence) do not hold for tables with cells
with low expected frequencies (see Agresti, 2003).

In conditional tests of independence one is interested in the hypergeometric
distribution on all tables with given margins. This corresponds to sampling
integer-weighted bipartite graphs conditional on vertex strengths. Diaconis
and Sturmfels (1998) proposed a simple ‘switch’ Markov chain to sample from
such tables. We describe this in more detail in Section 3. It suffers slow
mixing in sparse tables.

Diaconis and Sturmfels (1998) also proposed an algebraic algorithm to
construct a connected Markov chain in the context of incomplete tables. Other
MCMC methods proposed to sample incomplete tables also rely on computing a
Markov basis (Aoki and Takemura, 2005; Rapallo, 2006). The computational cost
of computing a Markov basis is exponential in the size of the table.
Additionally, the computation is example specific; i.e. a new basis must be
computed for each pattern of structural zeros considered.

Chen et al. (2005) introduced the first SIS method for uniform sampling of
contingency tables with given marginals. Chen (2007) extended this to
incomplete tables. Eisinger and Chen (2017) develop a sampler with improved
efficiency, particularly in sparse graphs. We will compare our sampler to SIS
approaches in our applications.

An R-package implementing the new algorithms is available at
https://github.com/jscott6/cgsampr.

##  2 Notation and Background

A graph $G\colonequals(V,E)$ is a pair with $V$ being the vertex set and $E$
the set of edges. Throughout we assume an unweighted graph has no multiple
edges, and note that the other case can be dealt with as an integer-weighted
graph. If the context requires clarification, we use $V(G)$ and $E(G)$ to
denote the sets belonging to $G$. We denote an edge from vertex $u$ to vertex
$v$ by $uv$. If the graph is undirected, $uv$ is equivalent to $vu$. An
integer-weighted graph is a triple $G\colonequals(V,E,c)$. The function
$c\colon V\times V\rightarrow\mathbb{N}_{0}$ assigns a positive integer valued
weight to each $uv\in E(G)$, and $0$ to each $uv\notin E$.

The in- and out-degrees of a vertex are the number of edges to and from the
vertex respectively. The in- and out-strengths of a vertex of a weighted graph
are the sum of the weights of edges to and from the vertex respectively. If
the graph is undirected, there is no distinction between in and out, so we
simply use the terms degree and strength of a vertex.

Two undirected graphs with the same vertex set have the same degree sequence
if every vertex has the same degree in both graphs. A similar terminology is
used for directed graphs, where both the in- and out-degrees have to be equal
for every vertex. We use analogous definitions for equivalent vertex strength
sequences.

##  3 State-Dependent Kernel Selection

Consider the following problem. Let $r\colonequals(r_{1},...,r_{I})$ and
$c\colonequals(c_{1},...,c_{J})$ be non-negative integer vectors, and let
$\mathcal{X}$ denote the set of all $I\times J$ non-negative integer matrices
such that the row and column marginals equal $r$ and $c$ respectively. Assume
$\mathcal{X}$ is non-empty. The task is to construct a Markov chain ergodic
with respect to the uniform distribution on $\mathcal{X}$.

Diaconis and Sturmfels (1998) describe a simple Markov chain for this purpose.
Given $X_{n}=x$, pick a pair of rows and a pair of columns uniformly at
random. The walk proceeds by sampling from the conditional distribution of the
delineated subtable given all other entries. An update takes the form

$\begin{matrix}+\Delta&-\Delta\\\ -\Delta&+\Delta\end{matrix}$

for $\Delta$ sampled uniformly from the admissible range: integers which do
not induce negative values in the subtable.

A Markov chain on $\mathcal{X}$ is completely characterized by its (Markov
transition) kernel $K$, a regular conditional distribution, where $K(x,A)$
represents the probability that the next state of the chain is in $A$
measurable given that the current state is $x$.

In this example, the kernel $Q$ of the chain can be viewed as randomly
selecting from a set of other kernels. Indeed, let $\mathcal{Z}$ be the
collection of indices of all $2\times 2$ sub-arrays of $I\times J$ tables. The
Gibbs update along each $z\in\mathcal{Z}$ defines a kernel $K_{z}$ on
$(\mathcal{X},\mathcal{B})$. We define a scan order as a method of choosing a
particular kernel from this collection at each iteration of the chain. The
aforementioned chain is an example of a random scan procedure; corresponding
to mixing the kernels $\\{K_{z}\\}$. The kernel of the chain is then
$Q\colonequals\sum_{z\in\mathcal{Z}}K_{z}/|\mathcal{Z}|$, where
$|\mathcal{Z}|$ represents the cardinality of $\mathcal{Z}$.

The chain suffers poor mixing in sparse matrices as $\Delta$ is often
degenerate at 0. State-dependent kernel selection asks whether we can improve
mixing by allowing the scan order to depend on the current state of the Markov
chain, whilst maintaining ergodicity with respect to the target distribution.

In the Diaconis-Sturmfels chain a useful strategy might attach density only on
the subset of $\mathcal{Z}$ for which the range of $\Delta$ would be non-zero.
However, this will not generally maintain invariance with respect to $\pi$. As
we will see, this can be overcome by either carefully constructing the scan
order, or by modifying the kernels $K_{z}$ themselves.

More generally, suppose we have a collection of kernels
$K\colonequals\\{K_{z}\colon z\in\mathcal{Z}\\}$ on a measurable space
$(\mathcal{X},\mathcal{B})$. Equip the index set $\mathcal{Z}$ with a sigma-
algebra $\mathcal{F}$. One defines a conditional distribution $\mathbb{V}$, on
$\mathcal{F}$ given the current state of the chain. For each
$x\in\mathcal{X}$, the measure $\mathbb{V}_{x}$ denotes the law of a random
variable representing the kernel from which to sample next. We will assume
throughout that the map $z\rightarrow K_{z}(x,A)$ is $\mathcal{F}$-measurable
for each $x\in\mathcal{X}$ and $A\in\mathcal{B}$.

The kernel of the state-dependent Markov chain $Q$ on
$(\mathcal{X},\mathcal{B})$ is defined by the integrals

$Q(x,\cdot)\colonequals\int
K_{z}(x,\cdot)\mathbb{V}_{x}(\mathop{}\\!\mathrm{d}z)\text{ for all
}x\in\mathcal{X},$ (1)

which are measures in the second argument. Suppose a Markov chain defined by
$Q$ has some current state $x$. The chain samples some kernel $K_{z}$
according to the measure $\mathbb{V}_{x}$, and proceeds to sample the next
state of the chain from the measure $K_{z}(x,\cdot)$.

###  3.1 Decomposing a Kernel

Any Markov kernel can be represented as a weighted average of a collection of
kernels. Definition 3.1 helps make this precise.

######  Definition 3.1.

$(K,\mathbb{V})$ is a decomposition of a kernel $Q$ on
$(\mathcal{X},\mathcal{B})$ if (1) holds.

Definition 3.1 implies that every kernel $Q$ has an ‘identity’ decomposition,
given by $(\\{K_{1}\\},\mathbb{V})$ with $K_{1}=Q$ and
$\mathbb{V}_{x}(\\{1\\})=1$ for all $x\in\mathcal{X}$. The decomposition of a
kernel is typically not unique.

######  Definition 3.2.

Let $(K,\mathbb{V})$ be a decomposition of $Q$. Suppose $\mathbb{V}_{x}$ is
dominated by a $\sigma$-finite measure $\mu$ for all $x\in\mathcal{X}$. Then,
the decomposition is symmetric if there exist densities $\\{f_{x}\\}$ (each
with respect to $\mu$) such that for each $z$ and $x$, $f_{x}(z)=f_{y}(z)$ for
$K_{z}(x,\cdot)$-almost every $y$.

Suppose $Q$ has a symmetric decomposition $(K,\mathbb{V})$. Definition 3.2
implies that, at any given iteration of the chain, once $K_{z}$ has been
sampled, the density of sampling $K_{z}$ under almost every possible new state
of the chain is unchanged. Any state-independent kernel selection is symmetric
(for example, random scan). The ‘identity’ decomposition of any kernel is
symmetric. Definition 3.2 generalises this concept to state-dependent mixing.

As an example, consider a three-state state space, as depicted in Figure 1.
The left figure defines three kernels on this space. Let $K_{1}$, $K_{2}$ and
$K_{3}$ be defined by the dotted, solid and dashed transitions in the figure
respectively. A naive chain might pick from these kernels randomly,
irrespective of the current state. If the chain is in state $i$, the kernel
$K_{i}$ cannot change the state. A better strategy might pick randomly from
the other two kernels, guaranteeing that the chain moves to a new state. This
state-dependent strategy is easily seen to lead to a symmetric decomposition.

$1$$2$$3$$1$$1$$1$$1$$1$$1$$1$$1$$1$

$1$$2$$3$$.5$$.5$$.5$$.5$$.5$$.5$

Figure 1: Kernels used as examples of a decomposition in Section 3.1. Left:
$K_{1}$ (dotted), $K_{2}$ (solid), $K_{3}$ (dashed). Right: Kernel $Q$

######  Lemma 3.3.

A kernel $Q$ is reversible with respect to a distribution $\pi$ if and only if
it has a symmetric decomposition $(K,\mathbb{V})$ where every $K_{z}\in K$ is
reversible with respect to $\pi$.

If both $\mathcal{X}$ and $\mathcal{Z}$ are countable, then it is easy and
instructive to illustrate detailed balance when the conditions of Lemma 3.3
hold. In particular, suppose $(K,\mathbb{V})$ is symmetric and each $K_{z}$ is
reversible with respect to the uniform distribution on $\mathcal{X}$. Then

$Q(x,y)=\sum_{z\in\mathcal{Z}}\mathbb{V}_{z}(x)K_{z}(x,y)=\sum_{z\in\mathcal{Z}}\mathbb{V}_{z}(y)K_{z}(y,z)=Q(y,x).$

We have used the fact that $\mathbb{V}_{x}(z)=\mathbb{V}_{y}(z)$ whenever
$K_{z}(x,y)>0$, and that $K_{z}(x,y)=K_{z}(y,x)$.

Recall the three state Markov chain, shown in Figure 1. Each of the three
kernels is reversible with respect to the uniform distribution. Additionally,
the decomposition under the state-dependent strategy previously suggested is
symmetric. Lemma 3.3 thus implies the overall chain (whose transitions are
shown in Figure 1, right) also conserves the uniform distribution.

Lemma 3.3 implies we can transform the problem of establishing reversibility
of $Q$ with respect to a distribution into a problem of finding a
decomposition $(K,\mathbb{V})$ of $Q$ satisfying particular conditions. This
will be our strategy in Section 4.

###  3.2 Kernel Selection as an Auxiliary Variable

We cannot usually expect a kernel decomposition $(K,\mathbb{V})$ to be
symmetric. We now outline a more general sampling strategy. The kernel $Q$ can
be interpreted as defining a Markov chain $(Z_{n},X_{n})_{n\geq 0}$ on the
extended space $\mathcal{Z}\times\mathcal{X}$. Let the current state of this
chain be $y:=(z,x)$. If $y^{\prime}$ represents the next state then

$\Pr(y^{\prime}\in A_{1}\times A_{2}\mid
y)=\int_{A_{1}}K_{z}(x,A_{2})\mathbb{V}_{x}(\mathop{}\\!\mathrm{d}z)$ (2)

for $A_{1}$ and $A_{2}$ measurable. If the chain defined by (2) is ergodic
with respect to the joint distribution $\mathbb{V}\otimes\pi$, then the sub-
chain $(X_{n})$ with kernel $Q$ is ergodic with respect to its marginal
distribution $\pi$. This suggests a two stage sampler, a step of which is
shown below.

  1. 1.

Sample $Z_{n}\sim\mathbb{V}_{X_{n-1}}$.

  2. 2.

Sample $X_{n}\sim\mathbb{V}\otimes\pi(\cdot\mid Z_{n},T_{Z_{n}}(X_{n-1}))$.

Each $T_{z}$ is a statistic on $\mathcal{X}$. Both steps are Gibbs kernels
reversible with respect to the joint distribution $\mathbb{V}\otimes\pi$, and
so their composition also conserves the joint. We use this strategy in Section
5.

##  4 Sampling Unweighted Graphs

Let $G_{0}$ be a given directed or undirected graph with a finite vertex set
$V$. Let $\mathcal{F}$ be a subset of possible edges of a graph with vertex
set $V$. Let $\mathcal{G}$ be the set of all graphs $G$ with the same vertex
set and degree sequence as $G_{0}$, and additionally satisfying
$E(G)\cap\mathcal{F}=E(G_{0})\cap\mathcal{F}$. Our goal is to sample from the
uniform distribution on $\mathcal{G}$.

Intuitively, the set $\mathcal{F}$ represents edges known by design to be
present or absent. Given vertices $u$ and $v$, if $uv$ belongs to
$\mathcal{F}$ then $uv$ is either present in all graphs in $\mathcal{G}$, or
in none. We stress by design because the constraints imposed by the degree
sequence and $\mathcal{F}$ may imply that further edges are present or absent
in all graphs of $\mathcal{G}$. We call this set $\tilde{\mathcal{F}}$ the set
of known edges, and formally define it as

$\tilde{\mathcal{F}}=\\{\text{possible edges }uv:uv\in
G_{0}\Leftrightarrow(uv\in G\text{ for all }G\in\mathcal{G})\\}$

We show in Section 4.2 a method of obtaining $\tilde{\mathcal{F}}$.

For every graph $G\in\mathcal{G}$, we will extensively use two ‘neighborhood’
sets associated to each vertex $u$. The set $N_{G}(u)$ are the in-neighbors of
$u$, excluding any vertex $v$ for which the edge $vu$ is known. $M_{G}(u)$ is
the set of all vertices $v$ which are not out-neighbors of $u$, and for which
the absence of $uv$ is not known. These are defined as

$N_{G}(u)\colonequals\\{v\in V\colon vu\in
E(G),vu\notin\tilde{\mathcal{F}}\\},$ $M_{G}(u)\colonequals\\{v\in V\colon
uv\notin E(G),uv\notin\tilde{\mathcal{F}}\\}.$

Here is a sketch of one iteration of the scheme. An initial vertex $W_{0}$ is
sampled uniformly from the set of all vertices $v$ for which $N_{G}(v)$ is
non-empty, and set $n=0$. Sample $W_{1}$ uniformly from $N_{G}(W_{0})$, then
sample $W_{2}$ uniformly from $M_{G}(W_{1})$. Replace $W_{1}W_{0}$ in $E(G)$
with $W_{1}W_{2}$. Letting $n=n+2$, iterate this procedure, however in each
subsequent step $W_{n+1}$ cannot be $W_{n-1}$; this prevents the sampler
adding the edge $W_{n+1}W_{n+2}$, and removing it in the next iteration, and
should improve state space exploration. Iterate until $W_{n}$ is $W_{0}$, at
which point all degrees have been maintained. Algorithm 1 makes this precise.
Figure 2 shows a straightforward example of one step of the sampler.

Algorithm 1 One Iteration of the Unweighted Graph Sampler (UGS)

1:$G$, $\tilde{\mathcal{F}}$

2:$W_{-1}\leftarrow*$

3:$W_{0}\sim U(\\{v\in V\colon N_{G}(v)\neq\emptyset\\})$

4:$n\leftarrow 0$

5:repeat

6: $W_{n+1}\sim U(N_{G}(W_{n})\setminus\\{W_{n-1}\\})$

7: $W_{n+2}\sim U(M_{G}(W_{n+1}))$

8: $E(G)\leftarrow E(G)\setminus\\{W_{n+1}W_{n}\\}$

9: $E(G)\leftarrow E(G)\cup\\{W_{n+1}W_{n+2}\\}$

10: $n\leftarrow n+2$

11:until $W_{n}=W_{0}$

return $G$

$w_{0}$$w_{2}$$w_{1}$

$w_{0}$$w_{2}$$w_{1}$$w_{3}$

$w_{0}$$w_{2}$$w_{1}$$w_{3}$

Figure 2: One iteration of Algorithm 1 with two iterations in the while loop.
Left and center show graphs and quantities immediately prior to the first and
second edge swaps respectively. Right shows the returned graph. This move
corresponds to the well known ‘switch’ step originating from Ryser (1963).

###  4.1 Properties of Algorithm 1

Let $a^{r}$ denote the reverse of a finite sequence $a$. Given a graph, let
$a_{1}b_{1}\leftrightarrow a_{2}b_{2}$ denote the operation of replacing the
edge $a_{1}b_{1}$ with the edge $a_{2}b_{2}$. We will refer to this operation
as an (edge) swap. We call $a_{1}b_{1}\leftrightarrow a_{2}b_{2}$ viable if
and only if $a_{1}b_{1}$ is an edge, $a_{2}b_{2}$ is not an edge and both
$a_{1}b_{1}$ and $a_{2}b_{2}$ are not in $\mathcal{F}$.

A single iteration of Algorithm 1 samples a (random) sequence of vertices $W$.
Proposition 4.1 implies that this sequence will be finite, so that $W$ takes
the form $w_{0}w_{1}...w_{k}w_{0}$ for some $k$ odd. Let $\mathcal{W}$ be the
collection of sequences taking this form.

######  Proposition 4.1.

For any input graph $G\in\mathcal{G}$ and any $\mathcal{F}$, the expected
length of the vertex sequence $W$ (formed by Algorithm 1) is finite.

Let two sequences be equivalent if and only if they are either identical or
they are each others’ reverse. We let $\mathcal{Z}$ be the quotient set of
$\mathcal{W}$ by this equivalence relation.

We will associate each equivalence class $z\in\mathcal{Z}$ with a kernel on
$\mathcal{G}$. Fix any $z$ and let $w$ be a representative of $z$. Consider
the following Markov chain on $\mathcal{G}$. From the current state, attempt
to iteratively perform the swaps $w_{1}w_{0}\leftrightarrow w_{1}w_{2}$,
$w_{3}w_{2}\leftrightarrow w_{3}w_{4}$, … , $w_{k}w_{k-1}\leftrightarrow
w_{k}w_{0}$ to obtain the next state. We say this move is viable if and only
if all of the swaps are viable when applied iteratively. We refer to this
sequence of swaps as the swaps corresponding to $w$. If the swaps are not
viable, attempt the swaps corresponding to $w^{r}$; i.e.
$w_{k}w_{0}\leftrightarrow w_{k}w_{k-1}$, …, $w_{1}w_{2}\leftrightarrow
w_{1}w_{0}$. If neither swap sequence is viable, then the next state of the
chain is unchanged. We define $K_{z}$ as the kernel of this chain. Remark 4.2
implies that $K_{z}$ is well-defined; specifically, the definition is
independent of the chosen representative of $z$.

######  Remark 4.2.

If the sequences $w$ and $w^{r}$ are distinct and the swaps corresponding to
$w$ are viable, then the swaps corresponding to $w^{r}$ are not viable.

Let $K$ be the collection of these kernels. The conditional distribution
$\mathbb{V}$ on $\mathcal{Z}$ is defined implicitly by the law of $W$.
Formally, the sampler selects a kernel $K_{z}$ by sampling a vertex sequence
$w\in z$. $K_{z}$ would also be selected if $w^{r}$ were sampled. Lemma 4.3
implies that $Q$ is reversible with respect to the uniform distribution on
$\mathcal{G}$.

######  Lemma 4.3.

$(K,\mathbb{V})$ is a symmetric decomposition of $Q$, and each $K_{z}\in K$ is
reversible with respect to the uniform distribution on $\mathcal{G}$.

Proposition 4.4 holds by additionally showing the chain is connected.

######  Proposition 4.4.

The Markov chain is ergodic with respect to the uniform distribution on
$\mathcal{G}$.

###  4.2 Identifying all Known Edges/Non-Edges

We show how to determine all structurally fixed edges/non-edges prior to
sampling; or in other words, how to determine $\tilde{\mathcal{F}}$ from
$\mathcal{F}$ and the degree sequence. Our method makes use of auxiliary
graphs, which we now define.

Given any graph $G\in\mathcal{G}$ we construct an auxiliary bipartite digraph
$B\colonequals(U,V,E)$ as follows. Let $U=\\{u_{i}\\}$ and $V=\\{v_{i}\\}$ for
$i=1,...,n$. Fix any vertex $i$ and any vertex $j$ in $V(G)$. We will add an
edge to $B$ if and only if $ij$ is not in $\mathcal{F}$. If additionally $ij$
is in $E(G)$, add $v_{j}u_{i}$ to $E(B)$, otherwise add $u_{i}v_{j}$. Figure 3
shows an example of one such graph.

$1$$2$$3$$4$

$u_{1}$$u_{2}$$u_{3}$$u_{4}$$v_{1}$$v_{2}$$v_{3}$$v_{4}$

Figure 3: Left: an undirected graph. Right: Auxiliary graph constructed from
it, as outlined in Section 4.2.

Let $B_{\mathcal{G}}$ denote the collection of all graphs generated this way
from the set $\mathcal{G}$. Proposition 4.5 shows that we can identify all
known edges/non-edges prior to sampling by identifying all strongly connected
components of any graph in $B_{\mathcal{G}}$. This can be done using a depth-
first search on $B_{G}$, followed by another depth-first search on the
transposed graph. The running time of this preprocessing procedure is
$\Theta(2n+n^{2})$ (Cormen, 2009, chap. 22).

######  Proposition 4.5.

Fix any graph $B\in B_{\mathcal{G}}$, and partition its vertex set into
strongly connected components $S_{1},...,S_{K}$. The vertex pair $ij$ belongs
to $\tilde{\mathcal{F}}$ if and only if either there is no edge incident to
$u_{i}\in S_{k}$ and $v_{j}\in S_{l}$, or if $k\neq l$.

##  5 Sampling Weighted Graphs

Let $G_{0}$ be a given integer-weighted directed or undirected graph with a
finite vertex set $V$. Let $\mathcal{F}$ be a subset of possible edges of a
graph with vertex set $V$. Let $\mathcal{G}$ be the set of all graphs $G$ with
the same vertex set and strength sequence as $G_{0}$, and additionally
satisfying $c_{G}(uv)=c_{G_{0}}(uv)$ for all $uv\in\mathcal{F}$. Our goal is
to sample from the uniform distribution on $\mathcal{G}$.

Define $\mathcal{Z}$ as in Section 4.1. Fix $z\in\mathcal{Z}$ and let $w$ be a
vertex sequence belonging to $z$. Let $n_{w}(uv)$ denote the difference
between the number of occurrences of the possible edge $uv$ in $w_{1}w_{0}$,
$w_{3}w_{2}$,…,$w_{k}w_{k-1}$ and the number of occurrences in $w_{1}w_{2}$,
$w_{3}w_{4}$,…,$w_{k}w_{0}$.

We will define $K_{z}$ as the kernel of the following Markov chain. From the
current state $G$, the chain can only move to graphs whose edge strengths take
the value $c_{G}(uv)+n_{w}(uv)\Delta$ for all $uv$, and for some $\Delta$ in
its admissible range: integers for which the resulting graph is in
$\mathcal{G}$. We will denote the range of $\Delta$ by
$[\Delta_{low},\Delta_{up}]$. Let $G_{\Delta^{*}}$ be the graph obtained by
sampling $\Delta=\Delta^{*}$. The next state of the chain is determined by
sampling $\Delta$ proportional to $\mathbb{V}_{G_{\Delta^{*}}}(Z_{n})$ for
each $\Delta^{*}\in[\Delta_{low},\Delta_{up}]$, where $\mathbb{V}$ is the (as
yet, undefined) law representing our kernel selection strategy.

Observe that there exists a statistic $T_{z}$ satisfying

$T_{z}^{-1}\circ
T_{z}(G)=\\{G_{\Delta^{*}}:\Delta^{*}\in[\Delta_{low},\Delta_{up}]\\}.$

Thus, sampling $\Delta$ as suggested is the same as the second stage of the
sampling scheme presented in Section 3.2, and the overall chain defined by $Q$
will will be reversible with respect to the uniform distribution on
$\mathcal{G}$.

###  5.1 Kernel Selection Strategy

$\Delta$ will be degenerate at $0$ if there exists $uv$ and
$u^{\prime}v^{\prime}$ such that $n_{w}(uv)$ is positive,
$n_{w}(u^{\prime}v^{\prime})$ is negative, and
$c_{G}(uv)=c_{G}(u^{\prime}v^{\prime})=0$. As the sparsity of $G$ increases,
the proportion of kernels in the collection $K$ which can move the chain to a
new state decreases. Thus, any state-independent strategy will suffer a
degradation in performance. Our state-dependent strategy avoids this.

We alternately samples from two sets associated with each vertex $u$. We
redefine $N_{G}(u)$ and $M_{G}(u)$ as

$N_{G}(u)\colonequals\\{v\in V\colon vu\in E(G),vu\notin\mathcal{F}\\},$
$M_{G}(u)\colonequals\\{v\in V\colon uv\notin\mathcal{F}\\}.$

The kernel sampling procedure is as follows. Let $W_{-1}=*$ and sample
$W_{0}\sim U(\\{v\in V\colon N_{G}(v)\neq\emptyset\\})$. Letting $n=0$, repeat
the following until termination.

  1. 1.

Sample $W_{n+1}\sim U(N_{G}(W_{n})\setminus\\{W_{n-1}\\})$, unless the set
from which to sample is empty, in which case return $K_{\text{id}}$.

  2. 2.

If $W_{0}\in M_{G}(W_{n+1})\setminus\\{W_{n}\\}$, then return $K_{[W]}$, where
$W\colonequals W_{0}...W_{n+1}W_{0}$. Otherwise sample $W_{n+2}\sim
U(M_{G}(W_{n+1})\setminus\\{W_{n}\\})$, unless the set from which to sample is
empty, in which case return $K_{\text{id}}$. Let $n=n+2$.

In practice, returning $K_{id}$ is rare and occurs primarily in cases of
extreme sparsity; specifically when some vertex has only one in-neighbour.
Figure 4 graphically represents several example sequences $W$.

$w_{0}$$w_{1}$$w_{2}$$w_{3}$$+\Delta$$-\Delta$$+\Delta$$-\Delta$

$w_{0}$$w_{1}$$w_{2}$$w_{3}$$w_{4}$$w_{5}$$+\Delta$$-\Delta$$+\Delta$$-\Delta$$+\Delta$$-\Delta$

$w_{0}$$w_{1}$$w_{2,6}$$w_{3,7}$$w_{4}$$w_{5}$$w_{8}$$w_{9}$$+\Delta$$-\Delta$$+2\Delta$$-\Delta$$-\Delta$$+\Delta$$-\Delta$$+\Delta$$-\Delta$

Figure 4: Sampled edge (black solid), sampled non-edge (black dashed), edge
(gray solid) and fixed edge (gray dotted). If there are no fixed edges/non-
edges, $W$ will have length 4 (left). Otherwise, longer sequences are
required, as shown in the center and right examples. The right example
demonstrates vertices may be revisited.

###  5.2 Conditional Distribution of $\Delta$

Suppose we sample a vertex sequence $w\in z$ using the strategy proposed in
Section 5.1. To ease notation, define $G_{low}$ and $G_{up}$ as the graphs
obtained at $\Delta_{low}$ and $\Delta_{up}$ respectively. It is often not
possible to sample $w$ from $G_{low}$, or to sample $w^{r}$ from $G_{up}$.
This is why we collapse $w$ and $w^{r}$ into $z$; doing so ensures we can
always sample $z$ from each graph in the support of $K_{z}(G,\cdot)$.

Suppose $(\Delta_{low},\Delta_{up})$ is non-empty. Let $G^{*}$ be any graph
obtained from sampling some $\Delta^{*}\in(\Delta_{low},\Delta_{up})$. Define
$\alpha_{int}$ as the probability of sampling $z$ from $G^{*}$. This quantity
is the same regardless of the chosen $G^{*}$. If on the other hand the range
$(\Delta_{low},\Delta_{up})$ is empty, let $\alpha_{int}$ take an arbitrary
finite value.

The random variable $\Delta$ then has conditional measure

$\frac{1}{d}\left(\mathbb{V}_{G_{low}}(z)\delta_{\Delta_{low}}+\mathbb{V}_{G_{up}}(z)\delta_{\Delta_{up}}+\alpha_{int}\sum_{\Delta_{low}<\tilde{\Delta}<\Delta_{up}}\delta_{\tilde{\Delta}}\right)$
(3)

where
$d:=\mathbb{V}_{G_{low}}(z)+\mathbb{V}_{G_{up}}(z)+\alpha_{int}\max(\Delta_{up}-\Delta_{low}-1,0)$.

$\mathbb{V}_{G_{low}}(z)$, $\mathbb{V}_{G_{up}}(z)$ and $\alpha_{int}$ are
easily computed by following the details of the kernel selection strategy in
Section 5.1. Algorithm 2 gives pseudo-code for one iteration of the sampler.

By construction, the chain is reversible with respect to the uniform
distribution. Proposition 5.1 holds by additionally showing the chain is
connected.

######  Proposition 5.1.

The Markov chain is ergodic with respect to the uniform distribution on
$\mathcal{G}$.

Algorithm 2 One Iteration of the Weighted Graph Sampler (WGS)

1:$G$, $\mathcal{F}$

2:$\text{visits}(uv)\leftarrow 0$; $\text{edges}\leftarrow\\{\\}$;
$(\Delta_{low},\Delta_{up})\leftarrow(-\infty,\infty)$; $W_{-1}\leftarrow*$;
$n\leftarrow 0$

3:$W_{0}\sim U(\\{v\in V\colon N_{G}(v)\neq\emptyset\\})$

4:repeat

5: if $N_{G}(W_{n})\setminus\\{W_{n-1}\\}=\emptyset$ then return $G$

6: else $W_{n+1}\sim U(N_{G}(W_{n})\setminus\\{W_{n-1}\\})$

7: if $W_{0}\in M_{G}(W_{n+1})\setminus\\{W_{n}\\}$ then $W_{n+2}\leftarrow
W_{0}$

8: else if $M_{G}(W_{n+1})\setminus\\{W_{n}\\}=\emptyset$ then return $G$

9: else $W_{n+2}\sim U(M_{G}(W_{n+1})\setminus\\{W_{n}\\})$

10: $\text{edges}\leftarrow\text{edges}\cup\\{W_{n+1}W_{n},W_{n+1}W_{n+2}\\}$

11: $\text{visits}(W_{n+1}W_{n})\leftarrow\text{visits}(W_{n+1}W_{n})+1$

12: $\text{visits}(W_{n+1}W_{n+2})\leftarrow\text{visits}(W_{n+1}W_{n+2})-1$

13: $n\leftarrow n+2$

14:until $W_{n}=W_{0}$

15:for each $\text{edge}\in\text{edges}$ do

16: if $\text{visits}(\text{edge})>0$ then
$\Delta_{low}\leftarrow\max\left(\Delta_{low},-\lfloor\frac{c_{G}(\text{edge})}{\text{visits}(\text{edge})}\rfloor\right)$

17: if $\text{visits}(\text{edge})<0$ then
$\Delta_{up}\leftarrow\min\left(\Delta_{up},\lfloor\frac{c_{G}(\text{edge})}{\text{visits}(\text{edge})}\rfloor\right)$

18:Sample $\Delta$ according to (3)

19:for each $\text{edge}\in\text{edges}$ do

20: $c_{G}(\text{edge})\leftarrow
c_{G}(\text{edge})+\text{visits}(\text{edge})\Delta$

return G

##  6 Simulation Study

Figure 5: Comparative performance of WGS and DS in simulation study of Section
6.

Methods used in this section, and in Section 7, were programmed in C/C++, and
run on an Intel Core i5-6360U 2GHz CPU. Here we investigate the comparative
performance of Algorithm 2 (WGS) and the Diaconis & Sturmfels chain (DS)
introduced in Section 3 in sparse contingency tables.

Let $A$ and $B$ be independent, uniform categorical variables with $L$ levels.
For each $L$ in $10,15,...100$, we simulate 200 bivariate samples and
construct the corresponding $L\times L$ contingency table. We then use the
samplers to approximately sample from the uniform distribution over all tables
with the same margins, and obtain 100,000 samples of the likelihood ratio
statistic computed from these tables. Thinning was set to $L$.

Define the (empirical) mixing rate of a Markov chain to be the observed
proportion of iterations that change the state of the chain. Figure 5 shows
plots this quantity for each sampler as we varied $L$. In addition Figure 5
plots the effective sample size per second (ESS/s) computed from the sample of
LR statistics.

The mixing rate of DG rapidly approaches zero, while for WGS it levels at
around 0.5. The effective sample size per second of DG rapidly falls rapidly,
in line with the mixing rate, while it reduces at a much slower rate for WGS.
When $L=100$, WGS appears to be over 1000 times more efficient than DG. Figure
5 provides trace plots of the likelihood ratio statistics.

Eisinger and Chen (2017) develop efficient SIS methods for sampling tables
from the uniform distribution over all tables with given margins. Their method
labeled $SIS-G$ (coded in C) took 227 seconds to produce 1000 samples of
$100\times 100$ tables with both margins equal to $(5,1,...,1)$, evaluated on
a laptop with a 2.2 GHz Intel Core i7 processor. The authors provide code for
a cell-by-cell SIS method SIS-G*. Using SIS-G* on the same example, we
estimated an ESS/s of around 7. Figure 5 shows that with $L=100$, the ESS/s
for the LR statistics produced by WGS was over 1400. Therefore, it appears WGS
can produce independent tables orders of magnitude faster than available SIS
methods in large, spare tables.

##  7 Applications

Reported standard errors were computed using spectral methods from R’s coda
package. These estimates were compared to those obtained using batch means,
and where feasible, bootstrapping. These latter estimates are not reported as
there was little discernible difference from those obtained by spectral
methods. Thinning used in each method was set to approximately equate the
resulting standard errors. We used burn-in equivalent to 20% of samples
obtained.

###  7.1 Ecological Networks

A food web encodes predator-prey relationships within a group of species, and
has a natural representation as a digraph. Each species is a node in the graph
and a link exists from species A to species B iff B consumes A.

Ecologists wish to identify and explain structural patterns in observed food
webs including motifs, diet contiguity, intervality, connectance and
compartmentalization. We will focus on assessing the tendency towards
compartmentalization in food webs. Compartmentalization describes the extent
to which species can be partitioned into distinct groups such that linkage
density within groups is greater than that between groups (Girvan and Newman,
2002; Krause et al., 2003). The level of compartmentalization in food webs is
an important determinant of the spread of ecological perturbations. If a food
web is highly compartmentalized, a perturbation should spread quicker within
groups than between groups. Thus, higher compartmentalization should reduce
systemic risks in the network and increase stability. Compartmentalization may
also determine the effect of a perturbation to the network. An effect reducing
compartmentalization, like the introduction of a generalist predator, may
reduce the stability of the network (Guimerà et al., 2010).

Figure 6: Food web of the Chesapeake bay ecosystem.

Figure 6 depicts the food web of 33 species in the Chesapeake bay in the
summer. The data was collected by Baird and Ulanowicz (1989) and is provided
in the R package Cheddar (Hudson et al., 2018).

Pimm and Lawton (1980) proposed a statistic $\bar{C}$ to measure the level of
compartmentalization in a food web. Here we describe a directed analogue of
this statistic. Let $G$ represent a food web of $n$ species, and $i$ and $j$
be two species. Let $c_{ij}$ be the number of shared predators of species $i$
and $j$ as a proportion of the total number of predators of $i$ and $j$.
$\bar{C}$ is then the mean of the off-diagonal elements of $(c_{ij})$.

$\bar{C}\colonequals\frac{1}{n(n-1)}\sum_{i=1}^{n}\sum_{j=1,j\neq
i}^{n}c_{ij}$

$\bar{C}$ takes values in $[0,1]$ and higher values are associated with
greater levels of compartmentalization.

We begin by testing whether the observed level of $\bar{C}_{0}=0.0260$ can be
considered high when compared to the set of all graphs with the same in-degree
and out-degree sequence as $G$ (Ulrich and Gotelli, 2007). With thinning of 5,
Algorithm 1 (UGS) took around 2 second to obtain 100,000 samples. The
estimated p-value was $0.0163\pm 4.3\mathrm{e}{-4}$, complementing previous
results suggesting food webs have high compartmentalization when compared to
random graphs where species have an equal probability of linking to each other
species (Krause et al., 2003; Rezende et al., 2009).

Chen (2007) develop a sequential importance sampling strategy SIS_CP1 for
sampling zero-one tables, but which is equivalently capable of uniformly
sampling unweighted digraphs with fixed degree sequence and an arbitrary
pattern of structural non-edges. SIS_CP1 took 33 seconds to obtain 100,000
samples, estimating a p-value of $0.0158\pm 4.3\mathrm{e}{-4}$.

Guimerà et al. (2010) find that compartmentalization observed in real food
webs is not unusual when compared to networks generated under niche models,
and conclude that ‘compartmentalization can be explained solely by the niche-
valued ranking of species’.

We attempt to test this hypothesis for the Chesapeake bay food web. We compute
the chain averaged trophic level (Williams and Martinez, 2004) for each
species, and assume any given species is forbidden from consuming other
species with a higher trophic level. The resulting forbidden links should help
to control for the food web’s trophic structure. The assumption induces 565
forbidden edges in the null distribution.

Again using thinning of 5, UGS took 2 seconds to obtain 100,000 samples. The
estimated effective sample size was over 95,000, giving an estimated p-value
of $0.0568\pm 7.5\mathrm{e}{-4}$. At a significance level of $\alpha=0.05$, we
can no longer conclude that the level of compartmentalization in the
Chesapeake food web is unusual under the null distribution. Our method of
determining trophic structure is relatively crude, and a closer analysis of
the food web is warranted before drawing any conclusions.

SIS_CP1, on the other hand, took 24 seconds to run and over 97% of the samples
produced were discarded as invalid, leaving only 3,069 to be used for
estimation. The estimated p-value was $0.0558\pm 6.3e\mathrm{-3}$. Using
alternative methods to calculate the species’ trophic levels gives rise to
other sets of forbidden edges. For some such patterns, SIS_CP1 was unable to
construct a single valid sample. SIS_CP1 cannot reliably sample graphs in the
face of arbitrary sets of fixed non-edges.

Forbidden links have a long history in ecological networks, and result from
factors including spatio-temporal uncoupling and morphological and
physiological-biochemical constraints (see Olesen et al., 2010). A realistic
specification of the null distribution of a hypothesis test must take into
account forbidden links, particularly as they can change the outcome of the
test. Our algorithm provides this flexibility, even in large sparse graphs.
Other MCMC methods generally require a Markov basis, which is computationally
infeasible for networks with more than a few nodes.

###  7.2 Incomplete Tables

A sample of size $N$ consisting of measurements on two categorical variables
can be represented as a two-way contingency table of dimensions $I\times J$,
where $I$ and $J$ represent the number of levels of the first and second
variables respectively.

A table with structurally fixed cells is referred to as incomplete. Incomplete
tables arise in several contexts. Particular combinations of the variables may
be impossible, forcing zero entries in the corresponding cells. Alternatively,
some observations may be missing. In some contexts, researchers may wish to
fit composite models by partitioning the cells into subsets, and fitting a
separate log-linear model for each group (Goodman, 1963, 1968; Fienberg,
1969)). See Bishop and Fienberg (1969) for extensive examples of incomplete
tables.

Pearson (1904) introduce a mobility table recording the occupations of 775
fathers and their sons. There are 14 separate job classifications. Given a
table $t$ and a log-linear model, define the chi-squared statistic $\chi^{2}$
as

$\chi^{2}=\sum_{i=1}^{r}\sum_{j=1}^{c}\frac{(t_{ij}-\hat{m}_{ij})^{2}}{\hat{m}_{ij}}$

where $\hat{m}$ is a matrix representing the MLEs of expected cell counts
under the log-linear model. For the independence model, we observe
$\chi^{2}=1005.45$. This appears high, indicating that the model of
independence is not appropriate.

Goodman (1965) propose that sons are liable to inherit their occupational
status from their fathers, but conditional on them changing occupation, their
occupation choice is considered independent from their father’s status. This
resulting model of ‘quasi-perfect mobility’ fixes the diagonals of the table.

Alternatively, Fienberg (1969) propose a graphical procedure to determine
cells with large interactions. They identify 14 such cells, namely $(1,1)$,
$(13,13)$, $(11,11)$, $(7,1)$, $(3,3)$, $(3,1)$, $(7,7)$, $(2,1)$, $(2,2)$,
$(4,4)$, $(12,12)$, $(10,10)$, $(6,8)$ and $(6,12)$. 10 of these are along the
diagonal.

We wish to test whether the remaining observed counts are plausible under the
independence model, conditional on the 14 cells being fixed. The model of
independence applied to tables with structural zeros is known as the quasi-
independence model (Goodman, 1968). The MLEs of the expected cell counts can
easily be computed using the IPFP procedure.

Diaconis and Efron (1985) propose the uniform distribution on the set of all
$I\times J$ tables with the same margins as an alternative to the independence
hypothesis. This can also be applied to incomplete tables. With obvious
modifications, our Markov chain can also be used to sample from the hyper-
geometric distribution to evaluate the hypothesis of quasi-independence.

Fixing the aforementioned 14 cells and applying the quasi-independence model
yields $\chi^{2}=345.21$. Using thinning of 50, WGS obtained 100,000 samples
in 3.5 seconds. The estimated p-value was $0.99987\pm 4.1e-5$ and the
estimated effective sample size was over 70,000. Conditional on fixing the 14
cells, the table appears to fit the quasi-independence hypothesis extremely
well. The oddly high p-value may perhaps be less surprising if you consider
that Feinberg’s graphical procedure has chosen the pattern of fixed cells in
order to remove large deviations from independence.

Diaconis and Sturmfels (1998) propose an alternative Markov chain based on
computational algebra, and requiring the computation of a Markov basis. Using
the software Macaulay 2 (Grayson and Stillman (n.d.)), we were unable to
compute a Markov basis for the support of the distribution.

Chen (2007) propose an SIS algorithm for performing the conditional volume
test on incomplete tables. They provide an implementation capable of handling
at most one structural zero in each column. Therefore, we can use their method
to test the quasi-perfect mobility model. The algorithm produced 100,000
samples in 1 second, however the variation in the importance weights was
extremely high, indicating that the sampling distribution is far from uniform
and giving a negligible effective sample size.

It appears current SIS methods for incomplete tables are not scalable in the
size of the table. This is also true for Markov chain methods using
computational algebra. Our sampler appears to largely overcome these
difficulties.

##  8 Discussion

This paper has developed a new class of MCMC samplers for two important
sampling problems. First, for sampling from the set of unweighted graphs
respecting prescribed vertex degrees. Second, for sampling from the set of
weighted graphs respecting prescribed vertex strengths.

The samplers appear more efficient than existing methods in sparse settings,
and also when there are arbitrary observed edges/non-edges. We have presented
examples where alternative MCMC methods are infeasible as they rely on
computing a Markov basis, and where existing SIS methods perform poorly. In
contrast, our methods do not require computing a Markov basis, and appear to
be orders of magnitude more efficient in these examples.

State-dependent mixing of Markov kernels is a general concept, and the
specific implementation of our samplers is not unique. The technique could be
used to develop alternative samplers specialized to particular settings. The
methods can be readily extended to sample from arbitrary distributions known
up to a normalization constant. Thus, the samplers can be adapted to carry out
Bayesian network tomography in the case of a star network topology. In
contrast, SIS methods are not readily adaptable to more general distributions.
A theoretical analysis of the mixing times of the new samplers is beyond the
scope of this paper. Future work could try to establish whether the chains are
rapidly mixing.

##  9 Acknowledgements

The work of the first author supported by an EPSRC Research Studentship.

## References

  * (1)
  * Agresti (2003) Agresti, A. (2003), Categorical data analysis, Vol. 482 John Wiley & Sons.
  * Aoki and Takemura (2005) Aoki, S., and Takemura, A. (2005), “Markov chain Monte Carlo exact tests for incomplete two-way contingency tables,” Journal of Statistical Computation and Simulation, 75(10), 787–812.
  * Baird and Ulanowicz (1989) Baird, D., and Ulanowicz, R. E. (1989), “The Seasonal Dynamics of the Chesapeake Bay Ecosystem,” Ecological monographs, 59(4), 329–364.
  * Bayati et al. (2010) Bayati, M., Kim, J. H., and Saberi, A. (2010), “A Sequential Algorithm for Generating Random Graphs,” Algorithmica, 58(4), 860–910.
  * Bezáková et al. (2012) Bezáková, I., Sinclair, A., Štefankovič, D., and Vigoda, E. (2012), “Negative Examples for Sequential Importance Sampling of Binary Contingency Tables,” Algorithmica, 64(4), 606–620.
  * Bishop and Fienberg (1969) Bishop, Y. M. M., and Fienberg, S. E. (1969), “Incomplete Two-Dimensional Contingency Tables,” Biometrics, pp. 119–128.
  * Blitzstein and Diaconis (2011) Blitzstein, J., and Diaconis, P. (2011), “A Sequential Importance Sampling Algorithm for Generating Random Graphs with Prescribed Degrees,” Internet Mathematics, 6(4), 489–522.
  * Chen (2007) Chen, Y. (2007), “Conditional Inference on Tables With Structural Zeros,” Journal of Computational and Graphical Statistics, 16(2), 445–467.
  * Chen et al. (2005) Chen, Y., Diaconis, P., Holmes, S. P., and Liu, J. S. (2005), “Sequential Monte Carlo Methods for Statistical Analysis of Tables,” Journal of the American Statistical Association, 100(469), 109–120.
  * Cormen (2009) Cormen, T. H. (2009), Introduction to algorithms MIT press.
  * Diaconis and Efron (1985) Diaconis, P., and Efron, B. (1985), “Testing for independence in a two-way table: new interpretations of the chi-square statistic,” The Annals of Statistics, pp. 845–874.
  * Diaconis and Sturmfels (1998) Diaconis, P., and Sturmfels, B. (1998), “Algebraic algorithms for sampling from conditional distributions,” Ann. Statist., 26(1), 363–397.
  * Eisinger and Chen (2017) Eisinger, R. D., and Chen, Y. (2017), “Sampling for Conditional Inference on Contingency Tables,” Journal of Computational and Graphical Statistics, 26(1), 79–87.
  * Fienberg (1969) Fienberg, S. E. (1969), “Preliminary Graphical Analysis and Quasi-Independence for two-way Contingency Tables,” Applied Statistics, pp. 153–168.
  * Girvan and Newman (2002) Girvan, M., and Newman, M. E. J. (2002), “Community Structure in Social and Biological Networks,” Proceedings of the national academy of sciences, 99(12), 7821–7826.
  * Goodman (1963) Goodman, L. A. (1963), “Statistical Methods for the Preliminary Analysis of Transaction Flows,” Econometrica: Journal of the Econometric Society, pp. 197–208.
  * Goodman (1965) Goodman, L. A. (1965), “On the Statistical Analysis of Mobility Tables,” American Journal of Sociology, 70(5), 564–585.
  * Goodman (1968) Goodman, L. A. (1968), “The Analysis of Cross-Classified Data: Independence, Quasi-Independence, and Interactions in Contingency Tables with or without Missing Entries: RA Fisher Memorial Lecture,” Journal of the American Statistical Association, 63(324), 1091–1131.
  * Grayson and Stillman (n.d.) Grayson, D. R., and Stillman, M. E. (n.d.), “Macaulay2, a Software System for Research in Algebraic Geometry,”, Available at http://www.math.uiuc.edu/Macaulay2/.
  * Guimerà et al. (2010) Guimerà, R., Stouffer, D. B., Sales-Pardo, M., Leicht, E. A., Newman, M. E. J., and Amaral, L. A. N. (2010), “Origin of compartmentalization in food webs,” Ecology, 91(10), 2941–2951.
  * Gustafsson (1980) Gustafsson, J.-E. (1980), “A solution of the conditional estimation problem for long tests in the Rasch model for dichotomous items,” Educational and Psychological Measurement, 40(2), 377–385.
  * Holland and Leinhardt (1981) Holland, P. W., and Leinhardt, S. (1981), “An Exponential Family of Probability Distributions for Directed Graphs,” Journal of the American Statistical Association, 76(373), 33–50.
  * Hudson et al. (2018) Hudson, L., Reuman, D., and Emerson, R. (2018), Cheddar: Analysis and Visualisation of Ecological Communities.  R package version 0.1-633.
  * Krause et al. (2003) Krause, A. E., Frank, K. A., Mason, D. M., Ulanowicz, R. E., and Taylor, W. W. (2003), “Compartments revealed in food-web structure,” Nature, 426(6964), 282.
  * McDonald et al. (2007) McDonald, J. W., Smith, P. W. F., and Forster, J. J. (2007), “Markov chain Monte Carlo exact inference for social networks,” Social Networks, 29(1), 127–136.
  * Milo et al. (2002) Milo, R., Shen-Orr, S., Itzkovitz, S., Kashtan, N., Chklovskii, D., and Alon, U. (2002), “Network Motifs: Simple Building Blocks of Complex Networks,” Science, 298(5594), 824–827.
  * Olesen et al. (2010) Olesen, J. M., Dupont, Y. L., O’Gorman, E., Ings, T. C., Layer, K., Melián, C. J., Trøjelsgaard, K., Pichler, D. E., Rasmussen, C., and Woodward, G. (2010), “From Broadstone to Zackenberg: Space, Time and Hierarchies in Ecological Networks,” in Advances in ecological research, Vol. 42 Elsevier, pp. 1–69.
  * Pearson (1904) Pearson, K. (1904), On the Theory of Contingency and its Relation to Association and Normal Correlation; On the General Theory of Skew Correlation and Non-Linear Regression Cambridge University Press.
  * Pimm and Lawton (1980) Pimm, S. L., and Lawton, J. H. (1980), “Are Food Webs Divided into Compartments?,” The Journal of Animal Ecology, pp. 879–898.
  * Rao et al. (1996) Rao, A. R., Jana, R., and Bandyopadhyay, S. (1996), “A Markov Chain Monte Carlo Method for Generating Random (0, 1) Matrices with Given Marginals,” Sankhya: The Indian Journal of Statistics, Series A (1961-2002), 58(2), 225–242.
  * Rapallo (2006) Rapallo, F. (2006), “Markov bases and structural zeros,” Journal of Symbolic Computation, 41(2), 164–172.
  * Rezende et al. (2009) Rezende, E. L., Albert, E. M., Fortuna, M. A., and Bascompte, J. (2009), “Compartments in a marine food web associated with phylogeny, body mass, and habitat structure,” Ecology Letters, 12(8), 779–788.
  * Roberts (2000) Roberts, J. M. (2000), “Simple methods for simulating sociomatrices with given marginal totals,” Social Networks, 22(3), 273–283.
  * Ryser (1963) Ryser, H. J. (1963), Combinatorial Mathematics Mathematical Association of America.
  * Snijders (1991) Snijders, T. A. B. (1991), “Enumeration and simulation methods for 0-1 matrices with given marginals,” Psychometrika, 56(3), 397–417.
  * Ulrich and Gotelli (2007) Ulrich, W., and Gotelli, N. J. (2007), “Null Model Analysis of Species Nestedness Patterns,” Ecology, 88(7), 1824–1831.
  * Verhelst (2008) Verhelst, N. D. (2008), “An Efficient MCMC Algorithm to Sample Binary Matrices with Fixed Marginals,” Psychometrika, 73(4), 705–728.
  * Wasserman and Faust (1994) Wasserman, S., and Faust, K. (1994), Social Network Analysis: Methods and Applications, Structural Analysis in the Social Sciences Cambridge University Press.
  * Williams and Martinez (2004) Williams, R. J., and Martinez, N. D. (2004), “Limits to trophic levels and omnivory in complex food webs: theory and data,” The American Naturalist, 163(3), 458–468.
  * Zhang and Chen (2013) Zhang, J., and Chen, Y. (2013), “Sampling for Conditional Inference on Network Data,” Journal of the American Statistical Association, 108(504), 1295–1307.

##  Appendix A Proofs

###### Proof of Lemma 3.3.

The identity decomposition immediately shows that the ‘only if’ part holds.
For the ‘if’ part, we must show reversibility of $Q$ with respect to $\pi$.
Without loss of generality (and for notational simplicity) assume $\pi$ and
$\mathbb{V}$ are dominated by one-dimensional Lebesgue measure. Reversibility
is then defined as

$\int_{A}\pi(x)Q(x,B)\mathop{}\\!\mathrm{d}x=\int_{B}\pi(x)Q(x,A)\mathop{}\\!\mathrm{d}x\quad\text{
for all }A,B\in\mathcal{\mathcal{B}}.$

Fix any $A$ and $B$ in $\mathcal{B}$, and define the densities $\\{f_{x}\\}$
as in Definition 3.2. Then

$\begin{split}\int_{A}\pi(x)Q(x,B)\mathop{}\\!\mathrm{d}x&=\int_{A}\int_{\mathcal{Z}}\int_{B}\pi(x)K_{z}(x,y)f_{x}(z)\mathop{}\\!\mathrm{d}y\mathop{}\\!\mathrm{d}z\mathop{}\\!\mathrm{d}x\\\
&=\int_{A}\int_{\mathcal{Z}}\int_{B}\pi(y)K_{z}(y,x)f_{y}(z)\mathop{}\\!\mathrm{d}y\mathop{}\\!\mathrm{d}z\mathop{}\\!\mathrm{d}x=\int_{B}\pi(x)Q(x,A)\mathop{}\\!\mathrm{d}x\end{split}$

as required. In the first step we have expressed the integral using densities.
In the second, we use reversibility of each $K_{z}$ with respect to $\pi$, and
the fact that for each $z\in\mathcal{Z}$ and $x\in\mathcal{X}$,
$f_{x}(z)=f_{y}(z)$ for $K_{z}(x,\cdot)$\- a.e. $y\in B$. A simple change of
variables then yields the result. ∎

###### Proof of Proposition 4.1.

Fix $G$ in $\mathcal{G}$. Define the second-order Markov chain $(Y_{n})_{n\geq
0}$, where $Y_{n}\colonequals(W_{n-1},W_{n},G_{n})$ and $G_{n}$ is defined as
follows. Let $G_{0}=G$, otherwise if $n$ is odd, let $G_{n}$ be the graph
obtained after $W_{n}W_{n-1}$ is removed from $E(G_{n-1})$. If $n$ is even,
$G_{n}$ is the graph obtained after $W_{n-1}W_{n}$ is added to $E(G_{n-1})$.
Define $\mathcal{Y}$ as the set of points reachable from $(*,W_{0},G)$ for
some $W_{0}$ in $\\{v:N_{G}(v)\neq\emptyset\\}$. Let
$D\colonequals(\mathcal{Y},E)$ be the digraph underlying this chain and define
$A$ as the subset of points $(v,u,G^{\prime})$ in $\mathcal{Y}$ for which
$G^{\prime}\in\mathcal{G}$. Let $T\colonequals\inf\\{n\geq 1\colon Y_{n}\in
A\\}$ be the first passage time of $A$. Proposition 4.1 is equivalent to
showing $\mathbb{E}(T)<\infty$. The following holds true, and will be shown at
the end of this proof.

$\text{From any }(v,u,G^{\prime})\in\mathcal{Y}\text{, there exists a simple
path to }A.$ (4)

We can bound the probability of traversing each edge in $D$ from below by some
constant $p>0$. Let $N$ denote the size of $\mathcal{Y}$. Suppose the chain is
at some state $y\notin A$. By (4), this implies the probability of hitting $A$
within the next $N$ steps is bounded from below by $p^{N}$. Hence,

$\begin{split}\mathbb{E}(T)&=\sum_{n=1}^{\infty}nP(T=n)\leq
N\sum_{k=1}^{\infty}kP(k-1<T/N\leq k)\leq
N\sum_{k=1}^{\infty}kp^{N}(1-p^{N})^{k-1}=Np^{-N}<\infty.\end{split}$

Thus it remains to show (4). We need the following observation repeatedly: for
any $(v_{1},u_{1},G_{1})(v_{2},u_{2},G_{2})\in E(D)$:

$(v^{*},u_{2},G_{2})(u_{2},u_{1},G_{1})\in E(D)\text{ if and only if
}(v^{*},u_{2},G_{2})\in\mathcal{Y}\text{ and }v^{*}\neq v_{2}.$ (5)

We now show (4). By definition, there exists a point $y_{0}:=(*,u_{0},G)$ and
a walk $y_{0}...y_{k}$ in $D$ such that $y_{k}=(v,u,G^{\prime})$. Given that
$V(D)$ is finite, continuing an arbitrary walk along $D$ from $y_{k}$ implies
we must eventually either return to $A$, or visit a graph already seen along
the walk. Denote the vertex visited at the $l$th step of this walk by
$y_{l}=(u_{l-1},u_{l},G_{l})$. If we revisit $A$ we are done, otherwise define

$n:=\inf\\{l>k:G_{l}=G_{m}\text{ for some }m<l\\}.$

The condition $G_{n}=G_{m}$ implies that $u_{n}=u_{m}$. Additionally
$u_{n-1}\neq u_{m-1}$, otherwise this would imply $G_{n-1}=G_{m-1}$, which
contradicts the definition of $n$. By (5),
$(u_{m-2},u_{m-1},G_{m-1})(u_{m-1},u_{m},G_{m})\in E(D)$ implies that
$(u_{n-1},u_{n},G_{n})(u_{m},W_{m-1},G_{m-1})\in E(D)$. Thus we can traverse
to $(u_{m},u_{m-1},G_{m-1})$. Iteratively applying (5) (which we can do as
$u_{l+1}\neq u_{l-1}$ for all $l\geq 0$) implies we can reach a state with
graph $G=G_{0}$, which must be in $A$, completing the proof of (4). ∎

###### Proof of Lemma 4.3.

We first show that $(K,\mathbb{V})$ is a symmetric decomposition. Fix any
$G\in\mathcal{G}$ and any $z\in\mathcal{Z}$, and let $w$ be a representative
of $z$. Let $p$ refer to the statement
‘$\mathbb{V}_{G}(z)=\mathbb{V}_{G^{*}}(z)$ for all $G^{*}$ for which
$K_{z}(G,G^{*})>0$’. It suffices to show that $p$ is true.

Consider a Markov chain with kernel $K_{z}$ and current state $G$. Suppose the
chain remains unchanged after one iteration of Algorithm 1. Then $p$ is true
trivially. Without loss of generality, suppose the swaps corresponding to $w$
are viable, and the chain moves to some $G^{*}\in\mathcal{G}$. Remark 4.2
implies swaps corresponding to $w^{r}$ are not viable. Since the swaps
corresponding to a sampled vertex sequence must be viable, $\mathbb{V}_{G}(z)$
is the probability $W=w$ given the chain is at $G$.

At $G^{*}$, the swaps corresponding to $w^{r}$ are viable. By an analogous
argument, it follows that $\mathbb{V}_{G^{*}}(z)$ is the probability $W=w^{r}$
given the chain is at $G^{*}$. One can deduce from Algorithm 1 that the
probability $W=w$ given the chain is at $G$ is equal to the probability
$W=w^{r}$ given the chain is at $G^{*}$. This holds because the degree
sequence is the same for either state.

We now show that each $K_{z}\in K$ is reversible with respect to the uniform
distribution. This is implied by detailed balance. Specifically, for each
$K_{z}\in K$ we show

$K_{z}(G,G^{*})=K_{z}(G^{*},G)\quad\text{ for all }G,G^{*}\in\mathcal{G}.$

Fix any $G$ and $G^{*}$. $K_{z}(G,G^{*})=1$ if and only if $K_{z}(G^{*},G)=1$,
because applying two iterations of a Markov chain with kernel $K_{z}$ from
some current state $G^{\prime}$, returns $G^{\prime}$. The result follows by
additionally observing that $K_{z}(G,G^{*})$ can only be zero or one. ∎

###### Proof of Proposition 4.4.

The proposition is implied by Lemma 4.3, and additionally showing the chain is
connected.

Fix any $G,G^{\prime}\in\mathcal{G}$, and suppose the current state of the
chain is $G$. Form a digraph $H$ as follows. For each vertex pair $uv$, if
$uv\in E(G)$ and $uv\notin E(G^{\prime})$, add a red edge $uv$ to $E(H)$. If
$uv\notin E(G)$ and $uv\in E(G^{\prime})$, add a blue edge $uv$ to $E(H)$.
Define an alternating cycle as a cycle whose edges are alternately red and
blue. $G$ and $G^{\prime}$ are equivalent if and only if $H$ has no edges.

Then $H$ is the union of a finite sequence of edge-disjoint alternating
cycles.

Fix any alternating cycle $v_{0}v_{1}...v_{k}v_{0}$ implied by this, ordered
so that the $v_{0}v_{1}$ is red. The Markov chain can sample
$W:=v_{0}v_{1}...v_{k}v_{0}$ with positive probability, yielding a new graph
$G^{\prime\prime}$, whilst removing all edges in $H$ corresponding to this
cycle. Iterate until $H$ has no more edges. ∎

###### Proof of Proposition 4.5.

For a given $\mathcal{F}$, the map from $\mathcal{G}$ to $B_{\mathcal{G}}$ is
injective, so the sampler can be thought of as a Markov chain ergodic with
respect to the uniform distribution on $B_{\mathcal{G}}$.

We briefly describe how to view the Markov chain as operating on
$B_{\mathcal{G}}$. An initial vertex $v_{j}$ is sampled from $V$. The chain
now samples $u_{i}$ from the out-neighborhood of $v_{j}$ and replaces the edge
$v_{j}u_{i}$ with $u_{i}v_{j}$. If $G$ is undirected, additionally switch
$v_{i}u_{j}$ with $u_{j}v_{i}$. Continue walking along the vertices of the
graph in this manner until the sampler returns to the initial vertex for the
first time.

Without loss of generality, suppose $G$ is directed. Fix $u_{i}\in S_{k}$ and
$v_{j}\in S_{l}$. If no edge is incident to $u_{i}$ and $v_{j}$ then
$ij\in\mathcal{F}\subseteq\tilde{\mathcal{F}}$. Otherwise if $k\neq l$, edges
between $S_{k}$ and $S_{l}$ are uniformly in one direction; say from $S_{k}$
to $S_{l}$. Suppose the Markov chain on $B$ traverses $u_{i}v_{j}$, replacing
it with $v_{j}u_{i}$. Returning to the initial vertex requires traversal of
$v_{j}u_{i}$. Hence, $u_{i}v_{j}$ can be flipped only an even number of times,
and the direction is unchanged. By Lemma 4.4, $ij\in\tilde{\mathcal{F}}$. If
$k=l$, $u_{i}v_{j}$ can be switched odd number of times, so
$ij\in\tilde{\mathcal{F}}$. The undirected case holds by an analogous
argument. ∎

###### Proof of Proposition 5.1.

It suffices to show connectedness. Define a metric
$d:\mathcal{G}\times\mathcal{G}\rightarrow\mathbb{N}_{0}$ on $\mathcal{G}$ by

$d(G,G^{\prime})\colonequals\sum_{u\in V}\sum_{v\in
V}|c_{G}(uv)-c_{G^{\prime}}(uv)|\quad\text{ for all
}G,G^{\prime}\in\mathcal{G}.$

Then $(\mathcal{G},d)$ is a metric space. Fix any two distinct graphs
$G,G^{\prime}\in\mathcal{G}$, and suppose the current state of the chain is
$G$. It suffices to show that one can construct a sampling step yielding a new
graph strictly closer to $G^{\prime}$ in this metric space.

Let $n_{uv}:=c_{G}(uv)-c_{G^{\prime}}(uv)$ for each vertex pair $uv$. We form
a multi-graph $H$ as follows. If $n_{uv}>0$, add $n_{uv}$ red copies of the
direction reversed edge $vu$ to $E(H)$, while if $n_{uv}<0$, add $-n_{uv}$
blue copies of $uv$ to $E(H)$. The graphs $G$ and $G^{\prime}$ are equivalent
if and only if $H$ has no edges. Define an alternating cycle in $H$ as a cycle
whose edges are alternately $red$ and $blue$.

It can be shown that $H$ can be expressed as the union of a finite number of
edge-disjoint alternating cycles.

Fix any such alternating cycle $v_{0}v_{1}...v_{l}v_{0}$ in $H$. Order the
cycle so that $v_{0}v_{1}$ is red. Letting $\mathbb{O}_{n}$ denote the set of
odd natural numbers less than or equal to $n$, we define

$k:=\inf\\{n\in\mathbb{O}_{l-2}:v_{n}v_{n+1}\notin\mathcal{F}\\}$

where we let $\inf\emptyset:=l$.

Under Algorithm 2, there is a positive probability of sampling the vertex
sequence $v_{0}v_{1}...v_{k}v_{0}$ given the chain is at $G$. Sampling
$\Delta=-1$ along this vertex sequence returns a new graph $G^{\prime\prime}$,
removing at least three edges from $H$ whilst adding at most one. Hence

$d(G^{\prime\prime},G^{\prime})\leq d(G,G^{\prime})-2.$

∎