© 2018-2022 by Alfred Galichon. Past and present support from NSF grant DMS-1716489, ERC grant CoG-866274 are acknowledged, as well as inputs from contributors listed here.
If you reuse material from this masterclass, please cite as:
Alfred Galichon, 'math+econ+code' masterclass on optimal transport and economic applications, January 2022. https://github.com/math-econ-code/mec_optim
Optimal assignment problem
Pairwise stability, Walrasian equilibrium
Computation
Galichon, Optimal Transport Methods in Economics, Ch. 3
Roth, Sotomayor(1990). Two-Sided Matching. Cambridge.
Koopmans and Beckmann (1957). "Assignment problems and the location of economic activities." Econometrica.
Shapley and Shubik (1972). The assignment game I: The core." IJGT.
Becker (1993). A Treatise of the Family. Harvard.
Gretsky, Ostroy, and Zame (1992). "The nonatomic assignment model." Economic Theory.
Burkard, Dell'Amico, and Martello (2012). Assignment Problems. SIAM.
Dupuy and Galichon (2014). "Personality traits and the marriage market." JPE.
Consider the problem of assigning a possibly infinite number of workers and firms.
Each worker should work for one firm, and each firm should hire one worker.
Workers and firms have heterogenous characteristics; let $x\in \mathcal{X}$ and $y\in\mathcal{Y}$ be the characteristics of workers and firms respectively.
Workers and firms are in equal mass, which is normalized to one. Both worker's types and firm's types have a distribution on $\mathcal{X}$ and $\mathcal{Y}$ which is assumed known.
It is assumed that if a worker $x$ matches with a firm $y$, the total output generated is $\Phi_{xy}$. The questions are then:
Optimality: what is the optimal assignment in the sense that it maximizes the overal output generated?
Equilibrium: what are the equilibrium assignment and the equilibrium wages?
Efficiency: do these two notions coincide?
The same tools have been used by Gary Becker to study the heterosexual marriage market, where $x$ is the man's characteristics, and $y$ is the woman's characteristics, and "wages" are replaced by "transfers".
We shall take a first look at marriage data. Dupuy and Galichon (JPE, 2014) study a marriage dataset where, in addition to usual socio-demographic features (such as education and age), other features that measure personality traits are reported.
The literature on quantitative psychology argues that one can capture relatively well an individual's personality along five dimensions, the "big 5" - consciousness, extraversion, agreableness, emotional stability, autonomy - assessed though a standardized questionaire.
In addition to this, we observed a (self-assessed) measure of health, risk-aversion, education, height and body mass index = weight in kg/ (height in m)$^2$. In total, the available characteristics $x_{i}$ of man $i$ and $y_{j}$ of woman $j$ are both $10$-dimensional vectors.
It is assumed that the surplus of interaction is given by $\Phi\left(x,y\right) =x ^{\intercal}Ay$, where $A$ is a given $10 \times 10$ matrix. (later in this course, we'll see how to estimate $A$ based on matched marital data).
Today, we solve a central planner's problem (a stylized version of the problem OKCupids would solve): given a population of men and a population of women, how do we mutually assign these in order to:
maximize matching surplus
attain a (hopefully) stable assignment.
# !python -m pip install -i https://pypi.gurobi.com gurobipy ## only if Gurobi not here
gurobipy?
Object `gurobipy` not found.
import pandas as pd
import numpy as np
import scipy.sparse as spr
import gurobipy as grb
from tabulate import tabulate
thepath = 'https://raw.githubusercontent.com/math-econ-code/mec_optim_2021-01/master/data_mec_optim/marriage_personality-traits/'
data_X = pd.read_csv(thepath + "Xvals.csv")
data_Y = pd.read_csv(thepath + "Yvals.csv")
data_X.head()
| educm | heightm | BMIm | healthm | consm | extram | agreem | emom | autom | riskym | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 2 | 186 | 28.905075 | 3 | -0.752877 | -0.360787 | -0.711276 | -0.291031 | 0.840217 | 0.479437 |
| 1 | 2 | 176 | 27.440599 | 3 | 0.345542 | -0.805524 | -0.251796 | -0.305475 | -0.064454 | 0.030303 |
| 2 | 3 | 187 | 23.163374 | 3 | -0.759678 | 0.898007 | -0.029462 | -0.672859 | -0.961691 | -0.556598 |
| 3 | 1 | 184 | 29.241493 | 2 | -0.455688 | -1.053375 | -0.041612 | 0.436133 | 0.121873 | 0.992084 |
| 4 | 1 | 174 | 23.781214 | 4 | -1.440239 | 1.163730 | 0.293750 | -0.538922 | 0.782285 | -1.401034 |
We standardize the data using:
sdX,sdY = data_X.std().values, data_Y.std().values
mX,mY = data_X.mean().values, data_Y.mean().values
feats_x_k, feats_y_l = ((data_X-mX)/sdX).values, ((data_Y-mY)/sdY).values
nbx,nbk = feats_x_k.shape
nby,nbl = feats_y_l.shape
We postulate that the form of the surplus function is
\begin{align*}
\Phi_{xy}=x^{\intercal} Ay
\end{align*}
where $x$ and $y$ are the 10-dimensional features of a man and a woman, and the form of $A$, a 10x10 matrix, is given and stored in the file affinitymatrix.csv. Again, we'll see later how to solve the econometrics problem of estimating $A$.
We load $A$:
aff_data = pd.read_csv(thepath + "affinitymatrix.csv")
A_k_l = aff_data.iloc[0:nbk,1:nbl+1].values
print(tabulate(A_k_l))
----- ----- ----- ----- ----- ----- ----- ----- ----- ----- 0.56 0.02 -0.08 0.02 -0.04 -0.01 -0.03 -0.04 0.05 -0.02 0.01 0.18 0.04 -0.01 -0.04 0.05 0.02 0.02 0.02 0.02 -0.05 0.05 0.21 0.01 0.06 0 -0.04 0.04 -0.01 0.01 -0.07 0 -0.06 0.14 -0.04 0.05 -0.04 0.04 0.02 0 -0.06 -0.03 0.07 0 0.14 0.07 0.04 0.06 -0.02 -0.01 0.01 -0.02 0.05 0.02 -0.06 0.02 -0.02 -0.01 -0.03 -0.05 0 0.01 -0.08 0.02 0.13 -0.14 0.02 0.11 -0.09 -0.04 0.03 0 0.12 0.04 0.21 0.05 -0.03 -0.04 0.08 0.01 0.02 0 0 0.01 -0.11 0.11 -0.04 0.03 -0.09 0.01 0 0.02 -0.03 0.02 0.01 -0.01 -0.01 -0.05 0.05 0.11 ----- ----- ----- ----- ----- ----- ----- ----- ----- -----
Next, we build the surplus matrix $\Phi_{xy}$ using:
Φ_x_y = feats_x_k @ A_k_l @ feats_y_l.T
Finally, we have one individuals of each type, so:
We are ready to build a class OTProblem which summarize the data of the problem.
class OTProblem():
def __init__(self,Φ_x_y,n_x=None,m_y=None):
self.Φ_a = Φ_x_y.flatten()
self.nbx , self.nby = Φ_x_y.shape # number of types on each parts
self.nbz = self.nbx + self.nby # total number of types
self.nba = self.nbx * self.nby # number of arcs
if n_x is None:
self.n_x=np.ones(self.nbx)
else:
self.n_x = n_x
if m_y is None:
self.m_y=np.ones(self.nby)
else:
self.m_y = m_y
return
marriageEx = OTProblem(Φ_x_y)
twoBytwoEx = OTProblem(np.eye(2))
Assume that the type spaces $\mathcal{X}$ and $\mathcal{Y}$ are finite.
The mass of workers of type $x$ is $n_{x}$; the mass of jobs of type $y$ is $m_{y}$, with $\sum_{x}n_{x}=\sum_{y}m_{y}$.
Let $\mu_{xy}$ be the mass of workers of type $x$ assigned to jobs of type $y$. Every worker is busy and every job is filled, thus
\begin{align*} \mu\in\mathcal{M}\left(n,m\right) :=\left\{ \mu_{xy}\geq 0, \sum_{y\in\mathcal{Y}}\mu_{xy}=n_{x}\text{ and }\sum_{x\in\mathcal{X}}\mu _{xy}=m_{y} \right\}. \end{align*}
(Note that this formulation allows for mixing, i.e. it allows for $\mu_{xy}>0$ and $\mu_{xy^{\prime}}>0$ to hold simultaneously with $y\neq y^{\prime}$.)
Assume the economic output created when assigning worker $x$ to job $y$ is $\Phi_{xy}$. Hence, the optimal assignment is
\begin{align*} \max_{\mu\geq0} & \sum_{xy}\mu_{xy}\Phi_{xy} \\ s.t.~ & \sum_{y\in\mathcal{Y}}\mu_{xy}=n_{x} \left[u_{x}\right] \\ & \sum_{x\in\mathcal{X}}\mu_{xy}=m_{y}~\left[v_{y}\right] \end{align*}
Theorem
The value of the primal problem coincides with the value of the dual problem
\begin{align*} \min_{u,v} & \sum_{x\in\mathcal{X}}n_{x}u_{x}+\sum_{y\in\mathcal{Y}} m_{y}v_{y}\\ s.t. & u_{x}+v_{y}\geq\Phi_{xy}~\left[\mu_{xy}\geq0\right] \end{align*}
Both the primal and the dual problems have optimal solutions. If $\mu$ is a solution to the primal problem and $\left(u,v\right)$ is a solution to the dual problem, then by complementary slackness,
\begin{align*} \mu_{xy}>0\text{ implies }u_{x}+v_{y}=\Phi_{xy}. \end{align*}
Proof. The proof follows from the min-cost flow duality result, but let us rewrite it anyway.
The value of the primal problem can be written as $\max_{\mu\geq0}\min_{u,v}S\left( \mu,u,v\right)$, where
\begin{align*} S\left( \mu,u,v\right) :=\sum_{xy}\mu_{xy}\Phi_{xy}+\sum_{x\in\mathcal{X}}u_{x}(n_{x}-\sum_{y\in\mathcal{Y}}\mu_{xy})+\sum_{y\in\mathcal{Y}}v_{y}(m_{y}-\sum_{x\in\mathcal{X}}\mu_{xy}) \end{align*}
but by the minmax theorem, this value is equal to $\min_{u,v}\max_{\mu\geq 0}S\left( \mu,u,v\right)$, which is the value of the dual problem.
Follows by noting that, for a primal solution $\mu$ and a dual solution $\left( u,v\right) $, then $S\left( \mu,u,v\right) =\sum_{xy}\mu_{xy} \Phi_{xy}$. $\blacksquare$
The following statements are equivalent:
$\mu$ is an optimal solution to the primal problem, and $\left( u,v\right) $ is an optimal solution to the dual problem, and
(i) $\mu\in \mathcal{M}\left( n,m\right) $
(ii) $u_{x}+v_{y}\geq\Phi_{xy}$
(iii) $\mu_{xy}>0$ implies $u_{x}+v_{y}\leq\Phi_{xy}$.
We saw the direct implication. But the converse is easy: take $\mu$ and $\left( u,v\right) $ satisfying (i)--(iii), Then one has
\begin{align*} dual\leq\sum_{x}n_{x}u_{x}+\sum_{y}m_{y}v_{y}=\sum_{xy}\mu_{xy}\left( u_{x}+v_{y}\right) \leq\sum_{xy}\mu_{xy}\Phi_{xy}\leq primal \end{align*}but by the MK duality theorem, both ends coincide. Thus $\mu$ is optimal for the primal and $\left( u,v\right) $ for the dual.
In a marriage context, an important concept is stability:
An outcome is a vector $\left( \mu,u,v\right) $, where $u_{x}$ and $v_{y}$ are $x$'s and $y$'s payoffs, and $\mu$ is a matching that is
\begin{align*} \mu\in\mathcal{M}\left( n,m\right). \end{align*}
A pair $xy$ is blocking if $x$ and $y$ can find a way of sharing their joint surplus $\Phi_{xy}$ in such a way that $x$ gets more than $u_{x}$ and $y$ gets more than $v_{y}$. Hence there is no blocking pair if and only if for every $x$ and $y$, one has
\begin{align*} u_{x}+v_{y}\geq\Phi_{xy}. \end{align*}
If $x$ and $y$ are actually matched, their utilities $u_{x}$ and $v_{y}$ need to be feasible, i.e. the above inequality should be saturated. Hence
\begin{align*} \mu_{xy}>0\text{ implies }u_{x}+v_{y}=\Phi_{xy}. \end{align*}
Definition:
A matching that satisfies primal feasbilitity, no blocking, and complementary slackness is called a stable matching.
As it turns out, these conditions are precisely the conditions that express complementarity slackness in the Monge-Kantorovich problem. Therefore, outcome$\left( \mu,u,v\right) $ is stable if and only if $\mu$ is a solution to the primal problem, and $\left( u,v\right) $ is a solution to the dual problem.
Back to the workers / firms interpretation and assume for now that workers are indifferent between any two firms that offer the same salary. We argue that $u_x $ can be interpreted as the equilibrium wage of worker $x$, while $v_y $ can be interpreted as the equilibrium profit of firm $y$. Indeed:
Proposition If $\left(u,v\right)$ is a solution to the dual of the Kantorovich problem, then
\begin{align*} u_{x} & =\max_{y\in\mathcal{Y}}\left( \Phi_{xy}-v_{y}\right) \label{conjug1}\\ v_{y} & =\max_{x\in\mathcal{X}}\left( \Phi_{xy}-u_{x}\right) . \label{conjug2} \end{align*}Therefore, $u_{x}$ can be interpreted as equilibrium wage of worker $x$, and $v_{y}$ as equilibrium profit of firm $y$. In this interpretation, all workers get the same wage at equilibrium.
Assume now that if a worker of type $x$ works for a firm of type $y$ for wage $w_{xy}$, then gets $\alpha_{xy}+w_{xy}$, where $\alpha_{xy}$ is the nonmonetary payoff associated with working with a firm of type $y$. The firm's profit is $\gamma_{xy}-w_{xy}$, where $\gamma_{xy}$ is the economic output.
If an employee of type $x$ matches with a firm of type $y$, they generate joint surplus $\Phi_{xy}$, given by%
\begin{align*} \Phi_{xy}=\underset{\text{employee's payoff}}{\underbrace{\alpha_{xy}+w_{xy}}}+\underset{\text{firm's payoff}}{\underbrace{\gamma_{xy}-w_{xy}}}=\alpha_{xy}+\gamma_{xy} \end{align*}which is independent from $w$.
Workers choose firms which maximize their utility, i.e. solve
\begin{align*} u_{x}=\max_{y}\left\{ \alpha_{xy}+w_{xy}\right\} \end{align*}
and $u_{x}=\alpha_{xy}+w_{xy}$ if $x$ and $y$ are matched. Similarly, the indirect payoff vector of firms is
\begin{align*} v_{y}=\max_{x}\left\{ \gamma_{xy}-w_{xy}\right\} \end{align*}
and, again, $v_{y}=\gamma_{xy}-w_{xy}$ if $x$ and $y$ are matched.
As a result,
\begin{align*} u_{x}+v_{y}\geq\alpha_{xy}+\gamma_{xy}=\Phi_{xy} \end{align*}and equality holds if $x$ and $y$ are matched. Thus, if $w_{xy}$ is an equilibrium wage, then the triple $\left( \mu,u,v\right)$ where $\mu$ is the corresponding matching, and $u_{x}$ and $v_{y}$ are defined by this and this defines a stable outcome.
Conversely, let $\left(\mu,u,v\right)$ be a stable outcome. Then let $\bar{w}_{xx}$ and \b{w}$_{xy}$ be defined by
\begin{align*} \bar{w}_{xy}=u_{x}-\alpha_{xy}\text{ and }w^l_{xy}=\gamma_{xy}-v_{y}. \end{align*}One has $\bar{w}_{xy}\geq b^l_{xy}$. Any $w_{xy}$ such that $\bar{w}_{xy}\geq w_{xy}\geq b^l_{xy}$ is an equilibrium wage. Indeed, $\mu_{xy}>0$ implies $\bar{w}_{xy} \geq b^l_{xy}$, thus this and this hold. Given $u$ and $v$, $w_{xy}$ is uniquely defined on the equilibrium path (ie. when $x$ and $y$ are such that $\mu_{xy}>0$), but there are multiple choices of $w$ outside the equilibrium path.
Note that all workers of the same type get the same indirect utility, but not necessarly the same wage.
This problem of computation of the Optimal Assignment Problem, more specifically of $\left(\mu,u,v\right)$, is arguably the most studied problem in Computer Science, and dozens, if not hundreds of algorithms exist, whose running time is polynomial in $\max\left(n,m\right)$, typically a power less than three of the latter.
Famous algorithms include: the Hungarian algorithm (Kuhn-Munkres); Bertsekas' auction algorithm; Goldberg and Kennedy's pseudoflow algorithm. For more on these, see the book by Burkard, Dell'Amico, and Martello, and a introductory presentation in http://www.assignmentproblems.com/doc/LSAPIntroduction.pdf.
Here, we will show how to solve the problem with the help of a Linear Programming solver used as a black box; our challenge here will be to carefully set up the constraint matrix as a sparse matrix in order to let a large scale Linear Programming solvers such as Gurobi recognize and exploit the sparsity of the problem.
Let $\mu$ and $\Phi$ be the matrices with typical elements $\left(\mu_{xy}\right) $ and $\left( \Phi_{xy}\right) $. We let $n$, $m$, $u$,$v$, and $1$ the column vectors with entries $\left( n_{x}\right)$, $\left( m_{y}\right) $, $\left( u_{x}\right) $, $\left( v_{y}\right)$, and $1$, respectively. The optimal assignment problem
\begin{align*} \max_{\mu\geq0} & \sum_{xy}\mu_{xy}\Phi_{xy} \\ s.t.~ & \sum_{y\in\mathcal{Y}}\mu_{xy}=n_{x}~\left[ u_{x}\right]\\ & \sum_{x\in\mathcal{X}}\mu_{xy}=m_{y}~\left[ v_{y}\right] \end{align*}Can be rewritten writes using matrix algebra as
\begin{align*} & \max_{\mu\geq0}Tr\left( \Pi^{\top}\mu\right) \\ & \mu 1_{\mathcal{Y}}=n\\ & 1_{\mathcal{X}}^{\top}\mu=m^{\top}. \end{align*}Recall that if $A$ is a $M\times P$ matrix and $B$ a $N\times Q$ matrix, and if $vec(A)$ vectorizes in the row-major order (i.e. concatenates the rows of $A$) then the Kronecker product $A\otimes B$ of $A$ and $B$ is a $MN\times PQ$ matrix such that
\begin{align*} vec\left( AXB^{\top}\right) =\left( A\otimes B\right) vec\left( X\right). \label{VecAndKronecker} \end{align*}In python, $A\otimes B$ is implemented by sparse.kron(A,B) of the library sparse of scipy.
The first constraint $I_{\mathcal{X}}\mu 1_{\mathcal{Y}}=n$, therefore vectorizes under the row-major order as
\begin{align*} \left( I_{\mathcal{X}} \otimes1_{\mathcal{Y}}^{\top}\right) vec\left( \mu\right) =vec\left( n\right), \end{align*}and similarly, the second constraint $1_{\mathcal{X}}^{\top}\mu I_{\mathcal{Y}}=m^{\top}$, vectorizes as
\begin{align*} \left( 1_{\mathcal{X}}^{\top} \otimes I_{\mathcal{Y}}\right) vec\left( \mu\right) =vec\left( m\right) . \end{align*}Note that the matrix $I_{\mathcal{X}} \otimes1_{\mathcal{Y}}^{\top}$ is of size $| \mathcal{X}| \times |\mathcal{X}| |\mathcal{Y}|$, and the matrix $ 1_{\mathcal{X}}^{\top} \otimes I_{\mathcal{Y}}$ is of size $ |\mathcal{Y}|\times |\mathcal{X}| |\mathcal{Y}|$; hence the full matrix involved in the left-hand side of the constraints is of size $ (|\mathcal{X}|+|\mathcal{Y}|)\times |\mathcal{X}| |\mathcal{Y}|$. In spite of its large size, this matrix is sparse. In python, the identity matrix $I_{N}$ is coded as scipy.sparse.identity(N).
We therefore introduce the margining-out matrix as:
$M=\binom{M_\mathcal{X}}{M_\mathcal{Y}}$
where $M_\mathcal{X}$ and $M_\mathcal{Y}$ are the partial margining-out matrices, defined as
$
\begin{eqnarray*}
M_\mathcal{X} := I_{\mathcal{X}}\otimes 1_{\mathcal{Y}}^{\top }\\
M_\mathcal{Y}=1_{\mathcal{X}}^{\top }\otimes I_{\mathcal{Y}}
\end{eqnarray*}
$
have the effect of taking the $x$ and $y$ margins, respectively.
We implement using:
def M_x_xy(nbx,nby):
return spr.kron(spr.identity(nbx), np.ones((1,nby)))
def M_y_xy(nbx,nby):
return spr.kron(np.ones((1,nbx)),spr.identity(nby))
def M_z_a(self):
return spr.vstack([M_x_xy(self.nbx,self.nby),M_y_xy(self.nbx,self.nby)])
OTProblem.M_z_a = M_z_a
For instance, let us see what the margining matrix is in the case of a 2x2 example:
twoBytwoEx.M_z_a().todense()
matrix([[1., 1., 0., 0.],
[0., 0., 1., 1.],
[1., 0., 1., 0.],
[0., 1., 0., 1.]])
We set $q=(n^\top,m^\top)^\top$, and implement using:
def q_z(self):
return(np.concatenate([self.n_x,self.m_y]))
OTProblem.q_z = q_z
The linear programming problem then becomes
$\begin{eqnarray}
& \max_{z\geq0}vec\left( \Phi\right) ^{\top} vec(\mu) \label{LPvectorized}\\
s.t.~ & M vec(\mu) = q
\end{eqnarray}
$
which is ready to be passed on to a linear programming solver such as Gurobi.
def solve_full_lp(self,OutputFlag=True):
m=grb.Model()
μ_a = m.addMVar(self.nba)
m.setObjective(self.Φ_a @ μ_a, grb.GRB.MAXIMIZE)
m.addConstr(self.M_z_a() @ μ_a == self.q_z())
m.setParam( 'OutputFlag', OutputFlag )
m.optimize()
if m.status == grb.GRB.Status.OPTIMAL:
μopt_a = np.array(m.getAttr('x')).reshape(self.nbx,self.nby)
popt_z = m.getAttr('pi')
uopt_x, vopt_y = popt_z[:self.nbx],popt_z[self.nbx:]
return (μopt_a.reshape((self.nbx,-1)),uopt_x,vopt_y)
OTProblem.solve_full_lp = solve_full_lp
(μ_x_y,u_x,v_y) = marriageEx.solve_full_lp()
Set parameter WLSAccessID
Set parameter WLSSecret
Set parameter LicenseID
Academic license - for non-commercial use only - registered to alfred.galichon@sciences-po.fr
Gurobi Optimizer version 9.5.0 build v9.5.0rc5 (linux64)
Thread count: 2 physical cores, 4 logical processors, using up to 4 threads
Academic license - for non-commercial use only - registered to alfred.galichon@sciences-po.fr
Optimize a model with 2316 rows, 1340964 columns and 2681928 nonzeros
Model fingerprint: 0x4ef07e60
Coefficient statistics:
Matrix range [1e+00, 1e+00]
Objective range [2e-07, 8e+00]
Bounds range [0e+00, 0e+00]
RHS range [1e+00, 1e+00]
Concurrent LP optimizer: dual simplex and barrier
Showing barrier log only...
Presolve removed 0 rows and 0 columns (presolve time = 5s) ...
Presolve removed 0 rows and 0 columns (presolve time = 12s) ...
Presolve time: 14.56s
Presolved: 2316 rows, 1340964 columns, 2681928 nonzeros
Ordering time: 0.02s
Barrier statistics:
AA' NZ : 1.341e+06
Factor NZ : 2.683e+06 (roughly 600 MB of memory)
Factor Ops : 4.144e+09 (less than 1 second per iteration)
Threads : 1
Objective Residual
Iter Primal Dual Primal Dual Compl Time
0 9.44954184e-13 -0.00000000e+00 1.16e+02 6.08e+01 1.36e+01 16s
1 -4.02310420e+01 1.92816327e+04 1.73e-12 9.10e-14 1.15e-01 17s
2 1.08229108e+02 7.94975758e+03 6.22e-15 8.51e-14 4.68e-02 17s
3 2.76235105e+02 4.89363140e+03 9.71e-14 8.79e-14 2.75e-02 18s
4 6.45576167e+02 2.68712610e+03 4.73e-14 5.33e-14 1.22e-02 19s
5 1.13888783e+03 2.23991775e+03 6.33e-15 5.42e-14 6.57e-03 20s
6 1.56194469e+03 2.11115875e+03 4.04e-14 5.77e-14 3.28e-03 21s
7 1.71990871e+03 2.07282241e+03 1.80e-14 6.07e-14 2.11e-03 21s
8 1.82885302e+03 2.03215701e+03 6.74e-14 8.53e-14 1.21e-03 22s
9 1.90821751e+03 2.00746350e+03 8.08e-13 6.39e-14 5.92e-04 23s
10 1.95046449e+03 1.98540910e+03 1.08e-12 6.39e-14 2.08e-04 23s
Barrier performed 10 iterations in 23.73 seconds (4.95 work units)
Barrier solve interrupted - model solved by another algorithm
Solved with dual simplex
Solved in 66451 iterations and 23.76 seconds (5.53 work units)
Optimal objective 1.973096540e+03
Who does man $0$ match with?
print('Man 0 matches with woman '+str(np.argwhere(μ_x_y[0,:] != 0)[0][0])+'.' )
Man 0 matches with woman 575.
A important variant of the problem exists with $\sum_{x\in\mathcal{X}}n_{x}\neq\sum_{y\in\mathcal{Y}}m_{y}$ and the primal constraints become inequality constraints. The duality then becomes
\begin{align*} \begin{array} [c]{rrr} \max_{\mu\geq0}\sum\mu_{xy}\Phi_{xy} & = & \min_{u,v}\sum_{x\in\mathcal{X} }n_{x}u_{x}+\sum_{y\in\mathcal{Y}}m_{y}v_{y}\\ s.t.~\sum_{y\in\mathcal{Y}}\mu_{xy}\leq n_{x} & & u\geq0,~v\geq0 \\ \sum_{x\in\mathcal{X}}\mu_{xy}\leq m_{y} & & u_{x}+v_{y}\geq\Phi_{xy} \end{array} \end{align*}def solve_partial_lp(self, OutputFlag = True):
m=grb.Model()
μ_a = m.addMVar(self.nba)
m.setObjective(self.Φ_a @ μ_a, grb.GRB.MAXIMIZE)
m.addConstr(self.M_z_a() @ μ_a <= self.q_z())
m.setParam( 'OutputFlag', OutputFlag )
m.optimize()
if m.status == grb.GRB.Status.OPTIMAL:
μopt_a = np.array(m.getAttr('x')).reshape(self.nbx,self.nby)
popt_z = np.array(m.getAttr('pi'))
uopt_x, vopt_y = popt_z[:self.nbx],popt_z[self.nbx:]
return (μopt_a.reshape((self.nbx,-1)),uopt_x,vopt_y)
OTProblem.solve_partial_lp = solve_partial_lp
(μ_x_y,u_x,v_y) = marriageEx.solve_partial_lp()
Gurobi Optimizer version 9.5.0 build v9.5.0rc5 (linux64)
Thread count: 2 physical cores, 4 logical processors, using up to 4 threads
Academic license - for non-commercial use only - registered to alfred.galichon@sciences-po.fr
Optimize a model with 2316 rows, 1340964 columns and 2681928 nonzeros
Model fingerprint: 0xbc69bfc8
Coefficient statistics:
Matrix range [1e+00, 1e+00]
Objective range [2e-07, 8e+00]
Bounds range [0e+00, 0e+00]
RHS range [1e+00, 1e+00]
Concurrent LP optimizer: dual simplex and barrier
Showing barrier log only...
Presolve removed 0 rows and 671378 columns
Presolve time: 2.35s
Presolved: 2316 rows, 669586 columns, 1339172 nonzeros
Ordering time: 0.03s
Barrier statistics:
AA' NZ : 6.696e+05
Factor NZ : 1.428e+06 (roughly 300 MB of memory)
Factor Ops : 1.073e+09 (less than 1 second per iteration)
Threads : 1
Objective Residual
Iter Primal Dual Primal Dual Compl Time
0 3.93757423e+04 6.74170472e+02 8.32e+01 4.46e+01 1.01e+01 4s
1 6.81803951e+02 1.45621301e+04 0.00e+00 6.35e-14 1.65e-01 4s
2 7.43476296e+02 6.31961841e+03 0.00e+00 4.66e-14 6.64e-02 5s
3 8.52296063e+02 3.58969923e+03 0.00e+00 3.18e-14 3.26e-02 5s
Barrier performed 3 iterations in 5.30 seconds (1.44 work units)
Barrier solve interrupted - model solved by another algorithm
Solved with dual simplex
Solved in 32828 iterations and 5.31 seconds (1.93 work units)
Optimal objective 1.973096540e+03
μ_x_y.sum()
1158.0
We saw that in the present example, even though everyone had the possibility of being unmatched, no individual actually took that option.
We can also study the payoff vectors.
u_x.min(),v_y.min()
(0.03885354583920321, 0.0)
We see that at least one woman gets a payoff of zero, which makes her indifferent between being matched and being unmatched. It looks like this solution is not the preferred one for women. Let's see if there is another one which is better for women.
Let's see if, among the solutions to the problem above, there is one which maximizes the welfare of women. Using https://www.gurobi.com/documentation/9.5/refman/py_model_setobjectiven.html we write:
def solve_dual_partial_multiobj_lp(self, favor_Y = True):
if favor_Y:
welfare_weights_z = np.concatenate([np.ones(self.nbx),-np.ones(self.nby) ]) # welfare weights for each individuals
else:
welfare_weights_z = np.concatenate([-np.ones(self.nbx),np.ones(self.nby) ]) # welfare weights for each individuals
m=grb.Model()
p_z = m.addMVar(self.nbz)
m.ModelSense = grb.GRB.MINIMIZE
m.setParam('OutputFlag', 0)
q_z = self.q_z()
m.setObjectiveN(q_z @ p_z, 0, priority = 2,weight = 1)
m.setObjectiveN(welfare_weights_z @ p_z,1,priority = 1, weight = 1)
m.addConstr(self.M_z_a().T @ p_z >= self.Φ_a)
m.optimize()
if m.status == grb.GRB.Status.OPTIMAL:
#μopt_a = np.array(m.getAttr('Pi')).reshape(self.nbx,self.nby)
popt_z = np.array(m.getAttr('x'))
uopt_x, vopt_y = popt_z[:self.nbx],popt_z[self.nbx:]
return (uopt_x,vopt_y) # μopt_a.reshape((self.nbx,-1))
OTProblem.solve_dual_partial_multiobj_lp = solve_dual_partial_multiobj_lp
u_x,v_y = marriageEx.solve_dual_partial_multiobj_lp(favor_Y=True)
u_x.min(),v_y.min()
Academic license - for non-commercial use only - registered to alfred.galichon@sciences-po.fr
(0.0, 0.038853545839203)
u_x,v_y = marriageEx.solve_dual_partial_multiobj_lp(favor_Y=False)
u_x.min(),v_y.min()
Academic license - for non-commercial use only - registered to alfred.galichon@sciences-po.fr
(0.03885354583920321, 0.0)
Exercise. Integral primal solutions. ([OTME], ex. 3.3) A doubly stochastic matrix is a $n \times n$ matrix $\mu$ which is a solution to
$M \mu = 1_{2n}$
where $M$ is the margining matrix seen above.
A permutation matrix is a doubly stochastic matrix whose entries are all integer.
The solution to $\max_{\mu \geq 0} \sum_{xy} \mu_{xy} \Phi_{xy} : M \mu = 1_{2n}$ is therefore a doubly stochastic matrix.
Write a program that takes as an input a doubly stochastic matrix solution to the optimal assignment problem and returns a permutation matrix solution to the same problem.
# answer:
import networkx as nx
def improve(μ_x_y):
nbx,nby = μ_x_y.shape
matched_pairs = [(x,nbx+y) for x in range(nbx) for y in range(nby) if ((μ_x_y[x,y] > 1e-9 ) and (μ_x_y[x,y] < 1-1e-9 ) ) ]
if len(matched_pairs)==0:
return False
g = nx.Graph()
g.add_nodes_from(list(range(nbx))+list(range(nbx,nby)))
g.add_edges_from(matched_pairs)
cycle = nx.find_cycle(g)
Δμ_x_y = np.zeros((nbx,nby ))
m = np.inf
for i,j in cycle:
if i < nbx:
Δμ_x_y[i,j-nbx] = 1
m = min(m,μ_x_y[i,j-nbx])
else:
Δμ_x_y[j,i-nbx] = -1
m = min(m,μ_x_y[j,i-nbx])
μ_x_y -= m* Δμ_x_y
return True
def purify(theμ_x_y):
μ_x_y = theμ_x_y.copy()
cont = True
while cont:
cont = improve(μ_x_y)
return(μ_x_y.round())
np.random.seed(107)
d = 10
theμ_x_y = 0.3*np.eye(d)[:,np.random.permutation(d)]+0.2*np.eye(d)[:,np.random.permutation(d)]+0.1*np.eye(d)[:,np.random.permutation(d)]+0.25*np.eye(d)[:,np.random.permutation(d)]+0.15*np.eye(d)[:,np.random.permutation(d)]
purify(theμ_x_y)
array([[0., 0., 0., 0., 0., 0., 0., 0., 0., 1.],
[0., 0., 0., 0., 1., 0., 0., 0., 0., 0.],
[0., 0., 1., 0., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 0., 1., 0., 0., 0.],
[0., 0., 0., 0., 0., 0., 0., 0., 1., 0.],
[0., 1., 0., 0., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 1., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 1., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 0., 0., 1., 0., 0.],
[1., 0., 0., 0., 0., 0., 0., 0., 0., 0.]])
Exercise. Integral dual solutions. ([OTME], ex. 3.2) Assume $c_\{xy\}∈{0,1}$ and consider the problem
$\begin{eqnarray}
V &=&\max_{u,v}\sum_{x}u_{x}-\sum_{y}v_{y}. \\
&&s.t.~u_{x}-v_{y}\leq c_{xy} \notag
\end{eqnarray}
$
(i) Show that this is a variant of the optimal transport problem seen in class.
(ii) Show that the value of problem is unchanged if one restricts the entries of u and v to be between 0 and 1.
(iii) Show that if $(u_{x},v_{y})\in [0,1]^2$, then
$ (u_{x},v_{y})=\int_0^1(1\{t\leq u_x\},1\{t\leq v_y\})dt.$
(iv) Assume $(u,v)$ is feasible for problem above (i.e. it satisfies the constraints). Show that for each $t\in[0,1]$, the vectors $(u^{t},v^{t})$ defined by $u_{x}^{t}=1\{t≤u_{x}\}$ and $v_{y}^{t}=1\{t≤v_{y}\}$ are also feasible.
(v) Write a program that in the above setting, takes a solution $(u_{x},v_{y})\in [0,1]^2$ and returns a solution $(u_{x},v_{y})\in \{0,1\}^2$.
# answer:
# t=0.5
# ut_x = np.int64(t <= u_x)
# vt_y = np.int64(t <= v_y)