“Almost-stable” matchings in the Hospitals / Residents problem with Couples

Abstract

The Hospitals / Residents problem with Couples (hrc) models the allocation of intending junior doctors to hospitals where couples are allowed to submit joint preference lists over pairs of (typically geographically close) hospitals. It is known that a stable matching need not exist, so we consider min bp hrc, the problem of finding a matching that admits the minimum number of blocking pairs (i.e., is “as stable as possible”). We show that this problem is NP-hard and difficult to approximate even in the highly restricted case that each couple finds only one hospital pair acceptable. However if we further assume that the preference list of each single resident and hospital is of length at most 2, we give a polynomial-time algorithm for this case. We then present the first Integer Programming (IP) and Constraint Programming (CP) models for min bp hrc. Finally, we discuss an empirical evaluation of these models applied to randomly-generated instances of min bp hrc. We find that on average, the CP model is about 1.15 times faster than the IP model, and when presolving is applied to the CP model, it is on average 8.14 times faster. We further observe that the number of blocking pairs admitted by a solution is very small, i.e., usually at most 1, and never more than 2, for the (28,000) instances considered.

Electronic supplementary material

The online version of this article (doi:10.1007/s10601-016-9249-7) contains supplementary material, which is available to authorized users.

1 Introduction

The Hospitals / Residents problem

The Hospitals / Residents problem (hr) [13] is a many-to-one allocation problem that models the assignment process involved in centralised matching schemes such as the National Resident Matching Program (NRMP) [42] which assigns graduating medical students to hospital posts in the USA. Analogous schemes exist in Canada [37] and Japan [39]. A similar process was used until recently to match medical graduates to Foundation Programme places in Scotland: the Scottish Foundation Allocation Scheme (SFAS) [19]. Moreover, similar matching schemes exist in the context of Higher Education admission in Hungary [4, 40], Spain [29], Turkey [3] and Ireland [38, 40]. The reader is referred to [40] for details of matching practices in a number of practical contexts throughout Europe.

An instance of hr consists of two sets of agents – a set \(R= \{r_{1} ,{\ldots } r_{n_{1}}\}\) containing residents and a set \(H= \{h_{1} ,{\ldots } h_{n_{2}}\}\) containing hospitals. Every resident expresses a linear preference over some subset of the hospitals, his preference list. The hospitals in a resident’s preference list are his acceptable partners; all other hospitals being unacceptable. Every hospital expresses a linear preference over those residents who find it acceptable. Further, each hospital hj∈H has a positive integral capacitycj, the maximum number of residents to which it may be assigned. A matchingM is a set of acceptable resident-hospital pairs such that each resident appears in at most one pair and each hospital hj belongs to at most cj pairs. If (ri,hj)∈M then ri is said to be assigned to hj, M(ri) denotes hj, and ri is an assignee of hj. Given ri∈R, if ri does not belong to any pair in M then ri is said to be unassigned. Given hj∈H, we let M(hj) denote the set of assignees of hj in M. Hospital hj is undersubscribed, full or oversubscribed according as |M(hj)| is less than, equal to, or larger than cj, respectively.

Roth [31] argued that a key property to be satisfied by any matching M in an instance I of hr is stability, which ensures that M admits no blocking pair in I. Informally, such a pair comprises a resident ri and a hospital hj, both of whom have an incentive to disregard their assignments (if any) and become matched to one another outside of M, undermining its integrity. A matching is stable if it admits no blocking pair. It is known that every instance of hr admits at least one stable matching, which can be found in time linear in the size of the instance [13].

The Hospitals / Residents problem with Couples

The Hospitals / Residents problem with Couples (hrc) is a generalisation of hr that is important in practical applications because it models the case where some of the residents may apply jointly in couples, so that they may be matched to hospitals that are geographically close to one another. In order to ensure this, a couple submits a joint preference list over pairs of hospitals, rather than individual hospitals. Matching schemes for junior doctors such as the NRMP [42] allow couples to apply jointly, as do assignment processes in the US Navy [28, 34, 36] (for which hrc is an appropriate problem model), for example.

Formally, an instance I of hrc consists of a set \(R=\{r_{1} ,{\ldots } r_{n_{1}}\}\) containing residents and a set \(H=\{h_{1} ,{\ldots } h_{n_{2}}\}\) containing hospitals. The residents in R are partitioned into two sets, S and S′. The set S consists of single residents and the set S′ consists of those residents involved in couples. There is a set C={(ri,rj):ri,rj∈S′} of couples such that each resident in S′ belongs to exactly one pair in C.

Each single resident ri∈S expresses a linear preference order over some subset of the hospitals, his acceptable hospitals; all other hospitals being unacceptable. Each couple (ri,rj)∈C expresses a joint linear preference order over a subset A of H×H where (hp,hq)∈A represents the simultaneous assignment of ri to hp and rj to hq. The hospital pairs in A represent those joint assignments that are acceptable to (ri,rj), all other joint assignments being unacceptable. Each hospital hj∈H expresses a linear preference order over those residents who find it acceptable, either as a single resident or as part of a couple, and as in the case of hr, each hospital hj∈H has a positive integral capacitycj.

