Introduction

matchingR is an R package which quickly computes the Gale-Shapley algorithm (Gale and Shapley 1962), Irving’s algorithm for the stable roommate problem (Irving 1985), and the top trading cycle algorithm (Shapley and Scarf 1973) for large matching markets. The package provides functions to compute the solutions to the stable marriage problem, the college admission problem, the stable roommates problem, and the house allocation problem.

The package may be useful when the number of market participants is large or when many matchings need to be computed (e.g., for simulation or estimation purposes). It has been used in practice to compute the Gale-Shapley stable matching with 30,000 participants on each side of the market.

Matching markets are common in practice and widely studied by economists. Popular examples include

• the National Resident Matching Program that matches graduates from medical school to residency programs at teaching hospitals throughout the United States
• the matching of students to schools including the New York City High School Match or the the Boston Public School Match (and many more)
• the matching of kidney donors to recipients in kidney exchanges.

Two-sided Matching Markets: Gale-Shapley Algorithm

Consider the following marriage market: There are N men and N women. Each man, m, receives utility uM(w, m) from a match with woman w, and similarly each woman receives a payoff of uW(m, w) from being matched with a man.

A matching assigns men to women such that each man is assigned to one woman and each woman is assigned to one man. A matching is stable if there is no man and woman who would jointly prefer to be matched to each other over their current spouses. In other words, a matching is stable if there are no pairs (m, w'), (m', w), such that m is matched with w', m' is matched with w, and both uW(m, w) > uW(m', w) and uM(m, w) > uM(m, w').

For example, we might have preferences for men given by

uM = matrix(c(1.0, 0.5, 0.0,
              0.5, 0.0, 0.5,
              0.0, 1.0, 1.0), nrow = 3, ncol = 3, byrow = TRUE)

##          cols
## rows      Man 1 Man 2 Man 3
##   Woman 1   1.0   0.5   0.0
##   Woman 2   0.5   0.0   0.5
##   Woman 3   0.0   1.0   1.0

In this example, man 1 receives a payoff of 1.0 from being matched to woman 1, a payoff of 0.5 from being matched to woman 2 and a payoff of 0.0 from being matched to woman 3 (same logic applies to men 2 and 3). Similarly, we might have preferences for women given by

uW = matrix(c(0.0, 1.0, 0.0,
              0.5, 0.0, 0.5,
              1.0, 0.5, 1.0), nrow = 3, ncol = 3, byrow = TRUE)

##        cols
## rows    Woman 1 Woman 2 Woman 3
##   Man 1     0.0     1.0     0.0
##   Man 2     0.5     0.0     0.5
##   Man 3     1.0     0.5     1.0

Here, columns in the matrix correspond to women, rows to men. In this example, woman 1 receives a payoff of 0.0 from being matched to man 1, a payoff of 0.5 from being matched to man 2 and a payoff of 1.0 from being matched to man 3 (same logic applies to women 2 and 3).

Instead of using payoff matrices, we can also represent preferences using preference orderings. The preference ordering that corresponds to uM is

prefM = matrix(c(1, 3, 3,
                 2, 1, 2,
                 3, 2, 1), nrow = 3, ncol = 3, byrow = TRUE)

##         cols
## rows     Man 1 Man 2 Man 3
##   Rank 1     1     3     3
##   Rank 2     2     1     2
##   Rank 3     3     2     1

prefM states that man 1 prefers woman 1 over woman 2 over woman 3, etc. The preference ordering that corresponds to uW is given by

prefW = matrix(c(3, 1, 3,
                 2, 3, 2,
                 1, 2, 1), nrow = 3, ncol = 3, byrow = TRUE)

##         cols
## rows     Woman 1 Woman 2 Woman 3
##   Rank 1       3       1       3
##   Rank 2       2       3       2
##   Rank 3       1       2       1

The matching algorithm discussed below can take either payoff matrices of the type uM and uW or preference orderings of the type prefM and prefW as arguments.

The Gale-Shapley algorithm works as follows: Single men (“the proposers”) sequentially make proposals to each of their most preferred available women (“the reviewers”). A woman can hold on to at most one proposal at a time. A single woman will accept any proposal that is made to her. A woman who already has a proposal will reject any proposal she values less than her current proposal in hand. If a woman receives a proposal from a man that she values more than her current proposal, she will accept the proposal and her previous match will rejoin the line of proposers. This process continues until all men are matched to all women.

