Recently my wife, who’s studying medicine, had to partake in what’s called an “internship lottery”. The basic idea is that there are 25 hospitals, each with a specific amount of vacancies, and students need to be assigned to internships.

Granted, some hospitals are more popular than others, so some students are bound to get disappointed. To keep things fair, there’s a lottery. Each student ranks all the hospitals, and then there’s an algorithm that takes into account the students’ preferences and the hospitals’ capacity, and assigns each student to a hospital.

Up till ~2015, the algorithm was an RSD (Random Serial Dictatorship, more on that soon), which eventually meant that allocations had a pretty random character. Around 2015 however, there has been an “upgrade” to the algorithm which was supposed to make students more satisfied (as a whole) with their assignments 1. Needless to say, many students were not satisfied with their assignments, which lead to a kind of folklore of strategies to maximize one’s odds of getting their desired assignment(s).

The researchers’ motivation for the new algorithm was to improve the overall satisfaction of students. And before going into details, this changed the aforementioned strategies. And this leads me to my motivations for looking into this thing.

First, as I mentioned, my wife, she’s going to participate in this lottery, and we wanted to test different strategies and measure their outcome. Secondary is my curiosity of figuring out how this thing works.

Before talking about the algorithm, I need to set up some notations.

There are $M$ hospitals, each has a capacity of $C_i$.

There are $N$ students, each has its own ranking of how they prefer the hospitals. Student $i$ ranks the hospitals in a list $R_i = \left(r_{i,1},…,r_{i,M}\right)$.

Back to the new algorithms. It’s actually composed of 3 steps:

  1. Run RSD a large amount of times. Then collect the results from all the runs to calculate the probability a student $i$ gets hospital $j$ ($p_{i,j}$).
  2. Successively trade probabilities between students to maximize a goal “happiness” (which depends on their preferences and calculated probabilities).
  3. Make assignments based on the improved probabilities. A sort of weighted lottery, if you will.

This is all very well, but to actually test the implementation, we need input data. This means we need the preferences of all the students. That’s 270 lists of 25 items. This information is not publicly available.

Generating the initial data set

What is available is a summary of preferences. It’s a table (actually a matrix $Q \in \mathbb{N}^{M \times M}$) that summarizes the amount of students which ranked hospital $j$ as their $i$’th choice.

What can we do with that? Well, remember the preference list $R_i$? It’s basically a vector, and it can be rewritten like this: $R_i = \left(r_{i,1},…,r_{i,M}\right) = \left(\rho_i(1), …, \rho_i(M)\right)$ where $\rho_i: [1,M] \rightarrow [1,M]$ is a permutation function that maps a hospital to a preference. This means we can write the preference vector as $R_i = \left[\rho_i\right] \cdot \left(1, …, M\right)^t = \left[\rho_i\right] \cdot H$.

If the student ranked hospital j at preference i, then the matrix $\left[\rho_i\right]$ will have a 1 in row/column i/j, and 0 at the rest of the row i and column i.

Which means, that if we sum all the permutation matrices of all students, and we look at the resulting matrix at row/colum i/i, we will see how many students made that choice!

\[Q = \sum_{i=1}^{N} \left[\rho_i\right]\]

Why is this insight important? Because if we normalize $Q$ by $N$, we get what’s called a Doubly Stochastic Matrix, which according to the Birkhoff–von Neumann theorem can be decomposed in the following manner:

\[\frac{1}{N} Q = \tilde Q = \sum_{i} \theta_i P_i\]

Where $P_i$ are permutation matrices and $0 \leq \theta_i \leq 1$.

We can apply the same theorem to decompose $Q$ directly:

\[Q = \sum_{i} n_i P_i\]

Where $P_i$ are (still) permutation matrices and $n_i$ will be the number of students that have that specific permutation. Of course $\sum n_i = N$.

The algorithm used to perform such a decomposition is the Birkhoff algorithm. Time to go off on a tangent about the Birkhoff algorithm.

Birkhoff algorithm

Starting with $Q_0 = Q$ the outline of the algorithm is as follows (for step i):

  1. Construct the positivity matrix $Q^{+}_i$ (has a 1 everywhere $Q_i$ is positive)
  2. Convert $Q^{+}_i$ to a bipartate graph $G_i$ with the rows being the left vertices, and the columns being the right ones.

  3. Find a perfect matching of $G_i$, which is a set of edges $\hat{G}_i$.

    A property of a matching is that no two edges have a common vertex.

    This corresponds to the equivalent row/column in the matrix containing only one 1. Combined with the matching being perfect (covering all vertices, i.e. having $M$ edges), corresponds to a matrix that’s doubly-stochastic.

    There are two questions here:

  4. Convert the matching back to a matrix $P_i$

  5. Find the smallest element of $Q_i$ when masked by $P_i$:

    \[n_i = \min_{\left(r,c\right) \in \mathrm{Matching}} \left[Q_i\right]_{r,c}\]
  6. We can now emit the $P_i$ matrix $n_i$ times

  7. $Q_{i+1} = Q_i - n_i P_i$

    When $Q_{i+1} = O$ we’re done.

Interactive example

In this example, we have the allocation matrix of 15 students competing for 4 hospitals (denoted by 4 colors).

Click next to cycle through the iteration of the algorithm and watch the middle to see the individual preferences as they are extracted (each preference line will consist of colored squares matching the hospitals):

Implementation

This wouldn’t really be a programming blog without code, so there (using numpy and networkx):

from typing import Iterator, Tuple

import numpy as np
import networkx as nx
from networkx import from_numpy_matrix
from networkx.algorithms.bipartite.matching import maximum_matching

def birkhoff(Q: np.matrix) -> Iterator[Tuple[int, np.matrix]]:
    M, _ = Q.shape

    while not np.all(Q == 0):
        # 1. Construct the positivity matrix of Q (Qp = Q-positive)
        Qp = np.zeroes_like(Q)
        Qp[Q.nonzero()] = 1

        # 2. Convert the positivity matrix to bipartate graph
        Z = np.zeros_like(G, dtype=int)
        # Exapnd to a bipartate matrix (rows and columns map to different nodes)
        B = np.block([[Z, G, G.T, Z])
        G = nx.from_numpy_matrix(B)

        # 3. Find the perfect matching edges (Gh = G-hat)
        Gh = maximum_matching(G, range(M))

        # 4. Convert the matching edges back to a matrix
        # We need to "fold" the vertices of the matching edges, as columns
        # are mapped to the range [M, 2M - 1]
        edges = set((x % M, y % M) for x, y in Gh.items())
        # Create a matrix from the edges
        P = np.zeroes_like(Q)
        P[list(edges)] = 1

        # 5. Find the coefficient
        n = min(Q[edges])

        # 6. Emit result
        yield n, P

        # 7. Prepare next iteration
        Q -= n * P

Now that we have a dataset to work on, we can continue with the implementation of the assignment algorithm. Stay tuned.

References

  1. Bronfman, S., Hassidim, A., Afek, A. et al. Assigning Israeli medical graduates to internships. Isr J Health Policy Res 4, 6 (2015). https://doi.org/10.1186/2045-4015-4-6 ↩