A matchingM in I is defined as in hr case, with the additional restriction that, for each couple (ri,rj)∈C, either both ri and rj appear in no pair of M, or else {(ri,hk),(rj,hl)}⊆M for some pair (hk,hl) that (ri,rj) find acceptable. In the former case, (ri,rj) are said to be unassigned, whilst in the latter case, (ri,rj) are said to be jointly assigned to (hk,hl). Given a resident ri∈R, the definitions of M(ri), assigned and unassigned are the same as for the hr case, whilst for a hospital hj∈H, the definitions of assignees, M(hj), undersubscribed, full and oversubscribed for hospitals are also the same as before.

We seek a stable matching, which guarantees that no resident and hospital, and no couple and pair of hospitals, have an incentive to deviate from their assignments and become assigned to each other outside of the matching. Roth [31] considered stability in the hrc context but did not define the concept explicitly. Whilst Gusfield and Irving [15] gave a formal definition of a blocking pair, it neglected to deal with the case that both members of a couple may wish to be assigned to the same hospital. A number of other stability definitions for hrc have since been given in the literature that address this issue (see [6] and [20, Section 5.3] for more details), including that of McDermid and Manlove [24], which we adopt in this paper. We repeat their definition again here for completeness.

Let I be an instance of hrc. A matching M is stable in I if none of the following holds:

1.

There is a single resident ri and a hospital hj, where ri finds hj acceptable, such that either ri is unassigned in M or prefers hj to M(ri), and either hj is undersubscribed in M or prefers ri to some member of M(hj).

2.

There is couple (ri,rj) and a hospital hk such that either

(a)

(ri,rj) prefers (hk,M(rj)) to (M(ri),M(rj)), and either hk is undersubscribed in M or prefers ri to some member of M(hk)∖{rj}or

(b)

(ri,rj) prefers (M(ri),hk) to (M(ri),M(rj)), and either hk is undersubscribed in M or prefers rj to some member of M(hk)∖{ri}.

3.

There is a couple (ri,rj) and a pair of (not necessarily distinct) hospitals hk≠M(ri), hl≠M(rj) such that (ri,rj) finds (hk,hl) acceptable, and either (ri,rj) is unassigned or prefers the joint assignment (hk,hl) to (M(ri),M(rj)), and either

(a)

hk≠hl, and hk (respectively hl) is either undersubscribed in M or prefers ri (respectively rj) to at least one of its assignees in M; or

(b)

hk = hl, and hk has two free posts in M, i.e., ck−|M(hk)|≥2; or

(c)

hk = hl, and hk has one free post in M, i.e., ck−|M(hk)|=1, and hk prefers at least one of ri,rj to some member of M(hk); or

(d)

hk = hl, hk is full in M, hk prefers ri to some rs∈M(hk), and hk prefers rj to some rt∈M(hk)∖{rs}.

A resident and hospital, or a couple and hospital pair, satisfying one of the above conditions, is called a blocking pair of M and is said to blockM.

Existing algorithmic results for hrc

An instance I of hrc need not admit a stable matching [31]. We call Isolvable if it admits a stable matching, and unsolvable otherwise. Also an instance of hrc may admit stable matchings of differing sizes [2]. Further, the problem of deciding whether a stable matching exists in an instance of hrc is NP-complete, even in the restricted case where there are no single residents and each hospital has capacity 1 [25, 30]. The decision problem is also W[1]-hard [22] when parameterized by the number of couples.

In many practical applications of hrc the residents’ preference lists are short. Let (α,β,γ)-hrc denote the restriction of hrc in which each single resident’s preference list contains at most α hospitals, each couple’s preference list contains at most β pairs of hospitals and each hospital’s preference list contains at most γ residents. Biró et al. [7] showed that deciding whether an instance of (0,2,2)-hrc admits a stable matching is NP-complete.

Heuristics for hrc were described and compared experimentally by Biró et al. [5]. As far as exact algorithms are concerned, Biró et al. [7] gave an Integer Programming (IP) formulation for finding a maximum cardinality stable matching (or reporting that none exists) in an arbitrary instance of hrc and presented an empirical evaluation of an implementation of their model, showing that their formulation was capable of solving instances of the magnitude of those arising in the SFAS application. Further algorithmic results for hrc are given in [6, 20, 23].

Most-stable matchings

Given that a stable matching need not exist in a given hrc instance I, a natural question to ask is whether there is some other matching that might be the best alternative amongst the matchings in I. Roth [32, 33] argued that instability in the outcome of an allocation process gives participants a greater incentive to circumvent formal procedures; it follows minimising the amount of instability might be a desirable objective. Eriksson and Häggström [11] suggested that the number of blocking pairs admitted by a matching is a meaningful way to measure its degree of instability.

Define bp(M) to be the set of blocking pairs relative to a matching M in I, and define a most-stable matching to be a matching M for which |bp(M)| is minimum, taken over all matchings in I. Clearly if I admits a stable matching M, then M is a most-stable matching in I. Let min bp hrc denote the problem of finding a most-stable matching, given an instance of hrc. Most-stable matchings have been studied from an algorithmic point of view in various matching problem contexts [1, 8, 9, 12, 16, 17] (see [20] for more details), including in humanitarian organisations [35]. Define (α,β,γ)-min bp hrc to be the restriction of min bp hrc to instances of (α,β,γ)-hrc.

Contribution of this work