For the preferences specified in uM and uW, we can compute the Gale-Shapley Algorithm by hand. Initially, all men are single.

1. • Man 1 proposes to woman 1, his most-preferred choice.
   • Unmatched men: 2, 3
2. • Man 2 proposes to woman 3, his most-preferred choice.
   • Unmatched men: 3
3. • Man 3 proposes to woman 3, his most-preferred choice.
   • Woman 3 now dumps man 2.
   • Unmatched men: 2
4. • Man 2 proposes to woman 1, his most-preferred available choice.
   • Woman 1 now dumps man 1.
   • Unmatched men: 1
5. • Man 1 proposes to woman 2, his most-preferred available choice.
   • All men are now matched.

The man-optimal stable matching is therefore:

The package computes the Gale-Shapley algorithm using the function galeShapley.marriageMarket:

matching = galeShapley.marriageMarket(uM, uW)

Note that we can obtain equivalent results when we use prefM and prefW as arguments:

matching = galeShapley.marriageMarket(proposerPref = prefM, reviewerPref = prefW)

The function galeShapley.marriageMarket returns a list matching that includes the vectors proposals and engagements with the final proposals and engagements, respectively. These two vectors contain the same information (i.e. they tell us who is matched with whom). For the example above, the vector of proposals contains

matching$proposals

##        cols
## rows    Proposed to Woman
##   Man 1                 2
##   Man 2                 1
##   Man 3                 3

The first element in the vector tells us that man 1 is matched with woman 2. Man 2 is matched to woman 1, and man 3 is matched to woman 3. The vector of engagement contains

matching$engagements

##          cols
## rows      Engaged to Man
##   Woman 1              2
##   Woman 2              1
##   Woman 3              3

The first element in the vector tells us that woman 1 is matched to man 2, woman 2 will be matched to man 1, and woman 3 will be matched to man 3.

We can then check if the computed matching is stable using the function checkStability. To check if a matching is stable, we check for each assignment (m,w) if there is some other woman w' that man m would rather be matched with and who would rather be matched to man m. This function will return true if the matching is stable and false otherwise.

galeShapley.checkStability(uM, uW, matching$proposals, matching$engagements)

## [1] TRUE

For the simple 3-by-3 example, we can perform this check by hand:

• Man 1 is matched to woman 2, his second-most preferred choice. His most preferred choice is woman 1. Woman 1 is matched with man 2 who she prefers over man 1. Thus man 1 cannot do better than woman 2.
• Man 2 is matched to woman 1, his second-most preferred choice. His most preferred woman is woman 3, who is matched with man 3. Since man 3 is her most-preferred choice, man 2 cannot do better than woman 1.
• Man 3 is matched to women 3, his most preferred choice, so he cannot do better.

Thus, this matching is stable.

The following examples illustrate some additional features of this package.

# set seed
set.seed(1)
# set number of men
nmen = 2500
# set number of women
nwomen = 2000
# generate preferences
uM = matrix(runif(nmen*nwomen), nrow = nwomen, ncol = nmen)
uW = matrix(runif(nmen*nwomen), nrow = nmen, ncol = nwomen)
# male-optimal matching
resultsM = galeShapley.marriageMarket(uM, uW)
str(resultsM)

## List of 4
##  $ proposals       : num [1:2500, 1] 927 1644 1965 1851 349 ...
##  $ engagements     : num [1:2000, 1] 386 471 1582 598 1657 ...
##  $ single.proposers: num [1:500] 8 11 17 19 26 30 34 37 49 53 ...
##  $ single.reviewers: num(0)

galeShapley.checkStability(uM, uW, resultsM$proposals, resultsM$engagements)

## [1] TRUE

# female-optimal matching
resultsW = galeShapley.marriageMarket(uW, uM)
str(resultsW)

## List of 4
##  $ proposals       : num [1:2000, 1] 386 471 1582 598 1657 ...
##  $ engagements     : num [1:2500, 1] 927 1644 1965 1851 349 ...
##  $ single.proposers: num(0) 
##  $ single.reviewers: num [1:500] 8 11 17 19 26 30 34 37 49 53 ...

galeShapley.checkStability(uW, uM, resultsW$proposals, resultsW$engagements)

## [1] TRUE

