Exact and sequential-importance-sampling fixed-margin samplers

Project.toml

@@ -2,7 +2,7 @@ name = "RandomBooleanMatrices"
 uuid = "9ae346a0-3d16-5633-ad70-ddb60ab77eac"
 license = "MIT"
 authors = ["mkborregaard <mkborregaard@snm.ku.dk"]
-version = "0.1.2"
+version = "0.2.0"
 
 [deps]
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
@@ -15,7 +15,7 @@ julia = "1"
 Random = "1"
 RandomNumbers = "1"
 SparseArrays = "1"
-StatsBase = "0.33, 0.34"
+StatsBase = "0.34"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

README.md

@@ -13,30 +13,110 @@ alternative approaches like the `quasiswap` algorithm in R's vegan package. Thes
 methods are neither efficient, nor are they guaranteed to sample the possible
 distribution of boolean vectors with a given row and column sum equally.
 
-Currently, the package only offers an implementation of the `curveball` algorithm
-of Strona et al. (2014). There are two forms: a `randomize_matrix!` function
-that will randomize a sparse boolean `Matrix` in-place, and a generator form:
+The package offers three algorithms, selected with the `method` keyword. There
+are two forms for each: a `randomize_matrix!` function that randomizes a sparse
+boolean `Matrix` in place, and a generator form:
 
 ```julia
 using SparseArrays, RandomBooleanMatrices
 
 # in-place
 m = sprand(Bool, 1000, 1000, 0.1)
-randomize_matrix!(m)
+randomize_matrix!(m)                    # curveball by default
+randomize_matrix!(m, method = sis)     # or sis / exact
 
 # using a Matrix generator object
 m = sprand(Bool, 1000, 1000, 0.1)
-rmg = matrixrandomizer(m)
+rmg = matrixrandomizer(m, method = sis)
 m1 = rand(rmg) # creates a new random matrix
 m2 = rand(rmg)
 
+# the curveball generator is warm-started from an independent draw, and `trades`
+# sets how many trades separate successive samples (raise it to thin correlation)
+rmg = matrixrandomizer(m, method = curveball, trades = 10_000)
+
 # You can also avoid copying by
 m3 = rand!(rmg)
 # but notice that this will not create a new copy of the Matrix, so generating multiple matrices at once with this is impossible
 ```
 
