Intro
Formulation
Naming
When I asked ChatGPT how different formulations of this problem should be called, it came up with these:
- “nonnegative integer matrices with fixed row and column sums”
- “contingency tables with fixed margins”
- “integer points in a transportation polytope”
- “matrix with prescribed margins”
- “bounded contingency tables”
- “fiber of a contingency table”
While I am going to leave these suggestions here for the entertainment of the reader—and maybe also the benefits of search engine optimization—I have elected to call this problem “Sampling Random Matrices under Fixed-Margins”.
Ideas
- Card Dealing: For each integer value of the fixed matrix mass, assign that mass to a random unfilled bin.
- Index Sampling: Create a function mapping between N and all of the possible worlds. Then sample from the space 0..num_worlds and simply apply the function to the sample.
- Noise Scaling: Generate a matrix of noise values, then scale rows and columns of that matrix to match the target margins. Afterwards, clamp to nearest integers and normalize to maintain margin invariants.
- Mass Mixing: Generate a valid world, and then perform “mass swaps” to traverse over the graph of possible worlds.
Card Dealing
The implementation of the card dealing algorithm looks something like this:
pub fn card_dealing<const M: usize, const N: usize, const T: usize>(
row_margins: [usize; M],
column_margins: [usize; N],
rng: &mut SmallRng,
) -> [[usize; N]; M] {
debug_assert!(T == M * N);
let mass_total: usize = row_margins.iter().sum();
debug_assert!(column_margins.iter().sum::<usize>() == mass_total);
let mut output = [[0; N]; M];
let mut remaining_row_mass = row_margins;
let mut remaining_col_mass = column_margins;
let mut active_rows = Vec::with_capacity(M);
let mut active_cols = Vec::with_capacity(N);
active_rows.extend((0..M).filter(|&row| remaining_row_mass[row] > 0));
active_cols.extend((0..N).filter(|&col| remaining_col_mass[col] > 0));
for _ in 0..mass_total {
let num_cols = active_cols.len();
let index = rng.random_range(0..active_rows.len() * num_cols);
let row_idx = index / num_cols;
let col_idx = index % num_cols;
let row = active_rows[row_idx];
let col = active_cols[col_idx];
output[row][col] += 1;
remaining_row_mass[row] -= 1;
remaining_col_mass[col] -= 1;
if remaining_row_mass[row] == 0 {
active_rows.swap_remove(row_idx);
}
if remaining_col_mass[col] == 0 {
active_cols.swap_remove(col_idx);
}
}
debug_assert!(output.iter().flatten().sum::<usize>() == mass_total);
return output;
}When testing the algorithm, I also experimented with avoiding the integer division and modulo by generating two numbers, one for the row, and another for the column.
To absolutely no-ones surprise, generating another random number was slower than my first version of the algorithm. What was interesting, however, was how much slower.
Here is the output for one of the regression tests of the function, which I ran on the alternative implementation (The regression/performance testing for the Monte Cardo project was done with criterion.rs).
Across the board, regression tests showed the function to be over 50% slower than the single rng version. While it’s nice we can optimize the function so much, the bigger point is this:
Since changing 1 rng call to 2 rng calls increased the runtime by half of the previous runtime, that suggests that around half of the runtime for this algorithm is spent generating random numbers.
If we can find a way to generate less random numbers for each world we sample, we may be able to find a much faster algorithm.
The output of this algorithm is not uniformly distributed over the possible worlds.