A Fast MCMC for the Uniform Sampling of Binary Matrices with Fixed Margins

TL;DR

This paper introduces the Rectangle Loop algorithm, a Markov chain Monte Carlo method that efficiently and uniformly samples binary matrices with fixed margins, outperforming traditional swap algorithms especially for large or sparse/dense matrices.

Contribution

The paper presents the Rectangle Loop algorithm, a novel MCMC approach that improves efficiency and theoretical properties over existing swap algorithms for binary matrix sampling.

Findings

Rectangle Loop algorithm is more efficient than swap algorithm.

Theoretically superior in Peskun's order.

Empirical results confirm improved performance for large or sparse/dense matrices.

Abstract

Uniform sampling of binary matrix with fixed margins is an important and difficult problem in statistics, computer science, ecology and so on. The well-known swap algorithm would be inefficient when the size of the matrix becomes large or when the matrix is too sparse/dense. Here we propose the Rectangle Loop algorithm, a Markov chain Monte Carlo algorithm to sample binary matrices with fixed margins uniformly. Theoretically the Rectangle Loop algorithm is better than the swap algorithm in Peskun's order. Empirically studies also demonstrates the Rectangle Loop algorithm is remarkablely more efficient than the swap algorithm.

Figures

Tables

Table 1. Table 1: Occurrence Matrix occurrence matrix of the finches on the Galapagos islands.

	Island
Finch	A	B	C	D	E	F	G	H	I	J	K	L	M	N	O	P	Q
1	0	0	1	1	1	1	1	1	1	1	0	1	1	1	1	1	1
2	1	1	1	1	1	1	1	1	1	1	0	1	0	1	1	0	0
3	1	1	1	1	1	1	1	1	1	1	1	1	0	1	1	0	0
4	0	0	1	1	1	0	0	1	0	1	0	1	1	0	1	1	1
5	1	1	1	0	1	1	1	1	1	1	0	1	0	1	1	0	0
6	0	0	0	0	0	0	0	0	0	0	1	0	1	0	0	0	0
7	0	0	1	1	1	1	1	1	1	0	0	1	0	1	1	0	0
8	0	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0
9	0	0	1	1	1	1	1	1	1	1	0	1	0	0	1	0	0
10	0	0	1	1	1	1	1	1	1	1	0	1	0	1	1	0	0
11	0	0	1	1	1	0	1	1	0	1	0	0	0	0	0	0	0
12	0	0	1	1	0	0	0	0	0	0	0	0	0	0	0	0	0
13	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1

Table 2. Table 2: All possible 3 × 3 3 3 3\times 3 binary matrices with [ 1 , 2 , 1 ] 1 2 1 [1,2,1] as both row and column sums

A	B	C	D	E

Table 3. Table 3: Transition matrices for swap algorithm (left), Curveball algorithm (middle), Rectangle Loop algorithm (right)

Swap	Curveball	Rectangle Loop

Table 4. Table 4: The comparison between swap algorithm and Rectangle Loop algorithm. Each algorithm is implemented 10000 10000 10000 iterations on 100 × 100 100 100 100\times 100 matrices with different filled portions. The third column records the number of successful swaps among the 10000 10000 10000 iterations, the last column records the average time per swap, respectively.

Method	Filled portion	Number of swaps	Time per swap (/s)
Rectangle Loop	$1\%$	$586$	$1.18\times {10}^{-5}$
Swap		$8$	$3.67\times {10}^{-4}$
Rectangle Loop	$5\%$	$977$	$5.30\times {10}^{-6}$
Swap		$42$	$3.52\times {10}^{-5}$
Rectangle Loop	$10\%$	$1838$	$3.23\times {10}^{-6}$
Swap		$156$	$1.25\times {10}^{-5}$
Rectangle Loop	$20\%$	$3271$	$2.64\times {10}^{-6}$
Swap		$509$	$5.68\times {10}^{-6}$
Rectangle Loop	$30\%$	$4222$	$2.10\times {10}^{-6}$
Swap		$803$	$5.06\times {10}^{-6}$
Rectangle Loop	$40\%$	$4794$	$1.27\times {10}^{-6}$
Swap		$1160$	$4.98\times {10}^{-6}$
Rectangle Loop	$50\%$	$5080$	$1.37\times {10}^{-6}$
Swap		$1271$	$5.36\times {10}^{-6}$

