From 0fd28f2059882d41e2caa5ac879761f30dc1c3a5 Mon Sep 17 00:00:00 2001 From: mkborregaard Date: Sat, 20 Jun 2026 21:54:17 +0200 Subject: [PATCH 1/7] Add exact and SIS fixed-margin samplers; update deps Add two independent-draw alternatives to the curveball sampler, both from Miller & Harrison, selected via the `method` keyword on `randomize_matrix!` and `matrixrandomizer`: - `exact`: exact uniform sampling via memoized counting of fixed-margin binary matrices (Miller & Harrison 2013). The generator counts once and reuses the result across draws. - `sis`: sequential importance sampling, a fast near-uniform proposal (Harrison & Miller 2013), suited to large matrices. Implements the uniform case; weighted sampling is left as a future extension. The dynamic programs are expressed as small dispatch-based functions over bespoke state types (ExactProblem/ExactCounts/Cursor, SISResidual/ SISWorkspace/ColumnScorer). Verified against brute-force counts and a uniformity check; the SIS proposal probability is validated by recovering the matrix count via importance sampling. Also bump StatsBase compat to 0.34 and the package version to 0.2.0. Co-Authored-By: Claude Opus 4.8 --- Project.toml | 4 +- README.md | 37 ++++- src/RandomBooleanMatrices.jl | 72 +++++++--- src/common.jl | 47 +++++++ src/exact.jl | 239 ++++++++++++++++++++++++++++++++ src/sis.jl | 259 +++++++++++++++++++++++++++++++++++ test/runtests.jl | 76 ++++++++++ 7 files changed, 707 insertions(+), 27 deletions(-) create mode 100644 src/common.jl create mode 100644 src/exact.jl create mode 100644 src/sis.jl diff --git a/Project.toml b/Project.toml index 3481647..770428b 100644 --- a/Project.toml +++ b/Project.toml @@ -2,7 +2,7 @@ name = "RandomBooleanMatrices" uuid = "9ae346a0-3d16-5633-ad70-ddb60ab77eac" license = "MIT" authors = ["mkborregaard big(1)) + prob = ExactProblem(p, length(rowsums), _pascal(length(colsums)), table, L) + c = copy(conjugate) + total = _count(prob, c, Cursor(1, 1, 0, p[2], c[1])) + table[copy(conjugate)] = total + ExactCounts(prob, rowperm, colsums, conjugate, total) +end + +# --------------------------------------------------------------------------- +# Sampling +# --------------------------------------------------------------------------- + +""" + _sample_plan!(plan, ex, rng) + +Draw a uniform fixed-margin matrix as a `plan`, where `plan[row, block]` is the +number of ones the sorted `row` places in column-block `block`. Each branch of +the recursion is taken in proportion to the number of completions it leaves, so +the resulting matrix is exactly uniform. +""" +function _sample_plan!(plan::Matrix{Int}, ex::ExactCounts, rng) + prob = ex.problem + c = copy(ex.conjugate) + fill!(plan, 0) + _sample!(plan, prob, c, Cursor(1, 1, 0, prob.p[2], c[1]), ex.total, rng) + plan +end + +function _sample!(plan::Matrix{Int}, prob::ExactProblem, c::Vector{Int}, cur::Cursor, total::BigInt, rng) + if cur.placed == prob.p[cur.row] + cur.row == prob.nrows + 1 && return + return _sample!(plan, prob, c, _next_row(prob, c, cur), total, rng) + end + width, range = _block_range(prob, c, cur) + target = rand(rng, zero(total):(total - one(total))) + cumulative = zero(BigInt) + for s in range + c[cur.block] -= s + completions = _count(prob, c, _next_block(prob, c, cur, s)) + cumulative += prob.binomial[width+1, s+1] * completions + if cumulative > target + plan[cur.row, cur.block] = s + _sample!(plan, prob, c, _next_block(prob, c, cur, s), completions, rng) + c[cur.block] += s + return + end + c[cur.block] += s + end + error("exact sampler reached an inconsistent state") +end + +# Turn a sampled plan into the actual matrix: for each row and block, pick which +# of the equal-residual columns receive the row's ones, uniformly at random. +function _plan_to_columns!(columns, plan::Matrix{Int}, ex::ExactCounts, rng) + residual = copy(ex.colsums) + for col in columns + empty!(col) + end + @inbounds for row in 1:ex.problem.nrows + for level in 1:ex.problem.L-1 + block = findall(==(level), residual) + for col in sample(rng, block, plan[row, level], replace = false) + push!(columns[col], ex.rowperm[row]) + residual[col] -= 1 + end + end + end + columns +end + +""" + _exact_sample!(columns, ex, rng) + +Draw one uniform fixed-margin matrix into `columns` (row indices per column). +""" +function _exact_sample!(columns, ex::ExactCounts, rng) + plan = Matrix{Int}(undef, ex.problem.nrows, ex.problem.L) + _sample_plan!(plan, ex, rng) + _plan_to_columns!(columns, plan, ex, rng) +end + +function _exact!(m::SparseMatrixCSC{Bool, Int}, rng = Random.GLOBAL_RNG) + rowsums, colsums = _margins(m) + columns = [Int[] for _ in colsums] + _exact_sample!(columns, _exact_counts(rowsums, colsums), rng) + _writecols!(m, columns) +end diff --git a/src/sis.jl b/src/sis.jl new file mode 100644 index 0000000..47e2a89 --- /dev/null +++ b/src/sis.jl @@ -0,0 +1,259 @@ +# Sequential importance sampling of fixed-margin binary matrices, after +# +# Harrison & Miller (2013), "Importance sampling for weighted binary random +# matrices with specified margins". +# +# The matrix is built one column at a time. Each column is a binary vector with a +# prescribed number of ones, drawn from a fast approximation of the conditional +# distribution that (a) respects the column total, (b) keeps the remaining +# problem completable (Gale–Ryser), and (c) prefers configurations in proportion +# to a combinatorial estimate of how many completions they admit. A dynamic +# program turns those per-row preferences and feasibility bounds into exact +# forward sampling probabilities for the column. This implements the uniform case +# (every fixed-margin matrix targeted with equal weight), the natural counterpart +# to the curveball sampler. + +# --------------------------------------------------------------------------- +# Combinatorial approximation of the per-row odds of a one +# --------------------------------------------------------------------------- + +""" + ColumnScorer + +A combinatorial approximation that scores how likely each row is to carry a one +in the column currently being sampled. Concrete scorers are callable: given a +row's residual sum they return that row's (unconstrained) probability of a one. +New approximations (e.g. Greenhill–McKay for sparse margins) can be added as +further subtypes without touching the dynamic program. +""" +abstract type ColumnScorer end + +""" + Canfield(nrows, ncols, total, weightA) + +The Canfield–Greenhill–McKay (2008) approximation used by Harrison & Miller for +near-regular margins. `ncols` and `total` describe the columns and ones that +remain *after* the current column, and `weightA` collects the curvature term of +the asymptotic enumeration so that scoring a row is a single `exp`. +""" +struct Canfield <: ColumnScorer + nrows::Int + ncols::Int + total::Int + weightA::Float64 +end + +function Canfield(nrows::Int, ncols::Int, total::Int, sumsq::Int) + weightA = if total == 0 || nrows * ncols == total + 0.0 + else + η = nrows * ncols / (total * (nrows * ncols - total)) + ν = η * (sumsq - total^2 / ncols) + η * (1 - ν) / 2 + end + Canfield(nrows, ncols, total, weightA) +end + +# Probability that a row with residual sum `r` carries a one in this column. +# Rows that are empty (r == 0) or saturated (r == ncols + 1) fall out at 0 and 1. +function (s::Canfield)(r::Int) + odds = r * exp(s.weightA * (1 - 2 * (r - s.total / s.nrows))) + odds / (s.ncols + 1 - r + odds) +end + +# --------------------------------------------------------------------------- +# The residual problem: margins still to be filled, evolving column by column +# --------------------------------------------------------------------------- + +""" + SISResidual + +The part of the problem still to be sampled. Rows are tracked by `residual` sum +and kept in descending `order`; columns are summarised by their `conjugate` +(for feasibility) together with the running `total` of ones, sum of squares +`sumsq` (for the scorer) and number `ncols` still to place. +""" +mutable struct SISResidual + residual::Vector{Int} # residual row sums, indexed by original row + order::Vector{Int} # original row indices, sorted by descending residual + conjugate::Vector{Int} # conjugate of the columns not yet sampled + total::Int # ones still to place + sumsq::Int # Σ colsum² over columns not yet sampled + ncols::Int # columns not yet sampled + nrows::Int +end + +function SISResidual(rowsums::Vector{Int}, colsums::Vector{Int}) + nrows = length(rowsums) + SISResidual(copy(rowsums), sortperm(rowsums, rev = true), + _conjugate(colsums, nrows), + sum(colsums), sum(abs2, colsums), length(colsums), nrows) +end + +# Remove a column of sum `k` from the column-side statistics, leaving the +# residual describing the *other* remaining columns (what a completion must fit). +function _advance!(res::SISResidual, k::Int) + @inbounds for level in 1:k + res.conjugate[level] -= 1 + end + res.total -= k + res.sumsq -= k^2 + res.ncols -= 1 + res +end + +# Decrement the rows that received a one and restore the descending order. +function _place!(res::SISResidual, rows) + @inbounds for r in rows + res.residual[r] -= 1 + end + sortperm!(res.order, res.residual, rev = true) + res +end + +# --------------------------------------------------------------------------- +# Per-column dynamic program +# --------------------------------------------------------------------------- + +""" + SISWorkspace + +Reusable scratch space for sampling one column: the per-row one-probabilities +`p`, the feasibility band `lo:hi` on the partial column sum after each row, the +forward transition probabilities `S`, and two rolling buffers for the backward +dynamic-programming weights. +""" +struct SISWorkspace + p::Vector{Float64} + lo::Vector{Int} + hi::Vector{Int} + S::Matrix{Float64} + gcur::Vector{Float64} + gnext::Vector{Float64} +end + +SISWorkspace(nrows::Int, maxcol::Int) = + SISWorkspace(Vector{Float64}(undef, nrows), Vector{Int}(undef, nrows), + Vector{Int}(undef, nrows), Matrix{Float64}(undef, maxcol + 1, nrows), + Vector{Float64}(undef, maxcol + 2), Vector{Float64}(undef, maxcol + 2)) + +# Score every row for the current column. +function _score!(ws::SISWorkspace, res::SISResidual, scorer::ColumnScorer) + @inbounds for i in 1:res.nrows + ws.p[i] = scorer(res.residual[res.order[i]]) + end + ws +end + +# Feasibility band for the partial column sum s_i = x_1 + … + x_i (rows in +# descending order). Gale–Ryser bounds it below; the column total and row count +# bound it above. `res.conjugate` is the conjugate of the *other* columns here. +function _band!(ws::SISWorkspace, res::SISResidual, k::Int) + m = res.nrows + rowsum = 0 # prefix of residual row sums + conjsum = 0 # prefix of the conjugate of the other columns + lo = 0 + @inbounds for i in 1:m + rowsum += res.residual[res.order[i]] + conjsum += res.conjugate[i] + lo = max(lo, rowsum - conjsum, k - (m - i)) + ws.lo[i] = lo + ws.hi[i] = min(k, i) + end + ws +end + +# Backward dynamic program (Harrison & Geman 2009): turn the per-row odds and the +# feasibility band into forward transition probabilities. `g[s]` is the total +# weight of feasible completions reaching column total `k` from partial sum `s`; +# `S[s+1, i]` is the resulting probability that row `i` is a one given partial sum +# `s` before it. +function _transitions!(ws::SISWorkspace, k::Int) + m = length(ws.p) + gnext, gcur, S = ws.gnext, ws.gcur, ws.S + fill!(gnext, 0.0) + gnext[k+1] = 1.0 # base case: all rows placed, sum == k + @inbounds for i in m:-1:1 + lo = i == 1 ? 0 : ws.lo[i-1] # band on the partial sum before row i + hi = i == 1 ? 0 : ws.hi[i-1] + p = ws.p[i] + fill!(gcur, 0.0) + total = 0.0 + for s in lo:hi + w1 = p * gnext[s+2] # row i is a one -> partial sum s+1 + w0 = (1 - p) * gnext[s+1] # row i is a zero -> partial sum s + weight = w0 + w1 + gcur[s+1] = weight + S[s+1, i] = weight > 0 ? w1 / weight : 0.0 + total += weight + end + total > 0 && (@views gcur[lo+1:hi+1] ./= total) # rescale to avoid underflow + gnext, gcur = gcur, gnext + end + ws +end + +# Forward pass: walk the rows, drawing each from its transition probability, and +# accumulate the log-probability of the column under the proposal. +function _draw_column!(rows, ws::SISWorkspace, order, k::Int, rng) + s = 0 + logq = 0.0 + @inbounds for i in eachindex(ws.p) + p = ws.S[s+1, i] + if rand(rng) < p + push!(rows, order[i]) + s += 1 + logq += log(p) + s == k && break + else + logq += log(1 - p) + end + end + logq +end + +""" + _sample_column!(rows, res, k, ws, rng) + +Sample one column of total `k` into `rows` (original row indices), update the +residual problem, and return the log proposal probability of the draw. +""" +function _sample_column!(rows, res::SISResidual, k::Int, ws::SISWorkspace, rng) + k == 0 && return 0.0 + _advance!(res, k) + _score!(ws, res, Canfield(res.nrows, res.ncols, res.total, res.sumsq)) + _band!(ws, res, k) + _transitions!(ws, k) + logq = _draw_column!(rows, ws, res.order, k, rng) + _place!(res, rows) + logq +end + +# --------------------------------------------------------------------------- +# Driver +# --------------------------------------------------------------------------- + +""" + _sis(rowsums, colsums, rng) + +Sample a fixed-margin binary matrix by sequential importance sampling, returning +the row indices of each column together with the log proposal probability of the +whole matrix. Columns are visited in order of decreasing sum, which the authors +find improves the approximation. +""" +function _sis(rowsums::Vector{Int}, colsums::Vector{Int}, rng) + res = SISResidual(rowsums, colsums) + ws = SISWorkspace(res.nrows, maximum(colsums, init = 0)) + columns = [Int[] for _ in colsums] + logq = 0.0 + for j in sortperm(colsums, rev = true) + logq += _sample_column!(columns[j], res, colsums[j], ws, rng) + end + columns, logq +end + +function _sis!(m::SparseMatrixCSC{Bool, Int}, rng = Random.GLOBAL_RNG) + rowsums, colsums = _margins(m) + columns, _ = _sis(rowsums, colsums, rng) + _writecols!(m, columns) +end diff --git a/test/runtests.jl b/test/runtests.jl index ea7c235..74942f6 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -5,6 +5,20 @@ using Test Random.seed!(1337) +# the row and column sums of a matrix, as plain vectors +margins(m) = (vec(sum(m, dims = 1)), vec(sum(m, dims = 2))) + +# brute-force count of the m×n binary matrices with the given margins +function countmatrices(rowsums, colsums) + m, n = length(rowsums), length(colsums) + count = 0 + for bits in 0:(2^(m * n) - 1) + A = [(bits >> ((j - 1) * m + (i - 1))) & 1 for i in 1:m, j in 1:n] + vec(sum(A, dims = 2)) == rowsums && vec(sum(A, dims = 1)) == colsums && (count += 1) + end + count +end + @testset "curveball" begin m = sprand(Bool, 8, 6, 0.2) m_old = copy(m) @@ -29,3 +43,65 @@ Random.seed!(1337) @test m3 != m4 end + +@testset "sis" begin + m = sprand(Bool, 8, 6, 0.3) + csm, rsm = margins(m) + randomize_matrix!(m, method = sis) + + @test margins(m) == (csm, rsm) + + m2 = rand(0:1, 6, 5) + csm, rsm = margins(m2) + rmg = matrixrandomizer(m2, method = sis) + m3 = rand(rmg) + m4 = rand(rmg) + + @test margins(m3) == (csm, rsm) # margins preserved + @test margins(m4) == (csm, rsm) + @test m3 != m4 # independent draws + + # a fully determined problem has a single solution, which must be returned + determined = sparse(Bool[1 1 1; 1 1 0; 1 0 0]) # unique matrix with margins (3,2,1) + @test countmatrices([3, 2, 1], [3, 2, 1]) == 1 + @test rand(matrixrandomizer(determined, method = sis)) == determined +end + +@testset "exact" begin + # counting agrees with brute-force enumeration + for (rs, cs) in [([2, 1, 1], [1, 2, 1]), ([2, 2, 1, 1], [2, 1, 2, 1]), + ([3, 1, 1, 1], [2, 2, 1, 1]), ([2, 2, 2], [2, 2, 2])] + @test RandomBooleanMatrices._exact_counts(rs, cs).total == countmatrices(rs, cs) + end + + m = sprand(Bool, 6, 5, 0.4) + csm, rsm = margins(m) + randomize_matrix!(m, method = exact) + @test margins(m) == (csm, rsm) + + m2 = rand(0:1, 5, 6) + csm, rsm = margins(m2) + rmg = matrixrandomizer(m2, method = exact) # counts once, reused below + m3 = rand(rmg) + m4 = rand(rmg) + + @test margins(m3) == (csm, rsm) + @test margins(m4) == (csm, rsm) + @test m3 != m4 + + # a fully determined problem returns its single solution + determined = sparse(Bool[1 1 1; 1 1 0; 1 0 0]) + @test rand(matrixrandomizer(determined, method = exact)) == determined + + # the exact sampler is uniform: every matrix with these margins appears, and + # frequencies are close to equal (seeded rng keeps this reproducible) + rs, cs = [2, 2, 1, 1], [2, 1, 2, 1] + seed = sparse(Bool[1 0 1 0; 1 1 0 0; 0 0 1 0; 0 0 0 1]) + @test margins(seed) == (cs, rs) + n = countmatrices(rs, cs) + rmg = matrixrandomizer(seed, MersenneTwister(20240620), method = exact) + samples = [rand(rmg) for _ in 1:200n] + counts = collect(values(Dict(s => count(==(s), samples) for s in unique(samples)))) + @test length(counts) == n # full support covered + @test maximum(counts) < 2 * minimum(counts) # roughly uniform +end From 9d9695f42ac84d3bb7d2781cf5b1b0923c79bb8c Mon Sep 17 00:00:00 2001 From: mkborregaard Date: Sat, 20 Jun 2026 22:10:02 +0200 Subject: [PATCH 2/7] Clarify curveball is sequential (correlated), exact/sis independent Co-Authored-By: Claude Opus 4.8 --- README.md | 27 ++++++++++++++++----------- 1 file changed, 16 insertions(+), 11 deletions(-) diff --git a/README.md b/README.md index a8cbfab..cc51508 100644 --- a/README.md +++ b/README.md @@ -39,20 +39,25 @@ m3 = rand!(rmg) ## 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; - repeated calls walk over the space of fixed-margin matrices. + fast in-place trade that performs a Markov step from the current matrix. It is + **sequential**: each call modifies the matrix in place and continues the chain + from where the previous call left off, so successive draws are *correlated*. + The chain's stationary distribution is the uniform one, so over a long run + (with a little burn-in / thinning) it samples fixed-margin matrices uniformly. * `sis` — sequential importance sampling (Harrison & Miller 2013). Each call - draws an *independent* matrix from a fast, near-uniform approximation of the + 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). Draws independent, - *exactly* uniform matrices, 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. - -The `exact` and `sis` methods are uniform over (resp. close to uniform over) all -binary matrices with the given margins. The `sis` proposal is the uniform special + * `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. The `sis` proposal is the uniform special case of the importance sampler of Harrison & Miller; non-uniform (weighted) sampling is a planned extension. From 4b88791c90c11b2c6307de0a86c1b4f209233aff Mon Sep 17 00:00:00 2001 From: mkborregaard Date: Sun, 21 Jun 2026 00:35:33 +0200 Subject: [PATCH 3/7] Speed up SIS: incremental row re-sort and band-local DP clearing Replace the per-column full re-sort with an O(ones-placed) repair of the descending row order (each placed row drops by one, so it only slides past the block holding its old value), tracked via an added inverse-permutation field. Also clear only the touched band of the DP buffer each row instead of the whole vector, and normalise with a reciprocal multiply. ~20% faster per draw with identical output; importance-weight count estimates unchanged. Co-Authored-By: Claude Opus 4.8 --- src/sis.jl | 48 +++++++++++++++++++++++++++++++++++++++--------- 1 file changed, 39 insertions(+), 9 deletions(-) diff --git a/src/sis.jl b/src/sis.jl index 47e2a89..ef3cb3f 100644 --- a/src/sis.jl +++ b/src/sis.jl @@ -76,6 +76,7 @@ and kept in descending `order`; columns are summarised by their `conjugate` mutable struct SISResidual residual::Vector{Int} # residual row sums, indexed by original row order::Vector{Int} # original row indices, sorted by descending residual + pos::Vector{Int} # inverse of `order`: where each row sits within it conjugate::Vector{Int} # conjugate of the columns not yet sampled total::Int # ones still to place sumsq::Int # Σ colsum² over columns not yet sampled @@ -85,8 +86,8 @@ end function SISResidual(rowsums::Vector{Int}, colsums::Vector{Int}) nrows = length(rowsums) - SISResidual(copy(rowsums), sortperm(rowsums, rev = true), - _conjugate(colsums, nrows), + order = sortperm(rowsums, rev = true) + SISResidual(copy(rowsums), order, invperm(order), _conjugate(colsums, nrows), sum(colsums), sum(abs2, colsums), length(colsums), nrows) end @@ -102,12 +103,29 @@ function _advance!(res::SISResidual, k::Int) res end -# Decrement the rows that received a one and restore the descending order. +# Decrement the rows that received a one and restore the descending order. Each +# of those rows dropped by exactly one, so it only has to slide past the block of +# rows still holding its previous value — an O(ones placed) repair rather than a +# full re-sort. Repairing the lowest-valued rows first keeps every swap local. function _place!(res::SISResidual, rows) + order, pos, residual = res.order, res.pos, res.residual @inbounds for r in rows - res.residual[r] -= 1 + residual[r] -= 1 + end + @inbounds for idx in lastindex(rows):-1:firstindex(rows) + r = rows[idx] + i = pos[r] + value = residual[r] + j = i + while j < res.nrows && residual[order[j+1]] > value + j += 1 + end + if j > i + displaced = order[j] + order[i], order[j] = displaced, r + pos[r], pos[displaced] = j, i + end end - sortperm!(res.order, res.residual, rev = true) res end @@ -172,23 +190,35 @@ function _transitions!(ws::SISWorkspace, k::Int) m = length(ws.p) gnext, gcur, S = ws.gnext, ws.gcur, ws.S fill!(gnext, 0.0) + fill!(gcur, 0.0) gnext[k+1] = 1.0 # base case: all rows placed, sum == k + lonext, hinext = k + 1, k + 1 # index range where gnext is nonzero + locur, hicur = 1, 0 # gcur is empty (all zero) @inbounds for i in m:-1:1 lo = i == 1 ? 0 : ws.lo[i-1] # band on the partial sum before row i hi = i == 1 ? 0 : ws.hi[i-1] p = ws.p[i] - fill!(gcur, 0.0) + for t in locur:hicur # clear only gcur's stale band + gcur[t] = 0.0 + end total = 0.0 for s in lo:hi - w1 = p * gnext[s+2] # row i is a one -> partial sum s+1 - w0 = (1 - p) * gnext[s+1] # row i is a zero -> partial sum s + w1 = p * gnext[s+2] # row i is a one -> partial sum s+1 + w0 = (1 - p) * gnext[s+1] # row i is a zero -> partial sum s weight = w0 + w1 gcur[s+1] = weight S[s+1, i] = weight > 0 ? w1 / weight : 0.0 total += weight end - total > 0 && (@views gcur[lo+1:hi+1] ./= total) # rescale to avoid underflow + if total > 0 # rescale to avoid underflow + invtotal = inv(total) + for s in lo:hi + gcur[s+1] *= invtotal + end + end + locur, hicur = lo + 1, hi + 1 gnext, gcur = gcur, gnext + lonext, hinext, locur, hicur = locur, hicur, lonext, hinext end ws end From 199ce5281bba064141bc2c1aa7eff82700695e93 Mon Sep 17 00:00:00 2001 From: mkborregaard Date: Sun, 21 Jun 2026 10:27:49 +0200 Subject: [PATCH 4/7] Curveball generator: SIS warm start + configurable trades (thinning) The curveball generator now warm-starts its chain from one independent SIS draw instead of the observed matrix, so the first sample is decorrelated from the input. A new `trades` keyword on randomize_matrix! and matrixrandomizer sets how many trades each draw performs (default 5*ncols), controlling thinning between successive samples. Refactor the generator to dispatch a single draw on a per-method sampler type (CurveballSampler/SISSampler/ExactCounts) rather than branching on the method enum, keeping each method's draw logic in its own file. Co-Authored-By: Claude Opus 4.8 --- README.md | 15 +++++--- src/RandomBooleanMatrices.jl | 68 ++++++++++++++++++++++-------------- src/curveball.jl | 19 ++++++++-- src/exact.jl | 8 +++++ src/sis.jl | 11 ++++++ test/runtests.jl | 9 +++++ 6 files changed, 98 insertions(+), 32 deletions(-) diff --git a/README.md b/README.md index cc51508..440e8f5 100644 --- a/README.md +++ b/README.md @@ -31,6 +31,10 @@ 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 @@ -40,10 +44,13 @@ m3 = rand!(rmg) * `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 modifies the matrix in place and continues the chain - from where the previous call left off, so successive draws are *correlated*. - The chain's stationary distribution is the uniform one, so over a long run - (with a little burn-in / thinning) it samples fixed-margin matrices uniformly. + **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 diff --git a/src/RandomBooleanMatrices.jl b/src/RandomBooleanMatrices.jl index 039cf91..bcd55d1 100644 --- a/src/RandomBooleanMatrices.jl +++ b/src/RandomBooleanMatrices.jl @@ -13,13 +13,14 @@ include("sis.jl") @enum matrixrandomizations curveball exact sis """ - randomize_matrix!(m [,rng]; method = curveball) + 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). A Markov step from - the current matrix; repeated calls explore the space of fixed-margin matrices. + * `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. @@ -27,36 +28,56 @@ column sums. The algorithm is chosen with `method`: 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) - method == curveball && return _curveball!(m, rng) +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 +# 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 - cache::Union{Nothing, ExactCounts} + 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") -# Precompute whatever a method reuses across draws. Only the exact method needs -# it: the (potentially expensive) matrix count is cached for repeated sampling. -_cache(m::SparseMatrixCSC{Bool, Int}, method::matrixrandomizations) = - method == exact ? _exact_counts(_margins(m)...) : nothing +# 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 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 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 ``` @@ -69,23 +90,18 @@ random2 = rand(rmg) """ matrixrandomizer(m, rng) = error("No matrixrandomizer defined for $(typeof(m))") matrixrandomizer(m) = error("No matrixrandomizer defined for $(typeof(m))") -function matrixrandomizer(m::AbstractMatrix, rng = Xoroshiro128Plus(); method::matrixrandomizations = curveball) +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, method, rng, _cache(sm, method)) + 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) +function matrixrandomizer(m::SparseMatrixCSC{Bool, Int}, rng = Xoroshiro128Plus(); method::matrixrandomizations = curveball, trades = nothing) sm = dropzeros(m) - MatrixGenerator{typeof(rng), SparseMatrixCSC{Bool, Int}}(sm, method, rng, _cache(sm, method)) + 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 exact -# method samples from the cached count; the others delegate to randomize_matrix!. -function _generate!(r::MatrixGenerator{R, SparseMatrixCSC{Bool, Int}}) where R - r.method == exact || return randomize_matrix!(r.m, r.rng, method = r.method) - columns = [Int[] for _ in 1:size(r.m, 2)] - _exact_sample!(columns, r.cache, r.rng) - _writecols!(r.m, columns) -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) diff --git a/src/curveball.jl b/src/curveball.jl index 1aff4f6..332cf4f 100644 --- a/src/curveball.jl +++ b/src/curveball.jl @@ -53,13 +53,13 @@ function _interdif!(v1, v2, inter, dif) ndiff, nshared end -function _curveball!(m::SparseMatrixCSC{Bool, Int}, rng = Random.GLOBAL_RNG) +function _curveball!(m::SparseMatrixCSC{Bool, Int}, rng = Random.GLOBAL_RNG, trades::Int = 5 * size(m, 2)) R, C = size(m) mcs = min(2maximum(diff(m.colptr)), size(m, 1)) not_shared, shared = Vector{Int}(undef, mcs), Vector{Int}(undef, mcs) newa, newb = Vector{Int}(undef, mcs), Vector{Int}(undef, mcs) - for rep ∈ 1:5C + for rep ∈ 1:trades A, B = rand(rng, 1:C,2) # use views directly into the sparse matrix to avoid copying @@ -82,3 +82,18 @@ function _curveball!(m::SparseMatrixCSC{Bool, Int}, rng = Random.GLOBAL_RNG) return m end + +""" + CurveballSampler(trades) + +Generator state for the curveball method: each draw advances the stored matrix by +`trades` curveball trades. The matrix is itself a running Markov chain, so the +generator warm-starts it from an independent SIS draw (see `matrixrandomizer`) to +decorrelate the first sample from the input, and `trades` controls how far the +chain moves between successive draws (thinning). +""" +struct CurveballSampler + trades::Int +end + +_draw!(m::SparseMatrixCSC{Bool, Int}, rng, s::CurveballSampler) = _curveball!(m, rng, s.trades) diff --git a/src/exact.jl b/src/exact.jl index ea3356a..0447dae 100644 --- a/src/exact.jl +++ b/src/exact.jl @@ -237,3 +237,11 @@ function _exact!(m::SparseMatrixCSC{Bool, Int}, rng = Random.GLOBAL_RNG) _exact_sample!(columns, _exact_counts(rowsums, colsums), rng) _writecols!(m, columns) end + +# Generator draw: the cached `ExactCounts` is the generator state, so repeated +# draws reuse the (expensive) count instead of recomputing it. +function _draw!(m::SparseMatrixCSC{Bool, Int}, rng, ex::ExactCounts) + columns = [Int[] for _ in 1:size(m, 2)] + _exact_sample!(columns, ex, rng) + _writecols!(m, columns) +end diff --git a/src/sis.jl b/src/sis.jl index ef3cb3f..6c546b2 100644 --- a/src/sis.jl +++ b/src/sis.jl @@ -287,3 +287,14 @@ function _sis!(m::SparseMatrixCSC{Bool, Int}, rng = Random.GLOBAL_RNG) columns, _ = _sis(rowsums, colsums, rng) _writecols!(m, columns) end + +""" + SISSampler() + +Generator state for the SIS method. SIS draws are independent, so the generator +keeps no precomputed state — each draw resamples the stored matrix from its +margins. +""" +struct SISSampler end + +_draw!(m::SparseMatrixCSC{Bool, Int}, rng, ::SISSampler) = _sis!(m, rng) diff --git a/test/runtests.jl b/test/runtests.jl index 74942f6..9065d5c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -42,6 +42,15 @@ end m4 = rand(rmg) @test m3 != m4 + + # the trade count is controllable; zero trades leaves the matrix untouched + m5 = sprand(Bool, 8, 6, 0.3) + @test randomize_matrix!(copy(m5), method = curveball, trades = 0) == m5 + @test matrixrandomizer(m5, method = curveball, trades = 7).state.trades == 7 + + # the generator warm-starts from an independent draw and preserves margins + csm5, rsm5 = margins(m5) + @test margins(rand(matrixrandomizer(m5, method = curveball, trades = 10))) == (csm5, rsm5) end @testset "sis" begin From 37990b47b0791ceee5710f56a58f62a0db21327b Mon Sep 17 00:00:00 2001 From: mkborregaard Date: Sun, 21 Jun 2026 10:34:32 +0200 Subject: [PATCH 5/7] Fix MatrixGenerator show: extend Base.show, not a shadowing local The bare `show` defined a new RandomBooleanMatrices.show instead of extending Base.show, so the custom display never fired. Qualify it and use print (no trailing newline) so it shows as intended. Co-Authored-By: Claude Opus 4.8 --- src/RandomBooleanMatrices.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/RandomBooleanMatrices.jl b/src/RandomBooleanMatrices.jl index bcd55d1..b5d2f2d 100644 --- a/src/RandomBooleanMatrices.jl +++ b/src/RandomBooleanMatrices.jl @@ -49,7 +49,7 @@ struct MatrixGenerator{R<:AbstractRNG, M} 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 From 1cd9ad60f8f99f9900cc607dbc58dbf0728dfcce Mon Sep 17 00:00:00 2001 From: mkborregaard Date: Sun, 21 Jun 2026 11:23:18 +0200 Subject: [PATCH 6/7] Add weighted importance sampling to SIS MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Generalise the SIS sampler to the Harrison & Miller (2013) weighted target P*(z) ∝ ∏ w[i,j]^z[i,j], of which the uniform sampler is the w≡1 special case. A WeightModel (UniformWeights / WeightMatrix) multiplies each row's odds of a one by a factor v_i built from the elementary symmetric polynomials of the row's remaining weights (carried in logs), leaving the dynamic program otherwise unchanged. Each draw now carries its log importance weight log(∏w^z / Q*), and the public importance_sampler(m, w) returns (matrix, logweight) pairs whose mean weight estimates the weighted count κ. Verified: the mean importance weight recovers brute-force κ within noise, the uniform case reduces to the matrix count, and cv² stays low. Scope: strictly positive weights; structural zeros, weight canonicalisation, and variance-based column ordering are deferred (they affect efficiency, not correctness). Co-Authored-By: Claude Opus 4.8 --- README.md | 42 +++++++++++++- src/RandomBooleanMatrices.jl | 59 +++++++++++++++++++- src/common.jl | 21 +++++++ src/sis.jl | 103 +++++++++++++++++++++++------------ src/weights.jl | 103 +++++++++++++++++++++++++++++++++++ test/runtests.jl | 41 ++++++++++++++ 6 files changed, 330 insertions(+), 39 deletions(-) create mode 100644 src/weights.jl diff --git a/README.md b/README.md index 440e8f5..82ee0e1 100644 --- a/README.md +++ b/README.md @@ -64,9 +64,45 @@ m3 = rand!(rmg) 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. The `sis` proposal is the uniform special -case of the importance sampler of Harrison & Miller; non-uniform (weighted) -sampling is a planned extension. +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 strictly positive 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 # strictly positive 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 +``` + +Structural zeros in `w` (forcing entries to zero), weight canonicalisation, and +the paper's variance-based column ordering are not yet implemented; the importance +weights are correct regardless, these only affect the sampler's efficiency. ## References diff --git a/src/RandomBooleanMatrices.jl b/src/RandomBooleanMatrices.jl index b5d2f2d..d3d7b7d 100644 --- a/src/RandomBooleanMatrices.jl +++ b/src/RandomBooleanMatrices.jl @@ -8,6 +8,7 @@ using StatsBase include("common.jl") include("curveball.jl") include("exact.jl") +include("weights.jl") include("sis.jl") @enum matrixrandomizations curveball exact sis @@ -106,7 +107,63 @@ _generate!(r::MatrixGenerator{R, SparseMatrixCSC{Bool, Int}}) where R = _draw!(r 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) -export randomize_matrix!, matrixrandomizer, matrixrandomizations +""" + 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 strictly positive 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*. + +Structural zeros in `w` are not yet supported. + +# 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()) + 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), 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, importance_sampler export curveball, exact, sis end diff --git a/src/common.jl b/src/common.jl index e841c81..3bd43d6 100644 --- a/src/common.jl +++ b/src/common.jl @@ -45,3 +45,24 @@ function _writecols!(m::SparseMatrixCSC{Bool, Int}, columns) 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 diff --git a/src/sis.jl b/src/sis.jl index 6c546b2..fc8dd66 100644 --- a/src/sis.jl +++ b/src/sis.jl @@ -136,13 +136,16 @@ end """ SISWorkspace -Reusable scratch space for sampling one column: the per-row one-probabilities -`p`, the feasibility band `lo:hi` on the partial column sum after each row, the -forward transition probabilities `S`, and two rolling buffers for the backward -dynamic-programming weights. +Reusable scratch space for sampling one column: the per-row probability `p` of a +one and the (possibly weighted) `one`-weight that biases it, the feasibility band +`lo:hi` on the partial column sum after each row, the forward transition +probabilities `S`, and two rolling buffers for the backward dynamic-programming +weights. For the uniform target `one == p`; a weight model scales `one` by its +per-row factor while leaving the zero-weight `1 - p` alone. """ struct SISWorkspace p::Vector{Float64} + one::Vector{Float64} lo::Vector{Int} hi::Vector{Int} S::Matrix{Float64} @@ -151,14 +154,29 @@ struct SISWorkspace end SISWorkspace(nrows::Int, maxcol::Int) = - SISWorkspace(Vector{Float64}(undef, nrows), Vector{Int}(undef, nrows), - Vector{Int}(undef, nrows), Matrix{Float64}(undef, maxcol + 1, nrows), + SISWorkspace(Vector{Float64}(undef, nrows), Vector{Float64}(undef, nrows), + Vector{Int}(undef, nrows), Vector{Int}(undef, nrows), + Matrix{Float64}(undef, maxcol + 1, nrows), Vector{Float64}(undef, maxcol + 2), Vector{Float64}(undef, maxcol + 2)) -# Score every row for the current column. -function _score!(ws::SISWorkspace, res::SISResidual, scorer::ColumnScorer) +# Score every row for the current column: `p` is the combinatorial probability of +# a one, and `one` is that probability scaled by the weight model's factor (so the +# odds of a one become p·v : 1-p). A weight factor of `Inf` flags a row the +# remaining weights force to a one, which we encode as p = 1. +function _score!(ws::SISWorkspace, res::SISResidual, scorer::ColumnScorer, + model::WeightModel, pos::Int) @inbounds for i in 1:res.nrows - ws.p[i] = scorer(res.residual[res.order[i]]) + row = res.order[i] + ρ = res.residual[row] + p = scorer(ρ) + v = _weightfactor(model, row, ρ, res.ncols, pos) + if isinf(v) + ws.p[i] = 1.0 + ws.one[i] = 1.0 + else + ws.p[i] = p + ws.one[i] = p * v + end end ws end @@ -197,14 +215,15 @@ function _transitions!(ws::SISWorkspace, k::Int) @inbounds for i in m:-1:1 lo = i == 1 ? 0 : ws.lo[i-1] # band on the partial sum before row i hi = i == 1 ? 0 : ws.hi[i-1] - p = ws.p[i] + one = ws.one[i] # weight of a one (odds one : zero = one : 1-p) + zero = 1 - ws.p[i] # weight of a zero for t in locur:hicur # clear only gcur's stale band gcur[t] = 0.0 end total = 0.0 for s in lo:hi - w1 = p * gnext[s+2] # row i is a one -> partial sum s+1 - w0 = (1 - p) * gnext[s+1] # row i is a zero -> partial sum s + w1 = one * gnext[s+2] # row i is a one -> partial sum s+1 + w0 = zero * gnext[s+1] # row i is a zero -> partial sum s weight = w0 + w1 gcur[s+1] = weight S[s+1, i] = weight > 0 ? w1 / weight : 0.0 @@ -223,40 +242,46 @@ function _transitions!(ws::SISWorkspace, k::Int) ws end -# Forward pass: walk the rows, drawing each from its transition probability, and -# accumulate the log-probability of the column under the proposal. -function _draw_column!(rows, ws::SISWorkspace, order, k::Int, rng) +# Forward pass: walk the rows, drawing each from its transition probability. +# Accumulate the log proposal probability `logq` of the column and, for a weighted +# target, the log target weight `logp = Σ log w` over the ones placed. +function _draw_column!(rows, ws::SISWorkspace, order, k::Int, model::WeightModel, label::Int, rng) s = 0 logq = 0.0 + logp = 0.0 @inbounds for i in eachindex(ws.p) p = ws.S[s+1, i] if rand(rng) < p - push!(rows, order[i]) + row = order[i] + push!(rows, row) s += 1 logq += log(p) + logp += _logweight(model, row, label) s == k && break else logq += log(1 - p) end end - logq + logq, logp end """ - _sample_column!(rows, res, k, ws, rng) + _sample_column!(rows, res, k, ws, model, pos, label, rng) -Sample one column of total `k` into `rows` (original row indices), update the -residual problem, and return the log proposal probability of the draw. +Sample the column with label `label` (at sampling position `pos`) and total `k` +into `rows` (original row indices), update the residual problem, and return the +log proposal probability and log target weight `(logq, logp)` of the draw. """ -function _sample_column!(rows, res::SISResidual, k::Int, ws::SISWorkspace, rng) - k == 0 && return 0.0 +function _sample_column!(rows, res::SISResidual, k::Int, ws::SISWorkspace, + model::WeightModel, pos::Int, label::Int, rng) + k == 0 && return 0.0, 0.0 _advance!(res, k) - _score!(ws, res, Canfield(res.nrows, res.ncols, res.total, res.sumsq)) + _score!(ws, res, Canfield(res.nrows, res.ncols, res.total, res.sumsq), model, pos) _band!(ws, res, k) _transitions!(ws, k) - logq = _draw_column!(rows, ws, res.order, k, rng) + logq, logp = _draw_column!(rows, ws, res.order, k, model, label, rng) _place!(res, rows) - logq + logq, logp end # --------------------------------------------------------------------------- @@ -264,27 +289,35 @@ end # --------------------------------------------------------------------------- """ - _sis(rowsums, colsums, rng) + _sis(rowsums, colsums, model, rng) -Sample a fixed-margin binary matrix by sequential importance sampling, returning -the row indices of each column together with the log proposal probability of the -whole matrix. Columns are visited in order of decreasing sum, which the authors -find improves the approximation. +Sample a fixed-margin binary matrix by sequential importance sampling under the +target described by `model` (uniform, or a [`WeightMatrix`]). Returns the row +indices of each column with the log proposal probability `logq` and log target +weight `logp` of the whole matrix; the importance weight of the draw is +`exp(logp - logq)`. Columns are visited in the order the model prescribes (by +decreasing sum), which the authors find improves the approximation. """ -function _sis(rowsums::Vector{Int}, colsums::Vector{Int}, rng) +function _sis(rowsums::Vector{Int}, colsums::Vector{Int}, model::WeightModel, rng) res = SISResidual(rowsums, colsums) ws = SISWorkspace(res.nrows, maximum(colsums, init = 0)) columns = [Int[] for _ in colsums] logq = 0.0 - for j in sortperm(colsums, rev = true) - logq += _sample_column!(columns[j], res, colsums[j], ws, rng) + logp = 0.0 + for (pos, label) in enumerate(_columnorder(model, colsums)) + dlogq, dlogp = _sample_column!(columns[label], res, colsums[label], ws, model, pos, label, rng) + logq += dlogq + logp += dlogp end - columns, logq + columns, logq, logp end +_sis(rowsums::Vector{Int}, colsums::Vector{Int}, rng) = + _sis(rowsums, colsums, UniformWeights(), rng) + function _sis!(m::SparseMatrixCSC{Bool, Int}, rng = Random.GLOBAL_RNG) rowsums, colsums = _margins(m) - columns, _ = _sis(rowsums, colsums, rng) + columns, _, _ = _sis(rowsums, colsums, UniformWeights(), rng) _writecols!(m, columns) end diff --git a/src/weights.jl b/src/weights.jl new file mode 100644 index 0000000..b2258e5 --- /dev/null +++ b/src/weights.jl @@ -0,0 +1,103 @@ +# Weighted targets for sequential importance sampling. +# +# The uniform sampler in sis.jl targets every fixed-margin matrix equally. +# Harrison & Miller (2013) generalise this to the weighted distribution +# +# P*(z) ∝ ∏_{ij} w[i,j]^z[i,j] over fixed-margin binary z, +# +# of which the uniform case (w ≡ 1) is a special case. The generalisation only +# changes the per-row odds of a one in the column being sampled: those odds are +# multiplied by a factor vᵢ built from the row's remaining weights (eq. 15). The +# factor uses the elementary symmetric polynomials of the weights over the +# not-yet-sampled columns, all carried in logs to keep their range in check. +# +# A draw is then no longer (near-)uniform, so it carries an importance weight +# f(z) = ∏ w^z / Q*(z); callers reweight Monte-Carlo estimates by it. Structural +# zeros (w[i,j] == 0), weight canonicalisation, and the variance-based column +# ordering of the paper are not handled here yet (see README); none affect the +# correctness of the importance weights, only the efficiency of the sampler. + +""" + WeightModel + +How the SIS proposal biases a one beyond the combinatorial [`ColumnScorer`]. +`UniformWeights` adds nothing (the uniform target); `WeightMatrix` carries an +entrywise weight matrix and the factors needed to bias the proposal toward it. +New targets dispatch [`_weightfactor`](@ref) and [`_logweight`](@ref) here +without touching the dynamic program. +""" +abstract type WeightModel end + +struct UniformWeights <: WeightModel end + +""" + WeightMatrix(w, rowsums, colsums) + +The target P*(z) ∝ ∏ w[i,j]^z[i,j]. Stores the log weights (for the importance +weight), and, with columns put in sampling order, the log elementary symmetric +polynomials `loge[k+1, i, pos] = log eₖ(w[i, columns from pos onward])` from which +the per-row factor vᵢ is read. Requires strictly positive weights. +""" +struct WeightMatrix <: WeightModel + logw::Matrix{Float64} # log original weights [row, col] (the importance weight) + logwbar::Matrix{Float64} # log weights [row, position] (columns in sampling order) + order::Vector{Int} # column sampling order (column labels) + loge::Array{Float64, 3} # loge[k+1, row, pos] = log eₖ(weights at positions pos:end) +end + +function WeightMatrix(w::AbstractMatrix, rowsums::Vector{Int}, colsums::Vector{Int}) + all(>(0), w) || throw(ArgumentError("weighted SIS currently requires strictly positive weights")) + logw = Matrix{Float64}(log.(float.(w))) + order = sortperm(colsums, rev = true) + logwbar = logw[:, order] + loge = _logesym(logwbar, maximum(rowsums, init = 0)) + WeightMatrix(logw, logwbar, order, loge) +end + +# log(exp(a) + exp(b)), with -Inf absorbing as the log of zero. +function _logaddexp(a::Float64, b::Float64) + a == -Inf && return b + b == -Inf && return a + m = max(a, b) + m + log1p(exp(-abs(a - b))) +end + +# Elementary symmetric polynomials of each row's weights over column suffixes, in +# logs. Built from the empty suffix backwards using eₖ(j:n) = eₖ(j+1:n) + +# w[j]·eₖ₋₁(j+1:n); `loge[k+1, i, pos]` covers columns at positions pos:n. +function _logesym(logwbar::Matrix{Float64}, maxk::Int) + m, n = size(logwbar) + loge = fill(-Inf, maxk + 1, m, n + 1) + loge[1, :, :] .= 0.0 # e₀ ≡ 1 over any suffix + @inbounds for pos in n:-1:1, i in 1:m + lw = logwbar[i, pos] + for k in 1:maxk + loge[k+1, i, pos] = _logaddexp(loge[k+1, i, pos+1], lw + loge[k, i, pos+1]) + end + end + loge +end + +# --------------------------------------------------------------------------- +# Interface consumed by the sampler (sis.jl). Plain-Int arguments keep these +# independent of the residual/workspace types. +# --------------------------------------------------------------------------- + +# Order in which to visit the columns (by decreasing sum). +_columnorder(::UniformWeights, colsums) = sortperm(colsums, rev = true) +_columnorder(model::WeightMatrix, colsums) = model.order + +# Factor multiplying a row's odds of a one in the column at sampling position +# `pos`, given the row's residual sum `ρ` and the number of columns `after` it. +_weightfactor(::UniformWeights, row, ρ, after, pos) = 1.0 +function _weightfactor(model::WeightMatrix, row, ρ, after, pos) + ρ == 0 && return 1.0 + logden = model.loge[ρ+1, row, pos+1] # log eρ over the columns after pos + logden == -Inf && return Inf # eρ == 0: the row is forced to a one here + lognum = model.loge[ρ, row, pos+1] # log eρ₋₁ + (after - ρ + 1) / ρ * exp(model.logwbar[row, pos] + lognum - logden) +end + +# Contribution of a placed one at (row, column) to log ∏ w^z. +_logweight(::UniformWeights, row, label) = 0.0 +_logweight(model::WeightMatrix, row, label) = model.logw[row, label] diff --git a/test/runtests.jl b/test/runtests.jl index 9065d5c..35f886f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,6 +3,8 @@ using SparseArrays using Random using Test +mean(x) = sum(x) / length(x) + Random.seed!(1337) # the row and column sums of a matrix, as plain vectors @@ -19,6 +21,19 @@ function countmatrices(rowsums, colsums) count end +# brute-force weighted count κ = Σ_z ∏ w^z over the binary matrices with the margins +function weightedcount(rowsums, colsums, w) + m, n = length(rowsums), length(colsums) + κ = 0.0 + for bits in 0:(2^(m * n) - 1) + A = [(bits >> ((j - 1) * m + (i - 1))) & 1 for i in 1:m, j in 1:n] + if vec(sum(A, dims = 2)) == rowsums && vec(sum(A, dims = 1)) == colsums + κ += prod(w[i, j]^A[i, j] for i in 1:m, j in 1:n) + end + end + κ +end + @testset "curveball" begin m = sprand(Bool, 8, 6, 0.2) m_old = copy(m) @@ -114,3 +129,29 @@ end @test length(counts) == n # full support covered @test maximum(counts) < 2 * minimum(counts) # roughly uniform end + +@testset "weighted sis" begin + m0 = sparse(Bool[1 0 1; 1 1 0; 0 1 0]) # margins (2,2,1) / (2,2,1) + cs, rs = margins(m0) + w = [1.0 2.0 0.5; 1.5 1.0 2.0; 0.5 1.0 1.0] + + sampler = importance_sampler(m0, w, MersenneTwister(11)) + draws = [rand(sampler) for _ in 1:100_000] + + # every draw is a fixed-margin matrix + @test all(d -> margins(d.matrix) == (cs, rs), draws) + @test draws[1].matrix != draws[2].matrix # independent draws + + # the mean importance weight is an unbiased estimate of the weighted count κ + weights = [exp(d.logweight) for d in draws] + @test isapprox(mean(weights), weightedcount(rs, cs, w), rtol = 0.03) + + # uniform weights reduce to the plain (unweighted) matrix count + uniform = importance_sampler(m0, ones(size(m0)), MersenneTwister(11)) + uweights = [exp(rand(uniform).logweight) for _ in 1:50_000] + @test isapprox(mean(uweights), countmatrices(rs, cs), rtol = 0.02) + + # input validation + @test_throws DimensionMismatch importance_sampler(m0, ones(2, 2)) + @test_throws ArgumentError importance_sampler(m0, [1.0 0.0 1.0; 1.0 1.0 1.0; 1.0 1.0 1.0]) +end From 8a1349e56c48b6b4ae74b668dc6953af96e64c7b Mon Sep 17 00:00:00 2001 From: mkborregaard Date: Sun, 21 Jun 2026 21:28:37 +0200 Subject: [PATCH 7/7] Weighted SIS: structural zeros and Sinkhorn canonicalization MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Allow nonnegative weights. A zero weight is a structural zero forbidding a one at that position: the log-space symmetric polynomials and the v_i factor drop those positions automatically, and a draw the zeros leave no way to complete is returned with logweight = -Inf (importance weight zero) rather than erroring. The v_i binomial term now uses the count of available (positive) positions in the row's remaining columns. Build the proposal from the Sinkhorn-canonical w̄ (eq. 16) by default (canonicalize=false to opt out); P* is invariant to it, so the importance weights are unchanged. Verified on zero-diagonal targets (κ recovered, 0% rejections) and that canonicalize on/off give the same estimate. Also make the exact m3!=m4 test deterministic and robust (it relied on an unseeded generator and a large enough support). Co-Authored-By: Claude Opus 4.8 --- README.md | 23 +++++--- src/RandomBooleanMatrices.jl | 16 ++++-- src/sis.jl | 34 ++++++----- src/weights.jl | 108 ++++++++++++++++++++++++++++------- test/runtests.jl | 29 +++++++--- 5 files changed, 152 insertions(+), 58 deletions(-) diff --git a/README.md b/README.md index 82ee0e1..8af15ca 100644 --- a/README.md +++ b/README.md @@ -76,17 +76,17 @@ P*(z) ∝ ∏ w[i,j] ^ z[i,j] ``` over binary matrices with the given margins. `importance_sampler(m, w)` builds a -sampler for strictly positive 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*`. +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 # strictly positive weights +w = rand(40, 30) .+ 0.5 # nonnegative weights sampler = importance_sampler(m, w) draw = rand(sampler) @@ -100,9 +100,14 @@ wts = exp.(getfield.(draws, :logweight)) Eh = sum(wts .* h.(getfield.(draws, :matrix))) / sum(wts) # for some statistic h ``` -Structural zeros in `w` (forcing entries to zero), weight canonicalisation, and -the paper's variance-based column ordering are not yet implemented; the importance -weights are correct regardless, these only affect the sampler's efficiency. +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 diff --git a/src/RandomBooleanMatrices.jl b/src/RandomBooleanMatrices.jl index d3d7b7d..913a607 100644 --- a/src/RandomBooleanMatrices.jl +++ b/src/RandomBooleanMatrices.jl @@ -132,14 +132,19 @@ 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 strictly positive weights the size of `m`; the uniform -special case is `w` all ones. Each `rand` returns a `(matrix, logweight)` pair: an +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*. -Structural zeros in `w` are not yet supported. +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 ``` @@ -151,11 +156,12 @@ draw.matrix # an independent fixed-margin sample draw.logweight # its log importance weight ``` """ -function importance_sampler(m::AbstractMatrix, w::AbstractMatrix, rng = Xoroshiro128Plus()) +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), rng) + WeightedSIS(rowsums, colsums, WeightMatrix(w, rowsums, colsums; canonicalize), rng) end function Random.rand(s::WeightedSIS) diff --git a/src/sis.jl b/src/sis.jl index fc8dd66..831cbe1 100644 --- a/src/sis.jl +++ b/src/sis.jl @@ -169,7 +169,7 @@ function _score!(ws::SISWorkspace, res::SISResidual, scorer::ColumnScorer, row = res.order[i] ρ = res.residual[row] p = scorer(ρ) - v = _weightfactor(model, row, ρ, res.ncols, pos) + v = _weightfactor(model, row, ρ, pos) if isinf(v) ws.p[i] = 1.0 ws.one[i] = 1.0 @@ -203,7 +203,8 @@ end # feasibility band into forward transition probabilities. `g[s]` is the total # weight of feasible completions reaching column total `k` from partial sum `s`; # `S[s+1, i]` is the resulting probability that row `i` is a one given partial sum -# `s` before it. +# `s` before it. Returns `false` if the column has no feasible completion — only +# possible when structural zeros forbid enough ones — so the caller can reject. function _transitions!(ws::SISWorkspace, k::Int) m = length(ws.p) gnext, gcur, S = ws.gnext, ws.gcur, ws.S @@ -229,17 +230,16 @@ function _transitions!(ws::SISWorkspace, k::Int) S[s+1, i] = weight > 0 ? w1 / weight : 0.0 total += weight end - if total > 0 # rescale to avoid underflow - invtotal = inv(total) - for s in lo:hi - gcur[s+1] *= invtotal - end + total > 0 || return false # no completion: a structural-zero dead-end + invtotal = inv(total) # rescale to avoid underflow + for s in lo:hi + gcur[s+1] *= invtotal end locur, hicur = lo + 1, hi + 1 gnext, gcur = gcur, gnext lonext, hinext, locur, hicur = locur, hicur, lonext, hinext end - ws + true end # Forward pass: walk the rows, drawing each from its transition probability. @@ -269,19 +269,20 @@ end _sample_column!(rows, res, k, ws, model, pos, label, rng) Sample the column with label `label` (at sampling position `pos`) and total `k` -into `rows` (original row indices), update the residual problem, and return the -log proposal probability and log target weight `(logq, logp)` of the draw. +into `rows` (original row indices), update the residual problem, and return +`(logq, logp, feasible)`: the log proposal probability, the log target weight, and +whether a completion existed (`false` is a structural-zero dead-end to reject). """ function _sample_column!(rows, res::SISResidual, k::Int, ws::SISWorkspace, model::WeightModel, pos::Int, label::Int, rng) - k == 0 && return 0.0, 0.0 + k == 0 && return 0.0, 0.0, true _advance!(res, k) _score!(ws, res, Canfield(res.nrows, res.ncols, res.total, res.sumsq), model, pos) _band!(ws, res, k) - _transitions!(ws, k) + _transitions!(ws, k) || return 0.0, 0.0, false logq, logp = _draw_column!(rows, ws, res.order, k, model, label, rng) _place!(res, rows) - logq, logp + logq, logp, true end # --------------------------------------------------------------------------- @@ -296,7 +297,9 @@ target described by `model` (uniform, or a [`WeightMatrix`]). Returns the row indices of each column with the log proposal probability `logq` and log target weight `logp` of the whole matrix; the importance weight of the draw is `exp(logp - logq)`. Columns are visited in the order the model prescribes (by -decreasing sum), which the authors find improves the approximation. +decreasing sum), which the authors find improves the approximation. With +structural zeros a draw can dead-end; it is returned with `logq = Inf` so that its +importance weight `exp(logp - logq)` is zero. """ function _sis(rowsums::Vector{Int}, colsums::Vector{Int}, model::WeightModel, rng) res = SISResidual(rowsums, colsums) @@ -305,7 +308,8 @@ function _sis(rowsums::Vector{Int}, colsums::Vector{Int}, model::WeightModel, rn logq = 0.0 logp = 0.0 for (pos, label) in enumerate(_columnorder(model, colsums)) - dlogq, dlogp = _sample_column!(columns[label], res, colsums[label], ws, model, pos, label, rng) + dlogq, dlogp, feasible = _sample_column!(columns[label], res, colsums[label], ws, model, pos, label, rng) + feasible || return columns, Inf, logp logq += dlogq logp += dlogp end diff --git a/src/weights.jl b/src/weights.jl index b2258e5..fb44298 100644 --- a/src/weights.jl +++ b/src/weights.jl @@ -12,10 +12,15 @@ # not-yet-sampled columns, all carried in logs to keep their range in check. # # A draw is then no longer (near-)uniform, so it carries an importance weight -# f(z) = ∏ w^z / Q*(z); callers reweight Monte-Carlo estimates by it. Structural -# zeros (w[i,j] == 0), weight canonicalisation, and the variance-based column -# ordering of the paper are not handled here yet (see README); none affect the -# correctness of the importance weights, only the efficiency of the sampler. +# f(z) = ∏ w^z / Q*(z); callers reweight Monte-Carlo estimates by it. +# +# Two refinements from the paper are included. Structural zeros (w[i,j] == 0, which +# forbid a one at that position) are handled: the symmetric polynomials and vᵢ drop +# those positions automatically (a zero weight is log -∞), and a draw that the +# zeros leave no way to complete is rejected with importance weight zero. And the +# weights are Sinkhorn-balanced into the canonical w̄ ∈ Λ(w) of eq. (16), which P* +# is invariant to but which makes the proposal scale-free and lower-variance. +# (Not yet done: the paper's variance-based column ordering.) """ WeightModel @@ -31,27 +36,36 @@ abstract type WeightModel end struct UniformWeights <: WeightModel end """ - WeightMatrix(w, rowsums, colsums) + WeightMatrix(w, rowsums, colsums; canonicalize = true) -The target P*(z) ∝ ∏ w[i,j]^z[i,j]. Stores the log weights (for the importance -weight), and, with columns put in sampling order, the log elementary symmetric -polynomials `loge[k+1, i, pos] = log eₖ(w[i, columns from pos onward])` from which -the per-row factor vᵢ is read. Requires strictly positive weights. +The target P*(z) ∝ ∏ w[i,j]^z[i,j], for nonnegative weights `w` (zeros forbid a +one at that position). Stores the log weights (for the importance weight) and, +with columns in sampling order, the per-row count of available (positive) +positions and the log elementary symmetric polynomials +`loge[k+1, i, pos] = log eₖ(positive weights of row i from position pos onward)`, +from which the per-row factor vᵢ is read. The proposal is built from the +Sinkhorn-canonical `w̄` unless `canonicalize = false`; the importance weight always +uses the original `w`. """ struct WeightMatrix <: WeightModel - logw::Matrix{Float64} # log original weights [row, col] (the importance weight) - logwbar::Matrix{Float64} # log weights [row, position] (columns in sampling order) + logw::Matrix{Float64} # log original weights [row, col] (the importance weight) + logwbar::Matrix{Float64} # log canonical weights [row, pos] (columns in sampling order) + navail::Matrix{Int} # navail[row, pos] = #positive weights at positions pos:end order::Vector{Int} # column sampling order (column labels) loge::Array{Float64, 3} # loge[k+1, row, pos] = log eₖ(weights at positions pos:end) end -function WeightMatrix(w::AbstractMatrix, rowsums::Vector{Int}, colsums::Vector{Int}) - all(>(0), w) || throw(ArgumentError("weighted SIS currently requires strictly positive weights")) - logw = Matrix{Float64}(log.(float.(w))) +function WeightMatrix(w::AbstractMatrix, rowsums::Vector{Int}, colsums::Vector{Int}; + canonicalize::Bool = true) + all(>=(0), w) || throw(ArgumentError("weights must be nonnegative")) + any(>(0), w) || throw(ArgumentError("weights must have a positive entry")) + logw = Matrix{Float64}(log.(float.(w))) # original weights, for the weight + wbar = canonicalize ? _canonical(w) : Matrix{Float64}(float.(w)) order = sortperm(colsums, rev = true) - logwbar = logw[:, order] + logwbar = log.(wbar)[:, order] # proposal weights, in sampling order + navail = _suffixcount(logwbar) loge = _logesym(logwbar, maximum(rowsums, init = 0)) - WeightMatrix(logw, logwbar, order, loge) + WeightMatrix(logw, logwbar, navail, order, loge) end # log(exp(a) + exp(b)), with -Inf absorbing as the log of zero. @@ -62,9 +76,21 @@ function _logaddexp(a::Float64, b::Float64) m + log1p(exp(-abs(a - b))) end +# Number of positive weights in each row over column suffixes: `navail[i, pos]` +# counts positions pos:n. A weight is positive exactly when its log is finite. +function _suffixcount(logwbar::Matrix{Float64}) + m, n = size(logwbar) + navail = zeros(Int, m, n + 1) + @inbounds for pos in n:-1:1, i in 1:m + navail[i, pos] = navail[i, pos+1] + (logwbar[i, pos] > -Inf) + end + navail +end + # Elementary symmetric polynomials of each row's weights over column suffixes, in # logs. Built from the empty suffix backwards using eₖ(j:n) = eₖ(j+1:n) + -# w[j]·eₖ₋₁(j+1:n); `loge[k+1, i, pos]` covers columns at positions pos:n. +# w[j]·eₖ₋₁(j+1:n); `loge[k+1, i, pos]` covers columns at positions pos:n. A zero +# weight (log -∞) contributes nothing, so structural zeros drop out naturally. function _logesym(logwbar::Matrix{Float64}, maxk::Int) m, n = size(logwbar) loge = fill(-Inf, maxk + 1, m, n + 1) @@ -78,6 +104,43 @@ function _logesym(logwbar::Matrix{Float64}, maxk::Int) loge end +""" + _canonical(w; tol = 1e-10, maxiter = 1000) + +Sinkhorn-balance `w` into the canonical `w̄ ∈ Λ(w)` of Harrison & Miller (eq. 16): +the scaling `w̄ = αᵢ βⱼ wᵢⱼ` whose row and column sums equal the number of positive +entries in that row/column. `P*` is invariant to this scaling, but it makes the +proposal scale-invariant and tends to reduce its variance. +""" +function _canonical(w::AbstractMatrix; tol::Float64 = 1e-10, maxiter::Int = 1000) + m, n = size(w) + wbar = Matrix{Float64}(float.(w)) + ni = [count(>(0), view(wbar, i, :)) for i in 1:m] # positive entries per row + mj = [count(>(0), view(wbar, :, j)) for j in 1:n] # positive entries per column + for _ in 1:maxiter + for i in 1:m # scale each row to sum nᵢ + s = sum(view(wbar, i, :)) + s > 0 || continue + f = ni[i] / s + @inbounds for j in 1:n + wbar[i, j] *= f + end + end + drift = 0.0 + for j in 1:n # scale each column to sum mⱼ + s = sum(view(wbar, :, j)) + s > 0 || continue + drift = max(drift, abs(s - mj[j])) + f = mj[j] / s + @inbounds for i in 1:m + wbar[i, j] *= f + end + end + drift < tol && break + end + wbar +end + # --------------------------------------------------------------------------- # Interface consumed by the sampler (sis.jl). Plain-Int arguments keep these # independent of the residual/workspace types. @@ -88,14 +151,17 @@ _columnorder(::UniformWeights, colsums) = sortperm(colsums, rev = true) _columnorder(model::WeightMatrix, colsums) = model.order # Factor multiplying a row's odds of a one in the column at sampling position -# `pos`, given the row's residual sum `ρ` and the number of columns `after` it. -_weightfactor(::UniformWeights, row, ρ, after, pos) = 1.0 -function _weightfactor(model::WeightMatrix, row, ρ, after, pos) +# `pos`, given the row's residual sum `ρ`. `0` forbids a one (structural zero); +# `Inf` forces one (the remaining positive weights leave no alternative). +_weightfactor(::UniformWeights, row, ρ, pos) = 1.0 +function _weightfactor(model::WeightMatrix, row, ρ, pos) ρ == 0 && return 1.0 + logw = model.logwbar[row, pos] + logw == -Inf && return 0.0 # structural zero: no one here logden = model.loge[ρ+1, row, pos+1] # log eρ over the columns after pos logden == -Inf && return Inf # eρ == 0: the row is forced to a one here lognum = model.loge[ρ, row, pos+1] # log eρ₋₁ - (after - ρ + 1) / ρ * exp(model.logwbar[row, pos] + lognum - logden) + (model.navail[row, pos+1] - ρ + 1) / ρ * exp(logw + lognum - logden) end # Contribution of a placed one at (row, column) to log ∏ w^z. diff --git a/test/runtests.jl b/test/runtests.jl index 35f886f..81399cc 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -103,15 +103,13 @@ end randomize_matrix!(m, method = exact) @test margins(m) == (csm, rsm) - m2 = rand(0:1, 5, 6) + m2 = rand(MersenneTwister(2718), 0:1, 5, 6) csm, rsm = margins(m2) - rmg = matrixrandomizer(m2, method = exact) # counts once, reused below - m3 = rand(rmg) - m4 = rand(rmg) + rmg = matrixrandomizer(m2, MersenneTwister(2719), method = exact) # counts once, reused + draws = [rand(rmg) for _ in 1:10] - @test margins(m3) == (csm, rsm) - @test margins(m4) == (csm, rsm) - @test m3 != m4 + @test all(d -> margins(d) == (csm, rsm), draws) + @test length(unique(draws)) > 1 # independent draws vary # a fully determined problem returns its single solution determined = sparse(Bool[1 1 1; 1 1 0; 1 0 0]) @@ -151,7 +149,22 @@ end uweights = [exp(rand(uniform).logweight) for _ in 1:50_000] @test isapprox(mean(uweights), countmatrices(rs, cs), rtol = 0.02) + # canonicalisation changes the proposal but not the (unbiased) estimate + off = importance_sampler(m0, w, MersenneTwister(11); canonicalize = false) + offweights = [exp(rand(off).logweight) for _ in 1:100_000] + @test isapprox(mean(offweights), weightedcount(rs, cs, w), rtol = 0.03) + + # structural zeros: a zero weight forbids a one there, and κ is still recovered + wz = [0.0 2.0 1.0; 1.5 0.0 2.0; 1.0 0.5 0.0] # zero diagonal + derange = sparse(Bool[0 1 0; 0 0 1; 1 0 0]) # margins (1,1,1), off-diagonal + zsampler = importance_sampler(derange, wz, MersenneTwister(3)) + zdraws = [rand(zsampler) for _ in 1:100_000] + @test all(d -> all(i -> !d.matrix[i, i], 1:3), zdraws) # never a one on a zero + zweights = [exp(d.logweight) for d in zdraws] + @test isapprox(mean(zweights), weightedcount([1, 1, 1], [1, 1, 1], wz), rtol = 0.05) + # input validation @test_throws DimensionMismatch importance_sampler(m0, ones(2, 2)) - @test_throws ArgumentError importance_sampler(m0, [1.0 0.0 1.0; 1.0 1.0 1.0; 1.0 1.0 1.0]) + @test_throws ArgumentError importance_sampler(m0, fill(-1.0, 3, 3)) # negative weight + @test_throws ArgumentError importance_sampler(m0, zeros(3, 3)) # no positive entry end