Smaller Universes for Uniform Sampling of 0,1-matrices with fixed row and column sums

TL;DR

This paper investigates Markov chain algorithms for sampling 0,1 matrices with fixed row and column sums, aiming to understand why practical convergence is faster than theoretical bounds suggest.

Contribution

It demonstrates that ignoring node labels in Markov chains leads to faster convergence, providing insights into the discrepancy between theory and practice.

Findings

Projected Markov chains converge faster than original chains

Ignoring node labels simplifies sampling without loss of generality

Experimental results show much fewer steps are needed for good samples

Abstract

An important problem arising in the study of complex networks, for instance in community detection and motif finding, is the sampling of graphs with fixed degree sequence. The equivalent problem of generating random 0,1 matrices with fixed row and column sums is frequently used as a quantitative tool in ecology. It has however proven very challenging to design sampling algorithms that are both fast and unbiased. This article focusses on Markov chain approaches for sampling, where a close-to-random graph is produced by applying a large number N of small changes to a given graph. Examples are the switch chain and Curveball chain, which are both commonly used by practitioners as they are easy to implement and known to sample unbiased when N is large enough. Within theoretical research, much effort has gone into proving bounds on N. However, existing theoretical bounds are impractically large for most applications while experiments suggest that much fewer steps are needed to obtain a good sample. The contribution of this article is twofold. Firstly it is a step towards better understanding of the discrepancy between experimental observations and theoretically proven bounds. In particular, we argue that while existing Markov chain algorithms run on the set of all labelled graphs with a given degree sequence, node labels are unimportant in practice and are usually ignored in determining experimental bounds. We prove that ignoring node labels corresponds to projecting a Markov chain onto equivalence classes of isomorphic graphs and that the resulting projected Markov chain converges to its stationary distribution at least as fast as the original Markov chain. Often convergence is much faster, as we show in examples, explaining part of the difference between theory and experiments...

Figures

Equations

$$P_{x[y]}=P_{x^{\prime }[y]}$$

$$d_V(\mu ,\nu )=\underset{A\subset \Omega }{\max }|\mu (A)-\nu (A)|$$

$$\tau (ϵ)=\underset{x\in \Omega }{\max }\underset{T}{\min }\{T|d_V(P_x^t,\pi )\le ϵ\text{mbox}forallt>T\}.$$

$$(\mu f^{-1})([x]):=\mu (f^{-1}([x]))=\sum _{z\in [x]}\mu (z).$$

$$\overline{\mu }\overline{P}([y])=\sum _{[x]\in \overline{\Omega }}\overline{\mu }([x]){\overline{P}}_{[x][y]}=\sum _{[x]\in \overline{\Omega }}\sum _{z\in [x]}\mu (z){\overline{P}}_{[x][y]}=\sum _{z\in \Omega }\mu (z)P_{z[y]}$$

$$\overline{\mu }\overline{P}([y])=\sum _{y\in [y]}\sum _{z\in \Omega }\mu (z)P_{zy}=\sum _{y\in [y]}(\mu P)(y)=(\mu P)f^{-1}([y])$$

$$\overline{1_x}=\{\begin{array}{ll}\frac{1}{|[x]|}& \text{mbox}whenx\in [x]\\ 0& \text{mbox}otherwise.\end{array}$$

$$\hat{\tau }(ϵ):=\underset{[x]\in \overline{\Omega }}{\max }\underset{T}{\min }\{T|d_V(P_{\overline{x}}^t,\pi )\le ϵ\text{mbox}forallt>T\}.$$

$$d_V(g\overline{\mu },g\overline{\nu })=\frac{1}{2}\sum _{x\in \Omega }\left|\frac{\overline{\mu }([x])}{|[x]|}-\frac{\overline{\nu }([x])}{|[x]|}\right|=\frac{1}{2}\sum _{[x]\in \overline{\Omega }}|\overline{\mu }([x])-\overline{\mu }([x])|=d_V(\overline{\mu },\overline{\nu }).$$

$$\mu P(y)=\sum _{x\in \Omega }\mu (x)P_{yx}=\sum _{[x]\in \overline{\Omega }}\sum _{x\in [x]}\frac{\overline{\mu }([x])}{|[x]|}{P}_{yx}=\sum _{[x]\in \overline{\Omega }}\frac{\overline{\mu }([x])}{|[x]|}{P}_{y[x]}=\sum _{[x]\in \overline{\Omega }}\frac{\overline{\mu }([x])}{|[x]|}{\overline{P}}_{[y][x]}$$

$$\mu P(y)=\sum _{[x]\in \overline{\Omega }}\frac{\overline{\mu }([x])}{|[y]|}{\overline{P}}_{[x][y]}=\frac{(\overline{\mu }\overline{P})([y])}{|[y]|}=g(\overline{\mu }\overline{P})([y]).$$

$$\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)\text{mbox}by\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right)\text{mbox}orviceversa.$$

$$\left(\begin{array}{ccccccc}1& 1& 0& 1& 0& 1& 0\\ 0& 0& 1& 1& 1& 0& 0\end{array}\right)\text{mbox}canbereplacedby\left(\begin{array}{ccccccc}0& 1& 1& 1& 1& 0& 0\\ 1& 0& 0& 1& 0& 1& 0\end{array}\right)\text{mbox}.$$

$$P_{G[H]}^S=\sum _{K\in [H]}{P}_{GK}^S=\sum _{K\in N_G\cap [H]}{P}_{GK}^S$$