# set seed
set.seed(1)
# set number of students
nstudents = 1000
# set number of colleges
ncolleges = 400
# generate preferences
uStudents = matrix(runif(ncolleges*nstudents), nrow = ncolleges, ncol = nstudents)
uColleges = matrix(runif(nstudents*ncolleges), nrow = nstudents, ncol = ncolleges)
# student-optimal matching
results = galeShapley.collegeAdmissions(studentUtils =  uStudents, collegeUtils =  uColleges, slots = 2)
str(results)

## List of 4
##  $ unmatched.students: num [1:200] 3 15 23 29 30 36 44 46 48 49 ...
##  $ unmatched.colleges: int(0) 
##  $ matched.colleges  : num [1:400, 1:2] 728 887 28 372 875 775 456 937 402 774 ...
##  $ matched.students  : int [1:1000, 1] 52 70 NA 210 155 170 238 16 371 391 ...

# check if matching is stable
galeShapley.checkStability(uStudents, uColleges, results$matched.students, results$matched.colleges)

## [1] TRUE

One-sided Matching Markets: Irving’s Algorithm

This package implements the algorithm by Irving (1985) for one-sided matching markets.

Consider the following example: A set of n potential roommates, each with ranked preferences over all the other potential roommates, are to be matched to rooms, two roommates per room. A matching is stable if there is no roommate r1 that would rather be matched to some other roommate d2 than to his current roommate r2 and the other roommate d2 would rather be matched to r1 than to his current roommate d1.

Preferences of potential roommates are summarized by an n − 1 × n dimensional matrix, e.g., if n = 6,

pref = matrix(c(3, 6, 2, 5, 3, 5,
                4, 5, 4, 2, 1, 1,
                2, 4, 5, 3, 2, 3,
                6, 1, 1, 6, 4, 4,
                5, 3, 6, 1, 6, 2), nrow = 5, ncol = 6, byrow = TRUE)

Column i represents the preferences of the ith roommate, and row j represents the ranking of the roommate whose index is encoded in that row. For example, in the above preference matrix, roommate 1 most prefers to be matched with roommate 3, followed by 4, followed by 2.

The function roommate.checkPreferences checks if a given preference order is complete, i.e. if all preferences are fully specified. If the preference order is complete, it will return the proper preference order using C++ style indexing (beginning at 0 instead of at 1). If the preference order is incomplete, the function will return an error.

roommate.checkPreferences(pref)

##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]    2    5    1    4    2    4
## [2,]    3    4    3    1    0    0
## [3,]    1    3    4    2    1    2
## [4,]    5    0    0    5    3    3
## [5,]    4    2    5    0    5    1

The algorithm proceeds in two phases.

In phase 1, potential roommates take turns sequentially proposing to the other roommates. Each roommate who is proposed to can accept or reject the proposal. A roommate accepts if he currently has made no better proposal which was accepted to another roommate. If a roommate has accepted a proposal, and then receives a better proposal, he rejects the old proposal and substitutes in the new proposal.

In the above example,

1. Roommate 1 begins by proposing to roommate 3, his most preferred roommate. 3, having no better offers, accepts.
2. 2 proposes to 6, who accepts.
3. 3 proposes to 2, who accepts.
4. 4 proposes to 5, who accepts.
5. 5 proposes to 3, who accepts. 3 cancels his proposal from 1.
6. 1, having no proposal, proposes to 4, who accepts.
7. 6 proposes to 5, who rejects, having a better proposal from 4.
8. 6 proposes to 1, who accepts.

In phase 2, we begin by eliminating all potential roommate matches which are worse than the current proposals held. For example, in the above example, 3 has a proposal from 5, and so we eliminate 1 and 6 from 3‘s column, and symmetrically eliminate 3 from 1 and 6’s column. This results in the following ’reduced’ preference listing:

   6, 2, 5, 3,
4, 5, 4, 2,    1
2, 4, 5, 3, 2,
6, 1,    6, 4, 4
   3,    1,    2

These preferences form what is called a ‘stable’ table, or, ‘s-table’. (‘Stable’ for short.) The defining characteristic of a stable table is that if i is the most preferred roommate on js list, then j is the least preferred roommate on is list. For example, 1 most prefers 4, but 4 least prefers 1.

