m+e+c case: matching with couples

Alfred Galichon (NYU & Sciences Po) and Antoine Jacquet (Sciences Po)

'math+econ+code' masterclass series

With python code examples

© 2018–2024 by Alfred Galichon. Research assistance from Alessandro Facchini is gratefully acknowledged. Past and present support from NSF grant DMS-1716489, ERC grant CoG-866274 are acknowledged, as well as inputs from contributors listed here.

If you reuse material from this masterclass, please cite as:
Alfred Galichon and Antoine Jacquet, 'Matching with couples', 'math+econ+code' masterclass series. https://www.math-econ-code.org/

References

• Scarf (1967). "The core of an N person game." Econometrica.
• Biró, Fleiner and Irving (2016). "Matching couples with Scarf's algorithm." Annals of Mathematics and Artificial Intelligence.
• Nguyen and Vohra (2018). "Near-feasible stable matchings with couples." American Economic Review.
• Galichon and Jacquet (2024). math+econ+code lecture on Scarf's ordinal basis algorithm https://www.math-econ-code.org/scarf-ordinal-basis.

Prerequisites

This lecture mainly builds upon the math+econ+code lecture on Scarf's ordinal basis algorithm https://www.math-econ-code.org/scarf-ordinal-basis. To import:

#!pip install mec --upgrade

from mec.gt import OrdinalBasis
import numpy as np

The setting

The following is based on Ngyen and Vohra's paper https://www.aeaweb.org/articles?id=10.1257/aer.20141188.

Primitives of the problem:

• set of single doctors $S$
• set of couples $C$ $\quad \implies \quad$ create set of doctors

$D = S \cup \big\{ d : (d,\cdot) \in C \text{ or } (\cdot, d) \in C \big\}$, we have $|D| = |S| + 2|C|$

• set of hospitals $H$
• set of regions $R$ (less regions than hospitals)

We build based on this:

• an assignment of hospitals to regions [random]
• the incidence matrix $M$ (matrix $A$ in Nguyen–Vohra and Scarf) [not random]
• the preference matrix $\Phi$ (matrix $C$ in Nguyen–Vohra and Scarf) [random]
• the capacity vector $q$ (vector $b$ in Nguyen–Vohra and Scarf) [random for hospitals]

We will use as illustration:

numSingles = 1  # 9
numCouples = 1
numHosps = 2
numRegions = 1

np.random.seed(7777)    # 0 by default

print('Total number of doctors:', numSingles + 2*numCouples)

Total number of doctors: 3

Regions

Hospitals are assigned randomly to regions, making sure that there is at least one hospital in each region.

def assign_regions(numHosps, numRegions):
    hospitalRegions = np.concatenate((np.arange(numRegions), # at least one hospital per region
                                      np.random.randint(0, numRegions, numHosps - numRegions))) 
    return np.random.permutation(hospitalRegions)

hospitalRegions = assign_regions(numHosps, numRegions)
print(hospitalRegions)

[0 0]

Incidence matrix $M$

The incidence matrix has one row per agent:

• first $|S|$ rows: single doctors,
• next $|C|$ rows: couples,
• last $|H|$ rows: hospitals.

Its columns correspond to all possible arrangements. The left side of the matrix is the autarky arrangements: it is just the identity matrix. The rest can be written by blocks: \begin{equation} \begin{pmatrix} \text I_S & & & M_{11} & \\ & \text I_C & & & M_{22} \\ & & \text I_H & M_{31} & M_{32} \\ \end{pmatrix}. \end{equation}

The blocks $M_{11}$ and $M_{31}$ correspond to arrangements involving a single doctor and a hospital ($|S| \times |H|$ columns).
The blocks $M_{22}$ and $M_{32}$ correspond to arrangements involving a couple and a hospital pair of the form