In Section 2 we show that (∞,1,∞)-min bp hrc is NP-hard and not approximable within \(n_{1}^{1-\varepsilon }\), for any ε>0, unless P=NP (recall that n1 is the number of residents in a given instance). In this highly restricted case of min bp hrc, each couple finds only one hospital pair acceptable and each hospital has capacity 1 (∞ refers to preference lists of unbounded length). We also show that (∞,∞,1)-min bp hrc and (2,1,2)-min bp hrc are solvable in polynomial time. These results help to narrow down the search for the boundary between polynomial-time solvable and NP-hard restrictions of min bp hrc (recall that (0,2,2)-min bp hrc is NP-hard [7]).

In Section 3 we present the first IP model for min bp hrc; indeed this model can be used to find a most-stable matching of maximum cardinality. This formulation extends our earlier IP model for hrc, presented in [7]. Then in Section 4 we present data from an empirical evaluation of an implementation of the IP model for min bp hrc applied to randomly-generated instances. We measure the mean solution time, mean size of a most-stable matching and mean number of blocking pairs admitted by a most-stable matching when varying (i) the number of residents, (ii) the number of couples, (iii) the number of hospitals and (iv) the lengths of the residents’ preference lists. Our main finding is that, over the 28,000 instances considered, the number of blocking pairs admitted by a most-stable matching is very small: it is usually at most 1, and never more than 2. This suggests that in a given hrc instance in practice, even if a stable matching does not exist, we may be able to find a matching with only a very small amount of instability.

Finally, in Section 5 we present the first Constraint Programming (CP) model for min bp hrc and evaluate its performance compared to the IP model over the instances used for the empirical analysis in Section 4. We observe that on average, the CP model is about 1.15 times faster than the IP model, and when presolving is applied to the CP model, it is on average 8.14 times faster.

Related work

Drummond et al. [10] presented SAT and IP encodings of hrc and investigated empirically their performance, along with two earlier heuristics for hrc, on randomly-generated instances. Their main aim was to measure the time taken to find a stable matching or report that none exists, and the proportion of solvable instances. They found that the SAT encoding gave the fastest method and was generally able to resolve the solvability question for the highest proportion of instances. In another paper [27], the same authors conducted further empirical investigations on random instances using an extension of their SAT encoding to determine how many stable matchings were admitted, and whether a resident Pareto optimal stable matching existed. We remark that the results in [10, 27] are not directly comparable to ours, because the stability definition considered in those papers is slightly weaker than that given by Definition 1. See Section A of the online supplement for a discussion of this issue.

Hinder [18] presented an IP model for a general stable matching problem with contracts, which includes hrc as defined here, as a special case. He conducted an empircal study on randomly-generated instances, comparing the performance of the IP model, its LP relaxation and a previously-published heuristic. Hinder showed that the LP relaxation finds stable matchings (when they exist) with much higher probability than the heuristic, and with probability quite close to the true value given by the IP model. The IP model terminates surprisingly quickly when the number of residents belonging to a couple is 10 %, but it should be emphasised that in Hinder’s random instances, all hospitals have capacity 1. In such a case our IP/CP models would be much simpler and need not involve the constraints corresponding to stability criteria 3(b), 3(c) and 3(d) in Definition 1, thus our runtime results are not directly comparable to Hinder’s.

To the best of our knowledge there have been no previous CP models for hrc, though a CP model for hr was given in [21], extending an earlier CP model for the classical Stable Marriage problem, the 1-1 restriction of hr [14]. A detailed survey of CP models for stable matching problems is given in [20, Section 2.5].

Nguyen and Vohra [26] proved a remarkable result, namely that it is always possible to find a stable matching in an instance of hrc if the capacity of each hospital can be adjusted (up or down) by at most 4, with the total capacity of the hospitals increasing by at most 9.

2 Complexity results for min bp hrc

In this section we present complexity and approximability results for min bp hrc in the case that preference lists of some or all of the agents are of bounded length. We begin with (∞,1,∞)-min bp hrc, the restriction in which each couple lists only one hospital pair on their preference list. Even in this highly restricted case, the problem of finding a most-stable matching is NP-hard and difficult to approximate. The proof of this result, given in Section B of the online supplement, begins by showing that, given an instance of (∞,1,∞)-hrc, the problem of deciding whether a stable matching exists is NP-complete. Then a gap-introducing reduction is given from this problem to (∞,1,∞)-min bp hrc.

Theorem 2

(∞,1,∞)-min bp hrcis NP-hard and not approximable within a factor of\(n_{1}^{1-\varepsilon }\), for any ε>0, unless P=NP, where n1is the number of residents in a given instance. The result holds even if each hospital has capacity 1.

We now turn to the case that hospitals’ lists are of bounded length. It will be helpful to introduce the notion of a fixed assignment in a given hrc instance I. This involves either (i) a resident-hospital pair (ri,hj) such that hj is the first choice of ri, and ri is among the first cj choices of hj, or (ii) a pair comprising a couple (ri,rj) and a pair of hospitals (hp,hq) such that hp (resp. hq) is the first choice of ri (resp. rj), and ri (resp. rj) is among the first cp (resp. cq) choices of hp (resp. hq). Clearly any stable matching must contain all the fixed assignments in I. By eliminating the fixed assignments iteratively, we arrive at the following straightforward result for (∞,∞,1)-hrc (the proofs of all the results stated in this section from this point onwards can be found in Section C of the online supplement).