The algorithm proceeds by finding and eliminating ‘rotations’. A rotation is a sequence of pairs of roommates, such that there is a distinct roommate in the first position of each pair, the second roommate in each pair least prefers the roommate he is paired with, the first roommate in each pair most prefers the roommate he is paired with, and finally, the first roommate in each pair ranks the second roommate in the following pair second (modulo the number of pairs, that is, the first roommate in the last pair ranks the second roommate in the first pair second) Once a rotation has been identified, removing it results in another stable table.

For example, (1, 4), (3, 2) is a rotation in the above table, because 1 loves 4, 3 loves 2, 4 hates 1, 2 hates 3, 2 is second on 1s list, and 4 is second on 3’s list. Eliminating this rotation involves 2 rejecting 3, 4 rejecting 1, and then we remove every successive potential roommate as well to preserve the stable table property, resulting in

   6,    5, 3,
   5, 4, 2,    1
2, 4, 5, 3, 2,
6, 1,    6, 4, 4
               2

A further rotation is (1, 2), (2, 6), (4, 5). Eliminating it yields

            3,
   5, 4, 2,    1
   4, 5, 3, 2,
6,

A final rotation is (2, 5), (3, 4). Eliminating it yields

            3,
         2,    1
   4, 5,
6,

Therefore, a stable matching is for 1 and 6 to match, 2 and 4 to match, and 3 and 5 to match.

results = roommate(pref = pref)
results

##      [,1]
## [1,]    6
## [2,]    4
## [3,]    5
## [4,]    2
## [5,]    3
## [6,]    1

The function roommate.checkStability can be used to check if the resulting matching is stable:

roommate.checkStability(pref = pref, matching = results)

## [1] TRUE

# generate preferences
N = 10
u = matrix(runif(N^2),  nrow = N, ncol = N)
results = roommate(utils = u)
results

##       [,1]
##  [1,]    3
##  [2,]    4
##  [3,]    1
##  [4,]    2
##  [5,]    7
##  [6,]   10
##  [7,]    5
##  [8,]    9
##  [9,]    8
## [10,]    6

roommate.checkStability(utils = u, matching = results)

## [1] TRUE

set.seed(1)
N = 512
u = matrix(runif(N^2),  nrow = N, ncol = N)
results = roommate(utils = u)
print(results)

## NULL

Kidneys and Housing: The Top Trading Cycle Algorithm

This package implements the top trading cycle algorithm.

Consider the following problem: A set of n agents each currently own their own home, and have preferences over the homes of other agents. The problem is to trade the homes between the agents in such a way so that no two agents want to swap homes.

Preferences of agents are summarized by an n × n dimensional matrix, e.g., if n = 4,

pref = matrix(c(4, 4, 2, 4,
                2, 1, 1, 1,
                1, 2, 3, 3,
                3, 3, 4, 2), nrow = 4, ncol = 4, byrow = TRUE)

Column i represents the preferences of the ith agent, and row j represents the ranking of the roommate whose index is encoded in that row. For example, in the above preference matrix, agent 1 most prefers the home of agent 4, followed by 2, followed by 1, followed by 3.

Roughly speaking, the top trading cycle proceeds by identifying cycles of agents, then eliminating those cycles until no agents remain. A cycle is a sequence of agents such that each agent most prefers the next agent’s home (out of the remaining unmatched agents), and the last agent in the sequence most prefers the first agent in the sequence’s home.

4,  4,  2,  4
2,  1,  1,  1
1,  2,  3,  3
3,  3,  4,  2

For example, for the above preference matrix, when all the agents are unmatched, the only rotation is {4}, representing the fact that agent 4 most prefers his own house. Therefore, the algorithm begins by matching agent 4 to himself, and then removing him from the pool:

        2
2,  1,  1
1,  2,  3
3,  3,

Now, a rotation is {1, 2}, because 1 most prefers 2s house, and 2 most prefers 1s house. So agents 1 and 2 will swap homes, leaving agent 3 all by himself.

        3

Therefore, the final matching is that agent 1 swaps with agent 2, and agents 3 and 4 keep their own homes.

results = toptrading(pref = pref)
results

##      [,1]
## [1,]    2
## [2,]    1
## [3,]    3
## [4,]    4

The function toptrading.checkStability can be used to check if the resulting allocation is in fact stable:

toptrading.checkStability(pref = pref, matchings = results)

## [1] TRUE