Equations

$$m_{i,j}=\{\begin{array}{l}1,\text{if there is an edge connecting}{u}_i\text{and}{v}_j,\\ 0,\text{otherwise}.\end{array}$$

$$\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right),\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right).$$

$$\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)\text{or}\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right),$$

$$S_{a-b}\doteq \{k:A_t(a,k)=1,A_t(b,k)=0\}$$

$$S_{b-a}\doteq \{l:A_t(a,l)=0,A_t(b,l)=1\},$$

$$A_{t+1}(a,l)=\{\begin{array}{l}1\text{if}l\in V\\ 0\text{if}l\in S_{a-b}\cup S_{b-a}\setminus V\\ A(a,l)\text{otherwise}\end{array}$$

$$A_{t+1}(b,k)=\{\begin{array}{l}1\text{if}k\in S_{a-b}\cup S_{b-a}\setminus V\\ 0\text{if}k\in V\\ A(b,k)\text{otherwise}\end{array}$$

$$\mathbb{P}(A,B)=\{\begin{array}{ll}\frac{1}{\left(\genfrac{}{}{0ex}{}{m}{2}\right)\left(\genfrac{}{}{0ex}{}{n}{2}\right)}& \text{If A and B are swappable},\\ 1-\frac{s(A)}{\left(\genfrac{}{}{0ex}{}{m}{2}\right)\left(\genfrac{}{}{0ex}{}{n}{2}\right)}& \text{If B=A},s(A)\doteq \#\{\text{C: C and A are swappable}\}\\ 0& \text{Otherwise}\end{array}$$

$$\frac{1}{\left(\genfrac{}{}{0ex}{}{m}{2}\right)\left(\genfrac{}{}{0ex}{}{n}{2}\right)}.$$

$$\min \{1,\frac{\pi (B)\mathbb{P}(B,A)}{\pi (A)\mathbb{P}(A,B)},\}$$

$$\mathbb{P}(A,B)=\mathbb{P}(B,A),$$

$${\mathbb{P}}_1(x,A)\ge {\mathbb{P}}_2(x,A)$$

$$\frac{1}{mn}\cdot \frac{1}{n-r_{i_1}}\cdot \frac{1}{c_{j_2}}.$$

$$\mathbb{P}_{r}(A,B)=\frac{1}{mn}\bigg{(}\frac{1}{n-r_{i_{1}}}\cdot\frac{1}{c_{j_{2}}}+\frac{1}{c_{j_{2}}}\cdot\frac{1}{n-r_{i_{2}}}+\frac{1}{n-r_{i_{2}}}\cdot\frac{1}{c_{j_{1}}}+\frac{1}{c_{j_{1}}}\cdot\frac{1}{n-r_{i_{1}}}\bigg{)}.$$

$$\mathbb{P}_{r}(B,A)=\frac{1}{mn}\bigg{(}\frac{1}{c_{j_{1}}}\cdot\frac{1}{n-r_{i_{2}}}+\frac{1}{n-r_{i_{2}}}\cdot\frac{1}{c_{j_{2}}}+\frac{1}{c_{j_{2}}}\cdot\frac{1}{n-r_{i_{1}}}+\frac{1}{n-r_{i_{1}}}\cdot\frac{1}{c_{j_{1}}}\bigg{)}.$$

$${\mathbb{P}}_r(A,B)=$$

$$\frac{1}{mn}\cdot \frac{1}{n-r_{i_2}}\cdot \frac{1}{c_{j_1}}+\frac{1}{mn}\cdot \frac{1}{c_{j_1}}\cdot \frac{1}{n-r_{i_1}}.$$

$$\underset{A\in {\Sigma }_{\text{r},\text{c}}}{\max }\Vert {\mathbb{P}}^k(A,\cdot )-\pi {\Vert }_{\text{TV}}=\frac{1}{2}\underset{A\in {\Sigma }_{\text{r},\text{c}}}{\max }\sum _{B\in {\Sigma }_{\text{r},\text{c}}}|{\mathbb{P}}^k(A,B)-\pi (B)|,$$

$${\bar{S}}^2(A)=\frac{1}{m(m-1)}\sum _{i\ne j}{s}_{ij}^2,$$

Full text

1 Introduction

The problem of sampling binary matrices with fixed row and column sums has attracted much attention in numerical ecology. In ecological studies, the binary matrix is called occurrence matrix. Rows usually corresponds to species, the columns, to locations. For example, the binary matrix shown on Table 1 is known as “Darwin’s Finch” dataset, which comes from Darwin’s studies of the finches on the Galapagos islands (an archipelago in the East Pacific). The matrix represents the presence/absence of 13 species of finches in 17 islands. A $``1"$ or $``0"$ in entry $(i,j)$ indicates the presence or absence of species $i$ at island $j$. It is clear from Table 1 that some pairs of species tend to occur together (for example, species 9 and 10) while some other pairs tend to be disjoint. Therefore, it is of our interest to investigate whether the cooperation/competition influences the distribution of species on islands, or the patterns found are just by chance.

Assuming different species have independent distributions on islands, then the observed binary matrix is simply a random sample from the uniform distribution of all the binary matrices with fixed margins. Table 2 gives an example of all configurations of $3\times 3$ binary matrices with $[1,2,1]$ as both row and column sums. Ideally, if we could list all the binary matrices with arbitrary size, then we could compare the pattern found in the observed matrix with others, to conclude whether the observed matrix is simply by chance. However, enumerating matrices with fixed margins is often impractical both theoretically and computationally for moderate size of matrices. Therefore, sampling such random matrices becomes the natural choice.

The problem of sampling such matrices also occurs in many other fields, with different names. For example, an equivalent formulation is uniformly sampling undirected bipartite graphs with given vertex degrees. A bipartite graph $G=(U,V,E)$ is a graph whose vertices are divided into two disjointed sets, denoted by $U=\{u_{1},\cdots,u_{m}\}$, $V=\{v_{1},\cdots,v_{n}\}$. $E$ is called the edge set where every edge connects one vertex in $U$ to one in $V$. The binary matrix $M=(m_{i,j})_{m\times n}$ is often called the bi-adjacency matrix of $G$ and is defined by

[TABLE]

Bipartite graphs are often used in network studies to model the interaction between two objects, for example, customers and products. It is often required to sample graphs with preserved degree sequence in network analysis uniformly. Throughout this paper, we will use the term ‘binary matrix’ instead of ‘bipartite graph’ to avoid confusion, although they are equivalent.

The algorithms of sampling binary matrices with fixed margins are divided into two classes. The first class of algorithms relies on the rejection sampling or importance sampling techniques, see [Sni91], [MP04], [CDHL05], [HM13] [HJ96] for examples. Importance sampling usually provides degree from non-uniform distribution, but it can be used to construct estimators to estimate the quantities of interest, such as the number of binary matrices with given margins. Chen et al. [CDHL05] introduced a sequential importance sampling (SIS) scheme to test the hypothesis we mentioned at the beginning of this paper on “Darwin’s Finch” dataset.

The second class falls into the Markov Chain Monte Carlo (MCMC) category and will be our main focus in this paper. The well-known “swap algorithm” has been used for decades. To the author’s best knowledge, it is first introduced by Besag and Clifford [BC89] in 1989 to solve a statistical testing problem. The swap algorithm has been formally proposed and analysed by Rao et al. [RJB96] and [KTV99] in the 1990s. A similar question is to sample matrices with non-negative integer entries, fixed row and column sums. Diaconis and Gangolli have proposed a random walk Metropolis algorithm [DG95]. Many variations and extensions of this algorithm are described by Diaconis and Sturmfels [DS*+*98].

The swap method attempts to make a single swap in each iteration, but when the matrix is large, or is mostly filled (or unfilled), the efficiency of swap algorithm can be relatively low. In 2008, Verhelst [Ver08] proposed a new MCMC algorithm based on the idea of performing multiple swaps per iteration. In 2014, Strona et al. [SNB*+*14] introduced the “Curveball algorithm”, which uses a ‘fair trade’ operation to replace the ‘swap’ operation in the swap algorithm, aiming for a faster mixing. The mathematical formulation of Curveball algorithm is equivalent to Verhelst’s algorithm, but with different implementation and reasoning. A nice survey and numerical comparisons of the existing algorithms can be found in a recent dissertation [Rec18]. The class of ‘multiple swaps’ algorithms tends to improve the mixing time empiricially. However, each step of the ‘multiple swaps’ algorithm is slower than the classical swap algorithm. Meanwhile, it is hard to compare the ‘multiple swaps’ algorithms and classical swap algorithm theoretically, as the corresponding Markov-chains have complicated behaviors and are therefore hard to analyze mathematically. The only existing result can be found in [CK18].

In this paper, we introduce a novel algorithm called Rectangle Loop algorithm. The algorithm is based on the classical swap algorithm, with a careful utilization of the matrix structure given by margins. We have also proved the resulting Markov Chain dominates the classical chain used in the swap algorithm in the sense of Peskun’s partial ordering [Pes73], and is easy to implement. Section 2 gives a review of swap algorithm and Curveball algorithm, including the details of both algorithms and a discussion. In Section 3 we introduce our new algorithm – Rectangle Loop algorithm. Section 4 proves the theoretical properties of Rectangle Loop algorithm. Section 5 gives numerical results.

2 Existing Methods

2.1 Swap Algorithm

The swap algorithm, or equivalently, swap chain is based on the idea of swapping checkerboard units. Here a checkerboard unit is a two by two matrix with one of the following forms:

[TABLE]

A swap means changing one checkerboard unit to the other.

Starting from an initial matrix, one chooses two rows and two columns uniformly at random among all rows and columns. If the resulting $2\times 2$ submatrix with entries in the intersection of these rows and columns is a checkerboard unit, it is swapped, otherwise, do nothing.

The swap algorithm is a Metropolis-type Markov chain Monte Carlo which converges to the uniform distribution.

2.2 Curveball Algorithm

The swap algorithm can often be inefficient, taking Darwin’s Finch data for example, there are ${13\choose 2}{17\choose 2}=10608$ submatrices with size $2\times 2$, however, only about $3\%$ of them are swappable. This means it requires a very large $T$ (the number of iterations) to ensure the generated degree is close to uniformly distributed. The Curveball algorithm provides another solution.

The Curveball algorithm uses ‘trade’ instead of ‘swap’ operation in each iteration. Steps 3-5 in Algorithm 2 gives an illustration of trading, it trades elements in column $V$ with elements in column $S_{a-b}$ for row $a$ and row $b$, preserving their row and column sums. Though seemingly complicated, there is a very intuitive explanation of the Curveball algorithm. We refer the readers to [SNB*+*14] for detailed illustrations.

3 Rectangle Loop Algorithm

The swap algorithm is proven to converge to uniform distribution. However, getting stuck at the same configuration is inefficient and thus the convergence could be very slow. For example, numerical experiments suggest that one would expect more than $30$ iterations before each successful swap using ‘Darwin’s Finch’ dataset. Assuming the randomly chosen row is mostly filled, such as row $1$ in Table 1, the two random chosen entries in this row would most likely be $[1,1]$ but the row of a ‘checkerboard unit’ has to be either $[1,0]$ or $[0,1]$. Therefore swapping rarely happens when the chosen row/column is mostly filled (or equivalently, mostly unfilled).

The Rectangle Loop algorithm is designed to increase the chance of swapping. The idea is illustrated in Figure $1$. In this example the target matrix is of size $5\times 5$, with row names $R_{1},\cdots,R_{5}$ and column names $C_{1},\cdots,C_{5}$. In each step, we choose one row and one column uniformly at random (Step A). Suppose $R_{2}$ and $C_{2}$ is chosen, with corresponding entry $1$ , the red number in the top middle plot of Figure 1. Then we randomly choose a [math] among all the [math]s in $R_{2}$ (Step B). Since there is only one [math] in $R_{2}$, which is at location $C_{4}$, this is our only choice. Again, we scan through all the entries in the same column with the [math] just chosen ($C_{4}$) and randomly choose a $1$ among all $1$s (Step C). In our example, the $1$s of $C_{4}$ are located at $R_{1}$ and $R_{4}$. Suppose we have chosen $(R_{4},C_{4})$. Now the three locations $(R_{2},C_{2}),(R_{2},C_{4}),(R_{4},C_{4})$ altogether give us the fourth one $(R_{4},C_{2})$, making the four entries a rectangle (Step D). If the fourth entry equals [math], then we swap the submatrix as we did in swap method (Step E). Otherewise the fourth entries equals $1$ and the original matrix is not changed. After Step A - E is iterated for many times, the resulting randomized matrices are used as representatives of uniformly distributed matrix with fixed margins.

The main difference between the Rectangle Loop algorithm and the swap algorithm is the sampling scheme. The Rectangle Loop algorithm is performing ‘conditional sampling’, making it more efficient than swap method, which is doing ‘unconditional choosing’. For example, suppose both the swap method and Rectangle Loop algorithm have chosen $R_{2}$ and $C_{2}$, an entry with value $1$. Then $R_{2}$ has only one [math] which is in column $4$. For the swap algorithm, the probability of correctly choosing $C_{4}$ is only $\frac{1}{4}$, as it is uniformly choosing among all columns. The Rectangle Loop algorithm, however, as the mechanism guarantees we only sample from the zero entries, chooses $C_{4}$ with probability $1$. Therefore it significantly increase the swapping probability, leading to a faster convergence than the swap chain.

The details of the Rectangle Loop algorithm is described in Algorithm 3. Noteworthy, when finding a [math] entry, we sample $1$ with the same column as the [math]. When finding a $1$, we sample [math] with the same row as the $1$. This ‘symmetric’ design ensures the algorithm that converges to the correct distribution, as will be proved in Section 4. The paths of sampled entries in each iteration always form a rectangle, that is where the name ‘Rectangle Loop algorithm’ comes from.

4 Theoretical Results

Given row sums $\textbf{r}=(r_{1},r_{2},\cdots,r_{m})$ and column sums $\textbf{c}=(c_{1},c_{2},\cdots,c_{n})$, we define $\Sigma_{\textbf{r},\textbf{c}}$ be the set of all matrices with row sums r and column sums c. The suffcient and necessary condition for $\Sigma_{\textbf{r},\textbf{c}}$ not being zero is given by Gale [G*+*57], Ryser [Rys09] in 1957. We call two matrices $A$, $B$ ‘swappable’ if one can transform to the other via one step swap algorithm. Equivalently, $A$ and $B$ only differs in a $2\times 2$ ‘checkerboard unit’. For the sake of simplicity, we assume henceforth $0<r_{i}<n,0<c_{j}<m$ for any $1\leq i\leq m,1\leq j\leq n$, as otherwise we could simply delete that degenerate row/column. The following theorem characterizes the limit distribution and transition probability of the swap chain.

Theorem 1**.**

Given $\textbf{r},\textbf{c}$ and an initial matrix $A_{0}\in\Sigma_{\textbf{r},\textbf{c}}$, the swap algorithm defines a Metropolis-type Markov chain with stationary distribution $\text{Unif}(\Sigma_{\textbf{r},\textbf{c}})$, transition kernel:

[TABLE]

and acceptance probability $1$.

Proof.

When $A$ and $B$ are swappable, there exists two rows $i_{1},i_{2}$ and two columns $j_{1},j_{2}$ such that $A$ and $B$ only differs in the $2\times 2$ submatrix extracted by row $i_{1},i_{2}$ and column $j_{1},j_{2}$. Therefore the probability of swapping $A$ to $B$ equals

[TABLE]

Recall that in the setting of Metropolis-Hastings, the acceptance probability from $A$ to $B$ is given by

[TABLE]

notice here the stationary distribution $\pi$ is designed to be $\text{Unif}(\Sigma_{\textbf{r},\textbf{c}})$, $\mathbb{P}(A,B)=\mathbb{P}(B,A)=\frac{1}{{m\choose 2}{n\choose 2}}$. Hence it is clear that the swap algorithm is a Metropolis-type Markov chain with $\text{Unif}(\Sigma_{\textbf{r},\textbf{c}})$ as stationary distribtuion and acceptance probability $1$, which justifies the correctness of the swap algorithm. ∎

The key point in swap algorithm is symmetry. When two different states $A,B$ are swappable, the associated transition probability is symmetric, i.e.,

[TABLE]

this ensures the chain has acceptance probability $1$.

In order to compare the efficiency of different Markov kernels with the same distribution, Peskun [Pes73] first introduced the following partial-ordering.

Let $\mathbb{P}_{1}$, $\mathbb{P}_{2}$ be two Markov transition kernels on the same state space $\mathcal{S}$ with same stationary distribution $\pi$, then $\mathbb{P}_{1}$ dominates $\mathbb{P}_{2}$ off the diagonal, $\mathbb{P}_{1}\succeq\mathbb{P}_{2}$, if

[TABLE]

for all $x\in\mathcal{S}$ and $A$ measurable with $x\notin A$.

When the state space is finite, as in our case, $\mathbb{P}_{1}\succeq\mathbb{P}_{2}$ iff all the off-diagonal entries of $\mathbb{P}_{1}$ are greater than or equal to the corresponding off-diagonal entries of $\mathbb{P}_{2}$. This indicates $\mathbb{P}_{1}$ has lower probability to get stuck in the same state, and is exploring the state space in a more efficient way. The following theorem shows the rectangale loop algorithm also has uniform distribution as stationary distribution, and the corresponding chain dominates the swap chain off the diagonal. For simplicity, we will use $\mathbb{P}_{s}$ and $\mathbb{P}_{r}$ to denote transition kernel for the swap chain and rectangle loop chain, respectively, omitting its dependency on $\textbf{r},\textbf{c},\text{Unif}(\Sigma_{\textbf{r},\textbf{c}}).$

Theorem 2**.**

Given $\textbf{r},\textbf{c}$ and an initial matrix $A_{0}\in\Sigma_{\textbf{r},\textbf{c}}$, the Rectangale loop algorithm defines a Metropolis-type Markov chain with stationary distribution $\text{Unif}(\Sigma_{\textbf{r},\textbf{c}})$. The transition kernel $\mathbb{P}_{r}$ dominates $\mathbb{P}_{s}$ off the diagonal.

Proof.

Given any two swappable configurations $A$ and $B$, we are aiming to show $\mathbb{P}_{r}(A,B)=\mathbb{P}_{r}(B,A)\geq\mathbb{P}_{s}(A,B)$.

As $A,B$ are swappable, there exists two rows $i_{1},i_{2}$ and two columns $j_{1},j_{2}$ such that $A$ and $B$ only differs in the $2\times 2$ submatrix extracted by row $i_{1},i_{2}$ and column $j_{1},j_{2}$. Without loss of generality, we assume the checkerboard unit corresponding to $A$ has the form $\begin{bmatrix}1&0\\ 0&1\end{bmatrix}$, as shown in Figure 4. Notice that there are four vertices of the $2\times 2$ submatrix and the Rectangle Loop algorithm chooses one arbitray row and column at its first step. This suggests the probability of transforming $A$ to $B$ equals the summation of four probabilities, each one corresponds to choosing one specific vertex of the ‘checkerboard unit’ .

Figure 2 illustrates the calculation of one path, starting from the vertex $(i_{1},j_{1})$. The possibility of choosing row $i_{1}$ and column $j_{1}$ is $\frac{1}{mn}$. Then one chooses a [math] among all the [math]s in row $i_{1}$, and there are $n-r_{i_{1}}$ of them. Therefore the possibility of choosing column $j_{2}$ equals $\frac{1}{n-r_{i_{1}}}$. Similarly, after choosing $j_{2}$, one chooses a $1$ among all $1$s in column $j_{2}$, and there are $c_{j_{2}}$ of them. Hence the possibility of choosing row $i_{2}$ equals $\frac{1}{c_{j_{2}}}$. The fourth entry is fixed after determining the first three entries thus the last step has probability $1$. Multipling the possibilities above altogether, the possibility of transforming $A$ to $B$, starting with $(i_{1},j_{1})$, equals

[TABLE]

The probability of transforming $A$ to $B$ with other starting vertex can be calculated accordingly. It turns out $\mathbb{P}_{r}(A,B)$ can be written as the following summation:

[TABLE]

Following the same strategy, $\mathbb{P}_{r}(B,A)$ can also be calculated below. Figure 3 illustrates the calculation of $\mathbb{P}_{r}(B,A)$ starting from $(i_{1},j_{1})$.

[TABLE]

After matching all the terms of $\mathbb{P}_{r}(B,A)$ with $\mathbb{P}_{r}(A,B)$, we conclude that $\mathbb{P}_{r}(B,A)=\mathbb{P}_{r}(A,B)$, which justifies the Rectangle Loop algorithm has Unif($\Sigma_{\textbf{r},\textbf{c}}$) as stationary distribution.

To show $\mathbb{P}_{r}(A,B)\geq\mathbb{P}_{s}(A,B)$, notice that $\mathbb{P}_{s}(A,B)=\frac{1}{{m\choose 2}{n\choose 2}}=\frac{4}{m(m-1)n(n-1)}$. and

[TABLE]

It is clear that $\mathbb{P}_{r}(A,B)$ can be written as the summation of four terms. Each term in the summation is greater than or equal to $\frac{1}{m(m-1)n(n-1)}$. Therefore we conclude that $\mathbb{P}_{r}(A,B)\geq\mathbb{P}_{s}(A,B)$ for any swappable $A,B$. This indicates $\mathbb{P}_{r}\succeq\mathbb{P}_{s}$, the Rectangle Loop algorithm is exploring the state space in a more efficient way than the swap algorithm. ∎

5 Simulation Results and Applications

5.1 A Concrete Example

Now we use the example used in [SNB*+*14] and [MP04] to compare the existing algorithms and the Rectangle Loop algorithm. The example below is concrete. The transition matrix can be calculated explicitly and convergence can be assessed analytically.

The five matrices shown in Table 2 are all possible configurations of $3\times 3$ binary matrices with $[1,2,1]$ as both row and column sums. The transition matrices for swap algorithm, Curveball algorithm and Rectangale loop algorithm is shown in Table 3.

Figure 4 shows the comparison of the three algorithms. Here we measure the distance between transition kernel $\mathbb{P}$ and the stationary distribution $\pi$ by total variation distance:

[TABLE]

where $k$ denotes the power. It is clear that all the algorithms converge exponentially to uniform distribution, but the Rectangle Loop algorithm converges faster than swap algorithm and Curveball algorithm.

5.2 Experiments on Empirical Mixing Time

For larger matrices, it is infeasible to calculate the transition matrix theoretically. To provide empirical justification for the advantage of the Rectangle Loop algorithm. We have designed the experiment as follows. For each $p=0.01,0.05,0.1,0.2,0.3,0.4,0.5$, a $100\times 100$ binary matrix is generated for which each entry has probability $p$ to be $1$. We ran each algorithm for $10000$ iterations and collected the corresponding number of swaps, as shown in Table 4. When the filled portion $p$ is small, the Rectangle Loop algorithm is extremely efficient, producing more than $73$ times more swaps than the swap algorithm. For large $p$, the advantage of Rectangle Loop algorithm is reduced, but still very significant. For $p=0.5$, the Rectangle Loop still produces $4$ times more swaps than the swap algorithm. Noteworthy, the zeros and ones play the symmetric rule in a binary matrix, therefore it is not necessary to generate the random matrix for $p>0.5$.

The result above justifies our theoretical result that the Rectangle Loop algorithm converges faster than the swap algorithm. However, the above experiments did not consider the running time for each iteration. In fact, one iteration of Rectangle Loop algorithm is computationally more expensive than that of swap algorithm. To investigate this issue, we also record the time per swap for both algorithms, as shown in the last column of Table 4. It turns out that the Rectangle Loop algorithm still has a significant advantage than the swap algorithm after the running time issue is taken into account. For $p=0.01$, the Rectangle Loop is about $31$ times more efficient than the swap algorithm. Even for $p=0.5$, Rectangle Loop algorithm is still about $4$ times more efficient than the swap algorithm.

We have also used the pertubation score suggested by Strona et.al. [SNB*+*14] to access convergence for both algorithm. Pertubation score of a matrix is defined by the fraction of cells differing from the corresponding ones of the initial matrix. It takes several iterations for each algorithm to stabilize around its expectation. It is shown in Figure 5 that it takes less iterations and less time for the Rectangle Loop algorithm to stabilize, suggesting a faster mixing than the swap algorithm.

5.3 Finch Data Applications

Going back to ‘Darwin’s Finch’ dataset, we use the test statistics $\bar{S}^{2}$ suggested by Roberts and Stone [RS90] to compare the three algorithms. $\bar{S}^{2}$ is defined by

[TABLE]

where $m$ is the number of rows of matrix $A$, $S=(s_{ij})=AA^{T}$. For the finch data, $\bar{S}^{2}=45.03$. Suppose this number is too large or too small, comparing with its expectation over all the matrices having the same margins as finch data. We would like to conclude that the cooperation/competition do influence the distribution of species. To investigate this, we implemented the swap algorithm, Curveball algorithm and Rectangle Loop algorithm on the same data for $20000$ iterations, using its average as an estimator for $\mathbb{E}({\bar{S}}^{2})$. The results are shown in Figure 6. After $20000$ iterations, Rectangle Loop algorithm gives an estimate of $42.135$ with standard deviation $0.537$, swap algorithm gives an estimate of $42.126$ with standard deviation $0.509$, Curveball algorithm gives an estimate of $42.191$, with standard deviation $0.590$. Therefore the observed data falls outside the three standard deviation boundaries for all three algorithms, suggesting strong evidence that the observed occurence matrix is not just by chance. Meanwhile, both swap algorithm and Rectangle algorithm gives similar estimations and lower standard deviations, which seem to be more accurate than the Curveball algorithm. Lastly, there is a significant pattern in Figure 6 that for both $\bar{S}^{2}$ and standard deviation estimation, the Rectangle Loop algorithm becomes stabilized much earilier than the swap algorithm, indicating a faster mixing.

6 Discussion

There is a growing tendency to study the behavior of binary matrices with fixed margins in numerous scientific fields, ranging from mathematics to natural science to social science. For example, mathematicians and computer scientists are interested in the total number of configurations of given margin sums. Ecologists use the so-called occurence matrix to model the presence/absence of species in different locations. Biologists use the binary matrix to model neuronal networks. Social scientists use the binary matrix for studying social network features.

One of the central and difficult problems is uniformly sampling binary matrices with given margins. In this article we have developed the Rectangle Loop algorithm which is efficient, intuitive and easy to implement. Theoretically, the algorithm is superior to the classical swap algorithm in Peskun’s order. In practice, the Rectangle Loop algorithm is notably more efficient than the swap approach. For a fixed number of iterations, Rectangle Loop algorithm produces $4-73$ times more successful swaps than the swap algorithm. For a fixed amount of time, Rectangle Loop algorithm still produces $4-31$ times more successful swaps than the swap algorithm. This suggests the Rectangle Loop algorithm is efficient both statistically and computationally.

There are many other problems that remain. From a theoretical point of view, it is important to give sharp bounds on the convergence speed of a given Markov chain. However, giving a useful running time estimate is often challenging in practical problems. It would be very interesting if the swap algorithm, Curveball algorithm and the Rectangle Loop algorithm can be investigated analytically. From an applied point of view, there are many factors that influence the performance of algorithms, such as running time per step (swap algorithm is the fastest, while Curveball algorithm is the slowest), initialization of the matrix, size of the matrix, ratio between row number and column numbers, filled proportions. Our empirical studies suggest that all the factors have a significant impact on the convergence speed for all the algorithms. It would be beneficial if more numerical experiments are carried out, yielding a complete and comprehensive comparison between all the existing algorithms.