Migration as Submodular Optimization

Migration is one of the greatest societal challenges facing humanity in the 21st Century. In 2015 there were 244 million international migrants in the world, suggesting a larger increase in the rate of international migration than was previously anticipated (McAuliffe and Ruhs 2018). A key reason behind this increase is the widely understood fact that migration often has significant benefits to the migrants and their families, as well as the host country and the country of origin. For example, migrants working in the United States earn wages that are higher by a median factor of 4.11 than those they would have earned in their home countries (Clemens, Montenegro, and Prichett 2008); and the money they send back to their countries of origin is a reliable source of foreign currency. But the increase in the rate of migration is also driven by globalization, social disparities, and, in no small part, by a vast number of refugees — including, most recently, millions who have fled war in Syria, persecution in Myanmar, and economic calamity in Venezuela.

These events have fueled a surge of attempts to put migration, and, especially, refugee resettlement, on scientific footing. Alvin Roth, a Nobel laureate in economics, nicely summarizes the problem in a 2015 opinion piece (Roth 2015):

“Refugee relocation is what economists call a matching problem, in the sense that different refugees will thrive differently in different countries. Determining who should go where, and not just how many go to each country, should be a major goal of relocation policy.”

This observation underlies work in market design, which draws on the extensive literature on matching problems such as school choice. It focuses on refugee resettlement mechanisms that elicit refugees’ preferences over localities, and output a matching that is desirable with respect to these preferences (Moraga and Rapoport 2014; Jones and Teytelboym 2018; Delacrétaz, Kominers, and Teytelboym 2016).

By contrast, a recent paper by Bansak et al. 2018, published in Science, takes a markedly different, data-driven approach to refugee resettlement, which employs machine learning and optimization. Their goal is to maximize the expected number of refugees who find employment.¹¹ 1 To be precise, they consider families and maximize the expected number of families such that at least one member is employed, but this does not make a significant technical difference, so we focus on individuals for ease of exposition. Using historical data from the United States and Switzerland, they train predictors that estimate the probability $p_{i\ell}$ that a given refugee $i$ would find employment in a given locality $\ell$. The optimal solution is then an assignment of refugees to localities that maximizes the sum of probabilities, subject to capacity constraints. Assuming their predictions of employment probabilities are accurate, Bansak et al. 2018 demonstrate that their approach leads to a 40%–70% increase in the number of employed refugees, compared to the actual outcomes in the United States and Switzerland.

On a high level, we subscribe to the data-driven approach, and believe that the assumptions made by Bansak et al. 2018 are quite reasonable in the context of refugee resettlement — most notably, the implicit assumption that the probability $p_{i\ell}$ can be estimated based only on information about the refugee $i$ and the locality $\ell$, and does not depend on where other refugees are assigned. But as this approach gains traction, we envision it being deployed on a larger scale, and ultimately informing international migration policy more broadly, especially in the context of labor migration.

The key observation behind our work is that, at that scale, competition effects would invalidate the foregoing assumption. Indeed, the larger the number of, say, engineers, who settle in a specific area, the less likely it is that any particular engineer would find a job. The immigration of approximately 331,000 Jews from the former Soviet Union to Israel in 1990 and 1991 serves as a case in point. A disproportionate number of these highly-educated immigrants were engineers and doctors, leading to a saturation of the local job markets: “a state with an oversupply of doctors could not possibly double the number of its doctors in several years” (Smooha 2008). Our goal in this paper, therefore, is to enhance the data-driven approach to migration by explicitly modeling competition effects, and directly optimizing for them.

Let $N$ be a set of agents or migrants, and $L$ be a set of localities. Each locality $\ell$ has a capacity of $q_{\ell}\in\mathbb{N}$, limiting the number of migrants it is willing to accept. A matching is a set $S\subseteq N\times L$ of migrant-locality pairs, such that every agent is matched to at most one locality and every locality $\ell$ to at most $q_{\ell}$ migrants. Our general aim is to find matchings that maximize a given function $f:2^{N\times L}\to\mathbb{R}_{\geq 0}$, which assigns a utility to every matching.

In this paper, we will assume that $f$ is either submodular or approximately submodular. To be submodular, $f$ must satisfy for all $S\subseteq T\subseteq N\times L$ and $x\notin T$ that $f(S\cup\{x\})-f(S)\geq f(T\cup\{x\})-f(T)$. To prove submodularity, it is sufficient to show the above for $T$ being of the shape $S\cup\{y\}$ with $y\notin S$. A function $f$ is supermodular iff $(-f)$ is submodular, i.e., iff $f(S\cup\{x\})-f(S)\leq f(T\cup\{x\})-f(T)$ for appropriate $S$, $T$ and $x$. We will assume all submodular functions to be monotone and normalized, i.e., $f(\emptyset)=0$.

