Sampling Random Matrices under Fixed-Margins — Part 1/3

Introduction

“Who has what cards?”

If you have ever played a card game (and I assume most of you have), then you have asked yourself a version of this question. You may know how many cards each player has, and you may even know how many of each card type is left, but who has which cards?

Barring straight-up cheating, there usually isn’t an easy way to know that information.

But you can make some educated guesses about what they might have and—maybe even more usefully—you can know what they do not have.

As part of my ongoing upgrades to monte-cardo (my Information Set Monte Carlo Tree Search solver for ladder shedding games, read more about that here) I found myself needing to ask the same question.

For more details, motivated readers should check out the appropriate blog and project posts, but a quick explanation will suffice for our purposes:

The algorithm relies on generating many different “hypotheses” for what hands the other players might have. Each of these hypothetical game states is what we call a “possible world”, and we evaluate each of our available moves over a bunch of these possible worlds so that we can get a better idea of how good that move is.

The problem that I am addressing in this blog post is the generation of those possible worlds. We want to do it quickly, and we want to generate a wide range of possible worlds (more on what that means in a bit).

But how do we create these possible worlds? For our purposes there is only one unknown factor in the game, and that is other players’ hands. Luckily, a computer program is just about the best card counter out there, so we do know two things:

1. We know exactly how many cards each player currently has. Every player gets to see the cards dealt, and the cards played from the hand are public knowledge.
2. We know exactly how many of each rank of card are left. The deck that we are playing with is public knowledge, and—once again—the cards that players play are public knowledge.

The Simple Solution

The natural starting point for solving this problem is simply dealing the cards, while respecting the constraints imposed by the margins. I’ll call this algorithm the “card-dealing” approach.

While the most naive version of this approach is to create a list of cards, and then shuffle and deal them, creating and shuffling a list is quite slow and unnecessary. We waste time shuffling items that are for our purposes identical.

Instead, the process will be something like this: The number of cards we have to deal corresponds to the total mass. Each step of the algorithm, we simply select a row and a column which have remaining unallocated mass, and then add one to the entry at that row and column. Then, we update the remaining unallocated mass for the row and column, so that we can skip selecting columns and rows who are already “full”.

pub fn card_dealing<const M: usize, const N: usize>(
    row_margins: [usize; M],
    column_margins: [usize; N],
    rng: &mut SmallRng,
) -> [[usize; N]; M] {
    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;
}

Stars-and-Bars

Before we move on to the next algorithm, let’s take a look at a simpler problem, and the solution to that problem.

The related problem is often formulated as balls in urns, or cookies to children, but the idea is the same. Distribute some fixed number of items into a fixed number of bins, uniformly.

This is exactly the 1-dimensional case of our possible worlds problem. And luckily for us, there is a great algorithm for doing this uniformly.

Greedy Stars-and-Bars

So how do we use the stars-and-bars approach for our 2d problems?

While I did try in vain to figure out a way that I could replicate the simplicity of stars-and-bars, I wasn’t able to. The ideal would be if we could just sample $N\times M$ total numbers, and then sort them in a way that would preserve all of our margins.

While not the only problem, one of the largest problems that an approach like this faces is that not every selection of $N\times M$ integers which are less than the total mass forms a valid solution to the possible world problem.

Instead of trying to construct some impossible miracle solution, let’s try something a bit simpler, and more reasonable.

The general approach of the “Greedy Stars-and-Bars” algorithm is to fill in our matrix one row or column at a time. If we do that, then we can just use normal Stars-and-Bars to sample the values for that row or column.

The only thing we need to ensure is that we do not violate the constraints of the rows or columns that we are “crossing” (that are perpendicular to the row or column that we are currently filling).

All we need to do is to select the row or column which has the smallest remaining unassigned mass. No matter what values come out of the Stars-and-Bars, they cannot exceed the perpendicular constraints, because the constraint that we are currently honoring is the smallest constraint “on the board”, so to speak.

Let’s do a quick example. Consider the following fixed-margin problem.

