Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
69,921 characters · 14 sections · 42 citation commands
Who With Whom? Learning Optimal Matching Policies
Keywords: Policy learning, Two-sided matching, Regularized optimal transport, Regret bound.
Active job market programs commonly involve caseworkers providing advice and support to unemployed workers during their job search process. Empirical research such as dromundo2022 indicates that there is substantial heterogeneity in both job seekers and caseworkers, and the overall performance of these programs is driven by who matches with whom. This suggests an avenue for policy intervention. Instead of allocating caseworkers randomly, as is currently common, a planner can centrally assign caseworkers in a sophisticated manner, taking into account the complementarity or substitutability of the characteristics of caseworkers and job seekers. A similar opportunity for policy intervention exists in many other contexts including matching medical doctors and patients dahlstrand2024 or hospitals mourot2025, teachers and classrooms graham2023, attorneys and defendants spurr1987, and tax auditors and taxpayers bergeron2025 among others.
As emphasized by graham2014, graham2020, reallocating individuals through a change in matching procedure can be less costly than other policy interventions, such as training or hiring personnel or introducing a new program. Methods for learning optimal matching policies therefore have the potential to have a large real-world welfare impact in a wide scope of applications. However, methods for efficiently learning an optimal matching policy from data and their statistical performance are not well studied.
This paper develops a method to learn welfare optimal matching policies for two-sided matching problems in which a planner decides who should match with whom based on the observable characteristics of the two sides. We formulate the learning problem as an empirical optimal transport problem with a match cost function estimated from training data. In this formulation, we ultimately uncover the joint distribution which minimizes the average match cost and has marginal distributions that coincide with the empirical distributions of observable characteristics for the two sides. In the standard optimal transport formulation, an optimal matching policy can be obtained by linear programming. The dimension of choice variables in this linear program is $n^2$, where $n$ is the number of pairs created in the matching process.
To keep the estimation of an optimal matching policy feasible even for large $n$, we propose maximizing the entropy regularized transport cost criterion with the Sinkhorn algorithm (iterated projection fitting). This computationally attractive procedure is proposed in Cuturi2013 and Galichon_Salanie_2022, where its analytical properties are also studied. In the case of matching caseworkers and job seekers, the number of matches to be formed in a given month within each unemployment agency ranges from a few hundred to a few thousand. The corresponding unregularized optimal transportation problem is a linear program with 100,000 choice variables, which is computationally demanding. Solving the regularized problem with the Sinkhorn algorithm, on the other hand, involves iterations of arithmetic operations, which drastically reduces the complexity of computation.
A novel contribution of this article is to derive a welfare regret bound for the estimated matching policy and characterize its convergence, taking into account that the cost function is estimated. In contrast to the existing work on empirical optimal transport, e.g., Rigollet2022, in which the cost function is assumed to be known, our welfare regret bound accounts for estimation errors and sampling uncertainty. Our bound is non-asymptotic, and we derive it under a minimal regularity condition. It exhibits a bias-variance tradeoff with respect to a regularization parameter, which can guide the choice of a value for the regularization parameter.
We perform an extensive numerical study to assess the performance of our proposed method. A first simulation exercise with a simple data generating process illustrates how the welfare gains associated with our method depend on (i) the size of the training sample used to estimate (nonparametrically) the cost function and (ii) the choice of regularization parameter for the regularized optimal transport problem. In addition, to gauge the magnitude of the welfare gain of our proposal in a real-world scenario, we apply our method to the assignment of caseworkers to job seekers in a job search assistance program. We calibrate a data generating processes using estimates of key parameters reported in dromundo2022, and vary (i) the size of the training sample used for cost estimation, (ii) the regularization parameter, and (iii) the strength of complementarities in the job finding rate function. Our results indicate that under reasonable levels of complementarities, and for realistic training sample sizes given the administrative data available, our method could improve (at virtually no cost) job finding rates by about 1 percentage point---an order of magnitude similar to very costly job search counseling interventions (e.g., behaghel2014).
The statistical treatment choice literature initiated by Manski2004 and Dehejia2005 has been a growing area of econometric research. Given a planner's welfare objective, there are multiple approaches to learning treatment choice rules from data. These include solving the statistical decision problem in finite samples \citep*{stoye2009minimax,stoye2012minimax, tetenov2012statistical, ishihara2021, yata2021, montielolea2023decision, kitagawa2022treatment}, solving for asymptotically optimal rules in limit experiments \citep*{HiranoPorter2009, christensen2025optimaldecisionrulespayoffs, masten2023minimax}, and performing empirical welfare maximization \citep*{kitagawa2018,KT19,Athey2021,MT17,sun2021empirical,Viviano21,KST21,Sakaguchi21,KLQ2025regret_aversion}. These works generally assume that a planner allocates a finite number of treatment arms over infinitely many individuals belonging to a superpopulation.\footnote{Kallus_Zhou_2018and Ai_Fang_Xie2024 study policy learning for the allocation of a continuous treatment.} In statistical policy choice for two-sided matching problems, in contrast, the number of treatment arms (e.g., caseworkers) is as large as the number of individuals (e.g., job seekers) and each treatment arm has a tight capacity constraint (e.g., each caseworker can assist only a certain number of job seekers per month). In this setting, how can an optimal matching policy be learned from the data and its welfare regret assessed? This paper exploits recent advances in the empirical optimal transport literature to answer these questions.
The statistical properties of empirical optimal transport have been extensively studied in the case where the cost function is known. See Genevay2019, Rigollet2022, and the references therein. On the other hand, the unknown cost function case is less studied. Hundrieser2023 considers estimation and inference for the value of optimal transport with an estimated cost function, and obtains a distributional approximation for asymptotically valid inference. In contrast, we focus on a non-asymptotic regret bound for the estimated matching policy applied to a finite number of individuals drawn from a superpopulation. We obtain a regret bound under the minimal regularity condition that the cost function is bounded, and this bound is valid for any size of training sample and any number of individuals to be matched.
Optimal transport methods have been applied to many economic and econometric problems. See Galichon2016 for a recent monograph on the topic. These include the identification and estimation of Hedonic models \citep*{Ekeland_etal_2004}, the identification and estimation of a matching market surplus \citep*{Dupuy_Galichon_2014,Galichon_Salanie_2022}, partial identification \citep*{Galichon_Henry_2011,lei2025}, multivariate quantile analysis \citep*{Chernozhukov_etal_2017}, measurement errors and data combination \citep*{Schennach_Starck_2022, DHaultfoeuille_etal_2025}, and causal inference \citep*{Gunsilius_2023, pouliot2025}. The class of two-sided matching problems considered in this paper corresponds to the original Monge-Kantorovich optimal transport problem. bhattacharya2009 studies estimation of optimal matching policies with an estimated cost functions and discrete characteristics. graham2011 and graham2014 model the two-sided matching problem as choosing a copula between distributions of single indices aggregating multi-dimensional characteristics of each side, and perform estimation and inference for the optimal index coefficients as well as copula parameters. The optimal transportation approach has several advantages over the copula approach. First, the optimal transportation approach can accommodate multi-dimensional characteristics flexibly without reducing the characteristics of each side to a single index. Second, the optimal transportation approach is computationally attractive, as the estimation of an optimal matching policy can be reduced to convex optimization.
The analytical and computational tools of optimal transport have been utilized in the treatment choice literature. adjaho2022externally and Kido2025 model the difference between the sample and target populations by their Wasserstein distance, and study how treatment choice can be made robust to failures of external validity. Generalizing the analysis of BhattacharyaDupas2012 with a binary treatment to the multi-valued treatment setting, Sunada_Izumi_2025 formulates treatment assignment under capacity constraints as an optimal transportation problem and studies the local asymptotic optimality of an estimated optimal transport map in limit experiments.
Consider a planner who matches one group of individuals $\mathcal{I} := \{1, \dots, n \}$ (e.g. job seekers) with another $\mathcal{J} := \{ 1, \dots, n\}$ (e.g. caseworkers). Matching is one-to-one and individuals are indivisible. Each individual $i \in \mathcal{I}$ is matched with only one individual $j \in \mathcal{J}$ and vice versa. Without loss of generality, we let the two sides consist of the same number of individuals, $|\mathcal{I}| = |\mathcal{J}| = n$. For example, if $|\mathcal{I}| > |\mathcal{J}|$, we can construct $\tilde{\mathcal{J}}$ such that $|\mathcal{I}| = |\tilde{\mathcal{J}}|$ by adding dummy individuals to $\mathcal{J}$. Any individual $i \in \mathcal{I}$ paired with a dummy individual in $\tilde{\mathcal{J}}$ is then an individual with no match in $\mathcal{J}$.
Before matching $\mathcal{I}$ and $\mathcal{J}$, the planner observes the characteristics of all individuals on both sides. Let $\mathbf{X}^n \equiv (X_i: i \in \mathcal{I})$, $X_i \in \mathcal{X} \subset \mathbb{R}^{d_x}$, be the characteristics of the $n$ individuals in $\mathcal{I}$. For instance, if $\mathcal{I}$ is a pool of job seekers, $X_i$ could include $i$'s education level, previous earnings, work experience, an index of the risk of long-term unemployment, etc. We denote the empirical distribution on $\mathcal{X}$ constructed upon $\mathbf{X}^n$ by $\mu_n$. Similarly, let $\mathbf{W}^n \equiv (W_j: j \in \mathcal{J})$, $W_j \in \mathcal{W} \subset \mathbb{R}^{d_w}$ be the observable characteristics of the $n$ individuals in $\mathcal{J}$. If $\mathcal{J}$ is a pool of caseworkers providing job search assistance, $W_j$ could include caseworker $j$'s demographic characteristics and an estimated index of value-added. We denote the empirical distribution on $\mathcal{W}$ constructed upon $\mathbf{W}^n$ by $\nu_n$.
Define an assignment to be a permutation: $\sigma: \mathcal{I} \to \mathcal{J}$, yielding $n$ pairs $(X_i, W_{\sigma(i)}), \; i=1,\dots,n$. An assignment $\sigma$ corresponds to a planner's edict dictating who in $\mathcal{I}$ matches with whom in $\mathcal{J}$. Let $\pi_{\sigma}$ be a probability distribution over the set of permutations $\{ \sigma \}$.
Following the conventions of the optimal transport framework, the outcome when $i$ is matched with $j$ is the cost of the match, $Y_i(j) \in \mathbb{R_+}$. If the natural measure of the outcome of a match is output rather than cost, this can be transformed to a cost measured in terms of lost output relative to some maximum potential output. In the example of matching job seekers and caseworkers, we define $Y_i(j)$ to be an indicator that takes value $Y_i(j)=1$ if job seeker $i$ assisted by caseworker $j$ remains unemployed within 6 months of enrolling in the program and value $Y_i(j)=0$ otherwise.
We interpret the cost measure $Y_i(j), (i,j) \in \mathcal{I} \times \mathcal{J}$ as a set of potential outcomes in the sense that for any $i \in \mathcal{I}$, exogenously changing their match from $j$ to $j'$ causally shifts individual $i$'s match cost from $Y_i(j)$ to $Y_i(j')$. We assume that the potential outcomes $(Y_i(j): (i,j) \in \mathcal{I} \times \mathcal{J})$ are random variables defined on a superpopulation, and admit the following decomposition:
where $c(x,w) := E[Y_i(j)|X_i =x, W_j = w]$ is the average cost if an individual $i$ with observable characteristics $X_i =x$ is exogenously matched with an individual $j$ with observable characteristics $W_j = w$. The average cost function is anonymous in the sense that $c(\cdot, \cdot)$ depends on $i$ and $j$ only through the values of their observable characteristics $(X_i,W_j)$. The residual term $\epsilon_i(j)$ captures the effect of any unobservable characteristics of $i$ and $j$ on their match cost. By construction, $E[\epsilon_i(j)|X_i, W_j] = 0$. We assume no-interference in the following sense: for every $(i,j) \in \mathcal{I} \times \mathcal{J}$
That is, the match cost for $i$ and $j$ is statistically independent of the characteristics of all other individuals.\footnote{Another implicit no-interference assumption is embedded in the notation $Y_i(j)$ for denoting the (potential) outcome when $i$ is matched with $j$---as it states that this outcome does not depend on the assignment of other individuals.}
We define the planner's objective function to be the expected total cost of all $n$ matches given the profiles of the observable characterstics for both sides, $\mathbf{X}^n$ and $\mathbf{W}^n$. Under the no-interference assumption ((ref)), the expected total cost of assignment $\sigma$ is given by
If the planner generates an assignment $\sigma \sim \pi_{\sigma}$, the expected cost taking into account the randomness of assignment is given by
We define the planner’s optimal assignment rule to be:
where $\Pi_{\sigma}$ is the set of distributions over the permutations of $n$ elements.
Since the number of permutations is $n!$, the number of support points of $\pi_{\sigma}$ grows rapidly with $n$, and the optimization problem ((ref)) is infeasible for all but very low values of $n$. To overcome this computational challenge, we employ Kantorovich's relaxation and transform ((ref)) into the canonical Monge-Kantorovich optimal transport problem for the empirical distributions of $(\mathbf{X}^n,\mathbf{W}^n)$.
Let $\Pi(\mu_n,\nu_n)$ be the set of joint distributions on $\mathcal{X} \times \mathcal{W}$ whose marginal distribution on $\mathcal{X}$ coincides with $\mu_n$ and whose marginal distribution on $\mathcal{W}$ coincides with $\nu_n$. Notice that $\pi \in \Pi(\mu_n,\nu_n)$ corresponds to the probability masses $(\pi(X_i,W_j): (i,j) \in \mathcal{I} \times \mathcal{J})$ supported on $\mathbf{X}^n \times \mathbf{W}^n$.
Since the joint distribution on $\mathbf{X}^n \times \mathbf{W}^n$ representing permutation $\sigma$ belongs to $\Pi(\mu_n,\nu_n)$, any $\pi_{\sigma} \in \Pi_{\sigma}$ implies a unique $\pi \in \Pi(\mu_n, \nu_n)$. Conversely, by Birkhoff's theorem Birkhoff1946, for any $\pi \in \Pi(\mu_n, \nu_n)$, there exists $\pi_\sigma \in \Pi_{\sigma}$ such that the joint distribution of $(X_i,W_{\sigma(i)}), \sigma \sim \pi_{\sigma}$ follows $\pi$. This implies that the search for an optimal assignment rule ((ref)) can be reduced to the following discrete finite optimal transportation problem:
The planner's goal is to learn and implement the optimal matching rule $\pi_n^{\ast}$ defined in ((ref)). However, in many applications, the average cost function $c(x,w)$ is unknown. This presents a fundamental challenge to learning the optimal policy. We consider estimating the cost function using training data and learning a matching policy based on the estimated cost function. To reduce the sensitivity of the estimated matching rule to estimation errors, we solve a sample analogue of the optimization problem ((ref)) with a regularization term added to the objective function.
The planner has access to training data consisting of $N \geq 1$ observed matches $(X_{\ell}, W_{\ell}: \ell=1, \dots, N)$ and the realized cost for each match $Y_{\ell}$, $\ell=1, \dots, N$. We assume that matches in the training data are randomized in the sense that conditional on the observable characteristics $X_{\ell}$ of the $\ell$-th individual, their match is assigned independently of their unobservable heterogeneity, denoted by $\epsilon_{\ell}$. This corresponds to the exogeneity condition considered in graham2011, graham2014,
Under this exogeneity assumption, the cost function $c(x,w)$ is identified as the conditional expectation function of $Y_{\ell}$ given the characteristics of the matched individuals $X_{\ell}$ and $W_{\ell}$.
As a performance guarantee, we derive a welfare regret bound for an estimated matching policy which impose no assumptions on the cost function estimator. Let $\hat{c}(\cdot,\cdot)$ denote an estimator for $c(\cdot,\cdot)$. For the sake of later analysis, we define the $\mathcal{L}^1$ and $\mathcal{L}^2$ errors of $\hat{c}$ relative to the true cost function $c$ with respect to the product measure $\mu \otimes \nu$ (i.e., the independent coupling) of $\mu$ on $\mathcal{X}$ and $\nu$ on $ \mathcal{W}$:
Given the definition of the optimal matching policy ((ref)), a natural approach is to replace the true cost function $c$ with its estimated counterpart $\hat{c}$, and maximize $\pi_n(\hat{c})$ with respect to $\pi_n \in \Pi (\mu_n, \nu_n)$.
This plug-in approach faces several issues. First, since the optimization problem is a linear program, the solution will be an extreme point of the polyhedron constraint set $\Pi(\mu_n, \nu_n)$. As solutions of this type can be sensitive to perturbations in the objective function, an optimal matching policy estimated using a plug-in approach may be sensitive to estimation errors in $\hat{c}$. Second, although linear programming is a type of convex optimization, it scales poorly and this limits the scope of applications. For example, when matching $n=1000$ heterogeneous job seekers and case workers, the dimension of $\pi_n$ and the number of constraints are of the order $1000^2$, which is challenging even for modern linear programming solvers.
Entropy regularized optimal transport (ROT) with iterative projection fitting (the Sinkhorn algorithm), as proposed in Cuturi2013 and Galichon_Salanie_2022, mitigates these concerns. Given a regularization parameter $1/\eta \geq 0$ and $(\mathbf{X}^n,\mathbf{W}^n)$, an empirical version of ROT considers the following optimization problem:
where $KL(\pi_n \| \mu_n \otimes \nu_n)$ is the Kullback-Leibler divergence of $\pi_n$ from the product empirical measure $\mu_n \otimes \nu_n$,
Since $KL(\pi_n \| \mu_n \otimes \nu_n)$ increases as $\pi_n$ diverges from the independent coupling $\mu_n \otimes \nu_n$, the additional term penalizes $\pi_n$'s that concentrate at the extreme points of $\Pi(\mu_n, \nu_n)$. This pushes $\hat{\pi}_n^{ROT}$ away from deterministic policies corresponding to permutations in $\Pi_{\sigma}$ and towards stochastic matching policies. The regularization parameter $1/\eta$ governs how close $\hat{\pi}_n^{ROT}$ is to pure random matching. A smaller $\eta$ corresponds to a larger penalty, and forces $\hat{\pi}_n^{ROT}$ closer to random matching. As such, regularization makes $\hat{\pi}_n^{ROT}$ less sensitive to estimation errors in $\hat{c}$. Conversely, as $\eta \to \infty$, the regularized solution $\hat{\pi}_n^{ROT}$ converges to the solution of the non-regularized empirical optimal transport problem. See, e.g., Proposition 4.1 in Peyre2020.
A computational advantage of ROT stems from the dual form of ((ref)). Let $f: \mathbf{X}^n \to \mathbb{R}$ and $g: \mathbf{W}^n \to \mathbb{R}$ be Lagrange multipliers for the equality constraints $\pi_n \in \Pi(\mu_n,\nu_n)$, and consider the Lagrangian of the primal problem:
Reversing the order of $\min_{\pi_n}$ and $\max_{f,g}$, we take the first-order conditions in $\pi_n$ given $(f,g)$ and obtain:
Plugging this first-order condition into the Lagrangean yields the following dual form:
where $\hat{f}_n: \mathbf{X}^n \to \mathbb{R}$ and $\hat{g}_n: \mathbf{W}^n \to \mathbb{R}$ are optimal Lagrange multipliers in the dual problem, and $(\mu_n \otimes \nu_n)(\cdot)$ denotes integration with respect to $\mu_n \otimes \nu_n$.
The dual problem is an unconstrained optimization problem with a concave objective function and two choice variables, $f$ and $g$, each of dimension $n$. The overall dimension of the choice variables in the dual problem is thus $2n$, whereas the dimension of the primal problem ((ref)) is $n^2$. Note that the dual objective function is invariant to $(f+a,g-a)$ for any translation vector $a \in \mathbb{R}^n$. Hence, with a location normalization of one of the vectors of choice variables such as
the solution of the dual problem can be made unique. The main computational advantage of the formulation (ref) is that it can be cast as a matrix scaling problem to which an iterative projection algorithm called the Sinkhorn algorithm can be applied. See, Proposition 4.3 in Peyre2020.
With the solution $(\hat{f}_n,\hat{g}_n)$ in hand, we can express the optimal plan as, for each $i \in \mathcal{I}$ and $j \in \mathcal{J}$
Note that $\hat{\pi}_n^{ROT}$ itself does not give pairings between $\mathcal{I}$ and $\mathcal{J}$. The last step is applying the Birkhoff–von Neumann algorithm to transform $\hat{\pi}_n^{ROT}$ into a distribution over permutations $\hat{\pi}_{\sigma}$. Matches are then formed by drawing $\sigma \sim \hat{\pi}_\sigma$. The average cost of the resulting pairs is $\hat{\pi}_n^{ROT}(c)$.
Given $(\mu_n, \nu_n)$, we have built a recommended matching policy $\hat{\pi}_n^{ROT}$. This section assesses the statistical behavior of $\hat{\pi}_n^{ROT} (c)$ in terms of its average cost performance of with respect to the sampling distribution of the training data $(Y_{\ell}, X_{\ell}, W_{\ell}: \ell=1, \dots, N) \sim P^N$ and the characteristics of the individuals to be matched $(\mathbf{X}^n, \mathbf{W}^n)$. Specifically, we derive a non-asymptotic upper bound for the regret of the average cost of $\hat{\pi}_n^{ROT}$ and characterize its convergence rate.
For this goal, define the oracle ROT matching policy given $(\mu_n,\nu_n)$ with knowledge of the true cost function $c(x,y)$ to be
Let $\pi_n^{\ast}(c)$ be the global minimum average cost. This corresponds to the average cost attained by $\pi^{\ast}_n \in \arg \min_{\pi_n \in \Pi(\mu_n, \nu_n)} \left\{ \pi_n(c) \right\}$, the oracle optimal matching policy without regularization. Let the average cost attained by a feasible regularized matching policy be $\hat{\pi}_n^{ROT}(c)$. We can then define the expected regret of implementing policy $\hat{\pi}_n^{ROT}$ given $(\mu_n,\nu_n)$ relative to the oracle unregularized policy $\pi_n^{\ast}$ as
In what follows, we obtain a uniform upper bound for $Re(\hat{\pi}_n^{ROT})$ and discuss how its average (over samples) depends on the size of the training sample, the regularization parameter $\eta$, and the size of the matching market $n$.
Since the regularization term is nonnegative, we have the following inequality for $Re(\hat{\pi}_n^{ROT})$:
where $ Re^{ROT}(\hat{\pi}_n^{ROT})$ is regret defined in terms of the regularized objective function:
There are two sources of randomness in regret and its upper bound ((ref)).
To bound the welfare regret, consider the following steps:
Note that, depending on the application, it may not be necessary to consider the randomness of the characteristics of the sample of individuals to be matched. In this case, expected regret must account only for the randomness of the training sample, and step 3 can be omitted.
Based on the decomposition of ((ref)), we bound the regret of the regularized cost $Re^{ROT}(\hat{\pi}_n^{ROT}) $ and the regularization bias separately.
Consider bounding $Re^{ROT}(\hat{\pi}_n^{ROT})$ as follows:
where
is the dual objective function of the regularized OT problem, $(\hat{f}_n,\hat{g}_n) = \arg \sup_{f,g} \Phi_n(\hat{c},f,g)$, and $(f_n,g_n) = \arg \sup_{f,g} \Phi_n(c,f,g)$. Note that the inequality ((ref)) holds as $(f_n,g_n)$ maximizes $\Phi_n(c,f,g)$, which implies that $\Phi_n(c,f_n,g_n) \geq \Phi_n(c,\hat{f}_n,\hat{g}_n)$.
We obtain an upper bound for $Re^{ROT}(\hat{\pi}_n^{ROT})$ by bounding each of the two difference terms in ((ref)).
We impose the following regularity condition.
In our application to a job search assistance program, the cost measure is the probability of unemployment, which satisfies this assumption with $\bar{c} = 1$.
We first present the following lemma, which is borrowed from Rigollet2022.
For completeness, we provide a proof of this lemma in Appendix (ref). Using this lemma, we obtain the following proposition, which gives finite sample upper bounds for $Re^{ROT}(\hat{\pi}_n^{ROT})$.
This proposition shows that the bound ((ref)) can be expressed as the sum of $L^1$-error and the squared $L^2$-error of the cost function estimator multiplied by the constant $e^{2 \eta \bar{c}}$. This constant comes from a uniform upper bound for $\hat{p}_n(x,y)$ implied by Lemma (ref). If $\hat{c}$ is consistent, the squared-$L^2$ error term vanishes faster than the $L^1$ error term, and the leading term of the convergence rate is governed by the $L^1$ error of the cost function estimator. For instance, if $c(x,y)$ is parametric and the parameters can be estimated at $N^{-1/2}$-rate, then $E_{P^N} \left[ \|\hat{c}-c\|_{L^1(\mu \otimes \nu)} \right] = O(N^{-1/2})$ and $E_{P^N} \left[ \|\hat{c}-c\|^2_{L^2(\mu \otimes \nu)} \right] = O(N^{-1})$.
The next proposition provides a uniform upper bound for the regularization bias defined in ((ref)).
Combining Propositions (ref) and (ref), the decomposition in ((ref)) yields the following bounds:
The first term in these bounds captures the variance component of regret. This grows exponentially with $\eta$ (i.e., as the regularization parameter $1/\eta$ decreases), but decreases as the size of the training sample $N$ increases and the integrated $L^1$ and $L^2$-errors of $\hat{c}$ shrink.
The second term in these bounds captures the regularization bias, which decreases as $\eta$ increases. Hence, the choice of $\eta$ exhibits a bias-variance trade-off.
Suppose that the convergence rate of the integrated $L^1$ error is $N^{-\alpha}$, e.g., the parametric rate of $\alpha = 1/2$, the nonparametric rate of $\alpha = -1/3$, etc. Then, if $\eta$ is chosen to be proportional to $\frac{\tilde{\alpha}}{2 \bar{c}} \log (N)$ with $\tilde{\alpha} < \alpha$, the upper bounds in ((ref)) and ((ref)) converge to zero at the rate $\log n / \log N$.
This first numerical example is designed to be a simple data generating process within which we can explore the performance of ROT and unregularized OT plans learned using an estimated cost function. We evaluate how the performance---in terms of regret relative to the unregularized oracle optimal policy---of the ROT plan derived from $\hat{c}$ depends on (i) the size of the training sample and (ii) the value of the regularization parameter $1/\eta$.
\paragraph{Data generating process.} The characteristics of job seekers $X_{\ell}$ and caseworkers $W_{\ell}$ both consist of a single variable drawn independently from the uniform distribution on $[0,1]$. \[ X_{\ell}, W_{\ell} \sim \mathcal{U}[0,1]. \] In the calibrated exercise of the next subsection, the job seekers' characteristic will correspond to a covariate-based prediction of their risk of remaining unemployed 6 months after entering unemployment, and the caseworkers' characteristic will correspond to their ability to place individuals in employment within 6 months of their registration as unemployed job seekers. The true job-finding probability, as a function of the job seeker's characteristic $x$ and their assigned caseworker's characteristic $w$. is given by: \[ p(x, w) = \frac{1}{3}x^2 + \frac{1}{3}w^2 + \frac{1}{3}xw. \] Since $\frac{\partial^2 p}{\partial x \partial w} > 0$, the optimal plan associated with this “production” function is positive assortative matching (PAM). The cost function used in the computation of the OT plans is $c(x,w) = 1-p(x,w) \in \left[0,1\right]$. The vector of outcomes for all pairs---a vector whose $\ell$th entry is a binary indicator taking value 1 if job seeker $\ell$ found a job within 6 months---is drawn as: \[ Y_{\ell} \sim \mathcal{B}(1-c(X_{\ell}, W_{\ell})), \] where $\mathcal{B}(p)$ denotes the Bernoulli distribution with success probability $p$.
The training samples used to estimate the cost function are samples of size $N \in \{500,$ $5\ 000, 50\ 000, 500\ 000\}$ drawn from the DGP described above. The “main” sample of job seekers and caseworkers to be matched one-to-one consists of 100 job seekers and 100 caseworkers (slots). This is a rough approximation of the number of matches that have to be made every month in a given unemployment agency.\footnote{See footnote 6 in the next subsection for a back-of-the-envelop computation of the monthly number of matches to be made in a given unemployment agency.} In this example we do not explore sampling uncertainty coming from the observations to be matched. The characteristics of this sample are equally spaced over the support of $X$ and $W$, i.e., $[0,1]$.
\paragraph{Cost Function Estimation}
To obtain an estimated cost function $\hat{c}(x, w)$, we train an XGBoost model on the training sample, regressing $Y_\ell$ on $(X_\ell, W_\ell)$. Figure (ref) shows example estimated cost functions for different values of the training sample size $N$.
\paragraph{Welfare evaluation of the different transport plans}
We compare the average (over training samples) welfare generated by the following matching policies:
The ROT plans are computed using the Sinkhorn algorithm in the log domain for numerical stability. We vary both the training sample size $N$ and the regularization parameter $1/\eta$.
We report the welfare gain (or loss) of each ROT plan relative to random matching (the default policy in the current French Public Employment Services system), rescaled by the welfare gap between random matching and PAM:
\[ \frac{(\mu_n \otimes \nu_n)(c) - \pi_n^{\mathrm{ROT}}(c)}{(\mu_n \otimes \nu_n)(c) - \pi_n^*(c)} \quad \text{and} \quad \frac{(\mu_n \otimes \nu_n)(c) - \hat{\pi}_n^{\mathrm{ROT}}(c)}{(\mu_n \otimes \nu_n)(c) - \pi_n^*(c)}. \]
This measure takes a negative value if the ROT plan performs worse than random matching, a value between 0 and 1 if the ROT plan outperforms random matching, and takes value 1 if the ROT plan coincides with unregularized OT plan.
We choose to report rescaled welfare gains/losses as this numerical example is not calibrated in any way---as opposed to the one in the following subsection. Due to the computational burden associated with the estimation of the cost function, the quantities reported are averages over only 30 repetitions of the (i)-training sampling/(ii)-cost estimation/(iii)-transport plan computation procedure.
Figure (ref) reports the average welfare gain of the oracle ROT plans learnt using the true cost function. As expected, as the regularization parameter goes to $0$, the average relative gain approaches that of the oracle unregularized plan.
Figure (ref) reports the welfare gains associated with ROT plans learnt using an estimated cost function. Interestingly, for the two smallest sample sizes, stronger regularization is beneficial in terms of welfare.
To better understand the performance of our procedure in the context of the French Public Employment Services data,\footnote{We are currently constructing the dataset necessary for our analysis, under a data access agreement with the French Public Employment Services.} we calibrate our simulations using results on caseworker value added reported in dromundo2022.
\paragraph{Data generating process.} A job seeker's characteristics consist of a single variable $X_{\ell}$ drawn from the $\operatorname{Beta}(\alpha = 3.7939, \beta=8.8634)$ distribution. $\alpha$ and $\beta$ are calibrated to match the moments of the distribution of the predicted 6-month job finding probability reported in Table 3 of dromundo2022.\footnote{Specifically, we match the above and below median averages (respectively 0.20 and 0.40) of this predicted probability in the sample of job seekers studied in dromundo2022.} Figure (ref) shows the resulting distribution. A caseworker's characteristics consists of a single variable $W_{\ell}$ drawn from the $\mathcal{N}(0,1)$ distribution. In Figure 1 of dromundo2022, estimates of caseworkers' value added appear to follow a normal distribution.
The true job-finding probability as a function of the job seeker's characteristic $x$ and their assigned caseworker's caracteristic $w$ is given by a standard logistic function: \[ p(x, w) = \frac{\exp(a + b \cdot x + c \cdot w + d \cdot x \cdot w)}{1 + \exp(a + b \cdot x + c \cdot w + d \cdot x \cdot w)} \] where the parameters $a$, $b$, $c$ and $d$ are calibrated to yield a given degree of complementarity. In particular, we define a gap $\gamma \in \{0.02, 0.06, 0.10\}$ in the effect of a one standard deviation increase in $w$ (from 0 to 1) on the job finding rate of job seekers with $x=0.20$ compared to job seekers with $x = 0.40$. We use the values $x = 0.20$ and $x = 0.40$ here as they are, respectively, the average values of $x$ for job seekers with below and above median characteristics. More precisely, we calibrate $(a,b,c,d)$ so that the function $p(x,w)$ matches
The cost function is given by $1-p(x,w) \in \left[0,1\right]$.
The vector of outcomes ---a vector whose $\ell$th entry is a binary indicator taking value 1 if job seeker $\ell$ found a job within 6 months---is drawn as before as \[ Y_{\ell} \sim \mathcal{B}(p(X_{\ell}, W_{\ell})), \] where $\mathcal{B}(p)$ denotes the Bernoulli distribution with success probability $p$.
Training samples---used to estimate the cost function---are of size $N \in \{500,$ $5\ 000, 50\ 000, 500\ 000, 5\ 000\ 000\}$ and follow the DGP described above. The 5 million observation sample is added to the sample sizes considered in Section (ref) as we expect the size of our training sample to be of this order of magnitude (several million observations).\footnote{Using the ILO definition, the stock of job seekers in France at any given point in time is roughly 5 million. We may have to filter our data before estimating the cost function, but we will have access to several years of data.} The “main” sample of job seekers and caseworkers to be matched one-to-one consists of 200 job seekers and 200 caseworkers (slots). This is a rough approximation of the number of matches that have to be made every month in a given unemployment agency.\footnote{dromundo2022 document that unemployment agencies in the Paris area have an average of 30 caseworkers, and on average a caseworker is assigned 30 job seekers each quarter. A “caseworker” in our DGP can be interpreted as one caseworker slot, hence the number of job seekers determines the number of matches to be formed. In a typical Parisian agency, there are $30 \times 30 = 900$ matches a quarter, i.e., about $300$ per month. Parisian agencies are likely to be larger than the typical French agency. Therefore, in our numerical simulations, we match 100 or 200 job seekers and caseworkers (slots).} We do not explore the effect of sampling uncertainty coming from the observations to be matched. Instead, the job seekers to be matched are given characteristics corresponding to 200 quantiles of the calibrated Beta distribution, while caseworkers are given characteristics corresponding to 200 quantiles of a $N(0,1)$.
\paragraph{Cost Function Estimation}
To obtain an estimated cost function $\hat{c}(x, w)$, we train an XGBoost model on the training sample, regressing $Y_\ell$ on $(X_\ell, W_\ell, X_\ell \times W_\ell)$. Figures (ref) and (ref) allow us to compare the true and estimated cost functions as the training sample size $N$ varies. In these figures the difference in the marginal effect of $w$ is held fixed at $\gamma = 0.06$.
\paragraph{Welfare evaluation of the different transport plans}
We compare the average (over training samples) welfare generated by the following matching policies:
The ROT plans are computed using the Sinkhorn algorithm in the log domain for numerical stability. The unregularized OT plan (computed with $1/\eta = 0$) is obtained using the Hungarian assignment algorithm. Simulations are conducted for different values of the training sample size $N$, regularization parameter $1/\eta$, and “complementarity” parameter $\gamma$.
We report the welfare gain (or loss) of each ROT plan in comparison to random matching (the default policy in the current French Public Employment Services system) in both (i) absolute terms (the pp. change in the job finding rate) and (ii) relative terms (the difference in welfare scaled by the welfare gap between random matching and the oracle unregularized OT plan). \[ \mu_n \otimes \nu_n)(c) - \hat\pi_n^{\mathrm{ROT}}(c) \quad \text{and} \quad \frac{(\mu_n \otimes \nu_n)(c) - \hat{\pi}_n^{\mathrm{ROT}}(c)}{(\mu_n \otimes \nu_n)(c) - \pi_n^*(c)}. \] To approximate the expectation over sampling uncertainty with respect to the training sample, we report the average values of the above quantities over 200 repetitions of the (i)-training set sampling/(ii)-cost estimation/(iii)-transport plan computation procedure.
Before considering feasible/learnable policies, we first discuss (as in the previous numerical example) the performance of the oracle policies---i.e., the unregularized and regularized OT plans learned using the (unknown) true cost function $c(x,w)$, \[ \mu_n \otimes \nu_n)(c) - \pi_n^{\mathrm{ROT}}(c) \quad \text{and} \quad \frac{(\mu_n \otimes \nu_n)(c) - \pi_n^{\mathrm{ROT}}(c)}{(\mu_n \otimes \nu_n)(c) - \pi_n^*(c)}. \] These are presented in Figure (ref). An interesting pattern in subfigure (ref) is that the cost of regularization is higher (in relative terms) when complementarities, as “proxied” by the parameter $\gamma$, are stronger and the gains from reallocation are higher. In addition, the results shown in subfigure (ref) indicate that even with $\gamma$ as low as $0.02$, the potential gains from reallocation (measured by the absolute gain associated with the unregularized OT plan, $1 / \eta=0$) approach $1$ percentage point. Given the absence of any reallocation cost, this is a sizeable economic gain. In comparison, behaghel2014 report intention-to-treat estimates for intensive job search counseling programs in France of between 1.5 pp. (se.: 0.8) for private providers and 2.8 pp. (se.: 1.2) for public providers, and the sizeable costs per program participant lead the authors to adopt a cautious stance on the cost-effectiveness of these programs---at least when privately provided.
Figure (ref) presents heatmaps describing the joint distribution of $(X,W)$ implied by the oracle transport plans. As expected given the complementarities featured in the true cost function, these oracle plans place most probability mass on the diagonal of the product space. This is optimal as the marginal effect of increasing the caseworker characteristic $w$ on the job finding probability is increasing with the jobs seekers' characteristic $x$ over most of its support, ---i.e., $\frac{\partial^2}{\partial x\partial w} p(x,w) > 0$ is increasing in $x$ over most of the support of $X$. However, the oracle unregularized OT plan is not pure positive assortative matching---as can be seen in panel (b) of Figure (ref). This is because we set the job finding probability function to be a logistic function.\footnote{Job seekers with $x \approx 0.5$ have a job finding probability equal to $0.5$ when matched with a caseworker whose characteristic is equal to the mean of the $W$ distribution (i.e. $w=0$). Since the derivative of the logistic function with respect to its index is highest at index $= \operatorname{logistic}^{-1}(0.5)$, in this setting the unregularized optimal plan assigns the top caseworkers to job seekers with $x \approx 0.5$.}
Now, we can turn to the performance of policies learned using estimated cost functions. Figure (ref) presents similar graphs to Figure (ref) for these feasible policies for one iteration of the (i)-training set sampling/(ii)-cost estimation/(iii)-transport plan computation procedure---with the cost function estimated using a training sample of $5,000,000$ observations. Reassuringly, when the cost function is estimated from a training dataset of this size, the resulting plans are quite similar to their oracle counterparts in Figure (ref). That said, the unregularized OT plan shown in panel (b) is noisier than its oracle counterpart, and appears to place some mass on a comparable support to the regularized plan shown in panel (c).
Figure (ref) reports similar quantities to Figure (ref) for feasible (R)OT plans learned using an estimated cost function. Overall, Figure (ref) parallels Figure (ref). Unlike the first numerical example, there is no clear pattern in terms of the benefits (or costs) of regularization for smaller training samples.\footnote{Although some amount of regularization appears to be, if anything, costless for the smallest training sample sizes.} For larger training sample sizes, the welfare gain from learned policies approaches that of the oracle unregularized optimal plan.
This paper studies policy learning for two-sided matching policies. We propose the use of entropy regularized empirical optimal transport and assess its welfare performance when matching caseworkers and job seekers in the context of a data generating processes calibrated to match French administrative data.
We leave several questions for future research. First, the performance of the estimated regularized policy can be sensitive to the choice of penalty parameter, and there is currently no method to tune it in a data-driven manner. Second, in our welfare regret analysis we obtain only an upper bound. A lower bound for regret remains to be investigated. Third, we set the base measure of the entropy penalty to independent coupling. Depending on the economic interpretation attached to the penalty term, it may be sensible to choose a different base. Investigation of the economic justification for the choice of penalty, in terms of both the base and the divergence measure, is left for future work. Future versions of this paper will also attempt to apply the method to real-world administrative data on the matching of job seekers and caseworkers to go beyond the simulations currently presented.