Proposition 3

An instance I of (∞,∞,1)-hrcadmits exactly one stable matching, which can be found in polynomial time.

We now consider the (2,1,2)-hrc case. The process of satisfying a fixed assignment involves matching together the resident(s) and hospital(s) involved, deleting the agents themselves (and removing them from the remaining preference lists). This may uncover further fixed assignments, which themselves can be satisfied. Once this process terminates, we say that all fixed assignments have been iteratively satisfied. Let I be the (2,1,2)-hrc instance that remains. It turns out that I has a special structure, as the following result indicates.

Lemma 4

An arbitrary instance of (2,1,2)-hrcinvolving at least one couple and in which all fixed assignments have been iteratively satisfied must be constructed from sub-instances of the form shown in Fig. 1, in which all of the hospitals have capacity 1.

An instance of (2,1,2)-hrc containing an arbitrary number of couples and an arbitrary number of residents that has no unsatisfied fixed assignments. Here residents with a subscript s are single residents, whilst those with a subscript c belong to couples. The structure of this instance is described in more detail in Section C of the online supplement

It is then straightforward to find a most-stable matching in each such sub-instance.

Lemma 5

Let I′be an instance of (2,1,2)-hrcof the form shown in Fig.1. If I′has an even number of couples then I′admits a stable matching M. Otherwise I′admits a matching M such that |bp(M)|=1 in I′.

Using Lemmas 4 and 5, it follows that we can find a most-stable matching in an instance I of (2,1,2)-hrc as follows. Assume that M0 is the matching in I in which all fixed assignments have been iteratively satisfied, and assume that the corresponding deletions have been made from the preference lists in I, yielding instance I′. Lemma 4 shows that I′ is a union of disjoint sub-instances I1,I2,…,It, where each Ij is of the form shown in shown in Fig. 1 (1≤j≤t). Let j (1≤j≤t) be given and let Nj be the number of couples in Ij. Lemma 5 implies that, if Nj is even, we may find a stable matching Mj in Ij, otherwise we may find a matching Mj in Ij such that |bp(Mj)|=1 in Ij. It follows that \(M=\cup _{j=0}^{t} M_{j}\) is a most-stable matching in I. This leads to the following result.

Theorem 6

(2,1,2)-min bp hrcis solvable in polynomial time.

It remains open to resolve the complexity of (p,1,q)-hrc for constant values of p and q where max{p,q}≥3.

3 An integer programming formulation for min bp hrc

In this section we describe our IP model for min bp hrc, which extends the earlier IP model for hrc presented in [7] (we discuss relationships between the two models at the end of this section). Let I be an instance of hrc where \(R=\{r_{1},r_{2},\dots ,r_{n_{1}}\}\) is the set of residents and \(H=\{h_{1},h_{2},\dots ,h_{n_{2}}\}\) is the set of hospitals; we will denote by J the IP model corresponding to I. To streamline the exposition we will only present some of the constraints in J; the full description of J is contained in Section D of the online supplement.

The IP model J is based on modelling the various types of blocking pairs that might arise according to Definition 1, and allowing them to be counted by imposing a series of linear inequalities. The variables are defined for each resident, whether single or a member of a couple, and for each element on his preference list (with the possibility of being unassigned). A further consistency constraint ensures that each member of a couple obtains hospitals from the same pair in their list, if assigned. A suitable objective function then enables the number of blocking pairs to be minimised. Subject to this, we may also maximise the size of the constructed matching.

Notation

We first define some required notation in I. Without loss of generality, suppose residents r1,r2…r2c are in couples. Thus \(r_{2c+1}, r_{2c+2}{\ldots } r_{n_{1}}\) comprise the single residents. Again, without loss of generality, suppose that the couples are (r2i−1,r2i)(1≤i≤c). A crucial component of the IP model is a mapping between the joint preference list of a couple \(\mathcal C_{i} = (r_{2i-1}, r_{2i})\) and individual preference lists for r2i−1 and r2i. Suppose that the joint preference list of \(\mathcal C_{i}\) is

From this list we say that \(h_{\alpha _{1}}, h_{\alpha _{2}}{\ldots } h_{\alpha _{l}}\) and \(h_{\beta _{1}}, h_{\beta _{2}}{\ldots } h_{\beta _{l}}\) are the individual preference lists for r2i−1 and r2i respectively. Note that a given hospital hj may appear more than once in the individual preference list of a resident belonging to a couple.

For a resident ri∈R (whether single or a member of a couple), let l(ri) denote the length of a resident ri’s individual preference list. Moreover let pref (ri,p) denote the hospital at position p of ri’s individual preference list.

For a hospital hj∈H, let l(hj) denote the length of hj’s preference list over individual residents. For an acceptable resident-hospital pair (ri,hj), let rank (hj,ri) = q denote the rank that hospital hj assigns resident ri, where 1≤q≤l(hj). Thus, rank (hj,ri) is equal to the number of residents that hj prefers to ri plus 1. Further, for each j(1≤j≤n2) and q(1≤q≤l(hj)), let the set R(hj,q) contain resident-position pairs (ri,p) such that ri∈R is assigned a rank of q by hj and hj is in position p(1≤p≤l(ri)) on ri’s individual list. Hence

Variables in the IP model

