Anytime Capacity Expansion in Medical Residency Match by Monte Carlo Tree Search

1 Introduction

This paper considers the capacity expansion problem in twosided matchings, where the policymaker is allowed to allocate some extra seats as well as the standard seats. The two-sided matchings have a lot of real applications such as medical residency match [Roth, 1991] and school choice [Abdulkadiro˘glu and S¨onmez, 2003]. The theory has been extensively developed across computer science and economics [Roth and Sotomayor, 1990; Manlove, 2013]. To establish our model and concepts, we use medical residency match as a running example. In this literature, each hospital (school) accepts a limited number of residents (students). Such capacity constraints (or quotas) are assumed to be known and fixed in advance.

However, in practice, even for hospitals, or its stakeholders, it is often uncertain how the capacity constraints are specified beforehand. Those could be flexible and variable more than the previous studies preclude. By pooling some extra funding, hospitals can accept a few more residents. Schools may be granted some budgets or dispatched teachers from their district according to the demand.

∗Atsushi Iwasaki was supported by JSPS KAKENHI Grant Numbers JP21H04890 and JP20K20752. A full version of this paper is found in https://arxiv.org/abs/2202.06570.

Such a situation involves a lot of resource allocation problems where capacities are predefined, e.g., how many jobs each machine takes in the jobshop scheduling problem. How we should choose the capacities has been paid less attention than it deserves, though they influence the performance generated by allocations.

A classical comparative statistic already analyzes the effect of such expansion. If the capacity of some hospitals is expanded, the welfare of every resident weakly improves [Roth and Sotomayor, 1990]. However, it is rarely investigated how capacities should be expanded to improve the total welfare. Bobbio et al. [2021] have initiated the question of how limited extra seats should be allocated, keeping resulting matchings stable. The stability is a key concept for the two-sided matching and if a matching is stable, no groups of residents and hospitals have profitable deviations. The celebrated deferred acceptance (DA) algorithm is known to find a stable matching [Gale and Shapley, 1962]. The capacity expansion with DA, i.e., finding the optimal allocation of extra seats in the welfare of residents, is formulated as an integer quadratic programming [Bobbio et al., 2021], which is computationally challenging to solve exactly. In addition to the exact method, they also developed some greedy heuristics, which run very effectively yet are suboptimal.

This paper proposes an alternative method to solve the capacity expansion problem where the upper confidence tree (UCT) searches the space of capacity expansions. 1 Each pattern of them has a resident-optimal stable assignment that DA finds. Not only does UCT obtain an optimal solution given a sufficiently large time, but also a policymaker can stop UCT anytime to obtain a reasonably good solution, which is important given the hardness of obtaining an optimal solution. We then characterize a good tree representation so that the tree search method efficiently exploits the structure of the capacity expansion problem. There have been a certain amount of studies on two-sided matchings in the AI community, although this literature has been established mainly in fields across algorithms and economics [Kamada and Kojima, 2015; Bir´o et al., 2010; Fragiadakis et al., 2016; Goto et al., 2014; Goto et al., 2016]. In many application domains, various distributional constraints are imposed on an outcome, e.g., regional maximum quotas are imposed on hospitals in urban areas to allocate more doctors to rural areas [Kamada and Kojima, 2015].

1 An implementation of our method is available at https://github. com/CyberAgentAILab/mcts-capacity-expansion.

2 Problem Setup

Let D = {d1, d2, . . . , dD} be the set of residents (doctors) and H = {h1, h2, . . . , hH} be the set of hospitals. We consider a bipartite graph where D and H are connected by edges E. Each edge represents that the resident is accepted by the hospital. Residents and hospitals have their own preferences. We denote h ≻d h′ if resident d prefers hospital h to hospital h′ . Similarly, we denote d ≻h d′ if hospital h prefers resident d over resident d′ . A vector q = (q1, q2, . . . , qH) represents the quota of the hospitals. An instance of the residency match problem is a tuple Γ = (D, H, ≻, q), where ≻ denotes all preferences of the residents and the hospitals.