+## Methods
+
+  * `curveball` (default) — the curveball algorithm of Strona et al. (2014). A
+    fast in-place trade that performs a Markov step from the current matrix. It is
+    **sequential**: each call continues the chain from where the previous one left
+    off, so successive draws are *correlated*. The chain's stationary distribution
+    is the uniform one, so over a run it samples fixed-margin matrices uniformly.
+    The `matrixrandomizer` generator warm-starts the chain from an independent
+    `sis` draw, so even the first sample is decorrelated from the input matrix, and
+    the `trades` keyword sets how many trades separate successive draws — raise it
+    to thin out the correlation between samples.
+  * `sis` — sequential importance sampling (Harrison & Miller 2013). Each call
+    draws an **independent** matrix from a fast, near-uniform approximation of the
+    fixed-margin distribution. Scales to large matrices and is the natural choice
+    when independent draws are wanted.
+  * `exact` — exact uniform sampling (Miller & Harrison 2013). Each call draws an
+    **independent**, *exactly* uniform matrix, and is feasible for small matrices
+    (roughly `rows + cols ≤ 100`, or very sparse). The matrix count is computed
+    once by `matrixrandomizer` and reused by every `rand`, so the generator form
+    is much cheaper than `randomize_matrix!(m, method = exact)` for repeated
+    sampling.
+
+In short: `curveball` is a sequential Markov chain (correlated draws, uniform in
+the limit), while `exact` and `sis` produce independent draws (exactly uniform,
+resp. close to uniform) on every call.
+
+## Weighted (importance) sampling
+
+The `sis` method above is the uniform special case of the importance sampler of
+Harrison & Miller (2013), which targets the weighted distribution
+
+```
+P*(z) ∝ ∏ w[i,j] ^ z[i,j]
+```
+
+over binary matrices with the given margins. `importance_sampler(m, w)` builds a
+sampler for nonnegative weights `w` (the same size as `m`); each `rand` returns an
+*independent* draw together with the log of its importance weight `∏ w^z / Q*(z)`.
+Reweighting by `exp(logweight)` gives Monte-Carlo estimates: the mean weight
+estimates the normalising constant `κ = Σ_z ∏ w^z`, and a weight-weighted average
+of any statistic estimates its expectation under `P*`.
+
+```julia
+using SparseArrays, RandomBooleanMatrices
+
+m = sprand(Bool, 40, 30, 0.3)
+w = rand(40, 30) .+ 0.5                 # nonnegative weights
+sampler = importance_sampler(m, w)
+
+draw = rand(sampler)
+draw.matrix                              # an independent fixed-margin sample
+draw.logweight                           # log( ∏ w^z / Q*(z) )
+
+# estimate κ = Σ_z ∏ w^z and E[h(Z)] under P*
+draws  = [rand(sampler) for _ in 1:1000]
+wts    = exp.(getfield.(draws, :logweight))
+κ̂      = sum(wts) / length(wts)
+Eh     = sum(wts .* h.(getfield.(draws, :matrix))) / sum(wts)   # for some statistic h
+```
+
+A zero weight `w[i,j] == 0` is a **structural zero**: it forbids a one at that
+position (useful e.g. for a zero diagonal, the uniform distribution over directed
+graphs). A draw that the zeros leave no way to complete comes back with
+`logweight == -Inf` (importance weight zero) and is simply discarded. The weights
+are also **Sinkhorn-canonicalised** before sampling (pass `canonicalize = false`
+to opt out); this never changes the importance weights, only the sampler's
+variance. The one remaining piece of the paper not yet implemented is its
+variance-based column ordering, which only affects efficiency.
+
 ## References
 
 Strona, G., Nappo, D., Boccacci, F., Fattorini, S. & San-Miguel-Ayanz, J. (2014)
 A fast and unbiased procedure to randomize ecological binary matrices with fixed row and column totals.
 Nature Communications, 5, 4114.
+
+Miller, J.W. & Harrison, M.T. (2013) Exact sampling and counting for fixed-margin matrices.
+The Annals of Statistics, 41, 1569-1592.
+
+Harrison, M.T. & Miller, J.W. (2013) Importance sampling for weighted binary random matrices with specified margins.
+arXiv:1301.3928.

src/RandomBooleanMatrices.jl

@@ -5,39 +5,80 @@ using RandomNumbers.Xorshifts
 using SparseArrays
 using StatsBase
 
+include("common.jl")
 include("curveball.jl")
+include("exact.jl")
+include("weights.jl")
+include("sis.jl")
 
-@enum matrixrandomizations curveball
+@enum matrixrandomizations curveball exact sis
 
 """
-    randomize_matrix!(m; method = curveball)
-
-Randomize the sparse boolean Matrix `m` while maintaining row and column sums
+    randomize_matrix!(m [,rng]; method = curveball, trades = 5 * size(m, 2))
+
+Randomize the sparse boolean Matrix `m` in place while maintaining its row and
+column sums. The algorithm is chosen with `method`:
+
+  * `curveball` — the curveball trade of Strona et al. (2014). Applies `trades`
+    trades, each a Markov step from the current matrix; repeated calls explore the
+    space of fixed-margin matrices.
+  * `sis` — sequential importance sampling (Harrison & Miller 2013). Draws an
+    independent matrix from a fast, near-uniform approximation. Suited to large
+    matrices.
+  * `exact` — exact uniform sampling (Miller & Harrison 2013). Draws an
+    independent, exactly uniform matrix; feasible for small matrices. The counting
+    step is repeated on every call, so prefer [`matrixrandomizer`](@ref), which
+    counts once and reuses the result.
+
+`trades` applies only to `curveball` and sets how far the chain moves per call.
 """