A function $f$ is $\epsilon$-approximately submodular for some $\epsilon>0$ if and only if there is a submodular function $\smash{\hat{f}}$ such that $(1-\epsilon)\,\smash{\hat{f}}(S)\leq f(S)\leq(1+\epsilon)\,\smash{\hat{f}}(S)$. Note that $f$ itself need not be monotone.

A (finite) matroid is a pair $(G,\mathcal{F})$ of a finite ground set $G$ and a set $\mathcal{F}\subseteq 2^{G}$ of independent sets such that (i) $\emptyset\in\mathcal{F}$, (ii) if $S\in\mathcal{F}$, then every subset $T$ of $S$ is also in $\mathcal{F}$, and (iii) if $S,T\in\mathcal{F}$ such that $|S|=|T|+1$, there is some $x\in S\setminus T$ such that $T\cup\{x\}\in\mathcal{F}$. For any set $S\subseteq G$, let its rank $r(S)$ denote the size of the largest independent subset of $S$, and let its span $\mathit{sp}(S)$ be the set of all elements $x\in G$ such that $r(S)=r(S\cup\{x\})$. An important class of matroids are the partition matroids that partition the ground set $G$ into disjoint subsets $G_{1},\dots,G_{k}$ and designate $S\subseteq G$ as independent iff $|S\cap G_{i}|$ is at most some cap $g_{i}$ for all $i$.

Both the constraint of matching each agent to at most one locality and the constraint induced by locality quotas easily translate into partition matroids over $N\times L$.³³ 3 We can even deal with migrant-locality incompatibilities, by removing incompatible pairs from the ground set. In the first matroid, we restrict a matching to select at most one edge out of the set of edges incident to each agent $i$; in the second matroid, to at most $q_{\ell}$ edges out of the set of edges incident to each locality $\ell$. Then, the matchings are exactly the sets that are independent in both matroids. Therefore, maximizing employment among matchings is an instance of maximizing a submodular function subject to multiple matroid constraints.

To optimize for employment, we must choose a way of predicting it under a given matching. For instance, we might use techniques from machine learning to generalize past observations about employment success. Alternatively, we might adopt a hand-crafted model — such as the ones presented in Section 4 — and set only its parameters based on data. The latter approach ensures that our predictor behaves reasonably on all inputs and requires less data for fitting.

However we choose to predict employment, the corresponding employment function will have a fairly complicated shape. As a result, we face the obstacle that optimizing the function exactly would be impractical from a computational viewpoint. Fortunately, Fisher, Nemhauser, and Wolsey 1978 found that a simple greedy algorithm gives an approximation guarantee of $\smash{\frac{1}{P+1}}$ when optimizing a monotone, submodular function $z$ subject to $P$ matroid constraints.⁴⁴ 4 Lee et al. 2010 report a polynomial-time algorithm that achieves an approximation ratio of $\frac{1}{P+\epsilon}$ for $P\geq 2$ partition matroids, where $\epsilon$ can be made arbitrarily small. However, the degree of the polynomial grows very fast as $\epsilon$ becomes smaller, making their algorithm impractical for our purposes. Furthermore, in practice, the same greedy algorithm has been found to perform much better than this theoretical guarantee, often giving results that are close to optimal (Kempe, Kleinberg, and Tardos 2003).

Call the matroids $\mathcal{M}_{p}=(N,\mathcal{F}_{p})$ for $1\leq p\leq P$, let $\mathit{sp}^{p}$ denote the span with respect to $\mathcal{M}_{p}$, and let the intersection of the matroids be $\mathcal{F}\coloneqq\smash{\bigcap_{p=1}^{P}}\mathcal{F}_{p}$. To refer to the marginal contribution of an element $j$ with respect to a set $S$, we set $\rho_{j}(S)\coloneqq z(S\cup\{j\})-z(S)$.

The greedy algorithm initializes $S^{0}\leftarrow\emptyset$, $N^{0}\leftarrow N$ and $t\leftarrow 1$. Then, it proceeds through the following steps, wherein $t$ denotes the number of the current iteration:

Step 0. If $N^{t-1}=\emptyset$, stop with $S^{t-1}$ as the greedy solution.

Step 1. Select $i(t)\in N^{t-1}$ for which $z(S^{t-1}\cup\{i(t)\})$ is maximal, with ties settled arbitrarily.

Step 2a. If $S^{t-1}\cup\{i(t)\}\notin\mathcal{F}$, set $N^{t-1}\leftarrow N^{t-1}\setminus\{i(t)\}$ and return to step 0.