A matching M ∈ M is a subset of E such that the following criteria are met. Each resident d ∈ D has at most one edge. Each hospital h ∈ H can have at most qh edges. We say that resident d is unassigned if there is no edge from d. We denote M(d) and M(h) to represent the hospital assigned to resident d and the set of the residents assigned to hospital h, respectively. Given a matching M, we say a pair (d, h) ∈ E is a blocking pair if (1) resident d is unassigned or prefers hospital h to M(d), and (2) hospital h is such that |M(h)| < qh or prefers resident d over at least one resident in M(h). Matching M is stable if there is no blocking pair.

Gale and Shapley [1962] showed that a stable matching always exists and proposed the deferred-acceptance (DA) algorithm that yields a stable matching. The resident-proposing DA finds the resident-optimal stable matching where no resident is better off among all stable matchings [Roth and Sotomayor, 1990]. That matching is also obtained by minimizing the total resident ranking under the stability constraints. A resident ranking is denoted by rankd(h) the rank of hospital h according to the resident d's preference. This is also referred to as the cost of the match (d, h). Given a matching M, the total resident ranking is defined as

DA gives a matching that minimizes Eq. (1) among all stable matchings [Bobbio et al., 2021].

This paper improves the resident utilities by allowing additional seats. Letting t ∈ NH be a non-negative vector, we consider an expanded matching problem in which the capacity of each hospital h is the sum of regular capacity qh and the extended capacity th. A reasonable extension of the capacity should not be very large and can accommodate the demand of each hospital. Given this, we consider the following optimization problem: Let

be the set of matchings. The capacity expansion problem is the following:

where Sdh = {h′ ∈ H : rankd(h′ ) ≤ rankd(h)} is the set of indices of hopsitals that resident d prefers at least as hospital h; similarly, Tdh = {d′ ∈ D : rankh(d′ ) < rankh(d)} is the set of indices of residents that hospital h prefers more than resident d. Namely, keeping stability, we maximize the residents' welfare subject to the capacity constraints of the hospitals, where the capacity can be relaxed up to B seats. Moreover, we assume that each hospital h has its own expansion limit bh.

Remark 1. (Complete preference) The optimization of Eq. (2) requires a complete preference of all residents and hospitals, which is often impractical. For example, in the real data of a residency matching program, an average resident only applies to four or five hospitals. In this case, we can add a dummy hospital of infinite capacity and rank all unpreferred hospitals after the dummy hospital.

Bobbio et al. [2021] proved the computational hardness of the problem of Eq. (2) without hospital-wise limit bh.

Proposition 1. (Hardness result [Bobbio et al., 2021]) The decision version of Eq. (2) is NP-complete even in the case where there is no hospital-wise limit bh.

The key idea of our algorithm for this problem is that, when we fix an expansion vector t ∈ Θ, the capacity expansion problem of Eq. (2) boils down to the standard matching instance Γ = (D, H, ≻, q + t) that we can solve efficiently by DA. 2 Therefore, a tree search on the space of Θ combined with DA is expected to give an efficient solution.

3 Method

This section proposes our algorithm to solve Eq. (2). We apply Upper Confidence Tree (UCT) [Kocsis and Szepesv´ari, 2006], which belongs to a class of Monte Carlo tree search methods [Browne et al., 2012]. One of the most successful examples of UCT is for abstract games such as the Game of Go [Wang and Gelly, 2007; Yoshizoe et al., 2011; Silver et al., 2016] where the goal is to find the best next move in a game tree. Another aspect of UCT is diverse recommendation method, which minimizes total regret, searches for different solutions of reasonable objective value [Bosc et al., 2018]. In this paper, we use UCT to find a global optimum node in a tree Tall.

2 The computational complexity of DA is O((D + H)2 ).