-function randomize_matrix!(m::SparseMatrixCSC{Bool, Int}, rng = Random.GLOBAL_RNG; method::matrixrandomizations = curveball)
-    if method == curveball
-        return _curveball!(m, rng)
-    end
+function randomize_matrix!(m::SparseMatrixCSC{Bool, Int}, rng = Random.GLOBAL_RNG; method::matrixrandomizations = curveball, trades::Int = 5 * size(m, 2))
+    method == curveball && return _curveball!(m, rng, trades)
+    method == sis       && return _sis!(m, rng)
+    method == exact     && return _exact!(m, rng)
     error("undefined method")
 end
 
-
-SparseMatrixCSC{Bool, Int}
+# A generator stores a method-specific sampler that both selects the draw
+# algorithm (by dispatch, via `_draw!`) and carries whatever that method reuses
+# across draws.
+const GeneratorState = Union{CurveballSampler, SISSampler, ExactCounts}
 
 struct MatrixGenerator{R<:AbstractRNG, M}
     m::M
-    method::matrixrandomizations
     rng::R
+    state::GeneratorState
 end
 
-show(io::IO, m::MatrixGenerator{R, SparseMatrixCSC{Bool, Int}}) where R = println(io, "Boolean MatrixGenerator with size $(size(m.m)) and $(nnz(m.m)) occurrences")
+Base.show(io::IO, m::MatrixGenerator{R, SparseMatrixCSC{Bool, Int}}) where R = print(io, "Boolean MatrixGenerator with size $(size(m.m)) and $(nnz(m.m)) occurrences")
+
+# Build the sampler for `method`, performing any one-off preparation of `sm`:
+# `exact` caches its (expensive) matrix count, and `curveball` warm-starts the
+# stored matrix from an independent SIS draw so the chain begins in the bulk of
+# the distribution rather than at the input, then records its per-draw `trades`.
+function _sampler!(method::matrixrandomizations, sm::SparseMatrixCSC{Bool, Int}, rng, trades)
+    method == sis   && return SISSampler()
+    method == exact && return _exact_counts(_margins(sm)...)
+    if method == curveball
+        _sis!(sm, rng)
+        return CurveballSampler(something(trades, 5 * size(sm, 2)))
+    end
+    error("undefined method")
+end
 
 """
-    matrixrandomizer(m [,rng]; method = curveball)
+    matrixrandomizer(m [,rng]; method = curveball, trades = 5 * size(m, 2))
 
-Create a matrix generator function that will return a random boolean matrix
-every time it is called, maintaining row and column sums. Non-boolean input
-matrix are interpreted as boolean, where values != 0 are `true`.
+Create a matrix generator that returns a random boolean matrix every time it is
+called, maintaining row and column sums. Non-boolean input matrices are
+interpreted as boolean, where values != 0 are `true`. See [`randomize_matrix!`](@ref)
+for the available `method`s.
+
+For `method = exact` the matrix count is computed once here and reused by every
+draw. For `method = curveball` the chain is warm-started from an independent `sis`
+draw, so even the first sample is decorrelated from `m`; `trades` then sets how
+many curveball trades separate successive draws — raise it to reduce correlation
+between samples.
 
 # Examples
 ```
@@ -50,15 +91,85 @@ random2 = rand(rmg)
 """
 matrixrandomizer(m, rng) = error("No matrixrandomizer defined for $(typeof(m))")
 matrixrandomizer(m) = error("No matrixrandomizer defined for $(typeof(m))")