$$=\sum _{\sigma (K)\in N_{G^{\prime }}\cap [H]}{P}_{G^{\prime }\sigma (K)}^S=\sum _{\sigma (K)\in [H]}{P}_{G^{\prime }\sigma (K)}^S=P_{G^{\prime }[H]}^S.$$

$$1_{\overline{A}}(B)=\{\begin{array}{ll}\frac{1}{|[A]|}& \text{mbox}ifB\in [A]\\ 0& \text{mbox}otherwise.\end{array}$$

Full text

This article focusses on Markov chain approaches for sampling, where a close-to-random graph is produced by applying a large number $N$ of small changes to a given graph. Examples are the switch chain and Curveball chain, which are both commonly used by practitioners as they are easy to implement and known to sample unbiased when $N$ is large enough. Within theoretical research, much effort has gone into proving bounds on $N$. However, existing theoretical bounds are impractically large for most applications while experiments suggest that much fewer steps are needed to obtain a good sample.

The contribution of this article is twofold. Firstly it is a step towards better understanding of the discrepancy between experimental observations and theoretically proven bounds. In particular, we argue that while existing Markov chain algorithms run on the set of all labelled graphs with a given degree sequence, node labels are unimportant in practice and are usually ignored in determining experimental bounds. We prove that ignoring node labels corresponds to projecting a Markov chain onto equivalence classes of isomorphic graphs and that the resulting projected Markov chain converges to its stationary distribution at least as fast as the original Markov chain. Often convergence is much faster, as we show in examples, explaining part of the difference between theory and experiments. The speed-up comes from the fact that the projected Markov chain starts from a better (more uniform) initial distribution and runs on a smaller state space.

Secondly, when node labels are important, we show that faster sampling from the space of labelled graphs can be achieved by combining the projected Markov chain with a simple linear-time preprocessing step. As a result of independent interest, we prove that this approach removes the need to introduce ‘hexagonal moves’ in the switch and Curveball chain when sampling labelled simple directed graphs.

1 Introduction

The uniform sampling of graphs with fixed degree sequence has attracted a large research effort in network science [20, 21, 2, 22]. Samples of random networks are used to determine the significance of properties of real-world networks. For instance, to study the clustering in social networks [25], to understand which subgraphs form the important building blocks of a network [19] or to find out if a network is expected to be connected given its degree sequence [21]. The equivalent problem of uniformly sampling a 0,1 matrix with fixed row and column sums is considered as one of the most useful ‘null model approaches’ in ecology [11, 28].

Markov chain algorithms, such as the switch chain and Curveball chain, are a popular approach to the above sampling problem [22, 15, 28]. Here a graph is randomised by repeatedly applying small degree preserving changes. These algorithms are known to converge to the uniform distribution on the set of labelled graphs with a given degree sequence [2, 22, 30, 5]. Even though Markov chain algorithms are easy to implement and provide a flexible sampling framework, they have a serious drawback: in general it is unknown how many changes need to be applied to obtain a close-to-random graph.

Two completely separate communities appear to work on the Markov chain approach to this problem. The first one provides us with important theoretical insights to understand the problem in greater depth. It finds graph classes for which the Markov chain can be proven to be efficient [7, 12, 18, 10, 13, 1]. Unfortunately the polynomial upper bounds for the running time possess such large exponents, that they can never be used to draw a graph at random in a real world scenario. Furthermore, the graph classes considered (matrices with identical sums for rows or columns) barely occur in the real world.

The second community uses implementations of these algorithms in practice [19, 2, 28, 21], i.e. they need them for their research to create null models. Applied researchers often stop the Markov chains after a fixed number of steps, using some assessment to judge that this number was large enough such that the sampling was done according to an almost uniform distribution.

Even though existing theoretical results give impractically large limits, we are optimistic about the speed of this class of algorithms. Several experiments [24, 28, 6, 23] indicate that both the well known switch chain [26, 29, 15, 8, 22] and the lesser known Curveball algorithm [30, 28, 5] are quite fast. The only problem is that we can not prove how good they really are, i.e. we have a lack of theory.

In this article we offer a partial explanation for the discrepancy between theoretical and experimental bounds. In most applications, the network statistic of interest only depends on the structure of the network, i.e. it is a topological property. In practice, convergence of the statistic of interest is used as an indicator that a Markov chain has converged to its stationary distribution [2, 23]. This approach ignores node labels when judging the convergence of the Markov chain. We formally show that ignoring node labels corresponds to projecting a Markov chain onto equivalence classes of isomorphic graphs. We prove that the projected chain converges at least as fast as the original chain and give several examples where convergence is much faster. The speed up is due to sampling from a (often much) smaller state space. In some applications node labels are important, for instance when determining the number of expected edges between certain individuals or communities. We show that faster sampling can be achieved by combining the projected Markov chain with a linear-time preprocessing step. The resulting improved run-time is of clear benefit to practitioners. Furthermore, our contribution opens new pathways for theoretical research, by reducing the size of the state space. As a result of independent interest, we prove that this combination of a projected Markov chain and preprocessing step results in an ergodic Markov chain for all directed graphs, that is, it removes the need for ‘hexagonal moves’[22] (Theorem 13).

The remainder of this article is organised as follows. We start with a general description of projected Markov chains and prove that the mixing time of these chains is smaller than or equal to the mixing time of the original chain. We then briefly discuss well-known Markov chains used for the sampling of graphs: the switch and Curveball chain. We show that these Markov chains can be projected onto isomorphism classes of graphs and that this results in faster mixing. Furthermore we introduce a preprocessing step which allows us to speed-up the switch and Curveball chain. We give several explicit examples of the method. Finally we discuss limitations and potential extensions to this framework.

2 Applications

The following examples illustrate the main idea behind our speed-up of the switch and Curveball Markov chains.

Example 1**.**

Motif finding is a popular tool in network analysis [19]. A motif is defined as a small subgraph which appears significantly more frequently in an observed (real-world) graph than in randomly generated graphs. The switch chain is often used to generate such random graphs. It samples a graph uniformly at random from the space of all graphs with a given degree sequence. As a small example, Figure 1 illustrates the six different graphs with degree sequence $(2,2,3,2,1)$. Note that the three graphs on the left have the same topology and the three graphs on the right have a second different topology. For motif finding, it is not necessary to generate a sample from all six labelled graphs. We only need to know the probability with which we sample each of these two classes $G$ and $H$, as this allows us to compute the expected number of occurrences of a given subgraph. For instance in this small example, we find a graph with topology $G$ with probability $\nicefrac{{1}}{{2}}$ and we find a graph with topology $H$ also with probability $\nicefrac{{1}}{{2}}$. Hence, the expected number of triangles equals $0.5$.

Example 2**.**

We are interested in the buying behaviour of customers. We may represent this as a bipartite graph $G$ where the primary nodes represent customers, the secondary nodes represent products and the edges indicate that a customer has bought an item. Say we have observed four customers who each bought two items. In total there are four items and each item has been bought twice. We want to know the probability that the customers can be divided into two groups (see Figure 2(a)) while fixing the number of items bought per customer and the number of times each item is bought. To do so, we may generate samples of bipartite graphs with the degree sequence $k=((2,2,2,2),(2,2,2,2))$ using the switch or the Curveball chain, and estimate this probability. Both algorithms sample a graph uniformly at random (provided we run them for long enough) from the set of $90$ distinct labelled bipartite graphs with degree sequence $k$. Only $18$ of these states correspond to the situation where we can split the customers into two types. Hence we will find a probability close to $0.2$ provided we take a large enough sample and run the chains for long enough.

But the property of interest, if the customers can be split into two groups, is a topological property and does not depend on the labelling of the nodes. So in fact we are interested in sampling from a much smaller state space, that of unlabelled bipartite graphs with the given degrees. When removing the node labels, we find that there are only two distinct graphs (as illustrated in Figure 2(b)). We will later see that we can obtain the probabilities of sampling either of these two topologies by projecting the switch or Curveball chain. After projecting, fewer switches and trades are required to get close to the stationary distribution, largely due to the reduced number of graphs we are sampling from. Figure 2(c) illustrates this by showing the Markov chain of the projected switch chain.

These two applications show that it is often unnecessary to sample from the set of labelled graphs. Instead only the topology of the sampled networks is important. The next example shows that even when the labels are important, sampling can be sped up by making use of the projection mechanism.

Example 3**.**

We are studying a social network and want to know what the probability of two specific individuals being connected is given the number of connections of all individuals in the network. In this case, the labels of the nodes, i.e. who they represent is important. However, we can still benefit from the speed-up of sampling isomorphic graphs. To see this, we return to the toy-example in Figure 1. If we generate a sample with the projected switch chain we obtain a given representative of class $G$, say $G_{1}$, roughly half the time and a representative of $H$, say $H_{1}$, the rest of the time. To obtain a uniform sample from the six labelled graphs we can use this sample and apply a random node relabelling to all of the sampled graphs. Note that this relabelling has to preserve the degrees of the nodes. Hence in this example we simply choose a random permutation of the node labels $v_{1},v_{2}$ and $v_{4}$. Now we obtain a random sample from all labelled graphs with degree sequence $(2,2,3,2,1)$.

3 Projected Markov chains

We now introduce the framework of projected Markov chains. We show that the projection of a Markov chain has two equivalent interpretations. Firstly we can think of the projection as a Markov chain on equivalence classes. That is the chain $X_{0},X_{1},\dots,X_{t}$ becomes $\overline{X_{0}},\overline{X_{1}},\dots,\overline{X_{t}}$ and the state space is reduced in size: $\overline{\Omega}:=\Omega\backslash\sim$. Secondly, we can interpret the projected chain as running the original chain with an alternative starting distribution: instead of starting in a single state the chain starts from the uniform distribution on states in a single equivalent class.

In this article we consider discrete time Markov chains $\mathcal{M}=(\Omega,P)$ with finite state space $\Omega$ and transition matrix $P$. It is well-known that such a Markov chain converges to a unique stationary distribution $\pi$ if it is ergodic. If the Markov chain is time-reversible then this distribution satisfies $\pi_{i}P_{ij}=\pi_{j}P_{ji}$ for all $X_{i},X_{j}\in\Omega$. All Markov chains discussed in this article are finite, ergodic and time-reversible.

We will be interested in Markov chains that satisfy condition (1) as defined below, and their projected or lumped chain. We use the following result of Wilson [17, Theorem 2.5].

Lemma 4** (Projected Markov chains).**

Let $\mathcal{M}=(\Omega,P)$ be a Markov chain and let $\sim$ be an equivalence relation on $\Omega$ with equivalence classes $[x]\in\bar{\Omega}$. Assume that $P$ satisfies

[TABLE]

whenever $x\sim x^{\prime}$, and where $P_{x[y]}:=\sum_{z\in[y]}P_{xz}$. Then $\bar{\mathcal{M}}=(\bar{\Omega},\bar{P})$ with $\bar{P}_{[x][y]}:=P_{x[y]}$ is a Markov chain. $\bar{\mathcal{M}}$ is called the projected chain.

The stationary distribution of the projected Markov chain is proportional to the sizes of the equivalence classes. The following lemma was proved by Grone et al. [14, Propostion 3].

Lemma 5**.**

The projected chain $\overline{\mathcal{M}}=(\overline{\Omega},\overline{P})$ satisfies $\pi_{[x]}P_{[x][y]}=\pi_{[y]}P_{[y][x]}$ where $\pi_{[x]}=\sum_{x\in[x]}\pi(x)$ and hence has stationary distribution $\overline{\pi}=(\pi_{[x_{1}]},\dots,\pi_{[x_{n}]})$.

The mixing time of a Markov chain quantifies how quickly the chain approaches its stationary distribution. It is defined in terms of the variation distance between distributions.

Definition 6**.**

Let $\mu,\nu:\Omega\rightarrow[0,1]$ be probability distributions on $\Omega$. Their variation distance is defined as

[TABLE]

and can be shown to equal $\nicefrac{{1}}{{2}}\sum_{x\in\Omega}|\mu(x)-\nu(x)|$.

Let $P^{t}_{x}$ be the distribution of the Markov chain at time $t$ when started from state $x$, that is when started from the one-point distribution $\mathbf{1}_{x}$. This distribution $\mathbf{1}_{x}(y)$ equals $1$ when $x=y$ and [math] otherwise. When the complete transition matrix $P$ is known, the distribution $P^{t}_{x}$ can be computed by $t$ times right multiplying $\mathbf{1}_{x}$ with $P$, i.e. $P^{t}_{x}=\mathbf{1}_{x}P^{t}$. The mixing time of a Markov chain is defined as

[TABLE]

Informally, the mixing time is the maximum number of steps needed to get within distance $\epsilon$ of the stationary distribution regardless of the starting state. We now show that the mixing time $\overline{\tau}(\epsilon)$ of a projected chain is smaller or equal to the mixing time $\tau(\epsilon)$ of the original chain.

Lemma 7**.**

Let $\mathcal{M}=(\Omega,P)$ be a finite, ergodic Markov chain with stationary distribution $\pi$ and satisfying (1). Then $\tau(\epsilon)\geq\overline{\tau}(\epsilon)$.

Proof.

We will show that for any $x\in\Omega$, $t\in\mathbb{N}$ we have $d_{V}(P^{t}_{x},\pi)\geq d_{V}(\overline{P}^{t}_{[x]},\overline{\pi})$ which gives the result. Let $f:\Omega\rightarrow\overline{\Omega}$ be the function that maps a state $x$ to its equivalence class $[x]$. Furthermore, for a probability distribution $\mu$ on $\Omega$ let $\mu f^{-1}$ be the probability distribution on $\overline{\Omega}$ given by:

[TABLE]

Notice that $\overline{\pi}([x])=\sum_{x\in[x]}\pi(x)=\pi f^{-1}([x])$, i.e. $\overline{\pi}$ equals $\pi f^{-1}$. Furthermore the ‘one-point’ starting distribution $\overline{P}_{[x]}^{0}=\mathbf{1}_{[x]}$ equals $\mathbf{1}_{x}f^{-1}$ for any $x\in[x]$, that is $\overline{P}_{[x]}^{0}=P_{x}^{0}f^{-1}$. We now show that if $\overline{\mu}=\mu f^{-1}$ then also $(\overline{\mu}\overline{P})=(\mu P)f^{-1}$. We evaluate $(\overline{\mu}\overline{P})$ on a class $[y]$ and find

[TABLE]

where the first equality comes from writing out the matrix multiplication. The second equality uses $\overline{\mu}=\mu f^{-1}$ and the last equality uses the definition of $\overline{P}$. Next we obtain

[TABLE]

by using $P_{z[y]}=\sum_{y\in[y]}P_{zy}$ and again recognizing the matrix multiplication. Thus we now know that $\overline{P}_{[x]}^{t}=P_{x}^{t}f^{-1}$ for all $t$ and $x$. The proof now follows from [17, Lemma 7.9] where it is shown that $d_{V}(\mu,\nu)\geq d_{V}(\mu f^{-1},\nu f^{-1})$ for any $\mu$ and $\nu$. ∎

We may think of the projected chain $\overline{\mathcal{M}}$ as the original chain $\mathcal{M}$ started from the uniform distribution on an equivalence class $[x]$. That is with starting distribution

[TABLE]

We will denote by $P^{t}_{\overline{x}}$ the probability distribution of $\mathcal{M}$ at time $t$ with starting distribution $\overline{\mathbf{1}_{x}}$. We now show that the ‘mixing time’, $\hat{\tau}(\epsilon)$, of $\mathcal{M}$ when started from $\overline{\mathbf{1}_{x}}$ is exactly equal to that of $\overline{\mathcal{M}}$, i.e. starting $\mathcal{M}$ from $\overline{\mathbf{1}_{x}}$ is at least as fast as starting it from $\mathbf{1}_{x}$ (Lemma 3.3.). To do so we define

[TABLE]

Lemma 8**.**

Let $\mathcal{M}=(\Omega,P)$ be a finite, ergodic Markov chain with stationary distribution $\pi$ the uniform distribution and satisfying (1). Then $\hat{\tau}(\epsilon)=\overline{\tau}(\epsilon)$.

Proof.

To prove this statement we will show that $d_{V}(P^{t}_{\overline{x}},\pi)=d_{V}(\overline{P}^{t}_{[x]},\overline{\pi})$.

Let $\overline{\mu}:\overline{\Omega}\rightarrow[0,1]$ be a probability distribution on $\overline{\Omega}$, we define a probability distribution $g\overline{\mu}$ on $\Omega$ by $g\overline{\mu}(x):=\nicefrac{{\overline{\mu}([x])}}{{|[x]|}}$. Now clearly

[TABLE]

Furthermore $P^{0}_{\overline{x}}=\overline{\mathbf{1}_{x}}=g\mathbf{1}_{[x]}=gP^{0}_{[x]}$. We next show that if a probability distribution $\mu$ on $\Omega$ can be written as $\mu=g\overline{\mu}$ then $\mu P=g(\overline{\mu}\overline{P})$. We derive

[TABLE]

Now using detailed balance for the the projected chain $\frac{|[x]|}{|\Omega|}\overline{P}_{[x][y]}=\frac{|[y]|}{|\Omega|}\overline{P}_{[y][x]}$ we obtain

[TABLE]

This implies $P^{t}_{\overline{x}}=g\overline{P}^{t}_{[x]}$ for all $t$ and $x$ and thus $d_{V}(P^{t}_{\overline{x}},\pi)=d_{V}(g\overline{P}^{t}_{[x]},g\overline{\pi})=d_{V}(\overline{P}^{t}_{[x]},\overline{\pi})$. ∎

In this framework we are able to compare the mixing time of Markov chains directly in terms of the variation distance. In the theoretical literature the mixing time of a Markov chain is often bounded using the spectral gap upper and lower bound. The spectral gap of an ergodic finite Markov chain $\mathcal{M}=(\Omega,P)$ is defined as follows. Let us denote the left eigenvalues of the transition matrix $P$ by $\lambda_{N}\leq\dots\leq\lambda_{1}$. It is a classical result that $-1<\lambda_{N},\lambda_{1}=1$ and $\lambda_{2}<1$ (see for example [17, Lemma 12.1.].) Let us denote by $\lambda^{*}:=\max\{|\lambda_{i}|:1<i\leq N\}$ the eigenvalue with largest absolute value smaller than 1. The value $1-\lambda^{*}$ is often called spectral gap. The main effort for bounding the mixing time of a Markov chain often goes into finding an expression or a bound for the spectral gap. We will later show that our method can sometimes be used to find an explicit expression for the spectral gap of the projected chain (i.e. Example 14).

4 Markov chains for sampling graphs

There are two commonly used Markov chain algorithms designed for the sampling of graphs; we briefly discuss the switch chain $\mathcal{M}^{S}$ and the Curveball chain $\mathcal{M}^{C}$. Both exist in several flavours: for the sampling of bipartite graphs, simple directed graphs, directed graphs and simple undirected graphs. We describe the algorithms in terms of bipartite graphs, in some sense the most general case. The bi-adjacency matrix of a bipartite graph is an $n\times n^{\prime}$ matrix where $n$ is the number of primary nodes and $n^{\prime}$ the number of secondary nodes. The $(i,j)$-th entry of this matrix equals one if there is an edge between primary node $p_{i}$ and secondary node $s_{j}$, otherwise it equals zero. We will describe the switch chain and Curveball chain as algorithms that randomise the bi-adjacency matrix of a bipartite graph while keeping its row and column sums fixed. Note that this corresponds exactly to sampling a bipartite graph with fixed degrees.

A Markov chain $\mathcal{M}=(\Omega,P)$ is described by its state space $\Omega$ and its transition probabilities $P$. Given a binary matrix $A$ with row and columns sums $k=((r_{1},\dots,r_{n}),(c_{1},\dots,c_{n^{\prime}}))$, both $\mathcal{M}^{S}$ and $\mathcal{M}^{C}$ have as their state space the set of all binary $n\times n^{\prime}$ matrices with row and column sums $k$. We denote this state space by $\Omega_{k}$.

In practice, both Markov chains are started from a specific state $X_{0}$ (a binary matrix) and each transition corresponds to making a small change to the current state $X_{i}$. The switch chain applies switches: replacing a submatrix

[TABLE]

The Curveball algorithm applies trades: a trade randomly exchanges ‘ones’ between two selected rows. For instance, the rows

[TABLE]

The top row had ‘tradeable ones’, i.e. where the bottom row equals zero, in columns 1,2 and 6 before the trade and the bottom row had tradeable ones in columns 3 and 5. A trade corresponds to randomly selecting three columns from these five available columns. In this example columns 2,3 and 5 were chosen. Notice that also columns 1,2 and 3 could be chosen which corresponds to applying a switch. One trade can apply multiple switches at once.

There are several versions of the switch chain known in the literature [15], here we use the chain which randomly selects two non-zero matrix entries, and applies a switch if possible, i.e. if the $(2\times 2)$-submatrix formed by the rows and columns corresponding to these entries are as in equation (2). The Curveball chain that we use in this article proceeds by selecting a pair of rows, $i$ and $j$ at random (with probability $\binom{n}{2}^{-1}$) and applying a trade with probability $\binom{s_{i}+s_{j}}{s_{i}}^{-1}$ where $s_{i}$ is the number of columns where row $A_{i}$ equals $1$ and row $A_{j}$ equals [math] and $s_{j}$ is the number of columns where row $A_{i}$ equals [math] and row $A_{j}$ equals $1$, see [30, 28] for more information.

To define the projected switch and Curveball chain we use the framework discussed in Section 3. The equivalence relation that we will use is that of graph isomorphism, that is graphs are equivalent when they have the same topology. Formally, two (bipartite) graphs $G$ and $H$ are isomorphic if there is a bijective map $\sigma$ between their nodes such that edges are preserved, i.e. such that the edge $(\sigma(u),\sigma(v))$ is present in $H$ if and only if the edge $(u,v)$ is present in $G$. Note that, graph isomorphisms are degree preserving maps. In general it can be hard to decide if two graphs are isomorphic, but creating isomorphic graphs is simple: pick a random labelling of the nodes of a given graph $G$ to obtain an isomorphic graph $H$.

We say two bipartite graphs $G\sim H\in\Omega$ are equivalent if and only if they are isomorphic as bipartite graphs, i.e. if there exists a graph isomorphism $\sigma=(\sigma_{1},\sigma_{2})$, where $\sigma_{1}$ ($\sigma_{2}$) maps the primary (secondary) nodes of $G$ to the primary (secondary) nodes of $H$. We need to show that the transition matrices $P^{S}$ and $P^{C}$ of the switch chain and the Curveball chain respectively are of form (1) under this equivalence relation. Intuitively, this holds because applying a specific switch or trade to isomorphic graphs will lead to isomorphic graphs. This implies that the probability of ending up in a given equivalence class is equal for graphs that are in the same equivalence class. In Figure 3 we illustrate this for the switch chain of a simple directed graph on five nodes. The following lemma is a formal statement of the above argument in the setting of bipartite graphs, but can easily be generalized to other graph classes.

Lemma 9**.**

Let $\mathcal{M}^{S}=(\Omega,P^{S})$, $\mathcal{M}^{C}=(\Omega,P^{C})$ be the switch and Curveball chain for a given bipartite degree sequence. Both $P^{S}$ and $P^{C}$ are of the form (1).

Proof.

Let $G=(P,S,E)$ and $G^{\prime}=(P,S,E^{\prime})$ be bipartite graphs in $\Omega$ such that $G\sim G^{\prime}$. That is, there exist degree preserving isomorphisms $\sigma=(\sigma_{1},\sigma_{2})$, with $\sigma_{1}:P\rightarrow P$ and $\sigma_{2}:S\rightarrow S$ such that for any $e=\{p,s\}\in E$ we have $\{\sigma(p),\sigma(s)\}\in E^{\prime}$. We will write $\sigma(G):=G^{\prime}$. Let $H$ be switch-adjacent to $G$ with regard to a specific switch: i.e. $H=(V,E\backslash\{\{p_{i},s_{j}\},\{p_{k},s_{l}\}\}\cup\{\{p_{i},s_{l}\},\{p_{k},s_{j}\}\}$. Then the graph $\sigma(H)=(V,E^{\prime}$ $\backslash\{\{\sigma(p_{i}),\sigma(s_{j})\},\{\sigma(p_{k}),\sigma(s_{l})\}\}$ $\cup\{\{\sigma(p_{i}),\sigma(s_{l})\},\{\sigma(p_{k}),\sigma(s_{j})\}\})$ is switch-adjacent to $\sigma(G)$ and furthermore $H\sim\sigma(H)$ by definition. That is, all graphs that are switch-adjacent to $G$ are isomorphic to graphs that are switch adjacent to $\sigma(G)$ under the graph isomorphism $\sigma$ and $P^{S}_{GH}=P^{S}_{\sigma(G)\sigma(H)}$. We now write $N_{G}$ for the set of graphs that is switch-adjacent to $G$.

For any equivalence class $[H]$ and $G\sim G^{\prime}$ with graph isomorphism $\sigma$ we obtain

[TABLE]

[TABLE]

Since any trade equals a sequence of $k$ switches the result immediately follows for $P^{C}$. ∎

In particular, the projections $\overline{\mathcal{M}}^{S}$ and $\overline{\mathcal{M}}^{C}$ of the switch chain and the Curveball chain with respect to $\sim$ are well-defined. Lemma 5 now tells us that the stationary distribution of these projected chains is proportional to the size of the equivalence classes. Hence, if we generate a sample using the projected chains we obtain each topology with the correct probability: the probability of sampling a graph in $\Omega$ with the given topology. In practice, when we are running experiments where we are only interested in the topology of the sampled networks, we could argue that we already use the projected chain. To illustrate this, we elaborate on Example 2 with respect to the switch chain.

Example 10**.**

In Example 2 we wanted to know the probability that a bipartite graph with degrees $k=((2,2,2,2),(2,2,2,2))$ is disconnected. Due to the small size of $\Omega$ and $\overline{\Omega}$ we can explicitly compute $P^{S}$ and $\overline{P}^{S}$ (see Figure 2(c)) and hence determine the mixing time for a given $\epsilon$. For $\epsilon=0.001$ we find $\tau(\epsilon)=28$ and $\overline{\tau}(\epsilon)=6$. This means that after $28$ switches, the probability of obtaining any specific graph $G$ of the $90$ distinct labelled graphs with degrees $k$ is roughly $\frac{1}{90}$. However, the probability of obtaining a graph with topology $G$ is already roughly $\frac{1}{5}$ after $6$ steps in $\overline{P}^{s}$.

In general, if we know the mixing time of the switch chain theoretically and run the chain $N$ times for $\tau(\epsilon)$ steps to obtain a sample of size $N$, we could be taking much longer than necessary because the property of interest (and any other topological property) already converges after $\overline{\tau}(\epsilon)$ steps.

In fact, for any property of interest (motifs, number of connected components, diameter) we may try to project the chain to an even smaller state space. To see this, let the property be given as a function $f:\Omega\rightarrow\mathbb{R}$. We say two states $s,s^{\prime}$ are equivalent if and only if $f(s)=f(s^{\prime})$. Hence, $f$ decomposes $\Omega$ in equivalence classes. If this equivalence relation satisfies 3.1 we obtain a projected chain with smaller mixing time than the original chain. Moreover, since many properties of interest are topological measurements, we know that we can always project the chain down to the isomorphism classes (even if the equivalence relation given by $f$ does not satisfy 3.1). In practice the convergence of the switch chain is often estimated through the convergence of the property of interest. That is, the estimated convergence is the convergence of a projected chain. This can explain part of the difference between theoretically proven bounds and experimentally observed bounds.

In certain applications, e.g. Example 3, node labels are important. We now show how we can speed up sampling by adding a preprocessing step to the switch or Curveball chain. We define the preprocessing step in terms of the bi-adjacency matrix.

Definition 11** (Preprocessing step).**

Let $A$ be a binary $(n\times n^{\prime})$-matrix in $\Omega_{k}$ for fixed row and column sums $k$. Let $R$ be the set of all permutations $\rho:$ $(1,\dots,n)$ $\rightarrow$ $(1,\dots,n)$ such that the row sums $\sum_{j=1}^{n^{\prime}}A_{\rho(i)j}=r_{i}$ for all $i$ and let $S$ be the set of all permutations $\sigma:$ $\{1,\dots,n^{\prime}\}$ $\rightarrow$ $\{1,\dots,n^{\prime}\}$ such that the column sums $\sum_{i=1}^{n}A_{i\sigma(j)}=c_{j}$. We define a preprocessing step which randomly selects a $\rho\in R$ and $\sigma\in S$ to form the matrix $B_{ij}=A_{\rho(i)\sigma(j)}$. This preprocessing step can be implemented by choosing a random order for sets of rows with equal row sum and a random order for sets of columns with equal column sum in $O(n+n^{\prime})$ [9].

We can think of including the preprocessing step as starting the Markov chain from the uniform distribution on a single equivalence class. That is with starting distribution

[TABLE]

Importantly, adding the preprocessing step to the Curveball chain or the switch chain, does not change their convergence to the uniform distribution since the stationary distribution is independent of the starting distribution of the chains. However, it may speed up the convergence.

For undirected and directed graphs we can similarly introduce a preprocessing step.

Definition 12** (Preprocessing step).**

Let $G$ be a (directed) graph with $n$ nodes and adjacency matrix $A$. Let $P$ be a partition of the nodes with equal degree (equal in- and out-degree in case of directed graphs). Let $R$ be the set of all permutations $\rho:$ $\{1,\dots,n\}$ $\rightarrow$ $\{1,\dots,n\}$ which respect the partition $P$, i.e. nodes are permuted within each partition separately. We define the following preprocessing step: randomly select a permutation $\rho\in R$ and apply it to the nodes of $G$, then return the resulting adjacency matrix.

This finishes our discussion of how the framework of projected Markov chains may improve the mixing time in applications. We finish this section by proving that the inclusion of a preprocessing step removes the need for ‘hexagonal moves’ in the switch and Curveball chain for directed graphs as proved by Rao et al [22]. The intuition behind this proof is simple: the reason that the state graph of $M^{S}$ and $M^{C}$ is disconnected for some directed graphs is that the direction of certain directed triangles can not be reversed [3], which can be achieved with our preprocessing step.

Theorem 13**.**

The switch chain $M^{S}$ and Curveball chain $M^{C}$ have connected state space when including the preprocessing step in Definition 12.

Proof.

It has been shown that the above chains sample uniformly when including a pre-sampling step which assigns orientations to each ‘induced cycle set’ randomly [3]. Our suggested preprocessing step does this too, since it permutes nodes with equal in and out degree, and it was shown that three nodes that form an induced cycle set always have the same in- and out-degree [16]. Hence, our preprocessing also re-orients all induced cycle sets. ∎

In the next section we will discuss several examples where preprocessing improves the mixing time. Furthermore, we show that for certain families of graphs, the size of the state space $\overline{\Omega}$ is constant whereas the size of $\Omega$ grows quadratically in the number of rows of the matrices.

5 Examples: Smaller Universes and Faster Sampling

In this section we discuss several examples where the state space is reduced significantly by using the framework of projected Markov chains. We start with an example where a quadratically growing state space is reduced to constant size. Furthermore, we show that our method makes it possible to explicitly compute the spectral gap of the projected chain, which would be very complicated for the original chain since its state space is growing in size. The spectral gap is often used to bound the mixing time of Markov chains [17, Theorem 12.3].

Example 14**.**

Let $r_{n}$ be a vector of length $n$ where all entries are equal to two, i.e. $r_{n}=(2,\dots,2)$ and let $c_{n}$ be the vector $(n-1,n-1,1,1)$. Then $k_{n}=(r_{n},c_{n})$ are valid degrees for a bipartite graph with $n$ primary nodes and four secondary nodes $(s_{1},s_{2},s_{3},s_{4})$ whenever $n\geq 2$. Let $G$ be the bipartite graph consisting of $K_{n-1,2}$ and $K_{1,2}$, i.e. $s_{1}$ and $s_{2}$ are connected to the same $n-1$ primary nodes and $s_{3}$ and $s_{4}$ are connected to the same primary node. This graph belongs to the state space $\Omega_{k_{n}}$. The equivalence class $[G]$ of graph $G$ has size $n$ since there are $n$ choices for the label of the primary node in the disconnected $K_{1,2}$. There is one other equivalence class, the class $[H]$, the graphs in this class are connected, $s_{1}$ and $s_{2}$ share $n-2$ neighbours and are each connected to a single additional primary node. One of these two nodes is furthermore connected to $s_{3}$ and the other to $s_{4}$. The class $[H]$ has size $2n(n-1)$. Hence the size of the state space $\Omega_{k^{n}}$ equals $n(2n-1)$ and grows quadratically in terms of $n$. On the other hand, the size of the state space of the projected chain is independent of $n$ and always equal to two.

We are able to explicitly compute the spectral gap for the projected switch chain. The transition probability $P^{S}_{[G][H]}$ equals the probability of selecting an edge in $K_{n-1,2}$ and an edge in $K_{1,2}$ and hence equals $\binom{m}{2}^{-1}4(n-1)$ with $m=2n$, the total number of edges. The transition probability $P^{S}_{[H][G]}$ equals $\binom{2n}{2}^{-1}2$ which can be seen by inspecting a specific graph $H\in[H]$. Let $H$ be the graph where $s_{1}$ is connected to $p_{1},\dots,p_{n-1}$ and $s_{2}$ is connected to $p_{1},\dots,p_{n-2},p_{n}$. Furthermore $p_{n-1}$ is connected to $s_{3}$ and $p_{n}$ is connected to $s_{4}$. The only switches that will give us a graph in $[G]$ are $(p_{n-1},s_{1})$ with $(p_{n},s_{4})$ and $(p_{n-1},s_{3})$ with $(p_{n},s_{4})$. The eigenvalues of $\overline{P}^{S}$ can be symbolically computed and equal $\lambda_{1}=1$, $\lambda_{2}=1-\nicefrac{{2}}{{n}}$. For $n\geq 2$, the spectral gap is given by $1-\lambda_{2}=\nicefrac{{2}}{{n}}$, leading to an upper bound of the mixing time of $O(n\log(n/\epsilon))$ [27].

In the next example we show that for a family of matrices, using our preprocessing step reduces the size of the state space impressively: from growing exponentially in the number of columns to always consisting of a single state.

Example 15**.**

Let $G_{l}$ be the bipartite graph with two primary nodes with degree $l$ and $2l$ secondary nodes with degree $1$. The size of the state space $\Omega_{l}$ of graphs with these degrees is $\binom{2l}{l}$ and grows exponentially in $l$. For the Curveball chain, a trade from $G_{l}$ reaches all realisations of the degree sequence with probability $\binom{2l}{l}^{-1}$. Now consider the graphs $G_{l}$ and $H_{l}$ where the neighbours of node $p_{1}$ in $G_{l}$ are the neighbours of node $p_{2}$ in $H_{l}$ and vice versa. To go from state $G_{l}$ to state $H_{l}$ with the switch chain, at least $l$ switches are needed. In other words, at least $l$ steps of the switch chain are needed to reach every state with positive probability. This is a clear example where the Curveball chain is the better choice in terms of mixing time. Finally, if we use our preprocessing step and project either chain with respect to the equivalence relation $\sim$, we find that only a single state remains since all states are isomorphic. Hence we only need to apply the preprocessing step and are left with a uniformly sampled labelled graph.

6 Discussion and conclusion

In this article we introduce a projected version of the switch and Curveball Markov chain where only the topology of the resulting graph is used. In many applications this is the main feature of interest and projecting can significantly reduce the size of the state space and hence improve the mixing time of a Markov chain. We furthermore introduce a preprocessing step that can be used in combination with the projected chain to obtain a random sample from the set of labelled graphs.

Clearly we can find examples where projecting does not alter the size of the state space, that is, any graph where each node has a unique degree leads to equivalence classes of size one and hence no reduction in the size of the state space. However, such degree sequences have smaller state spaces to begin with, exactly due to the absence of this redundant symmetry. In [4] an interesting relation is discussed between majorization of degree sequences and the size of the state space. This could turn out to be the key to showing all state spaces are small after projection.

Most theoretical bounds on the mixing time of the switch chain give a bound on the spectral gap. The spectral gap of a projected chain is smaller than or equal to the spectral gap of the original chain since it has been proven that the eigenvalues of a projected chain are a subset of the eigenvalues from the original chain [17, Theorem 12.8.(ii)]. Focusing the analysis of spectral gap bounds for projected chains would make an interesting area of future research.

Acknowledgement

The authors would like to thank Pieter Kleer for insightful comments and discussions on the content of this work.