The Snake Algorithm: A Rejection-Free Sampler for Binary Matrices with Fixed Margins

Abstract: We study uniform sampling of binary matrices with fixed row and column sums, a recurring problem in ecological null models, Rasch-model testing, network analysis, and combinatorics. We propose the Snake algorithm, a rejection-free Markov chain Monte Carlo sampler that grows an alternating path until its first self-intersection and flips the resulting loop. The chain is reversible and irreducible on the fixed-margin state space, hence has the uniform stationary distribution. We prove that one step flips on the order of $\sqrt{n}$ entries in sparse and balanced square regimes, give upper bounds on the per-step path length, and show that the resulting work per flipped entry is rate optimal in sparse and balanced regimes and near-optimal up to a polylogarithmic factor under a one-sided half-balanced condition. A Markov-chain comparison, combined with the recently established universal spectral-gap bound for the swap chain, proves that the lazy Snake chain is rapidly mixing for every feasible pair of margins; in the permutation-matrix case, the raw chain has the sharp total-variation mixing time $Θ(n \log n)$. We also describe a directed-graph extension and an equal-margin label-shuffling variant. Numerical experiments against Swap, Rectangle Loop, Curveball, sequential importance sampling, and a directed edge-swap algorithm show consistent gains in move size, wall-clock convergence, and sampling efficiency.

• The paper introduces a rejection-free, reversible sampler that flips uniquely closed alternating loops while preserving row and column sums, enabling every iteration to succeed.
• The algorithm flips Θ(√n) entries in sparse and balanced regimes at constant work per flipped entry, compared with severe rejection rates for local swap methods.
• The paper proves Θ(√n log n) mixing for permutation matrices, establishes rapid mixing for all feasible margins, and reports substantial empirical gains for matrix and directed-graph sampling.

Problem and motivation

Uniform sampling of binary matrices with prescribed row and column sums underlies null-model testing in ecology, goodness-of-fit testing for the Rasch model in psychometrics, network reconstruction from degree sequences, and several problems in theoretical computer science. The state space $\Sigma(\mathbf{r}, \mathbf{c})$ is typically exponential in the matrix dimensions, so enumeration is infeasible and MCMC methods targeting the uniform distribution are the standard approach. The dominant method, the swap algorithm of Besag–Clifford and Diaconis–Sturmfels, flips a $2\times 2$ checkerboard per accepted step. Its two structural weaknesses are a bounded move size (at most four entries) and a high rejection rate: on sparse matrices, the expected number of attempts per successful swap scales with $mn$. On permutation matrices, swap reduces to a lazy random transposition walk with mixing time $\Theta(n^3\log n)$ — over $10^9$ steps to mix a $1000\times 1000$ matrix.

The algorithm

The Snake algorithm grows an alternating path through the current matrix. Starting from a uniformly chosen cell, it alternates between row moves (from a one-entry to a uniformly chosen zero in the same row) and column moves (from a zero-entry to a uniformly chosen one in the same column), until the first self-intersection closes an alternating loop; all loop entries are then flipped. Because row and column coordinates are finite, termination occurs after at most $m+n+1$ visited cells, so every step succeeds — the chain is rejection-free — and the closing loop is unique. Flipping an alternating loop preserves all margins by construction.

The paper also develops two extensions. Directed-Snake samples simple directed graphs with fixed in- and out-degree sequences by representing the graph as a bipartite adjacency matrix with label-diagonal cells structurally fixed at zero; trajectories never touch these forbidden cells. Correctness holds under the standard connectivity condition via directed 2-switches and directed 3-cycle reorientations (Rao–Jana–Bandyopadhyay / Berger–Müller-Hannemann). Snake+ interleaves five ordinary transitions with uniform permutations within equal-margin row and column classes; this block kernel is generally non-reversible but preserves the uniform law.

Correctness

The main correctness theorem establishes that the transition kernel is irreducible and reversible with respect to $U(\Sigma(\mathbf{r},\mathbf{c}))$, hence has the uniform stationary distribution. Irreducibility follows from Ryser's theorem that checkerboard swaps connect $\Sigma(\mathbf{r},\mathbf{c})$, since each swap is realizable as a length-four loop with positive probability. Reversibility is proved by constructing an involution on first-closing histories that maps paths from $A$ to $2\times 2$0 onto equal-probability reversed paths from $2\times 2$1 to $2\times 2$2; the key observation is that flipping a balanced loop leaves the number of zeros in each touched row and ones in each touched column unchanged, so the denominators of the transition probabilities match.

Move size and computational cost

The efficiency analysis centers on two quantities: $2\times 2$3, the number of flipped entries, and $2\times 2$4, the inner-loop length, with expectations conditional on any current state and uniform over states. The results are:

