Efficient uniform sampling of certain Boolean matrices

Question

I want to uniformly sample $(0,1)$-matrices with certain constraints. In particular, I want sum of each row to be rowsum and the sums of columns to equal columnsum. If uniform sampling wouldn't be a requirement this could be solved using SatisfiabilityInstances, but to my understanding there's no way to guarantee sampling properties with it.

I have implemented the following, both ugly and slow version of a function to solve this problem:

ClearAll@matrixSample;
matrixSample[rows_, columns_, rowsum_, columnsum_] :=
 Enclose[
  (* Check consistency of arguments *)
  ConfirmAssert[rows rowsum == Total@columnsum, 
   "Conflicting rows, rowsum and columnsum"];
  With[
    {bitpositions =
      (* List of individual bit positions sets to perform over rows *)
      Flatten@MapIndexed[Table[First@#2 - 1, #1] &, columnsum]},
    (* Retry until successful *)
    NestWhile[
     Enclose[
       (* Add a bit at a time *)
       Fold[#1 +
          Outer[Times,
           (* Add to a random row which has not reached rowsum bit limit
              and doesn't have this bit set *)
           UnitVector[rows,
            (* Confirm that such a row is found,
               otherwise start from scratch *)
            First@RandomChoice@ConfirmBy[
               Position[#1, {l_, v_} /; 
                 l < rowsum && BitGet[v, #2] == 0], # != {} &]],
           {1, 2^#2}] &,
        (* Initial {bitcount, bitvector} values *)
        ConstantArray[{0, 0}, rows],
        bitpositions]] &,
     Failure[x, x], FailureQ]]
   (* Decode integers to rows *)
   // Map[Reverse@IntegerDigits[Last@#, 2, columns] &]]

It does work:

(m = matrixSample[8, 10, 7, {7, 6, 6, 6, 6, 6, 5, 5, 5, 4}]) // MatrixForm

$\left( \begin{array}{cccccccccc} 1 & 1 & 1 & 0 & 1 & 1 & 1 & 1 & 0 & 0 \\ 1 & 1 & 1 & 0 & 1 & 1 & 0 & 1 & 1 & 0 \\ 1 & 1 & 1 & 1 & 1 & 1 & 1 & 0 & 0 & 0 \\ 1 & 1 & 1 & 1 & 0 & 0 & 1 & 0 & 1 & 1 \\ 1 & 0 & 1 & 1 & 0 & 1 & 0 & 1 & 1 & 1 \\ 1 & 1 & 0 & 1 & 1 & 1 & 1 & 1 & 0 & 0 \\ 1 & 1 & 0 & 1 & 1 & 0 & 0 & 1 & 1 & 1 \\ 0 & 0 & 1 & 1 & 1 & 1 & 1 & 0 & 1 & 1 \\ \end{array} \right)$

Results fulfill specified constraints:

{Plus @@ Transpose@m, Plus @@ m}

{{7, 7, 7, 7, 7, 7, 7, 7}, {7, 6, 6, 6, 6, 6, 5, 5, 5, 4}}

But performance is abysmal (I'd want to generate tens of millions of samples of larger matrices):

RepeatedTiming[
 matrixSample[8, 10, 7, {7, 6, 6, 6, 6, 6, 5, 5, 5, 4}];, 10]

{0.0141238, Null}

How to improve or rewrite this function?

EDIT: Some performance improvement can be acquired by penalising sampling of rows with have lots of bits set when choosing a row to add a new bit. This lowers chances of failure:

ClearAll@matrixSample;
matrixSample[rows_, columns_, rowsum_, columnsum_] :=
 Enclose[
  (* Check consistency of arguments *)
  ConfirmAssert[rows rowsum == Total@columnsum, 
   "Conflicting rows, rowsum and columnsum"];
  With[
    {bitpositions =
      (* List of individual bit positions sets to perform over rows *)
      Flatten@MapIndexed[Table[First@#2 - 1, #1] &, columnsum]},
    (* Retry until successful *)
    NestWhile[
     Enclose[
       (* Add a bit at a time *)
       Fold[#1 +
          Outer[Times,
           (* Add to a random row which has not reached rowsum bit limit
              and doesn't have this bit set *)
           UnitVector[rows,
            (* Confirm that such a row is found,
               otherwise start from scratch *)
            RandomChoice[
             (* Penalise rows with lots of entries
                on candidate sampling *)
             Rule @@ Transpose@With[{bps = #1},
                {(rowsum - bps[[First@#, 1]])/rowsum, First@#} & /@
                ConfirmBy[
                 Position[#1, {l_, v_} /; 
                   l < rowsum && BitGet[v, #2] == 0], # != {} &]]]],
           {1, 2^#2}] &,
        (* Initial {bitcount, bitvector} values *)
        ConstantArray[{0, 0}, rows],
        bitpositions]] &,
     Failure[x, x], FailureQ]]
   (* Decode integers to rows *)
   // Map[Reverse@IntegerDigits[Last@#, 2, columns] &]]

But does it skew sampling? At least it speeds it up:

RepeatedTiming[
 matrixSample[8, 10, 7, {7, 6, 6, 6, 6, 6, 5, 5, 5, 4}];, 10]

{0.00621864, Null}

EDIT 2: Further performance improvement can be achieved by adding bits to each column in a single operation:

ClearAll@matrixSample;
matrixSample[rows_, columns_, rowsum_, columnsum_] :=
 Enclose[
  (* Check consistency of arguments *)
  ConfirmAssert[rows rowsum == Total@columnsum, 
   "Conflicting rows, rowsum and columnsum"];
  With[
    {bitcounts =
      (* {bitposition, bitcount} list *)
      Transpose@{Range@columns - 1, columnsum}},
    (* Retry until successful *)
    NestWhile[
     Enclose[
       (* Add a column of bits at a time *)
       Fold[#1 +
          Outer[Times,
           (* Add to random rows which have not reached
              rowsum bit limit *)
           ReplacePart[ConstantArray[0, rows],
            (* Confirm that enough such rows are found, 
               otherwise start from scratch *)
            RandomSample[
              (* Penalise rows with lots of entries on
                 candidate sampling *)
              Rule @@ Transpose@With[{bps = #1},
                 {(rowsum - bps[[First@#, 1]])/rowsum, #} & /@
                 With[{count = Last@#2},
                  ConfirmBy[
                    Position[#1, {l_, _} /; l < rowsum],
                    Length@# >= count &]]], Last@#2] -> 1],
           {1, 2^First@#2}] &,
        (* Initial {bitcount, bitvector} values *)
        ConstantArray[{0, 0}, rows],
        bitcounts]] &,
     Failure[x, x], FailureQ]]
   (* Decode integers to rows *)
   // Map[Reverse@IntegerDigits[Last@#, 2, columns] &]]

This is roughly 8 times faster than the original:

RepeatedTiming[
 matrixSample[8, 10, 7, {7, 6, 6, 6, 6, 6, 5, 5, 5, 4}];, 10]

{0.00177051, Null}

Still, it's godawful code. Would there be any better approach?

Answer 1

This is graph theory in disguise: you are asking to uniformly sample bipartite graphs with given degrees.

There are two major approaches to similar sampling problems with prescribed degrees: the configuration model, and MCMC sampling based on degree-preserving edge switches. The latter is a bit more involved, but the first one is easily implemented.

I will use my IGraph/M package for this.

Needs["IGraphM`"]

IGBigraphicalQ tests if such a graph exists at all based on the Gale-Ryser theorem.

Let us create two degree sequences, i.e. row sums and column sums, in such a way that there is actually a graph that realizes them (i.e. the pair is bigraphical). We can do this by generating a random bipartite graph and taking its degrees.

n1 = 5; n2 = 8; m = 20;
{d1, d2} = TakeDrop[VertexDegree@IGBipartiteGameGNM[n1, n2, m], n1]

n1 and n2 are row and column numbers, m is the total number of matrix elements, d1 and d2 is the row and column sums we obtained.

The idea is to imagine each node with as many dangling connections, stubs, as its degree. Then we randomly connect up these stubs, and hope to not get a multigraph. If we do, we discard it and re-start.

Let us represent a stub with the index of the vertex it belongs to:

makeStubs[d_] := Join @@ MapIndexed[ConstantArray[First[#2], #1] &, d]

Then keep sampling until we get a simple graph (i.e. no matrix elements are greater than 1). Wanting to do such a thing is so common that I even included a convenience function for it in IGraph/M: IGTryUntil keeps trying until a condition is satisfied.

edges = IGTryUntil[DuplicateFreeQ]@
  Transpose@{RandomSample@makeStubs[d1], RandomSample@makeStubs[d2]}

Now can build the matrix:

mat = SparseArray[edges -> ConstantArray[1, Length[edges]], {n1, n2}];
MatrixPlot[mat]

Check the constraints:

{Total /@ mat == d1, Total[mat] == d2}
(* {True, True} *)

---

Note that this approach does sample exactly uniformly. The problem with it is the many retries that might be needed. If m is large, we may in practice never randomly stumble upon a realization that does not have multi-edges (even if one exists), so the algorithm will never stop in practice.

In this case, you can do one of three things:

• Give up uniformity and try to avoid multi-edges in some ad-hoc way.
• Use an MCMC sampler based on edge switches. I never looked into this for bipartite graphs, but one would need to check that all realizations are reachable with the chosen kind of edge swap, and construct the Markov chain carefully in order to ensure uniformity. In principle, even then the question remains: does the chain mix fast enough? Note that for this you would first need to construct a single realization, which should be doable with an adapted Havel-Hakimi type construction.
• There is a third type of approach, where stubs are connected up in such a way as to avoid multi-edges, and uniformity is broken, but the actual (non-uniform) distribution can be derived and corrected for. It is a type of importance sampling. I think our recent paper gives an intuitive description and has the references to the original work on this type of approach. But implementing this for the bipartite case would be quite involved.

Answer 2

This answer will construct a single matrix and base on that, populate some other matrices. Sampling is not involved.

You might find a faster and easier solution in the world of mathematics. Here is my solution (a programmer mindset):

f::invalid = "Can't produce a result.";

f[row_, col_] := 
 Block[{temp = ConstantArray[0, Length@col], tcol = col, result = {}},
  
  result = 
   Rest@FoldList[
      ReplacePart[temp, 1, List /@ Ordering[tcol = tcol - #1, -#2]] &,
       temp, row] // Quiet;
  
  (*validate the output*)
  If[Total[result] == col && Total[Transpose@result] == row, 
   Return[result], Message[f::invalid]]]

You just need to provide the sum of rows as the first argument and the sum of columns as the second:

f[{1, 2, 1}, {2, 1, 1}]

Output (in TableForm with heading):

Speed

f[ConstantArray[7, 8], {7, 6, 6, 6, 6, 6, 5, 5, 5, 4}]; // RepeatedTiming

Output:

{0.0000981148, Null}

Result:

How it works?

Let's solve this example step by step:

row = {1, 2, 1};
col = {2, 1, 1};

Create a list of zero as a template:

temp = ConstantArray[0, 3];

The first row needs one cell to be 1, so we'll replace one of the zeros in the template in the position of the highest column given by Ordering:

row1 = ReplacePart[temp, Rule[#, 1] & /@ Ordering[col, -1]]

(*Out: {1, 0, 0} *)

The second row needs two cell, we'll replace 2 zero in the highest column after subtraction with the first calculated row:

row2 = ReplacePart[temp, Rule[#, 1] & /@ Ordering[col - row1, -2]]

(*Out: {0, 1, 1} *)

The third row needs one cell and we'll use the highest position in (col - row1 - row2):

r3 = ReplacePart[temp, Rule[#, 1] & /@ Ordering[col - r1 - r2, -1]]

(*Out: {1, 0, 0} *)

Basically in each step, we try to fill the column that has a larger number. if we had {3,2,1} for columns and {1,2,3} for rows. first row would be {1,0,0} because Ordering[{3, 2, 1}, -1] gives {1} the position of the column to fill.

Since in each step, we do the same calculation along with the previous answer, FoldList becomes useful. Also, the final code uses another form of ReplacePart which is faster.

Sometimes these steps produce inaccurate results, so to prevent that, I put validation at the end.

Where they're inaccurate?

Because of the validation at the end, you'll be warned with Null as the output but if you remove that section for a moment and try:

(* Modified version: without validation *)

f[{3, 2, 1}, {3, 0, 0}]

Result:

You'll see because of using Ordering, there is no condition to check whether filling the columns are allowed or not and that, may produce unexpected results.

Populating some other matrices

By rotating 1 and 0 in a rectangular shape, we could create a new matrix with satisfying the conditions. For example for the code f[{1, 2, 1}, {2, 1, 1}] we have:

Considering this rectangle:

we could swap values:

Code:

populateSwaps[table_List, n_ : Infinity] :=
 Tuples /@
    Take[Flatten[#, 1], UpTo[n]] &@
  (Tuples /@
    MapAt[
     Take[Permutations[Flatten@Position[#, 1], {2}], UpTo[n]] &
     ,
     Cases[#, {_, x_List} /; Count[x, 1] > 1, {1}, UpTo[n]] &@
      Map[{List@#, Plus @@ Extract[table, List /@ #]} &, 
       Permutations[Range[Length@table], {2}]]
     , {All, 2}])


swap[data_List, swap_List] := Block[{tdata = data},
  ({tdata[[Sequence @@ #2]], 
       tdata[[Sequence @@ #1]]} = {data[[Sequence @@ #1]], 
       data[[Sequence @@ #2]]}) & @@@
   Transpose[{swap[[;; ;; 2]], swap[[2 ;; ;; 2]]}];
  Return[tdata]]

swapValues[data_, n_ : Infinity] := 
 swap[data, #] & /@ populateSwaps[data, n]

Use swapValues to populate all the possible swaps from the base matrix:

Length@swapValues@f[{1, 2, 1}, {2, 1, 1}] // RepeatedTiming

(*Out: {0.00114618, 24} *)

You could also limit the number of swaps with the second argument:

swapValues[f[{1, 2, 1}, {2, 1, 1}], 1]; // RepeatedTiming

(*Out: {0.000180499, Null} *)

Result:

Also, I recommend posting this question on Code Golf, that's where you'll find super creative ways in different languages.