The currently most-constrained row or column is $r_0$, with $3$ remaining unassigned mass. We will start by performing stars-and-bars across that row. The sampler gives us $[0, 1, 2]$ as the row.

Notice that we couldn’t violate any constraints, because all the values sampled would be $\leq 3$. The next most constrained row or column is $c_2$. Here, we will sample only for the entries which are not already assigned, so we only need to distribute the $3$ remaining mass among two bins.

You simply continue that process until the entire matrix is filled. It will always generate a valid solution to the fixed-margin problem.

What is nice about this algorithm is that it does not depend on the total mass. The card-dealing algorithm, for example, is $\mathcal{O}(T)$, where $T$ is the total mass of the problem.

The greedy stars-and-bars algorithm is $\mathcal{O}(N + M)$, instead. It scales with the size of the matrix.

fn stars_and_bars<const N: usize>(
    mass: usize,
    slots: usize,
    output: &mut [usize; N],
    rng: &mut SmallRng,
) {
    debug_assert!(slots <= N);
    debug_assert!(slots > 0);

    if slots == 1 {
        output[0] = mass;
    }

    let num_bars = slots - 1;
    let mut bar_positions = [0usize; N];

    for i in 0..num_bars {
        bar_positions[i] = rng.random_range(0..=mass);
    }

    bar_positions[..num_bars].sort_unstable();

    let mut previous_bar = 0;

    for idx in 0..num_bars {
        let bar = bar_positions[idx];
        output[idx] = bar - previous_bar;
        previous_bar = bar;
    }

    output[num_bars] = mass - previous_bar;

    debug_assert_eq!(output[..slots].iter().sum::<usize>(), mass);
}

pub fn greedy_stars_and_bars<const M: usize, const N: usize>(
    row_margins: [usize; M],
    column_margins: [usize; N],
    rng: &mut SmallRng,
) -> [[usize; N]; M] {
    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);

    let mut row_allocations = [0; N];
    let mut col_allocations = [0; M];

    active_rows.extend((0..M).filter(|&row| row_margins[row] > 0));
    active_cols.extend((0..N).filter(|&col| column_margins[col] > 0));

    while !active_rows.is_empty() && !active_cols.is_empty() {
        // First, find our candidates
        let mut best_row_index = 0;
        let mut best_row_index_location = 0;
        let mut best_row_value = usize::MAX;
        for (row_location, row_index) in active_rows.iter().enumerate() {
            if remaining_row_mass[*row_index] < best_row_value {
                best_row_index = *row_index;
                best_row_value = remaining_row_mass[*row_index];
                best_row_index_location = row_location;
            }
        }

        let mut best_col_index = 0;
        let mut best_col_index_location = 0;
        let mut best_col_value = usize::MAX;
        for (col_location, col_index) in active_cols.iter().enumerate() {
            if remaining_col_mass[*col_index] < best_col_value {
                best_col_index = *col_index;
                best_col_value = remaining_col_mass[*col_index];
                best_col_index_location = col_location;
            }
        }

        if best_row_value < best_col_value {
            stars_and_bars(best_row_value, active_cols.len(), &mut row_allocations, rng);
            for (col, allocation) in active_cols.iter().zip(row_allocations) {
                output[best_row_index][*col] = allocation;
                remaining_col_mass[*col] -= allocation;
            }
            active_rows.swap_remove(best_row_index_location);
        } else {
            stars_and_bars(best_col_value, active_rows.len(), &mut col_allocations, rng);
            for (row, allocation) in active_rows.iter().zip(col_allocations) {
                output[*row][best_col_index] = allocation;
                remaining_row_mass[*row] -= allocation;
            }
            active_cols.swap_remove(best_col_index_location);
        }
    }

    return output;
}

Next Steps

Our whole goal was to create a uniform sampler over all possible worlds, and neither of these algorithms is even that!

Seems like if we want an unbiased sampler, I am going have to get ~creative~…

In the next post in this series, I will explore some more ideas about how we can build a truly uniform sampler, and explain why it’s a harder problem than you think (The scaling laws of combinatorics always seem to be getting in our way).