For each i(1≤i≤n1) and p(1≤p≤l(ri)), J has a variable xi,p∈{0,1} such that xi,p=1 if and only if ri is assigned to his pth-choice hospital. Also, for each i (1≤i≤n1) and p = l(ri)+1, J has a variable xi,p∈{0,1} such that xi,p=1 if and only if ri is unassigned. Let X={xi,p:1≤i≤n1∧1≤p≤l(ri)+1}.

J also contains variables 𝜃i,p∈{0,1} for each i(1≤i≤n1) and p(1≤p≤l(ri)). The intuitive meaning of a variable 𝜃i,p is that 𝜃i,p=1 if and only if resident ri is involved in a blocking pair with the hospital at position p on his individual preference list, either as a single resident or as part of a couple.

Constraints in the IP model

We firstly add constraints to J which force every variable to be binary valued. Next we ensure that matching constraints are satisfied, as follows. As each resident ri∈R is assigned to exactly one hospital or is unassigned (but not both), \({\sum }_{p=1}^{l(r_{i})+1} x_{i,p}=1\) must hold for all i(1≤i≤n1). Similarly, since a hospital hj may be assigned at most cj residents, xi,p=1 where pref (ri,p) = hj for at most cj residents, and hence for all j(1≤j≤n2), \({\sum }_{i=1}^{n_{1}}{\sum }_{p=1}^{l(r_{i})} \{x_{i,p} \in X :\,\, \)pref (ri,p) = hj}≤cj must hold.

For each couple (r2i−1,r2i), r2i−1 is unassigned if and only if r2i is unassigned, and r2i−1 is assigned to the hospital in position p in their individual list if and only if r2i is assigned to the hospital in position p in their individual list. Thus for all i(1≤i≤c) and p(1≤p≤l(r2i−1)+1), x2i−1,p = x2i,p must hold,

The remaining constraints in J allow the number of blocking pairs of a given matching to be counted. Each such constraint deals with a specific type of blocking pair that satisfies a given part of Definition 1. It allows a blocking pair to exist involving either (i) a single resident ri with the hospital at some position p on his list, or (ii) a couple (r2i−1,r2i) with the hospital pair at some position p on their joint list, if and only if 𝜃i,p=1. We illustrate the construction of J by giving the constraint corresponding to so-called “Type 1” blocking pairs, involving involve single residents, where Condition 1 of Definition 1 is satisfied. The other constraints may be dealt with in a similar fashion – see Section D of the online supplement for further details.

Type 1 blocking pairs

In a matching M in I, if a single resident ri∈R is unassigned or has a worse partner than some hospital hj∈H where pref (ri,p) = hj and rank (hj,ri) = q then hj must be fully subscribed with better partners than ri, for otherwise (ri,hj) blocks M. Hence if ri is unassigned or has worse partner than hj, i.e., \(\sum \limits _{p^{\prime }=p+1}^{l(r_{i})+1} x_{i,p^{\prime }}=1\), and hj is not fully subscribed with better partners than ri, i.e., \(\sum \limits _{q^{\prime }=1}^{q-1} \{x_{i^{\prime },p^{\prime \prime }} \in X : (r_{i^{\prime }}, p^{\prime \prime }) \in R(h_{j}, q^{\prime })\} < c_{j}\), then we require 𝜃i,p=1 to count this blocking pair. Thus, for each i(2c+1≤i≤n1) and p(1≤p≤l(ri)) we obtain the following constraint where pref (ri,p) = hj and rank (hj,ri) = q:

Objective functions in the IP model

A maximum cardinality most-stable matching M is a matching of maximum cardinality, taken over all most-stable matchings in I. To compute a maximum most-stable matching in J, we apply two objective functions in sequence.

First we find an optimal solution in J that minimises the number of blocking pairs. To this end we apply the objective function \(\min \sum \limits _{i=1}^{n_{1}} \sum \limits _{p=1}^{l(r_{i})} \theta _{i,p}\).

The matching M corresponding to an optimal solution in J will be a most-stable matching in I. Let k=|bp(M)|. Now we seek a maximum cardinality matching in I with at most k blocking pairs. Thus we add the following constraint to J, which ensures that, when maximising on cardinality, any solution also has at most k blocking pairs: \(\sum \limits _{i=1}^{n_{1}} \sum \limits _{p=1}^{l(r_{i})} \theta _{i,p} \leq k\).

The final step is to maximise the size of the matching, subject to the matching being most-stable. This involves optimising for a second time, this time using the following objective function: \(\max \sum \limits _{i=1}^{n_{1}} \sum \limits _{p=1}^{l(r_{i})} x_{i,p}\).

The following result, which establishes the correctness of the IP formulation, is proved in Section D of the online supplement.

Theorem 7

Given an instance I ofmin bp hrc, let J be the corresponding IP model as defined above. A maximum cardinality most-stable matching in I is exactly equivalent to an optimal solution to J.

We remark that the IP model presented in this section develops the earlier model for hrc [7] with the addition of the 𝜃i,p variables. There are similarities between the constraints (with these variables omitted) when comparing the two models. However in the hrc model [7] essentially all stability constraints had to be satisfied, whereas in the min bp hrc model a blocking pair is allowed at the expense of a 𝜃i,p variable having value 1, which allows the number of blocking pairs to be counted. Suitable placement of the 𝜃i,p variables within the constraints from the hrc model allows this condition on the 𝜃i,p variables to be enforced.