UCT consists of a sequential process of traversing a tree structure where each tree node i is associated with its value vi. We use n ∈ {1, 2, . . . ,N} to denote the number of rounds. The entire tree Tall is typically very large, and UCT tries to develop a subtree of it. We denote the subtree of round n by T (n). Each node i in the current tree is associated with a tuple (Vi,Ni) ∈ R2 , where Vi is the sum of rewards and Ni the number of times at which node i is traversed. The value ˆµi := Vi/Ni is an estimator of vi, and its standard deviation is proportional to √1/Ni. Each round consists of selection, development, simulation, and backpropagation steps.

Selection: At each iteration, UCT traverses the tree from the root. At each node i, it chooses a child node c with the maximum UCB value:

where Cp > 0 is the parameter that determines how much exploration it attempts.

Development: When it reaches a node k that is out of the current tree, it adds the node to the current tree.

Simulation: Based on the reached node k, it conducts a random play to find l, which is a leaf that stems from k.

Backpropagation: Letting vl be the value of the leaf l, it updates the statistics of all nodes that we have traversed in this round: Vi ← Vi + vl,Ni ← Ni + 1 for each node i during the backpropagation.

The four steps above are quite standard in UCT. For a more algorithmic description, see the full version.

3.1 Tree Representation of Expansion Space

Unlike abstract games where the game tree is inherent, the relation between a tree and the capacity expansion in our problem is nontrivial. This section describes the mapping from a tree node into an expansion t ∈ Θ.

In the capacity expansion problem defined in Eq. (2), DA yields an optimal solution given a fixed expansion vector t. Therefore, we consider a tree where each node corresponds to an expansion vector t, and each leaf is an expansion of limit ∑h∈H th = B. The value vi at each node i is the objective value of DA for the corresponding expansion.

We consider the following criteria are important in constructing a good tree representation of an expansion.

1. Faithfulness: The tree preserves the inclusion relationship of the original expansion vector space: if node j is a descendant of node i, then t(j) - t(i) ≥ 0, where t(i), t(j) be corresponding expansions of nodes i, j and t(j) - t(i) ≥ 0 means all the features of the vector t(j) - t(i) are non-negative.
2. Nonredundancy: The mapping is one-to-one. Namely, for any expansion vector t ∈ Θ such that ∑h th = B, there exists only one leaf of tree Tall. Redundancy in a tree representation compromises the search efficiency.
3. Decisiveness: Assuming that a tree representation is faithful, the tree representation is decisive if branching in a shallow layer is more informative than the one in a deep layer. In our case, allocations related to important hospitals should appear in a shallow layer of the tree.

Considering those criteria, we propose the three tree representations in Figure 1. In all of the three representations, the root node r corresponds to the zero vector t = (0, 0, . . . , 0) that corresponds to no expansion. Table 1 illustrates which properties are satisfied in each of the representations. We will discuss the idea behind them.

The iterative representation, which is the most straightforward, builds an H-ary tree. 3 Each edge corresponds to allocating a seat to one of the hospitals. As a result, each path from the root to depth B represents an expansion of size B. Although the iterative representation is faithful, it is redundant. For example, consider the case of three hospitals. Let 1 → 2 → 2 to denote a path on the tree that sequentially expands hospitals 1, 2, 2. An expansion vector (2, 1, 0) corresponds to three different path on the tree (1 → 1 → 2, 1 → 2 → 1, and 2 → 1 → 1). As a result, UCT with iterative tree representation searches an unnecessary large representation space, which increases the computational burden.

To deal with this issue, we introduce an improved tree representation of the expansion, which we call the iterative priority tree (IPT). For each node i, it only allocates a node with its priority higher than the most prioritized hospital that has been allocated a seat. IPT solves the issue of redunduncy in the iterative tree:

Proposition 2. (Nonredundancy of IPT) For each expansion vector t ∈ NH, there exists a unique node in IPT.

Proof. Each expansion uniquely corresponds to a nonincreasing sequence of numbers. For example, expansion t = (1, 3, 1) corresponds to 3 → 2 → 2 → 2 → 1. It is easy to confirm that, in the iterative priority tree, this is the unique path that leads to the desired expansion.