• No assumptions: $2\times 2$5... more precisely, the universal upper bound is $2\times 2$6.
• Sparse ($2\times 2$7 or $2\times 2$8): $2\times 2$9.
• Balanced ($mn$0, $mn$1, $mn$2): $mn$3.
• Rate optimality: in both sparse ($mn$4) and balanced regimes, $mn$5, i.e., constant work per flipped entry.
• Half-balanced: if either $mn$6 or $mn$7 with $mn$8, then $mn$9 — near-optimal up to a polylogarithmic factor, without requiring $\Theta(n^3\log n)$0.

The contrast with local methods is stark: in the sparsest square case, swap achieves only $\Theta(n^3\log n)$1 flipped entries per unit work versus $\Theta(n^3\log n)$2 for Snake. With hash-table bookkeeping of seen coordinates, each inner iteration costs $\Theta(n^3\log n)$3, so total per-step cost is $\Theta(n^3\log n)$4.

Mixing time

Two mixing results anchor the theory. First, for permutation matrices ($\Theta(n^3\log n)$5), one Snake step is right multiplication by a random cycle whose length has birthday scale $\Theta(n^3\log n)$6. The paper proves the sharp bound

$\Theta(n^3\log n)$7

The lower bound uses a coupon-collector obstruction (untouched labels remain fixed points); the upper bound combines Hough's character-ratio estimates for random $\Theta(n^3\log n)$8-cycles with a parity-smoothing lemma showing the sign representation contributes only $\Theta(n^3\log n)$9 per step. This is exponentially faster than the transposition walk's $10^9$0.

Second, for arbitrary feasible margins, a Dirichlet-form comparison shows each lazy swap move can be realized by Snake with probability at least $10^9$1 where $10^9$2. Combined with the universal spectral-gap bound $10^9$3 of Fu, Qin, and Wang, this yields

$10^9$4

so the lazy Snake chain is rapidly mixing for every feasible pair of margins. The authors are explicit that this universal estimate is conservative — it exploits only single-checkerboard trajectories — and is not expected to be sharp outside the permutation case.

Numerical experiments

Experiments use compiled Rcpp implementations. On Bernoulli-initialized $10^9$5 matrices across fill levels from 1% to 50%, Snake flips roughly 18 entries per iteration throughout, while Rectangle Loop and Swap degrade sharply as fill decreases (Swap drops to $10^9$6 flips per attempt at 1% fill); Snake also achieves the highest flipped-entries-per-second rate in every setting. In the ultrasparse permutation regime, mean flips per attempt grow from $10^9$7 at $10^9$8 to $10^9$9 at $1000\times 1000$0, consistent with the $1000\times 1000$1 law, while Swap averages fewer than $1000\times 1000$2 flips per attempt for $1000\times 1000$3.

On convergence diagnostics for $1000\times 1000$4 banded matrices with margin sums 10, Snake+ reaches the perturbation threshold $1000\times 1000$5 in $1000\times 1000$6 seconds versus $1000\times 1000$7 s for Curveball and $1000\times 1000$8 s for Rectangle Loop; Swap fails to reach any displayed threshold within the window. Against conditional-Poisson sequential importance sampling on exactly calibrated Rasch examples, both methods match exact hypergeometric tail probabilities within Monte Carlo error, but Snake is $1000\times 1000$9 ($m+n+1$0) and $m+n+1$1 ($m+n+1$2) more efficient by the variance–time product criterion. For directed graphs, D-Snake achieves roughly $m+n+1$3 higher effective sample throughput than directed edge-swap on a reciprocity null test (4965.7 vs. 311.7 effective samples per second), with both samplers agreeing on the null mean and placing the observed 480-dyad reciprocal ring far outside the fixed-degree null.

Limitations and open questions

Several caveats deserve emphasis. The half-balanced efficiency guarantee carries a $m+n+1$4 factor, which the authors do not claim to be tight. The universal rapid-mixing bound is explicitly conservative and known to be suboptimal even in the permutation case; sharp Snake mixing times for broader margin families remain open. The D-Snake irreducibility result is conditional on the standard connectivity assumption for directed degree sequences, and no general mixing-time analysis with diagonal constraints exists. The move-size results are diagnostics rather than wall-clock guarantees — large steps do not automatically imply short mixing for every statistic — though the mixing theorems partially address this concern. The SIS comparison relies on a specific conditional-Poisson proposal; other importance-sampling designs could narrow the gap.

Conclusion

This paper introduces a rejection-free MCMC sampler for binary matrices with fixed margins that provably makes nonlocal moves of order $m+n+1$5 in sparse and balanced regimes at rate-optimal work per flipped entry. It establishes rapid mixing for all feasible margins via comparison with the swap chain and a sharp $m+n+1$6 mixing time on permutation matrices, supported by consistent empirical gains over Swap, Rectangle Loop, Curveball, and SIS. The principal open problem is closing the gap between the conservative universal mixing bound and margin-specific behavior.