4 Empirical results from the IP model for min bp hrc

In this section we present data from an empirical evaluation of an implementation of the IP model for finding a maximum cardinality most-stable matching in an instance of min bp hrc. We considered the following properties for randomly-generated hrc instances: the time taken to find a maximum cardinality most-stable matching, the size of a maximum cardinality most-stable matching and the number of blocking pairs admitted by a most-stable matching. We show how these properties varied as we modified the number of residents, the percentage of residents involved in couples, the number of hospitals and the lengths of residents’ preference lists in the constructed instances.

Methodology

We ran all the experiments on an implementation of the IP model using the CPLEX 12.4 Java Concert API applied to randomly-generated instances of hrc.1 In these instances, the preference lists of residents and hospitals were constructed to take into account of the fact that, in reality, some hospitals and residents are more popular than others, respectively. Typically, the most popular hospital in the SFAS context had 5-6 times as many applicants as the least popular, and the numbers of applicants to the other hospitals were fairly uniformly distributed between the two extremes. Our constructed instances reflected this real-world behaviour. For more details about the construction of the instances and the correctness testing methodology, the reader is referred to [23, Chapters 6,7].

All experiments were carried out on a desktop PC with an Intel i5-2400 3.1Ghz processor with 8Gb of memory running Windows 7. To find a most-stable matching in an instance I of hrc we applied the following procedure. We first used the hrc IP implementation presented in [7] to find a maximum cardinality stable matching M in I if one exists. Clearly, if I is solvable then M is a maximum cardinality most-stable matching. However, if I was found to be unsolvable, we applied the min bp hrc IP model to I. In this case we applied a lower bound of 1 to the number of blocking pairs in a most-stable matching in I since we knew that no stable matching existed. All instances were allowed to run to completion. We remark that the min bp hrc model appears to be much more difficult to solve than the hrc model presented in [7], and thus the largest instances sizes considered here are smaller than the largest ones generated in the experimental evaluation in [7].

Experiment 1

In the first experiment we increased the number of residents while maintaining a constant ratio of couples, hospitals and posts to residents. For various values of x(50≤x≤150) in increments of 20, 1000 randomly generated instances were created containing x residents, 0.1x couples (and hence 0.8x single residents) and 0.1x hospitals with x available posts that were randomly distributed amongst the hospitals. Each resident’s preference list contained a minimum of 3 and a maximum of 5 hospitals. Figure 2 (and indeed all the figures in this section) shows the mean time taken to find a maximum cardinality most-stable matching, the mean size of a maximum cardinality most-stable solution (in each case over both solvable and unsolvable instances), and the mean and maximum number of blocking pairs admitted by most-stable matchings.

The results show that the time taken to find an optimal solution increases with x, with the min bp hrc formulation being more difficult to solve in general than the hrc formulation. The mean size of an optimal solution increases with x for both solvable and unsolvable instances (it is around 95 % of x for x=50, decreasing to around 93 % of x for x=150, with the optimal matching size for unsolvable instances being very slightly larger than that for solvable instances). Perhaps most interestingly, the maximum number of blocking pairs was 1, with the mean at most 0.1, and the mean number of unsolvable instances being 77.

Experiment 2

In our second experiment we increased the percentage of residents involved in couples while maintaining the same numbers of residents, hospitals and posts. For various values of x(0≤x≤30) in increments of 5, 1000 randomly generated instances were created containing 100 residents, x couples (and hence 100−2x single residents) and 10 hospitals with 100 available posts that were unevenly distributed amongst the hospitals. Each resident’s preference list contained a minimum of 3 and a maximum of 5 hospitals. The results for all values of x are displayed in Fig. 3.

The results show that the time taken to find an optimal solution increases with x; again the min bp hrc formulation is more difficult to solve in general than the hrc formulation. The mean size of an optimal solution decreases with x for both solvable and unsolvable instances; again the optimal matching size for unsolvable instances is slightly larger than that for solvable instances. As for Experiment 1, the maximum number of blocking pairs was 1, with the number of unsolvable instances increasing from 50 for x=5 to 224 for x=30.

Experiment 3

In our third experiment we increased the number of hospitals in the instance while maintaining the same numbers of residents, couples and posts. For various values of x(10≤x≤100) in increments of 10, 1000 randomly generated instances were created containing 100 residents, 10 couples (and hence 80 single residents) and x hospitals with 100 available posts that were unevenly distributed amongst the hospitals. Each resident’s preference list contained a minimum of 3 and a maximum of 5 hospitals. The results for all values of x are displayed in Fig. 4.

The results show that the time taken to find an optimal solution decreases with x; again the min bp hrc model solution time is slower than that for the hrc model. Clearly the problem is becoming less constrained as the number of hospitals increases. Also the mean size of an optimal solution decreases with x for both solvable and unsolvable instances; again the optimal matching size for unsolvable instances is slightly larger than that for solvable instances. This time the maximum number of blocking pairs was 2, with the mean number of blocking pairs decreasing from 0.08 for x=20 to 0.04 for x=100.

Experiment 4