Step 2b. If $S^{t-1}\cup\{i(t)\}\in\mathcal{F}$, set $\rho_{t-1}\leftarrow\rho_{i(t)}(S^{t-1})$, $S^{t}\leftarrow S^{t-1}\cup\{i(t)\}$, and $N^{t}\leftarrow N^{t-1}\setminus\{i(t)\}$.

Step 3. Set $t\leftarrow t+1$ and continue from step 0.

Although the greedy algorithm is an effective way of maximizing submodular functions, we are unlikely to have direct access to the submodular function we are interested in, even for the value queries required for the greedy algorithm. If we adopt a machine-learning approach, our predictor might get reasonably close to the real employment function but is unlikely to be exactly monotone and submodular. The same problem presents itself even for the hand-crafted models: While we prove that each model induces a perfectly submodular function, this function is defined as the expectation over a random process. For scenarios of nontrivial size, we can only estimate this function through repeated sampling, and the introduced noise may invalidate the monotonicity and submodularity of the function. In both cases, we can only expect our functions to be approximately submodular.

Fortunately, the greedy algorithm still performs well if $z$ is only approximately submodular. We adapt the classic result by Fisher, Nemhauser, and Wolsey 1978 to show the following approximation guarantee.

To the best of our knowledge, we are the first to give approximation bounds for optimizing approximately submodular functions subject to multiple matroid constraints. For a single matroid, Horel and Singer 2016 give a slightly sharper bound of $(2+\smash{\frac{2\,\epsilon}{1-\epsilon}}\,k)\,\smash{\hat{z}}(S)\geq\max_{S^{\prime}}\smash{\hat{z}}(S^{\prime})$, but their proof makes use of strong properties that only hold if $P=1$. By contrast, in order to capture our matching constraints we need to take the intersection of two matroids, that is, for our application $P=2$, and the corresponding worst-case approximation guarantee is roughly $1/3$ — assuming $\epsilon$ is sufficiently small. This is a realistic assumption, especially in our hand-crafted models, as standard concentration inequalities imply that sampling leads to exponentially fast convergence to the true function value.

Turning to the proof, let $U^{t}$ be the set of elements considered in the first $t+1$ iterations except for the element added to $S$ in iteration $t+1$, i.e., $U^{t}\coloneqq N\setminus N^{t}$. We make use of the following known results.

We are now ready for the theorem’s proof.

We propose three models for the submodular effects of competition between migrants. Each of them makes different assumptions about how migrants find their jobs and, thus, about how they compete with each other. In each model we are interested in the expected employment function that it induces, i.e., the function that takes as input agent-locality pairs, and outputs the expected number of agents that find employment (in each case this function depends on parameters of the model).

The approximation results of Section 3 provide a way of respecting competition effects when matching migrants to localities. But when is it worth adopting such an approach? After all, besides the potential benefits, embracing submodularity also entails an increase in modeling complexity, and it forces us to abandon the efficient optimization tools that are applicable to additive optimization.

If one chooses not to deal with these drawbacks, one can always just ignore submodularity: For each agent and locality, one may estimate the probability of her finding employment at this place in the absence of all other agents. One can then optimize the additive function induced by the aformentioned probabilities, say, by formulating the problem as an integer linear program, and hope that the result would do well under the true objective function. The approach of Bansak et al. 2018 can be understood as doing exactly this. Our goal is to show that — when competition effects are indeed present, and the objective function is submodular — accounting for submodularity is almost always the preferable choice and can lead to significantly better outcomes than the additive approach described above. To this end, we empirically evaluate the two approaches on all three models from Section 4.

Our simulation code is written in Python; we use Gurobi for additive optimization and IGraph for computing maximum bipartite matchings. All code is open source and available for reproduction at https://github.com/pgoelz/migration. We provide additional simulation results in Appendix C, where we also link to IPython notebooks with the exact setup of all of our experiments.

For increased performance, we reuse estimations of expected employment. When a model is queried with a matching that puts the same agents of profession $\pi$ — or just the same agents in case of the coordination model — in the same locality $\ell$ as a previous matching, the previous estimation of expected employment is used. Since all these estimations are obtained by averaging many random simulations, this should not significantly influence the experiments. Furthermore, since the additive algorithm is not affected, and since all final utility values are computed without memoization, this only disadvantages the greedy algorithm, strengthening our findings.

Due to the high number of parameters in all of our models, we adopt a standard setting that is applicable across all of them. Let $N$ be a set of 100 agents, split equally between two professions, 1 and 2. While we vary the number of localities, we distribute 100 jobs — 50 per profession — randomly over the localities, ensuring that each locality has at least one job. Furthermore, we set the capacity of each locality to exactly its number of jobs. Finally, we average 1 000 simulations to estimate expected employment whenever our models are queried with a matching.