While the IPT is nonredundant, there still remains some space for improvement. The batch tree (BT) representation fully exploits the priority information to allocate more than one seats to popular nodes, which reduces the depth of the optimal solution in the tree. Namely, each depth of BT corresponds to how many seats to allocate to each hospital. The following proposition states that the properties of BT:

Proposition 3. (Nonredundancy and decisiveness of BT) BT is nonredundant. It is also decisive, that is, if the optimal solution aligns with the expansion vector (i.e., t1 ≥ t2 ≥ t3 ≥ · · · ≥ tH), then BT minimizes the depth of the optimal node among all trees such that each edge allocates seats to one hospital.

Proof. The nonredundancy is trivial. In the following, we derive the decisiveness of the BT. Assume that the optimal expansion vector is such that t1 ≥ t2 ≥ t3 ≥ · · · ≥ tH, and let h be the first index such that th = 0. BT includes the path of depth h - 1 that leads to this allocation.

3.2 Ordering Hospitals

We introduce two ideas for ordering hospitals that we use in IPT and BT. The first idea is to order hospitals in terms of their popularity

3 Note that H is the number of the hospitals.

We use the term "popularity" because it is the total rank of hospital h in view of residents. Note that this value is also referred to as the Borda count in the context of social choice [Moulin, 1994]. A smaller value represents a more popular hospital.

The second idea is utilizing potential envies that may not be justified so that we could construct a better decision tree of the hospitals. Namely, let

be the potential envy that the residents have toward the ones matched to hospital h in matching M. The score indicates the number of residents who prefer hospital h over their matched hospital in M. For calculating M, we run DA with no expansion before running UCT to determine the ordering.

Avoiding multiple entries. Since the evaluation of each node (i.e., DA with a given expansion vector) is deterministic, we can restrict the node selection so that a node is never evaluated twice. In the backpropagation step, we mark the node k as "evaluated" if it is a leaf of Tall. Moreover, for each ascendant node, if all children of the node are evaluated, then we mark that node as evaluated. In the selection and simulation phase, the evaluated nodes are never visited again.

3.3 Consistency of the Method

This section analyzes the proposed method. The first result is a trivial consequence of avoiding multiple entries:

Proposition 4. (Worst-case sample complexity) UCT with iterative tree finds an optimal solution in N = |Tall| rounds.

Proof. Since the development step adds a node to T (n) at each round n, it follows that it obtains the optimal solution if N reaches the number of the nodes |Tall|.

We further propose a nontrivial analysis by following the analysis of the original UCT by [Kocsis and Szepesv´ari, 2006]. The largest difference is that Kocsis and Szepesv´ari [2006] analyzed the average quality of the selected action (i.e., regret), whereas we are interested in the quality of the best path, which corresponds to the optimal solution of the capacity expansion problem.

Let ˆµi,m = Vi/m be the mean of vi over the first m visits, and Ui,m,n = ˆµi,m + Cp√ log n m be the corresponding UCB value. Let W be the maximum number of the children of a node and DT be the depth of the tree. 4 Let v∗ i be the maximum value of the leaves that are descendants of node i. Let v∗ = v∗ r be the global optimal value, where r is the root node. The following assumption is essentially the same as the one in [Kocsis and Szepesv´ari, 2006]:

4 That is W = H,DT = B for iterative and iterative priority tree, or W = B, DT = H for the batch tree.

Assumption 1. (Confidence bound) We assume

Definition 1. (∆-optimal tree) Let ∆ > 0 be arbitrary. The ∆-optimal subtree is defined recursively as follows. First, it includes the root of the original tree. For each node i in the subtree, add each children c such that v∗ i - v∗ c ≤ ∆. In other words, ∆-optimal subtree is a subtree where each edge is suboptimal at most ∆. Let S(∆) be the number of the nodes in the ∆-optimal subtree.

Theorem 5. If Assumption 1 holds for any node i, then with probability at least 1 - o(N -1 ), UCT finds at least one node in the ∆-suboptimal tree if N satisfies

The proof of Theorem 5 is in the full version.