In our last experiment, we increased the length of the individual preference lists for the residents in the instance while maintaining the same numbers of residents, couples, hospitals and posts. For various values of x(2≤x≤6), 1000 randomly generated instances were created containing 100 residents, 10 couples (and hence 80 single residents) and 10 hospitals with 100 available posts that were unevenly distributed amongst the hospitals. Each resident’s preference list contained exactly x hospitals. The results for all values of x are displayed in Fig. 5.

The results show that increasing the preference list length makes the problem harder to solve; again the min bp hrc model is slower to solve than the hrc model. Also the mean size of an optimal solution increases with x for both solvable and unsolvable instances as more options become available in the preference lists (from 86.4 for x=2 to 97.5 for x=6 in the case of unsolvable instances); again the optimal matching size for unsolvable instances is slightly larger than that for solvable instances. The maximum number of blocking pairs was 1, with the mean at most 0.1, and the mean number of unsolvable instances being 81.

Discussion

The results presented in this section suggest that, even as we increase the number of residents or hospitals, the percentage of residents involved in a couple or the length of the residents’ preference lists, the number of blocking pairs admitted by a most-stable matching is very low. For most of the 28,000 instances generated in our experimental evaluation, the most-stable matchings found admitted at most 1 blocking pair, and the maximum number of blocking pairs admitted by any most-stable matching was never more than 2. These findings are essentially consistent with the results of Nguyen and Vohra [26], who showed that an unsolvable hrc instance only requires a small amount of perturbation in order to become solvable. Further empirical investigation is required to determine whether this behaviour is replicated for larger hrc instance sizes.

5 A constraint programming model for min bp hrc

In addition to the IP model, we designed a Constraint Programming model for min bp hrc and implemented this using the MiniZinc constraint modelling language.

We assume that residents’ preference lists are given by integer variables rpref [i][j], which play a similar role to pref (ri,j) in the IP model, and that hospitals’ ranking arrays are given by integer variables hrank [h,i], which are analogous to rank (hj,ri) in the IP model. The lengths of the preference lists of a resident ri and a hospital hj are given by rpref_len [i] and hpref_len [j] respectively. The capacity of a hospital hj is given by hosp _cap[j].

For each single resident ri, the model includes an integer variable single_pos [i] with domain (1,…,l(ri)+1), where l(ri) is the value of rpref_len [i], which takes the value j if ri is assigned her jth-choice hospital, or l(ri)+1 if ri is unassigned. For each couple i, we include an integer variable coup _pos[i] with a similar interpretation.

Each single resident’s single_pos [i] variable is channelled to an array of l(ri) boolean variables single_assigned [i], such that single_assigned [i][j]= true if and only if single_pos [i] = j, and a variable single_unassigned [i], such that single_unassigned [i]= true if and only if single_pos [i] = l(ri)+1. Similarly, we have boolean coup_assigned and coup_unassigned variables for each couple.

For each hospital i, and each position j on hospital i’s preference list, we have a boolean variable hosp_assigned [i][j] which is true if and only if hospital i is assigned its jth-choice resident. We include a constraint to ensure that hosp_assigned [i][j]= true if and only if a corresponding single_assigned or coup_assigned variable is also true. Furthermore, each hospital has a linear inequality constraint to ensure that its capacity is not exceeded.

For each position on the preference list of a single resident or couple, we create a boolean variable single_bp [i][j] or coup_bp [i][j] indicating whether the resident or couple, along with their jth-choice hospital, constitutes a blocking pair. For each type of blocking pair, we define a set of constraints and then give some brief intuition.

Type 1 blocking pairs

The hosp_would_prefer predicate for a hospital h and a position q on the preference list of h takes the value true if and only if h has fewer than hosp_cap [h] assigned residents in positions strictly preferable to position q on its preference list. (Note the redundancy in this predicate: all we actually need is the first sum(… )(… )<hosp_cap [h] constraint; the sum(… )(… )>0 constraint improves propagation.)

The constraint for Type 1 blocking pairs thus sets single_bp [i,j] to true if and only if ri is unassigned or prefers h to his partner, and h is undersubscribed or prefers ri at least one of its assignees, where h = rpref [i,j].

Type 2a/b blocking pairs

The hosp_would_prefer_exc_partner predicate on inputs h1, h2, q1, q2 (where h1, h2 are hospitals and q1, q2 are positions on their preference lists respectively) takes the value true if and only if (a) h1 = h2, q1<q2 and the number of h1’s assignees that it prefers to its q1th choice is less than hosp _cap[h1]−1, or (b) h1≠h2 or q1>q2 and the number of h1’s assignees that it prefers to its q1th choice is less than hosp _cap[h1].

Type 3a blocking pairs

The constraint for Type 3a blocking pairs thus sets coup_bp [i,j] to true if and only if couple (r1,r2) are unassigned or prefer (h1,h2) to their joint assignment, whilst for each k∈{1,2}, hk is undersubscribed or prefers rk to at least one of its assignees, where (r1,r2) is the ith couple and (h1,h2) is the hospital pair at position j of their joint list.

Type 3b/c/d blocking pairs

The hosp_would_prefer2 predicate for a hospital h and a position q on the preference list of h takes the value true if and only if h has fewer than hosp_cap [h]−1 assigned residents in positions strictly preferable to position q on its preference list. (Note the redundancy in this predicate: all we actually need is the first sum(… )(… )<hosp_cap [h]−1 constraint; the sum(… )(… )>1 constraint improves propagation.)