\begin{equation} (h, h') \in (H \cup \{0\})^2 - \{(0,0)\} \overset{\text{def}}{=} P. \end{equation}

Such pairs may take the form $(h,0)$ or $(0,h)$ if one of the partners remains unassigned, or $(h, h)$ if both partners are assigned to the same hospital.
The number of such hospital pairs is $|P| = (|H|+1)^2 - 1$. Thus the blocks $M_{22}$ and $M_{32}$ correspond to the last $|C| \times |P|$ columns.

Let's write the right side of the matrix $M$ in extended form, just leaving the lower-right block $M_{32}$ as is for now:

\begin{equation} \left[ \quad \begin{matrix} 1 & \cdots & 1 & & & & & & & \\ & & & 1 & \cdots & 1 & & & & \\ & & & & & & \ddots & & & \\ & & & & & & & 1 & \cdots & 1 \\ & & & & & & & & & \\ & & & & & & & & & \\ & & & & & & & & & \\ & & & & & & & & & \\ 1 & & & 1 & & & & 1 & & \\ & \ddots & & & \ddots & & \cdots & & \ddots & \\ & & 1 & & & 1 & & & & 1 \\ \end{matrix} \qquad \begin{matrix} & & & & & & & & & \\ & & & & & & & & & \\ & & & & & & & & & \\ & & & & & & & & & \\ 1 & \cdots & 1 & & & & & & & \\ & & & 1 & \cdots & 1 & & & & \\ & & & & & & \ddots & & & \\ & & & & & & & 1 & \cdots & 1 \\ & & & & & & & & & \\ & & & & M_{32} & & & & & \\ & & & & & & & & & \\ \end{matrix} \quad \right]. \end{equation}

Arrangements with a single doctor. The left-hand blocks are the usual margining-out matrices. Using the Kronecker product $\otimes$,

\begin{equation} M_{11} = \text I_S \otimes 1_H^\top, \qquad M_{31} = 1_S^\top \otimes \text I_H. \end{equation}

Arrangements with a couple. On the right-hand side, the matrix $M_{22}$ is also a standard margining-out matrix: $M_{22} = \text I_C \otimes 1_P^\top$.

The block $M_{32}$ is a little trickier. For a given couple we can write two incidence matrices: one for partner 1, the other for partner 2.

We consider the (extended) list of hospital pairs: $(h_1,h_1), (h_1,h_2), \dots, (h_1,h_H), (h_1,0), (h_2,h_1), (h_2, h_2), \dots, (0,h_1), \dots, (0,h_H), (0,0)$.

In this list, partner 1 is assigned to a hospital according to

\begin{equation} \begin{bmatrix} 1 & 1 & \cdots & 1 & & & & & & & & & \\ & & & & 1 & 1 & \cdots & 1 & & & & & \\ & & & & & & & & \ddots & & & & \\ & & & & & & & & & 1 & 1 & \cdots & 1 \\ \end{bmatrix} \end{equation}

while partner 2 is assigned according to

\begin{equation} \begin{bmatrix} 1 & & & & 1 & & & & & 1 & & & \\ & 1 & & & & 1 & & & & & 1 & & \\ & & \ddots & & & & \ddots & & \cdots & & & \ddots & \\ & & & 1 & & & & 1 & & & & & 1 \\ \end{bmatrix}. \end{equation}

These are also standard margining-out matrices. We sum them up to obtain the incidence matrix of the couple:

\begin{equation} \begin{bmatrix} 2 & 1 & \cdots & 1 & 1 & & & & & 1 & & & \\ & 1 & & & 1 & 2 & \cdots & 1 & & & 1 & & \\ & & \ddots & & & & \ddots & & \cdots & & & \ddots & \\ & & & 1 & & & & 1 & & 1 & 1 & \cdots & 2 \\ \end{bmatrix}. \end{equation}

We still need to remove the last row (which corresponds to the empty spot 0) and the last column (which corresponds to both partners being unassigned).

We then obtain $M_{32}$ by repeating this matrix horizontally, once per each couple.

We generate the matrix $M$ using:

def create_M(numSingles, numCouples, numHosps):
    numHospPairs = (numHosps+1)**2 - 1
    M11 = np.kron(np.eye(numSingles), np.ones((1, numHosps)))
    M21 = np.zeros((numCouples, numSingles * numHosps))
    M31 = np.kron(np.ones((1, numSingles)), np.eye(numHosps))
    M12 = np.zeros((numSingles, numCouples * numHospPairs))
    M22 = np.kron(np.eye(numCouples), np.ones((1, numHospPairs)))
    M32_unit_m = np.kron(np.eye(numHosps+1), np.ones((1, numHosps+1)))[:-1,:-1]
    M32_unit_f = np.kron(np.ones((1, numHosps+1)), np.eye(numHosps+1))[:-1,:-1]
    M32 = np.kron(np.ones((1, numCouples)), M32_unit_m + M32_unit_f) 
    M_right = np.block([[M11, M12], [M21, M22], [M31, M32]])
    return np.hstack((np.eye(numSingles+numCouples+numHosps), M_right))

create_M(numSingles, numCouples, numHosps)[:,numSingles+numCouples+numHosps:]

array([[1., 1., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 1., 1., 1., 1., 1., 1., 1., 1.],
       [1., 0., 2., 1., 1., 1., 0., 0., 1., 0.],
       [0., 1., 0., 1., 0., 1., 2., 1., 0., 1.]])

Preference matrix $\Phi$

The preference matrix $\Phi$ will have a very similar shape to the incidence matrix $M$. We can write it by blocks as

\begin{equation} \begin{pmatrix} 0 & & & \Phi_{11} & \\ & \ddots & & & \Phi_{22} \\ & & 0 & \Phi_{31} & \Phi_{32} \\ \end{pmatrix} \end{equation}

where this time the spots left blank will be taken by the large number $K$.

Preferences of single doctors over hospitals.
The popularity of each hospital $h$ is determined according to the model

$p_h = .99 \, |D| \, (0.8)^{k_h} + .01 \, |H| \quad$ where $k_h$ is a random integer between 1 and $|H|$.

To generate a ranking of hospitals for each single doctor, we draw hospitals $h$ one by one without replacement with probability proportional to $p_h$.
We then transform this ranking into a vector of utilities ranging from 1 (the lowest-ranked hospital) to $|H|$ (the highest-ranked).
Stacking these vectors in a pattern similar to the margining-out-matrix $M_{11}$, we obtain the matrix $\Phi_{11}$.

Preferences of couples over hospitals.
First, couples are 'unemployment-averse': they always prefer an arrangement in which both of them have a job, rather than only one of them.
Thus the $2|H|$ hospital pairs of the form $(h,0)$ and $(0,h)$ are at the bottom of the preference ordering for all couples: we order them at random (uniformly) with utilities from 1 to $2|H|$.

Second, couples prefer to be assigned to hospitals in the same region. We generate popularity scores $\nu_{hh'}$ for hospital pairs $(h, h') \in H^2$ with

$\nu_{hh'} = \begin{cases} .7 \, p_h p_{h'} & \text{if $h$ and $h'$ are in the same region,} \\ .3 \, p_h p_{h'} & \text{otherwise.} \\ \end{cases}$

Then, we draw hospital pairs $hh'$ without replacement with probability proportional to $\nu_{hh'}$ to build a couple's ranking, and transform this ranking into utilities ranging from $2|H|+1$ to $2|H|+|H|^2$.