Remark 2. (Implication) The solution of the node i in ∆- Second, we generate the dataset, JRMP, based on Japan Residency Matching Program 2007 [Kamada and Kojima, 2015], which matches medical hospital students (residents) with residency training programs. We extracted 1, 287 residents in the Tokyo district who match with 50 + 1 (1 for a dummy) hospitals with a resident-side preference. We set B ∈ {10, 30} for the limit of capacity expansion. For each B, we average the results for 10 instances.

Tables 2 and 3 illustrate the quality of solutions and its runtime for Set 1, respectively. The results for Set 2, which is omitted due to page limitation, are found in the full version. Also, we place Tables 4 and 5 for the JRMP data. Note that Table 2 indicates the average percentage gap in the total resident ranking (TRR) defined in Eq. (1):

Agg-Lin always outputs an optimal solution upon the completion. However, as in Table 3, it does not run within one hour for large instances with positive correlation α. In that case, we use the temporal value that Gurobi outputs, which can be suboptimal and the percentage gap can be negative.

Among the algorithms we consider, the two greedy methods (Grdy and LPH) run very fast (< 10 seconds) and output a suboptimal solution. Both of them are outperformed by all variants of ours. In particular, BT outperformed iterative tree and IPT, and using the popularity and envy-based orderings outperformed the random one.

Figure 2 describes the objective value as a function of rounds, which indicates an early termination of BT often has a satisfying solution. In summary, BT with envy/popularitybased ordering yields the best results among our algorithms. Our algorithms run 2-20 times faster than Agg-Lin, and the quality of the solution is close to optimal.

5 Conclusion

This paper sheds light on the capacity expansion in the twosided matching to consider a flexible allocation of extra seats within a given limit beforehand. To handle this NP-Complete problem, we develop a UCT-based search method and verify that it outperforms the previous approaches. Future works include (1) extending it to matchings with constraints that need not admit stable matching and (2) utilizing other tree search algorithms such as Nested Monte-Carlo Search [Cazenave, 2009; Rosin, 2011].

suboptimal tree is at least

Therefore, letting ∆ be sufficiently small and N be sufficiently large, UCT finds an almost optimal solution. The required number of N depends on S(∆), which corresponds to the number of close-to-optimal children.

4 Evaluation

This section empirically evaluates our algorithm via synthetic and real datasets. We compare our algorithms with the ones by [Bobbio et al., 2021]: the greedy algorithm (Grdy), the linear programming-based heuristic (LPH), and the aggregated linearization (Agg-Lin).

Let us enumerate existing algorithms: First, Grdy is an algorithm that allocates B expansion seats iteratively such that each allocation maximizes the marginal cost reduction. Second, LPH first computes a minimum cost matching without stability constraints via minimum-cost flow to fix the expansion seats, and then uses DA to obtain a stable matching. Finally, Agg-Lin is an exact method that uses McCormick envelopes and solves a mixed integer programming, which is computationally intensive. We warm-up Agg-Lin as suggested by [Bobbio et al., 2021].

For our UCT method, we set N = B×103 in synthetic data experiments, N = B × 102 in real data experiments. The value of Cp, which determines the tradeoff between exploration and exploitation, is set to be √0.002. We consider this parameter is robust enough to cover all settings. We implement our simulation with Python 3 and solve the mathematical programming using Gurobi, which is in favor of Agg-Lin. We restrict the runtime for each method to one hour.

Let us describe the procedure of generating the datasets. First, we build two synthetic datasets. Set 1 involves B but each hospital does not have its expansion limit bh. Set 2 involves B as well as bh for each hospital. Note that the setting of Set 1 is similar to the experimental section in [Bobbio et al., 2021]. Regarding the preference among the residents, we follow the setting [Goto et al., 2016], which involves a correlation parameter α ≥ 0. Larger value of α implies a stronger correlation among the preferences. Note that α = 0 (no correlation) corresponds to the setting of [Bobbio et al., 2021]. We set D := |D| = 1, 000, and conduct experiments varying the parameters H := |H|, B, bh, and α. For each combination of parameters, we average the results for 10 instances.