We quickly sketch how this setting translates into the idiosyncrasies of our individual models: In the retroactive-correction model, for each agent $i$, we uniformly select a probability between 0 and 1, which is used as $p_{i\ell}$ for all localities $\ell$. We choose to keep these probabilities equal between localities to make sure that our samples will contain agents that have significantly different overall chances of being employed, an aspect that we would expect to see in any real dataset. As the correction function for a locality $\ell$ and profession $\pi$, we choose $x\mapsto\min(x,c)$, where $c$ is the number of available jobs in $\ell$ and $\pi$. Note that, while the cap ensures that the number of agents in a locality is at most the total number of jobs, the number of agents of a certain profession can exceed the number of jobs in that profession, which lets the cap $c$ kick in. Likewise, in the interview model, we choose $p_{i\ell}$ uniformly at random and equal over all localities. Finally, we induce the compatibility probabilities in the coordination model by the professions, i.e., each agent has a single compatibility probability chosen uniformly at random for all jobs of her profession and is incompatible with all other jobs.

As shown in Fig. 1, the greedy algorithm outperforms additive optimization in nearly all cases, frequently increasing employment by 10% and more. Given the wide range of scenarios that are randomly generated, it is striking that the greedy heuristic virtually never leads to worse employment than additive maximization. This trend persists over all settings that we simulated. While the advantage of the greedy algorithm does not seem as pronounced for small numbers of localities in the retroactive-correction model, it leads to strong improvements over all locality numbers in the other two models.

In the following, we set the number of localities to 10 and vary the number of agents $n$ — still equally split between professions — instead. Each locality has a capacity of $n/10$ and has $n/10$ jobs, where the $n/2$ jobs of each profession are randomly distributed within these constraints. Figure 2 shows that the greedy algorithm leads to impressive improvements across all simulated numbers of agents, again especially pronounced in the interview and coordination models. While the gains become smaller for high numbers of agents in the correction model, the reason is innocuous: Looking at the absolute utilities instead of the ratio, we see that in these scenarios both additive and greedy get very close to an employment of half of the agents, which is the expected number of qualifying agents and thus optimal (see Section C.2).

In Sections C.3, C.4 and C.5, we vary further parameters, namely the number of professions, the level of specialization of localities, and the availability of jobs with respect to fixed capacities. In general, we find that the improvements from the greedy optimization persist across these parameter settings.

Note that the results of the greedy algorithm provide a lower bound on what can be achieved when submodularity is taken into account. Other heuristics with access to the submodular function might be able to further improve employment while retaining the low query complexity of the greedy algorithm.

The simulations in the previous section revealed that employment in all of our models can be improved by accounting for submodularity. Since this effect appeared consistently between three significantly different models of how migrants find employment, we conjecture that similar gains can be obtained in practice. The next logical step is to measure competition effects on real migration data.

Alas, it seems to be very difficult to obtain datasets with the necessary amount of detail. For instance, the Database on Immigrants in OECD Countries (DIOC)⁵⁵ 5 https://www.oecd.org/els/mig/dioc.htm contains data on a vast number of migrants in different countries. However, while File D contains information on the level of educational attainment and on the migrants’ current occupations, there is no information about the field of education or pre-migration occupation. As a result, it is unclear whether an unemployed migrant was looking for employment at all and, if so, which job market she participated in. An additional problem is that DIOC does not track migrants who return to their home countries as a result of unemployment.

By contrast, longitudinal studies such as the German Socio-Economic Panel⁶⁶ 6 https://www.diw.de/en/diw_02.c.222517.en/data.html contain extensive information about the educational background of individuals. What is more, they follow individual migrants over an extended period of time and therefore track their success on the labor market in perfect detail. Unfortunately, the longitudinal studies that we saw observe too few individuals spread over too many localities. To directly observe competition, we need data on a large fraction of the migrants competing for jobs in a locality, from multiple localities.

Even the paper by Bansak et al. 2018, for which the authors had access to sensitive migration data, conspicuously does not use any features relating to the refugees’ field of occupation. Indeed, despite the predictive power of such features for employment success in a locality, at least one major refugee-resettlement agency in the US does not track this information at all. We have reason to believe that the same is true for other resettlement agencies as well. A clear policy recommendation, therefore, is to compile the records in question and to make them available for research. Hopefully, data on present migration will help to better direct flows of migration in the future — benefitting both the migrants and the societies they join.

This work was partially supported by the National Science Foundation under grants IIS-1350598, IIS-1714140, CCF-1525932, and CCF-1733556; by the Office of Naval Research under grants N00014-16-1-3075 and N00014-17-1-2428; and by a Sloan Research Fellowship and a Guggenheim Fellowship. We would like to thank Lucas Leopold for his helpful pointers on navigating migration datasets.