-matrixrandomizer(m::AbstractMatrix, rng = Xoroshiro128Plus(); method::matrixrandomizations = curveball) =
-    MatrixGenerator{typeof(rng), SparseMatrixCSC{Bool, Int}}(dropzeros!(sparse(m)), method, rng)
-matrixrandomizer(m::SparseMatrixCSC{Bool, Int}, rng = Xoroshiro128Plus(); method::matrixrandomizations = curveball) =
-    MatrixGenerator{typeof(rng), SparseMatrixCSC{Bool, Int}}(dropzeros(m), method, rng)
+function matrixrandomizer(m::AbstractMatrix, rng = Xoroshiro128Plus(); method::matrixrandomizations = curveball, trades = nothing)
+    sm = SparseMatrixCSC{Bool, Int}(dropzeros!(sparse(m)))
+    MatrixGenerator{typeof(rng), SparseMatrixCSC{Bool, Int}}(sm, rng, _sampler!(method, sm, rng, trades))
+end
+function matrixrandomizer(m::SparseMatrixCSC{Bool, Int}, rng = Xoroshiro128Plus(); method::matrixrandomizations = curveball, trades = nothing)
+    sm = dropzeros(m)
+    MatrixGenerator{typeof(rng), SparseMatrixCSC{Bool, Int}}(sm, rng, _sampler!(method, sm, rng, trades))
+end
+
+# A single draw from the generator, in place in its stored matrix; the algorithm
+# is selected by the type of the stored sampler state.
+_generate!(r::MatrixGenerator{R, SparseMatrixCSC{Bool, Int}}) where R = _draw!(r.m, r.rng, r.state)
+
+Random.rand(r::MatrixGenerator{R, SparseMatrixCSC{Bool, Int}}) where R = copy(_generate!(r))
+Random.rand!(r::MatrixGenerator{R, SparseMatrixCSC{Bool, Int}}) where R = _generate!(r)
+
+"""
+    WeightedSIS
+
+A sequential importance sampler for the weighted fixed-margin target
+P*(z) ∝ ∏ w[i,j]^z[i,j]. Built by [`importance_sampler`](@ref); each `rand` draws
+an independent matrix and returns it with the log of its importance weight.
+"""
+struct WeightedSIS{R<:AbstractRNG}
+    rowsums::Vector{Int}
+    colsums::Vector{Int}
+    model::WeightMatrix
+    rng::R
+end
+
+Base.show(io::IO, s::WeightedSIS) =
+    print(io, "WeightedSIS sampler for a $(length(s.rowsums))×$(length(s.colsums)) fixed-margin matrix")
+
+"""
+    importance_sampler(m, w [,rng])
+
+Create a sequential importance sampler (Harrison & Miller 2013) for the weighted
+distribution over binary matrices with the row and column sums of `m`,
+
+    P*(z) ∝ ∏ w[i,j]^z[i,j],
+
+where `w` is a matrix of nonnegative weights the size of `m`; the uniform special
+case is `w` all ones. Each `rand` returns a `(matrix, logweight)` pair: an
+independent draw and the log of its importance weight `∏ w^z / Q*(z)`. Monte-Carlo
+estimates reweight by `exp(logweight)` — for example `mean(exp(logweight))`
+estimates the normalising constant `κ = Σ_z ∏ w^z`, and
+`sum(exp(logweight) .* h) / sum(exp(logweight))` estimates `E[h(Z)]` under P*.
 
-Random.rand(r::MatrixGenerator{R, SparseMatrixCSC{Bool, Int}}; method::matrixrandomizations = curveball) where R = copy(randomize_matrix!(r.m, r.rng, method = r.method))
-Random.rand!(r::MatrixGenerator{R, SparseMatrixCSC{Bool, Int}}; method::matrixrandomizations = curveball) where R = randomize_matrix!(r.m, r.rng, method = r.method)
+A zero weight `w[i,j] == 0` is a structural zero: it forbids a one at that
+position. A draw the zeros leave no way to complete is returned with
+`logweight == -Inf` (importance weight zero) and should simply be discarded.
+The proposal is built from the Sinkhorn-canonical form of `w` unless
+`canonicalize = false`; this never changes the importance weights, only the
+sampler's efficiency.
+
+# Examples
+```
+m = sprand(Bool, 20, 15, 0.3)
+w = rand(20, 15) .+ 0.5
+sampler = importance_sampler(m, w)
+draw = rand(sampler)
+draw.matrix      # an independent fixed-margin sample
+draw.logweight   # its log importance weight
+```
+"""
+function importance_sampler(m::AbstractMatrix, w::AbstractMatrix, rng = Xoroshiro128Plus();
+                            canonicalize::Bool = true)
+    size(w) == size(m) || throw(DimensionMismatch("weights `w` must match the size of `m`"))
+    sm = SparseMatrixCSC{Bool, Int}(dropzeros!(sparse(m)))
+    rowsums, colsums = _margins(sm)
+    WeightedSIS(rowsums, colsums, WeightMatrix(w, rowsums, colsums; canonicalize), rng)
+end
+
+function Random.rand(s::WeightedSIS)
+    columns, logq, logp = _sis(s.rowsums, s.colsums, s.model, s.rng)
+    (matrix = _columnsmatrix(columns, length(s.rowsums)), logweight = logp - logq)
+end
 