The constraint for Type 3b/c/d blocking pairs thus sets coup_bp [i,j] to true if and only if couple (r1,r2) are unassigned or prefer (h,h) to their joint assignment, whilst h either has two free posts (Type 3b), or h has one free post and prefers one of r1 or r2 to at least one of its assignees (Type 3c), or h is full and and prefers r1 to some assignee rk, and prefers r2 to at least one of its assignees apart from rk (Type 3d), where (r1,r2) is the ith couple and (h,h) is the hospital pair at position j of their joint list.

Experiments

The CP model was solved using the lazy clause solver Chuffed [41] on the same machine that was used for the experiments on the IP model as reported in Section 4. All instances were allowed to run to completion. We present results on the runtime of the CP model both with and without presolving. The presolve step, when included, specifies in advance which set S of resident-hospital pairs will block the solution (in practice we try out values of k=0,1,2,… and generate all subsets S of size k until we reach a feasible solution) and then performs preference list deletions in the knowledge that the pairs in S will block. This allows large reductions in the model size, and works well because the number of blocking pairs admitted by a most-stable matching is generally very small, as we saw in Section 4. We did not use presolve with the IP model, but we note that it may be possible to solve the IP model more quickly by carrying out this step.

Figure 6 plots the mean run times for each of the four experiments for the IP model and for the CP models with and without presolving: each plot in the top row shows results for the solvable instances in one experiment, and each plot in the bottom row shows corresponding results for the unsolvable instances. Table 1 shows the actual mean and median runtimes for each model, taken over all 28,000 instances \(\mathcal I\) across all four experiments, those instances from \(\mathcal I\) that were solvable and those from \(\mathcal I\) that were unsolvable.

Comparison of run times using CP (with and without presolve) and MIP models

Table 1

Summary of mean and median runtimes over all experiments (all timings are in seconds)

Instance type

Mean

Median

IP model

CP model

IP model

CP model

No presolve

Presolve

No presolve

Presolve

All

2.568

2.237

0.315

0.031

0.430

0.129

Solvable

0.034

1.839

0.143

0.016

0.402

0.127

Unsolvable

30.781

6.669

1.276

8.948

3.240

1.395

The CP model without presolve generally performs unfavourably for solvable instances. Here, the IP model is faster than the CP model with presolve; this is likely to be due to the fact that for such instances, the earlier IP model for hrc [7] is used instead of the more complex IP formulation for min bp hrc. For unsolvable instances, the CP model (with or without presolve) is faster than the IP model. This is likely to be due to the fact that the CP model for min bp hrc is more compact than its IP counterpart, involving fewer variables and constraints. Comparing total run time summed across all 28,000 instances, the CP model was 1.15 times faster than the CP model without presolve, and the CP model with presolve was 8.14 times faster than the IP model.

When solving the CP model, the distribution of runtimes for the case without presolve had a very long right tail; 14 of the 28,000 instances accounted for over half of the total run time. The longest-running instance took 17,617 seconds, and surprisingly this was a solvable instance (generated for Experiment 2). For this reason, Table 1 shows median run times as well as mean run times; from this we can see that the median runtime for the IP model is lower than that for the CP models for all instances and for solvable instances. However for unsolvable instances, the median runtime for CP without presolve is 2.762 times faster than the median runtime for IP, and this factor increases to 6.414 for CP with presolve.

6 Concluding remarks

In this paper we have presented complexity and approximability results for min bp hrc, showing that the problem is NP-hard and very difficult to approximate even in highly restricted cases. We have then presented IP and CP models, together with empirical analyses of both models applied to randomly-generated hrc instances. Our main finding is that most-stable matchings admit a very small number of blocking pairs (in most cases at most 1, but never more than 2) on the instances we generated. We also showed that on average the CP model is faster than the IP model, with the performance of the CP model being enhanced if presolving was carried out. As far as future work is concerned, it would be interesting to determine the effect of presolving on the IP model, and more generally, to investigate further methods to enable the models to be scaled up to larger instances, such as column generation in the case of the IP model, and variable / value ordering heuristics in the case of the CP model.

Footnotes

Notes

Acknowledgments

David Manlove and James Trimble were supported by Engineering and Physical Sciences Research Council grants EP/K010042/1 and EP/N508792/1. We would like to thank anonymous reviewers of earlier versions of this paper for their valuable comments, including the suggestion of references [28, 34, 35, 36], which have helped to improve the presentation of this paper.

Hinder, O. (2015). The stable matching linear program and an approximate rural hospital theorem with couples. In Proceedings of WINE '15: the 11th Conference on Web and Internet Economics, volume 9470 of Lecture Notes in Computer Science (pp. 433). Springer. Full version available from http://stanford.edu/ohinder/stability-and-lp/working-paper.pdf.

19.

Irving, R.W. (1998). Matching medical students to pairs of hospitals: a new variation on a well-known theme. In Proceedings of ESA ’98: the 6th Annual European Symposium on Algorithms, volume 1461 of Lecture Notes in Computer Science (pp. 381–392). Springer.Google Scholar

Copyright information

Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

This article is published under an open access license.
Please check the 'Copyright Information' section for details of this license and
what re-use is permitted.
If your intended use exceeds what is permitted by the license or if
you are unable to locate the licence and re-use information,
please contact the Rights and Permissions team.