In this way, a couple gets a utility vector over all hospital pairs in $P$ (including those with an unemployed partner).
We stack these vectors of utilities in a pattern similar to $M_{22}$ to form the matrix $\Phi_{22}$.

Preferences of hospitals over doctors.
All hospitals have the same strict preference ordering over doctors.
This is enough to assign a value to single doctors, but not to couples.
The value of a couple is then set as the value of the lower-ranked partner in that couple.
This creates a number of ties: for instance, when partner 1 is the lower-ranked in couple $c$, then hospital $h$ is indifferent between all arrangements $(c, (h, h'))$ for $h' \in H \cup \{0\}$.
In this case, these ties are broken using the couples preferences: $h$ prefers $(c, (h, h'))$ over $(c, (h, h''))$ if and only if couple $c$ does as well.

Similarly, when partner 2 is the higher-ranked in couple $c$, then hospital $h$ is indifferent between all arrangements $(c, (h', h))$ for $h' \neq h$.
Ties are resolved in the same way: $h$ prefers $(c, (h', h))$ over $(c, (h'', h))$ if and only if couple $c$ does as well.

We generate the matrix $\Phi$ with:

def create_Φ(numSingles, numCouples, numHosps, hospitalRegions, K=np.iinfo(np.int32).max):
    numApplicants = numSingles + 2*numCouples
    numHospPairs = (numHosps+1)**2 - 1
    
    popularity_h = .99*numApplicants*(.8)**np.random.randint(1, numHosps+1, numHosps) + (1-.99)*numHosps
    same_region_h_k = np.equal.outer(hospitalRegions, hospitalRegions).flatten()
    score_h_k = .3 * np.kron(popularity_h, popularity_h) + .4 * same_region_h_k * np.kron(popularity_h, popularity_h)
    value_doctors = np.random.permutation(numApplicants) + 1
    
    Φ12 = np.full((numSingles, numCouples * numHospPairs), K)
    Φ21 = np.full((numCouples, numSingles * numHosps), K)
    
    # Φ11: preferences of single doctors over hospitals
    Φ11 = np.full((numSingles, numSingles * numHosps), K)
    for d in range(numSingles):
        hospital_ranking = np.random.choice(numHosps, numHosps, replace=False, p=popularity_h/popularity_h.sum())
        Φ11[d, numHosps*d:numHosps*(d+1)] = numHosps - np.argsort(hospital_ranking)
        
    # Φ22: preferences of couples over hospitals pairs
    Φ22 = np.full((numCouples, numCouples * numHospPairs), K)
    indices_pairswith0 = list(range(numHosps, numHosps*(numHosps+1), numHosps+1)) + list(range(numHospPairs - numHosps, numHospPairs))
    indices_pairswithout0 = [i for i in range(numHospPairs) if i not in indices_pairswith0]
    for c in range(numCouples):
        utilities = np.ones(numHospPairs)
        utilities[indices_pairswith0] = np.random.permutation(2*numHosps)+1
        pair_ranking = np.random.choice(numHosps**2, numHosps**2, replace=False, p=score_h_k/score_h_k.sum())
        utilities[indices_pairswithout0] = numHospPairs - np.argsort(pair_ranking)
        Φ22[c, numHospPairs*c:numHospPairs*(c+1)] = utilities
        
    # Φ31: preferences of hospitals over single doctors
    Φ31 = np.kron(value_doctors[:numSingles], np.eye(numHosps))
    
    # Φ32: preferences of hospitals over couples
    Φ32 = np.full((numHosps, numCouples * numHospPairs), K, dtype=float) # float to break ties
    value_partner1 = value_doctors[numSingles:(numSingles+numCouples)]
    value_partner2 = value_doctors[-numCouples:]
    for c in range(numCouples):
        Φ32_unit_m = value_partner1[c] * np.kron(np.eye(numHosps+1), np.ones((1, numHosps+1)))[:-1,:-1]
        Φ32_unit_f = value_partner2[c] * np.kron(np.ones((1, numHosps+1)), np.eye(numHosps+1))[:-1,:-1]
        for h in range(numHosps):   # breaking ties
            columns_tied_m = numHospPairs*c + (numHosps+1)*h + np.arange(numHosps+1)
            Φ32_unit_m[h,(numHosps+1)*h:(numHosps+1)*(h+1)] += Φ22[c,columns_tied_m]/(Φ22[c,columns_tied_m].sum()+1)
            columns_tied_f = numHospPairs*c + h + (numHosps+1)*np.arange(numHosps+1)
            Φ32_unit_f[h,h+(numHosps+1)*np.arange(numHosps+1)] += Φ22[c,columns_tied_f]/(Φ22[c,columns_tied_f].sum()+1)
        Φ32_unit_m[Φ32_unit_m==0] = K
        Φ32_unit_f[Φ32_unit_f==0] = K
        Φ32_unit = np.minimum(Φ32_unit_m, Φ32_unit_f)
        Φ32[:,numHospPairs*c:numHospPairs*(c+1)] = Φ32_unit
    
    Φ_right = np.block([[Φ11, Φ12], [Φ21, Φ22], [Φ31, Φ32]])
    Φ_right[Φ_right == 0] = K
    Φ_left = np.full((numSingles+numCouples+numHosps, numSingles+numCouples+numHosps), K)
    np.fill_diagonal(Φ_left, 0)
    Φ = np.hstack((Φ_left, Φ_right))
    return Φ

create_Φ(numSingles, numCouples, numHosps, hospitalRegions, K=100)

array([[  0.        , 100.        , 100.        , 100.        ,
          1.        ,   2.        , 100.        , 100.        ,
        100.        , 100.        , 100.        , 100.        ,
        100.        , 100.        ],
       [100.        ,   0.        , 100.        , 100.        ,
        100.        , 100.        ,   7.        ,   8.        ,
          3.        ,   5.        ,   6.        ,   2.        ,
          1.        ,   4.        ],
       [100.        , 100.        ,   0.        , 100.        ,
          2.        , 100.        ,   1.5       ,   3.42105263,
          3.15789474,   1.35714286, 100.        , 100.        ,
          1.07142857, 100.        ],
       [100.        , 100.        , 100.        ,   0.        ,
        100.        ,   2.        , 100.        ,   1.42105263,
        100.        ,   3.35714286,   1.31578947,   3.14285714,
        100.        ,   1.21052632]])

Capacity vector $q$

For single doctors and couples, the capacity is of course 1.

For hospitals, we randomly generate seats to be allocated to hospitals, but making sure that (i) each hospital gets at least one seat, and (ii) the total number of seats is equal to the number of applicants.

To do this, we give one seat to each hospital (so $|H|$ seats in total), and we create $|S| + 2|C| - |H|$ extra seats to be allocated randomly.
Each of these extra seats is randomly allocated to a hospital, in proportion to the popularity of that hospital's region (popularity is itself chosen as the inverse of the hospital region's index, according to Zipf's law).

[Actually a small mistake here: seats should first be allocated to regions according to Zipf's law, and then seats within a region should be allocated uniformly to hospitals in that region. This is actually what the paper describes but the code doesn't seem to do that.]

We code this as follows:

def create_q(numSingles, numCouples, numHosps, hospitalRegions):
    numApplicants = numSingles + 2*numCouples
    popularity = 1/(hospitalRegions+1)    # Zipf's law
    extra_seats = np.random.choice(numHosps, numApplicants - numHosps, p=popularity/popularity.sum())
    q_h = np.bincount(extra_seats, minlength=numHosps) + 1    # at least 1 seat per hospital
    q = np.concatenate((np.ones(numSingles+numCouples), q_h))
    return q

print('Number of seats for each hospital:')
create_q(numSingles, numCouples, numHosps, hospitalRegions)[-numHosps:]

Number of seats for each hospital:

array([1., 2.])

We are now ready to create the ordinal bases problem with:

def create_nguyen_vohra(numSingles=270, numCouples=20, numHosps=18, numRegions=5): # default values from paper
    K = np.iinfo(np.int32).max
    hospitalRegions = assign_regions(numHosps, numRegions) 
    M_z_a = create_M(numSingles, numCouples, numHosps)
    Φ_z_a = create_Φ(numSingles, numCouples, numHosps, hospitalRegions, K)
    q_z = create_q(numSingles, numCouples, numHosps, hospitalRegions)
    print('In total', q_z[-numHosps:].sum(), 'seats.')
    return (Φ_z_a, M_z_a, q_z, K, 1e-5)


numSingles = 20
numCouples = 80
numHosps = 18
numRegions = 5
numHospPairs = (numHosps+1)**2 - 1

nguyen_vohra = OrdinalBasis(*create_nguyen_vohra(numSingles, numCouples, numHosps, numRegions))

In total 180.0 seats.

sol = nguyen_vohra.solve(verbose=1)
μ_a = sol['μ_a']

Let's first check if our solution $\mu$ is integral:

def is_integral(μ_a):
    int_a = μ_a == np.floor(μ_a)
    return np.all(int_a), np.where(~int_a)[0].tolist()

print(is_integral(μ_a))

(True, [])

If $\mu$ is not integral, we need a rounding procedure. Nguyen and Vohra suggest a rounding algorithm.

For a given solution $\mu$, there must be one hospital constraint

$\sum_{d \in D} \mu_{d, h} + \sum_{c \in C, h' \in H} (\mu_{c, h, h'} + \mu_{c, h', h}) \leq q_h$

which involves at most 2 fractional components of $\mu$. If not, it would mean all hospital constraints involve at least 3 fractional components of $\mu$, in which case...

We will sequentially solve the linear program:

\begin{align} \max_{\mu \geq 0} ~ & \sum_{d, h} \mu_{d,h} + \sum_{c,h,h'} 2 \mu_{c,h,h'} \\ \text{s.t.} ~& \mu_a = \mu_a^{t-1} \quad \text{if $\mu_a^{t-1}$ is 0 or 1} \\ & \textstyle\sum_{h \in H} \mu_{d, h} = 1 \quad \text{if that constraint binds with $\mu^{t-1}$} \\ & \textstyle\sum_{h \in H} \mu_{d, h} \leq 1 \quad \text{if that constraint does not bind with $\mu^{t-1}$} \\ & \textstyle\sum_{h, h' \in H} \mu_{c, h, h'} = 1 \quad \text{if that constraint binds with $\mu^{t-1}$} \\ & \textstyle\sum_{h, h' \in H} \mu_{c, h, h'} \leq 1 \quad \text{if that constraint does not bind with $\mu^{t-1}$} \\ & \textstyle\sum_{d \in D} \mu_{d, h} + \sum_{c \in C, h' \neq h} (\mu_{c, h, h'} + \mu_{c, h', h}) + \sum_{c \in C} 2\mu_{c, h, h} \leq q_h \quad \forall h \in H \\ & \textstyle\sum_{d, h} \mu_{d,h} + \sum_{c,h,h'} 2 \mu_{c,h,h'} \leq \sum_h q_h \end{align}

We impose the aggregate constraint

$\sum_{d, h} \mu_{d,h} + \sum_{c,h,h'} 2 \mu_{c,h,h'} \leq \sum_h q_h$.

To be completed.