-export randomize_matrix!, matrixrandomizer, matrixrandomizations
-export curveball
+export randomize_matrix!, matrixrandomizer, matrixrandomizations, importance_sampler
+export curveball, exact, sis
 
 end

src/common.jl

@@ -0,0 +1,68 @@
+# Helpers shared by the fixed-margin samplers (exact.jl and sis.jl).
+
+"""
+    _margins(m)
+
+Return the `(rowsums, colsums)` of a sparse boolean matrix, reading them directly
+off the `SparseMatrixCSC` structure: column sums are the gaps between column
+pointers, row sums are the occurrences of each row index.
+"""
+function _margins(m::SparseMatrixCSC{Bool, Int})
+   colsums = diff(m.colptr)
+   rowsums = zeros(Int, size(m, 1))
+   @inbounds for r in m.rowval
+      rowsums[r] += 1
+   end
+   rowsums, colsums
+end
+
+"""
+    _conjugate(sizes, len)
+
+The conjugate (Ferrers transpose) of a multiset of column sums, truncated to
+`len` levels: `conjugate[k]` is the number of columns whose sum is at least `k`.
+It encodes the column sums in the form the Gale–Ryser feasibility test needs.
+"""
+function _conjugate(sizes, len::Int)
+   conjugate = zeros(Int, len)
+   @inbounds for s in sizes, k in 1:min(s, len)
+      conjugate[k] += 1
+   end
+   conjugate
+end
+
+"""
+    _writecols!(m, columns)
+
+Overwrite the row indices of each column of `m` in place from `columns`, a vector
+of row-index lists. The column sums (and hence `m.colptr` and `m.nzval`) are
+unchanged, so only `m.rowval` is rewritten, with the rows of each column sorted.
+"""
+function _writecols!(m::SparseMatrixCSC{Bool, Int}, columns)
+   @inbounds for j in 1:size(m, 2)
+      rows = sort!(columns[j])
+      copyto!(view(m.rowval, m.colptr[j]:m.colptr[j+1]-1), rows)
+   end
+   m
+end
+
+"""
+    _columnsmatrix(columns, nrows)
+
+Assemble a fresh `nrows × length(columns)` sparse boolean matrix from `columns`, a
+per-column vector of row-index lists. Used by samplers that build a matrix from
+scratch rather than overwriting an existing one.
+"""
+function _columnsmatrix(columns, nrows::Int)
+   ncols = length(columns)
+   colptr = Vector{Int}(undef, ncols + 1)
+   colptr[1] = 1
+   @inbounds for j in 1:ncols
+      colptr[j+1] = colptr[j] + length(columns[j])
+   end
+   rowval = Vector{Int}(undef, colptr[end] - 1)
+   @inbounds for j in 1:ncols
+      copyto!(view(rowval, colptr[j]:colptr[j+1]-1), sort!(columns[j]))
+   end
+   SparseMatrixCSC(nrows, ncols, colptr, rowval, fill(true, length(rowval)))
+end