© 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 equilibrium transport and matching models in economics, June 2021. https://github.com/math-econ-code/mec_equil
[GS62] David Gale and Lloyd Shapley (1962). "College Admissions and the Stability of Marriage." American Mathematical Monthly 69 (1), pp. 9–14.
[A00] Hiroyuki Adachi (2000). "On a characterization of stable matchings." Economics Letters 68 pp. 43–49.
[GL21] Alfred Galichon, Flavien Léger (2021). "Matching algorithms as coordinate update". In progress.
[NY09] Muriel Niederle and Leeat Yariv (2009). "Decentralized Matching with Aligned Preferences." NBER Working Paper No. w14840.
[GGH21] Alfred Galichon, Octavia Ghelfi and Marc Henry (2021). "Stable and extremely unequal." Working paper.
[GH19] Alfred Galichon, and Yu-Wei Hsieh."Aggregate stable matching with money burning." SSRN id=2887732.
[GHS21] Alfred Galichon, Yu-Wei Hsieh and Maxime Sylvestre."Monotone comparative statics for submodular functions, with an application to aggregated deferred acceptance." In progress.
[GKW19] Alfred Galichon, Scott Kominers, and Simon Weber (2019). "Costly Concessions: An Empirical Framework for Matching with Imperfectly Transferable Utility." Journal of Political Economy 127 no. 6, pp. 2875-2925.
Assume strict preferences:
$\alpha_{xy} \neq \alpha_{xy^\prime}$ for $y\neq y^\prime$
$\gamma_{xy} \neq \gamma_{x^\prime y}$ for $x \neq x^\prime$.
The preferences are usually described in an ordinal way, by specifying a order relation for each agent on the other side of the market.
Let's first load the libraries we will need.
import numpy as np
import networkx as nx
import scipy.optimize as opt
We create the NTU_market class to encompass that information and store it in a convenient way:
class NTU_market:
def __init__(self,α_x_y,γ_x_y,n_x = np.array([]), m_y = np.array([])):
nbx,nby = α_x_y.shape
self.α_x_y = np.hstack((α_x_y,np.zeros((nbx,1) )))
self.γ_x_y = np.vstack((γ_x_y,np.zeros( (1,nby) )))
if n_x.size == 0:
n_x = np.ones(nbx)
if m_y.size == 0:
m_y = np.ones(nby)
self.n_x,self.m_y = n_x, m_y
self.largex,self.largey = nby+1, nbx+1
self.smallx,self.smally = -1, -1
self.nbx,self.nby = nbx, nby
self.αo_x_y = np.zeros((nbx,nby+1), dtype = 'int64')
self.γo_x_y = np.zeros((nbx+1,nby), dtype = 'int64')
self.prefslistα_x_y = np.zeros((nbx,nby+1), dtype = 'int64')
self.prefslistγ_x_y = np.zeros((nbx+1,nby), dtype = 'int64')
self.traceuo_x_t = np.array([])
self.traceu_x_t = np.array([])
for x in range(nbx):
thelistx = (- self.α_x_y )[x,:].argsort()
self.αo_x_y[x, thelistx] = nby - np.arange(nby+1)
self.prefslistα_x_y[x,:] = (thelistx) % (self.nby+1)
for y in range(nby):
thelisty = ( - self.γ_x_y)[:,y].argsort()
self.γo_x_y[thelisty, y] = nbx - np.arange(nbx+1)
self.prefslistγ_x_y[:,y] = (thelisty) % (self.nbx+1)
self.comp_nbsteps = -1
self.comp_time = -1
self.eq_μ_x_y = np.array([])
self.eq_u_x = np.array([])
self.eq_v_y = np.array([])
The following print_prefs method prints the ordinal preferences of the agents:
def list_prefs(x,xletter,yletter,αo_y):
nby = len(αo_y) - 1
xsprefs = xletter + str(x) +' : '
for y in range(nby):
ind = np.where(αo_y== (nby - y ))[0][0]
if (ind == nby):
break
xsprefs = xsprefs + yletter + str( ind ) + ' > '
return(xsprefs[:-3])
def print_prefs(self,xs=[],ys=[]):
if xs == [] and ys==[] :
xs = range(self.nbx)
ys = range(self.nby)
if xs != [] :
for x in xs :
print(list_prefs(x,'x','y',self.αo_x_y[x,:]))
if xs != [] and ys != [] :
print('===')
if ys != [] :
for y in ys :
print(list_prefs(y,'y','x',self.γo_x_y[:,y]))
NTU_market.print_prefs = print_prefs
We simulate such a market and print the corresponding preferences:
np.random.seed(seed=1000)
running_mkt = NTU_market(np.random.rand(3,2)-0.1,np.random.rand(3,2)-0.2)
running_mkt.print_prefs()
x0 : y0 > y1 x1 : y0 > y1 x2 : y0 > y1 === y0 : x1 > x2 y1 : x1 > x2 > x0
The following example is taken from [RS90], example 2.17
# example 2.17 in [RS90]
α_x_y = np.array([[0,1,2,3],[1,0,3,2],[2,3,0,1],[3,2,1,0]])
γ_x_y = np.array([[3,2,1,0],[2,3,0,1],[1,0,3,2],[0,1,2,3]])
rs_ex_2_17 = NTU_market(α_x_y,γ_x_y)
rs_ex_2_17.print_prefs()
x0 : y3 > y2 > y1 > y0 x1 : y2 > y3 > y0 > y1 x2 : y1 > y0 > y3 > y2 x3 : y0 > y1 > y2 > y3 === y0 : x0 > x1 > x2 > x3 y1 : x1 > x0 > x3 > x2 y2 : x2 > x3 > x0 > x1 y3 : x3 > x2 > x1 > x0
We shall see two notions of stable matchings with non-transferable utility (NTU):
Assume that there is one individual per type:
$n_x = 1$ and $m_y = 1$ for all $x\in\mathcal{X}$ and $y\in\mathcal{Y}$.
As before, $\mu_{xy}\in\{0,1\}$ will be a dummy variable equal to 1 iff $x$ is matched with $y$,
and $u_x$ and $v_y$ are the payoffs of worker $x$ and firm $y$ at equilibrium.
Definition. $(\mu,u,v)$ is a stable outcome in the Gale-Shapley sense if:
(i) Population constraints are satisfied
$\left\{
\begin{array}[l]
~\sum_{y} \mu_{xy} + \mu_{x0} = n_{x}=1 \\
\sum_{x}\mu_{xy} + \mu_{0y} = m_{y}=1
\end{array}
\right.$
(ii) There is no blocking pair and all individuals are rational:
$ \max \{ u_x - \alpha_{xy},v_y - \gamma_{xy} \} \geq 0$, and
$u_x \geq 0$ and $v_y \geq 0$
(iii) Strong complementarity:
$\mu_{xy}>0$ implies $u_x = \alpha_{xy}$ and $v_y = \gamma_{xy}$
$\mu_{x0}>0$ implies $u_x = 0$
$\mu_{0y}>0$ implies $v_y = 0$
In the Gale-Shapley definition:
$\mu_{xy}>0$ implies $u_x = \alpha_{xy}$ and $v_y = \gamma_{xy}$
This is stronger than assuming
$\mu_{xy}>0$ implies $\max \{ u_x - \alpha_{xy},v_y - \gamma_{xy} \} \geq 0$
Indeed, in the latter case, $(u_x,v_y)$ needs to be on the NTU frontier but doesn't have to be the efficient point $(\alpha_{xy},\gamma_{xy})$ as required in [GS62].
[GH19] and [GHS21] impose the latter condition as the definition of stability, which is weaker than Gale-Shapley. We shall see later why it may make sense not to impose the outcome to be Pareto efficient.
The following function detects stable matchings.
def is_GS_stable(self, μ_x_y = None, output=0 ):
if (min(self.n_x)<1) or (max(self.n_x)>1) or (min(self.m_y)<1) or (max(self.m_y)>1):
if output > 0 :
print('n_x or m_y do not only contain ones.')
return(False)
if (μ_x_y is None):
μ_x_y = self.eq_μ_x_y
μext_x_y0 = np.hstack([μ_x_y, 1-np.sum(μ_x_y,axis = 1).reshape(-1,1) ])
μext_x0_y = np.vstack([μ_x_y, 1-np.sum(μ_x_y,axis = 0).reshape(1,-1) ])
if (np.logical_and(μext_x_y0 != 0 , μext_x_y0 != 1 )).any() or (np.logical_and(μext_x0_y != 0 , μext_x0_y != 1 )).any():
if output > 0 :
print('The μ is not feasible.')
return(False)
uo_x = np.sum(μext_x_y0 * self.αo_x_y, axis = 1)
vo_y = np.sum(μext_x0_y * self.γo_x_y, axis = 0)
for x in range(self.nbx):
for y in range(self.nby):
if (self.αo_x_y[x,y] > uo_x[x]) and (self.γo_x_y[x,y] > vo_y[y]):
if output > 0 :
print('The matching is not stable. One blocking pair is: x'+str(x)+',y'+str(y)+'.')
return(False)
for x in range(self.nbx):
if (self.αo_x_y[x,self.nby] > uo_x[x]):
if output > 0 :
print('The matching is not stable. One blocking pair is: x'+str(x)+',y0.')
return(False)
for y in range(self.nby):
if (self.γo_x_y[self.nbx,y] > vo_y[y]):
if output > 0 :
print('The matching is not stable. One blocking pair is: x0,y'+str(y)+'.')
return(False)
if output > 0 :
print('The matching is stable.')
return (True)
NTU_market.is_GS_stable = is_GS_stable
Let's go back to example 2.17 in [RS90], and let's check the stability of the following matching:
μ_x_y = np.array([[0,0,0,1],[1,0,0,0],[0,1,0,0],[0,0,1,0]])
_ = rs_ex_2_17.is_GS_stable(μ_x_y,output=1)
The matching is not stable. One blocking pair is: x1,y3.
Let's try again with the following alternative matching:
μ_x_y = np.array([[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]])
_ = rs_ex_2_17.is_GS_stable(μ_x_y,output=1)
The matching is stable.
Let's construct some useful functions that will be useful later.
def μ_from_uo(self, uo_x):
return np.where( (self.αo_x_y == uo_x.reshape((-1,1))) , 1 , 0 )[:,0:-1]
NTU_market.μ_from_uo = μ_from_uo
def μ_from_vo(self, vo_y):
return np.where( self.γo_x_y == vo_y , 1 , 0 )[0:-1,:]
NTU_market.μ_from_vo = μ_from_vo
def uo_from_μ(self, μ_x_y = None):
if (μ_x_y is None):
μ_x_y = self.eq_μ_x_y
return np.sum(np.hstack([μ_x_y, (self.n_x-np.sum(μ_x_y,axis = 1)).reshape(-1,1) ]) * self.αo_x_y, axis = 1)
NTU_market.uo_from_μ = uo_from_μ
def vo_from_μ(self, μ_x_y = None):
if (μ_x_y is None):
μ_x_y = self.eq_μ_x_y
return np.sum(np.vstack([μ_x_y, (self.m_y-np.sum(μ_x_y,axis = 0)).reshape(1,-1) ]) * self.γo_x_y, axis = 0)
NTU_market.vo_from_μ = vo_from_μ
def u_from_μ(self, μ_x_y = None):
if (μ_x_y is None):
μ_x_y = self.eq_μ_x_y
return np.sum(np.hstack([μ_x_y, (self.n_x-np.sum(μ_x_y,axis = 1)).reshape(-1,1) ]) * self.α_x_y, axis = 1)
NTU_market.u_from_μ = u_from_μ
def v_from_μ(self, μ_x_y = None):
if (μ_x_y is None):
μ_x_y = self.eq_μ_x_y
return np.sum(np.vstack([μ_x_y, (self.m_y-np.sum(μ_x_y,axis = 0)).reshape(1,-1) ]) * self.γ_x_y, axis = 0)
NTU_market.v_from_μ = v_from_μ
The reference for this section is [GS62].
Principle: Workers make offers to firms that have not rejected them yet. Define:
$A^{t} \subseteq \mathcal{X} \times \mathcal{Y} =$set of available firms to worker at time $t$
$P^{t} \subseteq \mathcal{X} \times \mathcal{Y} =$set of proposals made by workers at time $t$
$E^{t}\subseteq \mathcal{X} \times \mathcal{Y} =$set of proposals kept by firms at the end of round $t$
For $B \subseteq \mathcal{X} \times \mathcal{Y}$, we introduce the following notations:
for $x \in \mathcal{X}$, denote $B(x)$ the set of $y \in \mathcal{Y}$ such that $xy \in B$.
for $y \in \mathcal{Y}$, denote $B(y)$ the set of $x \in \mathcal{X}$ such that $xy \in B$.
Deferred acceptance algorithm of Gale and Shapley.
At time $t=0$, all firms are avaiable to anyone $A^{0} =\mathcal{X} \times \mathcal{Y}$.
Iterate over $t$:
$ \left\{ \begin{array}{l} P^{t}=\left\{ xy\in \mathcal{X}\times \mathcal{Y}:y\in \arg \max_{y\in A^{t}\left( x\right) \cup \left\{ 0\right\} }\left\{ \alpha _{xy}\right\} \right\} \\ E^{t}=\left\{ xy\in \mathcal{X}\times \mathcal{Y}:x\in \arg \max_{x\in P^{t}\left( y\right) \cup \left\{ 0\right\} }\left\{ \gamma _{xy}\right\} \right\} \\ A^{t+1} =A^{t} \backslash \left\{ P^{t} \backslash E^{t} \right\} \end{array} \right. $
Repeat until $P^{t} =E^{t} $, i.e. no offer is rejected.
Let's implement as follows:
def solveGaleShapley(self,output=0, trace=False):
if (output>=2):
print("Offers made and kept are denoted +1; offers made and rejected are denoted -1.")
self.comp_nbsteps = 0
tracemax = self.nbx*self.nby
if trace:
self.traceuo_x_t = np.zeros((tracemax,self.nbx))
μA_x_y = np.ones((self.nbx+1, self.nby+1), dtype = 'int64') # initially all offers are non rejected
while True :
μP_x_y = np.zeros((self.nbx+1, self.nby+1), dtype = 'int64')
props_x = np.ma.masked_array(self.αo_x_y, μA_x_y[0:-1,:] ==0).argmax(axis = 1) # x's makes an offer to their favorite y
μP_x_y[range(self.nbx),props_x] = 1
μP_x_y[self.nbx,0:self.nby] = 1
μE_x_y = np.zeros((self.nbx+1, self.nby+1), dtype = 'int64') # y's retains their favorite offer among those made:
rets_y = np.ma.masked_array(self.γo_x_y, μP_x_y[:,0:-1] ==0).argmax(axis = 0)
μE_x_y[rets_y,range(self.nby)] = 1
rej_x_y = μP_x_y - μE_x_y # compute rejected offers
rej_x_y[self.nbx,:] = 0
rej_x_y[:,self.nby] = 0
μA_x_y = μA_x_y - rej_x_y # offers from x that have been rejected are no longer available to x
if output >= 2:
print("Round "+str(self.comp_nbsteps)+":\n" , ( (2*μE_x_y-1) * μP_x_y)[0:-1,0:-1])
if trace and self.comp_nbsteps<tracemax:
for x in range(self.nbx):
self.traceuo_x_t[self.comp_nbsteps,x] = np.sum(self.αo_x_y [x,:] * μP_x_y[x,:])
self.comp_nbsteps +=1
if np.max(np.abs(rej_x_y)) == 0:
if trace:
self.traceuo_x_t = self.traceuo_x_t[0:min(self.comp_nbsteps,tracemax),:]
break # if all offers have been accepted, then algorithm stops
self.eq_μ_x_y = μE_x_y[0:-1,0:-1]
return (0)
NTU_market.solveGaleShapley = solveGaleShapley
Test on our example:
running_mkt.solveGaleShapley(output=2,trace=True)
running_mkt.is_GS_stable()
Offers made and kept are denoted +1; offers made and rejected are denoted -1. Round 0: [[-1 0] [ 1 0] [-1 0]] Round 1: [[ 0 -1] [ 1 0] [ 0 1]] Round 2: [[0 0] [1 0] [0 1]]
True
We investigate a bit more systematically:
np.random.seed(seed=77)
nbx,nby=50,35
for i in range(10):
tstmkt = NTU_market(np.random.rand(nbx,nby)-0.3,np.random.rand(nbx,nby)-0.2)
tstmkt.solveGaleShapley()
print('Stable = '+str(tstmkt.is_GS_stable())+'; nb of steps = '+str(tstmkt.comp_nbsteps))
Stable = True; nb of steps = 71 Stable = True; nb of steps = 73 Stable = True; nb of steps = 83 Stable = True; nb of steps = 83 Stable = True; nb of steps = 93 Stable = True; nb of steps = 96 Stable = True; nb of steps = 78 Stable = True; nb of steps = 63 Stable = True; nb of steps = 69 Stable = True; nb of steps = 105
The reference for this section is [A00] and [GL21].
Adachi's algorithm.
Set initially $v^0_y$: say $v^0_y < \min_{x\in \mathcal{X}_0} \left\{ \gamma_{xy} \right\}$.
Iterate over t:
$\left\{ \begin{array}{l} u_{x}^{t+1}=\max \left\{ \max_{y\in \mathcal{Y}}\left\{ \alpha _{xy}:\gamma _{xy}\geq v_{y}^{t}\right\} ,\alpha _{x0}\right\} \\ v_{y}^{t+1}=\max \left\{ \max_{x\in \mathcal{X}}\left\{ \gamma _{xy}:\alpha _{xy}\geq u_{x}^{t+1}\right\} ,\gamma _{0y}\right\} \end{array} \right. $
until $u^{t+1} = u^{t}$.
We now implement Adachi's algorithm. It will be based on a pair of functions uo_from_vo and vo_from_uo which implement the maximization problem above, namely:
$\max \left\{ \max_{y\in \mathcal{Y}}\left\{ \alpha _{xy}:\gamma_{xy}\geq v_{y}\right\} ,\alpha _{x0}\right\}$, and
$\max \left\{ \max_{x\in \mathcal{X}}\left\{ \gamma _{xy}:\alpha
_{xy}\geq u_{x}\right\} ,\gamma _{0y}\right\}$.
def uo_from_vo(self,vo_y):
excluded = np.hstack([(self.γo_x_y < vo_y)[0:-1,:],np.zeros((self.nbx,1))])
return(np.ma.masked_array(self.αo_x_y, excluded ).max(axis = 1).data)
NTU_market.uo_from_vo = uo_from_vo
def vo_from_uo(self,uo_x):
excluded = np.vstack([(self.αo_x_y < uo_x.reshape((-1,1)))[:,0:-1],np.zeros((1,self.nby))])
return(np.ma.masked_array(self.γo_x_y, excluded ).max(axis = 0).data)
NTU_market.vo_from_uo = vo_from_uo
Adachi's algorithm iterates the loop above:
def solveAdachi(self,output=0,trace=False, startu_x = None, startv_y = None ):
self.comp_nbsteps = 0
tracemax = self.nbx*self.nby
if trace:
self.traceuo_x_t = np.zeros((tracemax,self.nbx))
if startu_x is None:
uo_x = self.largex * np.ones(self.nbx, dtype = 'int64') # x's utilities are highest
else:
uo_x = startu_x
if startv_y is None:
vo_y = self.smally * np.ones(self.nby, dtype = 'int64') # y's utilities are lowest
else:
vo_y = startv_y
while True :
uonew_x = self.uo_from_vo(vo_y) # each x proposes to favorite y among those willing to consider them
if (uonew_x == uo_x).all() :
if trace:
self.traceuo_x_t = self.traceuo_x_t[0:min(self.comp_nbsteps,tracemax),:]
break
uo_x = uonew_x
vo_y = self.vo_from_uo(uo_x) # each y proposes to favorite x among those willing to consider them
if output >= 2:
μP_x_y = self.μ_from_uo(uo_x)
μE_x_y = self.μ_from_vo(vo_y)
print("Round "+str(self.comp_nbsteps)+":\n","μ_P=\n" ,μP_x_y,"\n μ_E=\n", μE_x_y)
if trace and self.comp_nbsteps<tracemax:
self.traceuo_x_t[self.comp_nbsteps,:] = uo_x
self.comp_nbsteps += 1
self.eq_μ_x_y = self.μ_from_vo(vo_y)
return(0)
NTU_market.solveAdachi = solveAdachi
We try it on running_mkt:
running_mkt.solveAdachi()
running_mkt.is_GS_stable()
True
Next, we investigate a bit more Adachi's algorithm similarly to what we did above with Gale and Shapley's algorithm, and in comparison with the latter:
print("Times: Gale-Shapley; Adachi")
np.random.seed(seed=77)
nbx,nby=50,35
for i in range(10):
tstmkt = NTU_market(np.random.rand(nbx,nby)-0.3,np.random.rand(nbx,nby)-0.2)
res=tstmkt.solveGaleShapley()
if not tstmkt.is_GS_stable() :
raise Exception('Output of Gale Shapley is not stable')
i1 = tstmkt.comp_nbsteps
t1 = tstmkt.comp_time
mu1 = tstmkt.eq_μ_x_y
res=tstmkt.solveAdachi()
i2 = tstmkt.comp_nbsteps
t2 = tstmkt.comp_time
mu2 = tstmkt.eq_μ_x_y
if not tstmkt.is_GS_stable() :
raise Exception('Output of Adachi is not stable')
if np.max(np.abs(mu1 - mu2)) != 0:
raise Exception(seed,': results differ.')
print("GS steps = "+ str(i1)+"\t; A steps = "+str(i2))
print('All matching outcomes from Gale-Shapley and Adachi are stable and coincide.')
Times: Gale-Shapley; Adachi GS steps = 71 ; A steps = 10 GS steps = 73 ; A steps = 12 GS steps = 83 ; A steps = 10 GS steps = 83 ; A steps = 10 GS steps = 93 ; A steps = 12 GS steps = 96 ; A steps = 13 GS steps = 78 ; A steps = 11 GS steps = 63 ; A steps = 10 GS steps = 69 ; A steps = 11 GS steps = 105 ; A steps = 14 All matching outcomes from Gale-Shapley and Adachi are stable and coincide.
In this section, we show experimentally that Adachi's algorithm can equivalently be expressed by the following formulation, see [GL21]:
Adachi Bis algorithm.
Set initially $v^0_y$: say $v^0_y < \min_{x\in \mathcal{X}_0} \left\{ \gamma_{xy} \right\}$.
Iterate over t:
$\left\{ \begin{array}{l} u_{x}^{t+1}=\max \left\{ \max_{y\in \mathcal{Y}}\left\{ \alpha _{xy}:\gamma _{xy}\geq v_{y}^{t}\right\} ,\alpha _{x0}\right\} \\ v_{y}^{t+1}=\max \left\{ \max_{x\in \mathcal{X}}\left\{ \gamma _{xy}:\alpha _{xy}= u_{x}^{t+1}\right\} ,\gamma _{0y}\right\} \end{array} \right. $
until $u^{t+1} = u^{t}$.
def vo_from_uo_bis(self,uo_x):
excluded = np.vstack([(self.αo_x_y != uo_x.reshape((-1,1)))[:,0:-1],np.zeros((1,self.nby))])
return(np.ma.masked_array(self.γo_x_y, excluded ).max(axis = 0).data)
NTU_market.vo_from_uo_bis = vo_from_uo_bis
def solveAdachiBis(self,output=0,trace=False):
self.comp_nbsteps = 0
tracemax = self.nbx*self.nby
if trace:
self.traceuo_x_t = np.zeros((tracemax,self.nbx))
uo_x = self.largex * np.ones(self.nbx, dtype = 'int64') # x's utilities are highest
vo_y = self.smally * np.ones(self.nby, dtype = 'int64') # y's utilities are lowest
while True :
uonew_x = self.uo_from_vo(vo_y) # each x proposes to favorite y among those willing to consider them
if (uonew_x == uo_x).all() :
if trace:
self.traceuo_x_t = self.traceuo_x_t[0:min(self.comp_nbsteps,tracemax),:]
break
uo_x = uonew_x
vo_y = self.vo_from_uo_bis(uo_x) # each y proposes to favorite x among those willing to consider them
if trace and self.comp_nbsteps<tracemax:
self.traceuo_x_t[self.comp_nbsteps,:] = uo_x
self.comp_nbsteps += 1
self.eq_μ_x_y = self.μ_from_vo(vo_y)
return(0)
NTU_market.solveAdachiBis = solveAdachiBis
np.random.seed(seed=77)
nbx,nby=50,35
for i in range(10):
mkt_tst = NTU_market(np.random.rand(4,3)-0.2,np.random.rand(4,3)-0.3)
mkt_tst.solveAdachi(trace =True)
t1 = mkt_tst.traceuo_x_t
mkt_tst.solveAdachiBis(trace=True)
t2 = mkt_tst.traceuo_x_t
print("Discrepancy between Adachi and Adachi Bis=", np.max(np.abs(t2-t1)))
Discrepancy between Adachi and Adachi Bis= 0.0 Discrepancy between Adachi and Adachi Bis= 0.0 Discrepancy between Adachi and Adachi Bis= 0.0 Discrepancy between Adachi and Adachi Bis= 0.0 Discrepancy between Adachi and Adachi Bis= 0.0 Discrepancy between Adachi and Adachi Bis= 0.0 Discrepancy between Adachi and Adachi Bis= 0.0 Discrepancy between Adachi and Adachi Bis= 0.0 Discrepancy between Adachi and Adachi Bis= 0.0 Discrepancy between Adachi and Adachi Bis= 0.0
[GL21] reformulate Adachi's and Gale-Shapley algorithms as a variants of the (blockwise) Gauss-Seidel coordinate update algorithm.
As usual, change the sign of $u_{x}$ take $p=\left( u,-v\right) $
and define
$\left\{
\begin{array}{l}
Q_{x}\left( p\right) =p_{x}-\max_{y\in \mathcal{Y}}\left\{ \alpha
_{xy}:p_{y}\geq -\gamma _{xy},0\right\} ,x\in \mathcal{X} \\
Q_{y}\left( p\right) =p_{y}-\min_{x\in \mathcal{X}}\left\{ -\gamma
_{xy}:\alpha _{xy}\geq p_{x},0\right\} ,y\in \mathcal{Y}%
\end{array}%
\right.$
def Q_z(self,p_z):
uo_x = p_z[0:self.nbx]
vo_y = - p_z[self.nbx:(self.nbx+self.nby)]
uonew_x = self.uo_from_vo(vo_y)
vonew_y = self.vo_from_uo(uo_x)
return(np.append(uo_x - uonew_x,vonew_y - vo_y))
NTU_market.Q_z = Q_z
[GL21] reformulate Adachi's algorithm as a (blockwise) Gauss-Seidel algorithm.
Start by $p_{x}^{0}=\max_{y\in \mathcal{Y}_{0}}\left\{ \alpha _{xy}\right\} $ and $p_{y}^{0}=\max_{x\in \mathcal{X}_{0}}\left\{ -\gamma _{xy}\right\} $, which are such that $e\left( p\right) \geq 0$.
Adachi reinteprets as a blockwise Gauss-Seidel algorithm
$\left\{
\begin{array}{l}
p_{x}^{t+1}:Q_{x}\left( p_{x}^{t+1},\left( p_{y}^{t}\right) _{y}\right) =0,
\\
p_{y}^{t+1}:Q_{y}\left( p_{y}^{t+1},\left( p_{x}^{t+1}\right) _{x}\right) =0.%
\end{array}%
\right.$
def cux_z(self,p_z):
uo_x = p_z[0:self.nbx]
vo_y = - p_z[self.nbx:(self.nbx+self.nby)]
uonew_x = self.uo_from_vo(vo_y)
return (np.append(uonew_x , - vo_y))
NTU_market.cux_z = cux_z
def cuy_z(self,p_z):
uo_x = p_z[0:self.nbx]
vo_y = - p_z[self.nbx:(self.nbx+self.nby)]
excluded = np.vstack([(self.αo_x_y < np.repeat([uo_x],self.nby+1,axis = 1).reshape((self.nbx,-1)))[:,0:-1],np.repeat(False,self.nby).reshape((-1,self.nby))])
vonew_y = self.vo_from_uo(uo_x)
return( np.append(uo_x, - vonew_y ) )
NTU_market.cuy_z = cuy_z
def solveCU(self,output=0,trace=False):
self.comp_nbsteps = 0
tracemax = self.nbx*self.nby
if trace:
self.traceuo_x_t = np.zeros((tracemax,self.nbx))
p_z = np.append(self.largex * np.ones(self.nbx), - self.smally * np.ones(self.nby) )
while True :
pnew_z = self.cux_z(p_z)
pnew_z = self.cuy_z(pnew_z)
if (pnew_z == p_z).all() :
if trace:
self.traceuo_x_t = self.traceuo_x_t[0:min(self.comp_nbsteps,tracemax),:]
break
p_z = pnew_z
if trace and self.comp_nbsteps<tracemax:
self.traceuo_x_t[self.comp_nbsteps,:] = p_z[0:self.nbx]
self.comp_nbsteps += 1
self.eq_μ_x_y = self.μ_from_vo(- p_z[self.nbx:(self.nbx+self.nby)])
return(0)
NTU_market.solveCU = solveCU
We test numerically that Adachi and Gauss-Seidel coincide:
np.random.seed(1000)
mkt_tst = NTU_market(np.random.rand(80,120)-0.2,np.random.rand(80,120)-0.3)
mkt_tst.solveAdachi(trace =True)
t1 = mkt_tst.traceuo_x_t
mkt_tst.solveCU(trace=True)
t2 = mkt_tst.traceuo_x_t
print("Discrepancy between Adachi and Gauss-Seidel=", np.max(t2-t1))
Discrepancy between Adachi and Gauss-Seidel= 0.0
[GL21] reformulate Gale and Shapley's algorithm as a damped blockwise Gauss-Seidel algorithm.
Gale and Shapley's proposal step reformulates as
$\left\{
\begin{array}{l}
p_{x}^{t+1}=decr_{x}\left( p_{x}^{t}\right) \text{ if }Q_{x}\left(
p_{x}^{t+1},\left( p_{y}^{t}\right) _{y}\right) >0\text{, } \\
p_{x}^{t+1}=incr_{x}\left( p_{x}^{t}\right) \text{ if }Q_{x}\left(
p_{x}^{t+1},\left( p_{y}^{t}\right) _{y}\right) <0 \\
p_{x}^{t+1}=p_{x}^{t}\text{ if }Q_{x}\left( p_{x}^{t+1},\left(
p_{y}^{t}\right) _{y}\right) =0, \\
p_{y}^{t+1}:Q_{y}\left( p_{y}^{t+1},\left( p_{x}^{t+1}\right) _{x}\right) =0%
\end{array}%
\right.$
where
$decr_{x}\left( p\right) =\max_{y\in \mathcal{Y}_{0}}\left\{ \alpha
_{xy}:\alpha _{xy}<p\right\} $ is the next value below $p$ (in terms of the $\alpha _{xy}$'s), and
$incr_{x}\left(p\right) =\min_{x\in \mathcal{X}_{0}}\left\{ \alpha _{xy}:\alpha
_{xy}>p\right\} $ is the next value above $p$.
def damped_cux_z(self,p_z):
uo_x = p_z[0:self.nbx]
vo_y = - p_z[self.nbx:(self.nbx+self.nby)]
uonew_x = self.uo_from_vo(vo_y)
return (np.append(uo_x - np.where(uo_x > uonew_x, 1,0)+ np.where(uo_x < uonew_x, 1,0) , - vo_y))
NTU_market.damped_cux_z = damped_cux_z
def solveDampedCU(self,output=0,trace=False):
self.comp_nbsteps = 0
tracemax = self.nbx*self.nby
if trace:
self.traceuo_x_t = np.zeros((tracemax,self.nbx))
p_z = np.append(self.largex * np.ones(self.nbx), - self.smally * np.ones(self.nby) )
while True :
pnew_z = self.damped_cux_z(p_z) # each x proposes to favorite y among those willing to consider them:
pnew_z = self.cuy_z(pnew_z)
if (pnew_z == p_z).all() :
if trace:
self.traceuo_x_t = self.traceuo_x_t[0:min(self.comp_nbsteps,tracemax),:]
break
p_z = pnew_z
# each y proposes to favorite x among those willing to consider them:
if trace and self.comp_nbsteps<tracemax:
self.traceuo_x_t[self.comp_nbsteps,:] = p_z[0:self.nbx]
self.comp_nbsteps += 1
self.eq_μ_x_y = self.μ_from_vo(- p_z[self.nbx:(self.nbx+self.nby)])
return(0)
NTU_market.solveDampedCU = solveDampedCU
We test numerically that Gale and Shapley and damped Gauss-Seidel coincide:
np.random.seed(1000)
mkt_tst = NTU_market(np.random.rand(80,120)-0.2,np.random.rand(80,120)-0.3)
mkt_tst.solveGaleShapley(trace =True)
t1 = mkt_tst.traceuo_x_t
mkt_tst.solveDampedCU(trace=True)
t2 = mkt_tst.traceuo_x_t
print("Discrepancy between Gale and Shapley and coordinate update=", np.max(t2-t1))
Discrepancy between Gale and Shapley and coordinate update= 0.0
We build a wrapper that keeps track of the times.
from time import time
def solveDeferredAcceptance(self, algorithm = 'Adachi', output=0, trace=False):
start_time = time()
if algorithm == 'GS' :
self.solveGaleShapley(output,trace)
elif algorithm == 'Adachi' :
self.solveAdachi(output,trace)
elif algorithm == 'DARUM' :
self.solveDARUM(output,trace)
elif algorithm == 'CU' :
self.solveCU(output,trace)
else:
raise Exception("Algorithm " + algorithm + " is not implemented.")
self.comp_time = time() - start_time
if output >= 1:
print("Converged in ",self.comp_nbsteps," steps and ", self.comp_time, " seconds.")
if output == 1:
print("mu_x_y=",self.eq_μ_x_y)
return (0)
NTU_market.solveDeferredAcceptance = solveDeferredAcceptance
The following algorithm, described in [RS] p. 62, enumerates all stable matchings.
First, we implement a convenient method that returns a list of matched pairs.
def matched_pairs(self, μ_x_y = None):
if μ_x_y is None:
μ_x_y = self.eq_μ_x_y
nzx,nzy = np.nonzero(μ_x_y)
return [(nzx[i],nzy[i]) for i in range(len(nzx) )]
NTU_market.matched_pairs = matched_pairs
The following method returns the matchings (as designed by the payoff vectors) that are immediately above (or below) a given matching.
def next_us(self,uo_x,down = True):
μ_x_y = self.μ_from_uo(uo_x)
vo_y = self.vo_from_uo(uo_x)
if down:
incr= -1
else:
incr = 1
g = nx.DiGraph()
xs = list(range(self.nbx))
ys = list(range(self.nbx,self.nby))
g.add_nodes_from(xs+ys)
g.add_edges_from([(x,self.nbx+y) for x,y in self.matched_pairs(μ_x_y)])
#μprime_x_y = np.where( (self.αo_x_y.T == uo_x + incr).T , 1 , 0 )[:,0:-1]
μprime_x_y = np.where( (self.αo_x_y[:,0:-1].T == uo_x + incr).T & (self.γo_x_y[0:-1,:]*incr <= vo_y * incr), 1 , 0 )
compatible_pairs = self.matched_pairs(μprime_x_y)
g.add_edges_from([(self.nbx+j,i) for (i,j) in compatible_pairs])
cycles = list(nx.simple_cycles(g))
newus = []
for c in cycles:
xs_to_update = [x for x in c if (x <self.nbx)]
nextu_x = uo_x.copy()
nextu_x[xs_to_update] += incr
newus.append(nextu_x)
return(newus)
NTU_market.next_us = next_us
Let's run an example:
rs_ex_2_17.solveGaleShapley()
uo_x = np.sum(rs_ex_2_17.αo_x_y [:,0:-1] * rs_ex_2_17.eq_μ_x_y,axis = 1)
print(rs_ex_2_17.next_us(uo_x))
[array([4, 4, 3, 3], dtype=int64), array([3, 3, 4, 4], dtype=int64)]
The following method lists all stable matchings, starting with the lattice upper bound for side $\mathcal{X}$.
def enumerate_us(self):
self.solveGaleShapley()
uo_x = np.sum(self.αo_x_y [:,0:-1] * self.eq_μ_x_y,axis = 1)
us_list = [uo_x.tolist()]
i = 0
while (i < len(us_list)):
inds_alreadyhere = []
us_toadd = [arr.tolist() for arr in self.next_us(np.array(us_list[i],dtype='int64'))]
for u_x in us_toadd:
if u_x in us_list:
inds_alreadyhere.append(us_list.index(u_x))
us_toadd.remove(u_x)
inds_toadd = list(range(len(us_list),len(us_list)+len(us_toadd) ))
us_list = us_list + us_toadd
i += 1
return(np.array([u for u in us_list],dtype = 'int64') )
NTU_market.enumerate_us = enumerate_us
the_u_s = rs_ex_2_17.enumerate_us()
for u_x in the_u_s:
print('u_x = '+str(u_x)+ ' ; stable = '+str(rs_ex_2_17.is_GS_stable(μ_x_y)))
u_x = [4 4 4 4] ; stable = True u_x = [4 4 3 3] ; stable = True u_x = [3 3 4 4] ; stable = True u_x = [3 3 3 3] ; stable = True u_x = [3 2 2 3] ; stable = True u_x = [2 3 3 2] ; stable = True u_x = [2 2 2 2] ; stable = True u_x = [2 2 1 1] ; stable = True u_x = [1 1 2 2] ; stable = True u_x = [1 1 1 1] ; stable = True
Following [NY09], we define aligned preferences as the situation when $\alpha_{xy}=\gamma_{xy}=\varphi_{xy}$.
We setup a derived class for this called Aligned_NTU_market:
class Aligned_NTU_market(NTU_market):
def __init__(self,ϕ_x_y):
NTU_market.__init__(self,α_x_y = ϕ_x_y, γ_x_y = ϕ_x_y)
self.ϕ_x_y = ϕ_x_y
self.n,_ = ϕ_x_y.shape
ϕo_x_y = np.zeros(self.n* self.n)
ϕo_x_y[ϕ_x_y.flatten().argsort()] = 1+np.arange(self.n*self.n)
self.ϕo_x_y=ϕo_x_y.reshape((self.n,self.n))
The MaxMaxLex algorithm was proposed in [GGH21] to look for a stable matching in this type of markets. It looks for the pair $xy$ that maximizes $\varphi_{xy}$, matches $x$ and $y$, and removes them from the available individuals. We implement into:
def solveMaxMaxLex(self):
theϕo_x_y = self.ϕo_x_y.copy()
self.eq_μ_x_y = np.zeros((self.n,self.n), dtype = 'int64')
for k in range(self.n):
x,y = np.unravel_index(theϕo_x_y.argmax(),(self.n,self.n))
self.eq_μ_x_y[x,y]=1
theϕo_x_y[x,:]=0
theϕo_x_y[:,y]=0
return 0
Aligned_NTU_market.solveMaxMaxLex = solveMaxMaxLex
We compare the MaxMaxLex algorithm with Adachi:
np.random.seed(seed=77)
aligned_prefs_ex = Aligned_NTU_market(np.random.rand(20,20))
aligned_prefs_ex.solveMaxMaxLex()
print('MaxMaxLex:',aligned_prefs_ex.uo_from_μ(),'; stable = '+str(aligned_prefs_ex.is_GS_stable()))
aligned_prefs_ex.solveAdachi()
print('Adachi :',aligned_prefs_ex.uo_from_μ(),'; stable = '+str(aligned_prefs_ex.is_GS_stable()))
MaxMaxLex: [ 1. 20. 20. 20. 19. 16. 19. 20. 20. 20. 19. 20. 19. 14. 20. 20. 17. 20. 20. 20.] ; stable = True Adachi : [ 1. 20. 20. 20. 19. 16. 19. 20. 20. 20. 19. 20. 19. 14. 20. 20. 17. 20. 20. 20.] ; stable = True
The reference for this section are [GH19] and [GHS21].
We now allow for the possibility that there are more than one individual per type:
$n_x$ and $m_y$ don't need to be equal to one for all $x\in\mathcal{X}$ and $y\in\mathcal{Y}$.
As in the previous lecture, $\mu_{xy}$ is the mass of $xy$ pairs at equilibrium,
and $u_x$ and $v_y$ are the payoffs of worker $x$ and firm $y$ at equilibrium.
Definition. $(\mu,u,v)$ is an aggregate stable outcome if:
(i) Population constraints are satisfied
$\left\{
\begin{array}[l]
~\sum_{y} \mu_{xy} + \mu_{x0} = n_{x} \\
\sum_{x}\mu_{xy} + \mu_{0y} = m_{y}
\end{array}
\right.$
(ii) There is no blocking pair and all individuals are rational:
$ \max \{ u_x - \alpha_{xy},v_y - \gamma_{xy} \} \geq 0$, and
$u_x \geq 0$ and $v_y \geq 0$
(iii) Strong complementarity:
$\mu_{xy}>0$ implies $ \max \{ u_x - \alpha_{xy},v_y - \gamma_{xy} \} = 0$
$\mu_{x0}>0$ implies $u_x = 0$
$\mu_{0y}>0$ implies $v_y = 0$
We implement the distance-to-frontier function in this NTU matching setting as:
def DTF(self ,u_x,v_y):
return np.maximum(u_x.reshape((-1,1)) - self.α_x_y[:,0:-1], v_y - self.γ_x_y[0:-1,:])
NTU_market.DTF = DTF
The following function detects aggregate stable matchings:
def is_aggregate_stable(self, μ_x_y =None ,u_x=None,v_y =None, output=0 ):
if μ_x_y is None:
μ_x_y = self.eq_μ_x_y
if u_x is None:
u_x = self.eq_u_x
if v_y is None:
v_y = self.eq_v_y
if (np.min(μ_x_y) < 0) or (np.min(self.n_x-np.sum(μ_x_y,axis = 1))<0) or (np.min(self.m_y-np.sum(μ_x_y,axis = 0))<0) :
if output > 0 :
print('The μ is not feasible.')
return(False)
if (np.min(u_x)<0) or (np.min(v_y)<0):
if output > 0 :
print('u_x < 0 for some x or v_y < 0 for some y.')
return(False)
D_x_y = self.DTF(u_x,v_y)
if np.min(D_x_y)<0 :
if output > 0 :
print('There is a blocking pair.')
return(False)
if np.sum(μ_x_y * D_x_y)>0 :
if output > 0 :
print('Complementary slackness does not hold.')
return(False)
if output > 0 :
print('The matching is stable.')
return (True)
NTU_market.is_aggregate_stable = is_aggregate_stable
Consider the following example:
new_ex=NTU_market(np.array([[.2,.6],[.3,-.1],[.7,.5]]),np.array([[.7,.3],[.3,.6],[0.2,-.2]]) )
print('α_x_y=\n', new_ex.α_x_y[:,0:-1],'\nγ_x_y=\n',new_ex.γ_x_y[0:-1,:])
α_x_y= [[ 0.2 0.6] [ 0.3 -0.1] [ 0.7 0.5]] γ_x_y= [[ 0.7 0.3] [ 0.3 0.6] [ 0.2 -0.2]]
Solve this using Gale and Shapley into:
new_ex.solveGaleShapley()
theu,thev = new_ex.u_from_μ(),new_ex.v_from_μ()
print('\nu_x=',theu,'and v_y=',thev,'\n\nμ_x_y=\n',new_ex.eq_μ_x_y,'\n\n D_x_y=\n',new_ex.DTF(theu,thev))
u_x= [0.6 0.3 0. ] and v_y= [0.3 0.3] μ_x_y= [[0 1] [1 0] [0 0]] D_x_y= [[0.4 0. ] [0. 0.4] [0.1 0.5]]
It is interesting to note that while u_x= [0.6 0.3 0. ] and v_y= [0.3 0.3] are stable, one can decrease utility of some agents and remain stable (in the aggregate sense).
theu,thev = new_ex.u_from_μ(),new_ex.v_from_μ()
print("Is stable with full u and v :",new_ex.is_aggregate_stable(new_ex.eq_μ_x_y,theu,thev))
theu[0] = theu[0] - 0.3
theu[1] = theu[1] - 0.2
print("Is stable with altered u and v :",new_ex.is_aggregate_stable(new_ex.eq_μ_x_y,theu,thev))
Is stable with full u and v : True Is stable with altered u and v : True
In the next example we have one car and two passengers of the same type.
one_car_ex = NTU_market(np.array([[1]]),np.array([[1]]),n_x = np.array([1]), m_y = np.array([2]))
The sensible solution consists in matching the car to one of the passengers, and leave the other passenger unmatched. However, it is not a stable solution.
one_car_ex.is_GS_stable(np.array([[1]]))
False
However, the above solution is aggregate stable if $u=1$ for the car and $v=0$ for both passengers.
one_car_ex.is_aggregate_stable(np.array([[1]]),np.array([1]),np.array([0]))
True
The DARUM algorithm is an aggregate version of deferred acceptance. It is based on the specification of aggregate choice functions $C_{\mathcal{X}}$ and $C_{\mathcal{Y}}$ on each sides of the market:
$\mu \in C_{\mathcal{X}}\left( \overline{\mu }\right) $ means that $\mu $ is preferred (under $x$'s preferences) over all $\tilde{\mu%
}$ such that
$\tilde{\mu}\leq \overline{\mu }$ and $\sum_{y}\tilde{\mu}_{xy}\leq
n_{x}\forall x\in \mathcal{X}$.
$\mu \in C_{\mathcal{Y}}\left( \overline{\mu }\right) $ means that $\mu $ is preferred (under $y$'s preferences) over all $\tilde{\mu%
}$ such that
$\tilde{\mu}\leq \overline{\mu }$ and $\sum_{x}\tilde{\mu}_{xy}\leq m_{y}\forall y\in \mathcal{Y}$.
The DARUM algorithm is then:
Initialization: $\mu _{xy}^{A,0}=\min \left( n_{x},m_{y}\right) $
Repeat:
Until $| \mu ^{P,t}-\mu ^{D,t} |$ becomes below tolerance.
We implement the DARUM algorithm generically as follows -- we'll have to specify our aggregate choice functions later:
def solveDARUM(self,het1 = 'none',het2 = 'none',output=0,trace=False,tol = 1e-5):
self.comp_nbsteps = 0
tracemax = 100*self.nbx*self.nby
if (output>=2):
print("Offers made and kept are denoted +1; offers made and rejected are denoted -1.")
if trace:
self.traceu_x_t = np.zeros((tracemax,self.nbx))
self.traceuo_x_t = np.zeros((tracemax,self.nbx))
μA_x_y = np.array([[min(self.n_x[x],self.m_y[y]) for y in range(self.nby)] for x in range(self.nbx)])
# initially all offers are non rejected
while True :
μP_x_y,self.eq_u_x = self.aggregateChoice(0 , μA_x_y,het1) # the x's pick their preferred offers among those not rejected
μE_y_x,self.eq_v_y = self.aggregateChoice(1, μP_x_y[:,0:-1].T,het2)
μE_x_y = μE_y_x.T # the y's pick their preferred offers among those made
rej_x_y = μP_x_y[:,0:-1] - μE_x_y[0:-1,:] # compute rejected offers
μA_x_y = μA_x_y - rej_x_y # offers from x that have been rejected are no longer available to x
if output >= 2:
print("Round "+str(self.comp_nbsteps)+":\n" , ( (2*μE_x_y[0:-1,:]-1) * μP_x_y[:,0:-1]))
if trace and self.comp_nbsteps < tracemax:
self.traceu_x_t[self.comp_nbsteps,:] = self.eq_u_x
for x in range(self.nbx):
self.traceuo_x_t[self.comp_nbsteps,x] = np.sum(self.αo_x_y [x,:] * μP_x_y[x,:])
self.comp_nbsteps +=1
if np.max(np.abs(rej_x_y)) < tol:
if trace:
self.traceu_x_t = self.traceu_x_t[0:min(self.comp_nbsteps,tracemax),:]
self.traceuo_x_t = self.traceuo_x_t[0:min(self.comp_nbsteps,tracemax),:]
break # if all offers have been accepted (within tolerange), then algorithm stops
self.eq_μ_x_y = μE_x_y[0:-1,:]
return (0)
NTU_market.solveDARUM = solveDARUM
The function aggregateChoice above has not been defined yet; it will be a wrapper function to various aggregate choice functions.
def aggregateChoice(self,axis,μbar_x_y,heterogeneity = 'none'):
if heterogeneity == 'none':
μ_x_y,u_x = aggregateChoice_noHet(self,axis,μbar_x_y)
elif heterogeneity == 'logit':
μ_x_y,u_x = aggregateChoice_logit(self,axis,μbar_x_y)
else:
raise Exception("Heterogeneity " + heterogeneity + " is not supported.")
return (μ_x_y,u_x)
NTU_market.aggregateChoice = aggregateChoice
We shall implement the case with no heterogeneity (heterogeneity = 'none') first, and the logit case later.
We have in that case:
$\begin{array}{l}
C_{\mathcal{X}}\left( \overline{\mu }\right) =\arg \max_{\mu \geq 0}
&\left\{ \sum_{xy}\mu _{xy}\alpha _{xy}\right\} \\
s.t.~ &\sum_{y\in \mathcal{Y}}\mu _{xy}\leq n_{x} \\
~ &\mu _{xy}\leq \overline{\mu }_{xy}
\end{array}$
and
$\begin{array}{l}
C_{\mathcal{y}}\left( \overline{\mu }\right) =\arg \max_{\mu \geq 0}
&\left\{ \sum_{xy}\mu _{xy}\gamma _{xy}\right\} \\
s.t.~ &\sum_{x\in \mathcal{X}}\mu _{xy}\leq m_{y} \\
~ &\mu _{xy}\leq \overline{\mu }_{xy}
\end{array}$
We build this into:
def aggregateChoice_noHet(self,axis,μbar_x_y):
if axis == 0 : # if proposing side = x
n_x = self.n_x
nbx, nby = self.nbx, self.nby
prefs_x_y,α_x_y = self.prefslistα_x_y, self.α_x_y
else:
n_x = self.m_y
nbx, nby = self.nby, self.nbx
prefs_x_y, α_x_y = self.prefslistγ_x_y.T, self.γ_x_y.T
μ_x_y = np.zeros((nbx,nby+1))
u_x = np.zeros(nbx)
for x in range(nbx):
nxres = n_x[x]
for yind in prefs_x_y[x,]:
if yind == nby:
μ_x_y[x,yind] = nxres
break
if μbar_x_y[x,yind] > 0:
μ_x_y[x,yind] = min(nxres , μbar_x_y[x,yind])
nxres -= μ_x_y[x,yind]
if nxres == 0:
break
u_x[x]=α_x_y[x,yind]
return(μ_x_y,u_x)
We verify that in the case when there is one agent per type, DARUM coincides with Gale and Shapley at every step.
np.random.seed(seed=1000)
mkt_tst = NTU_market(np.random.rand(2,3)-0.2,np.random.rand(2,3)-0.3)
mkt_tst.solveDARUM(trace=True,output=0)
Darum_is_GS_stable = mkt_tst.is_GS_stable()
trace_Darum = mkt_tst.traceuo_x_t
mkt_tst.solveGaleShapley(trace=True,output=0)
GaleShapley_is_GS_stable = mkt_tst.is_GS_stable()
trace_GaleShapley = mkt_tst.traceuo_x_t
print("Both DARUM and Gale-Shapley outputs are stable=", (Darum_is_GS_stable & GaleShapley_is_GS_stable))
print("Discrepancy between DARUM and Gale-Shapley =", np.max(np.abs(trace_GaleShapley-trace_Darum)))
Both DARUM and Gale-Shapley outputs are stable= True Discrepancy between DARUM and Gale-Shapley = 0.0
one_car_ex.solveDARUM()
one_car_ex.eq_μ_x_y,one_car_ex.eq_u_x,one_car_ex.eq_v_y
(array([[1.]]), array([1.]), array([0.]))
We next implement the logit case, in which:
$\begin{array}{l}
C_{\mathcal{X}}\left( \overline{\mu }\right) =\arg \max_{\mu \geq 0}
&\left\{ \sum_{xy}\mu _{xy}\alpha _{xy} - \sum_{xy}\mu _{xy} \log \mu _{xy} \right\} \\
s.t.~ &\sum_{y\in \mathcal{Y}}\mu _{xy}\leq n_{x} \\
~ &\mu _{xy}\leq \overline{\mu }_{xy}
\end{array}$
The solution to the above problem has
$\mu _{xy}=\min \left\{ \mu _{x0}e^{\alpha _{xy}},\overline{\mu }% _{xy}\right\} $
where $\mu _{x0}$ solves
$\mu _{x0}+\sum_{y\in \mathcal{Y}}\min \left\{ \mu _{x0}e^{\alpha _{xy}},% \overline{\mu }_{xy}\right\} =n_{x}$
We build this into:
def aggregateChoice_logit(self,axis,μbar_x_y):
if axis == 0 : # if proposing side = x
n_x = self.n_x
nbx = self.nbx
nby = self.nby
α_x_y = self.α_x_y
else:
n_x = self.m_y
nbx = self.nby
nby = self.nbx
α_x_y = self.γ_x_y.T
μ_x_y = np.zeros((nbx,nby+1))
u_x = np.zeros(nbx)
for x in range(nbx):
thesolμ = opt.brentq(lambda theμ : (theμ+ np.minimum(theμ*np.exp(α_x_y[x,0:-1]),μbar_x_y[x,:]).sum() - n_x[x]) ,0,n_x[x])
μ_x_y[x,0:-1] = np.minimum(thesolμ*np.exp(α_x_y[x,0:-1]),μbar_x_y[x,:])
μ_x_y[x,nby] = thesolμ
u_x[x] = - np.log(thesolμ / n_x[x])
return(μ_x_y,u_x)
Run it on the one car example, for which we have:
$\left\{ \begin{array}{l}
\mu_{x0}+\min\{\mu_{x0},\mu_{0y} \} e = 1 \\
\mu_{0y}+\min\{\mu_{x0},\mu_{0y} \} e = 2
\end{array}\right.
$
thus $\mu_{0y}=\mu_{x0}+1$, so
$\mu_{x0} = (1+e)^{-1}$, $\mu_{0y} = (2+e)/(1+e)$
and $u_x = \log(1+e)$ and $v_y =\log(1+e) - \log(1+e/2)$.
one_car_ex.solveDARUM(het1='logit',het2='logit')
print('Solution from Darum: μ=',one_car_ex.eq_μ_x_y,', u=', one_car_ex.eq_u_x,', v=',one_car_ex.eq_v_y)
print('Predicted solution: μ=',np.exp(1)/(1+np.exp(1)),', u=', np.log(1+np.exp(1)),', v=',np.log(1+np.exp(1))-np.log(1+np.exp(1)/2))
Solution from Darum: μ= [[0.73105858]] , u= [1.31326169] , v= [0.45496415] Predicted solution: μ= 0.7310585786300049 , u= 1.3132616875182228 , v= 0.45496415414611713
As shown in the previous lecture, one can solve the problem using IPFP = Gauss-Seidel algorithm, as done in [GKW19].
This consists of iterately:
a. solving in $\mu_{x0}$ the equation
$\mu_{x0}+\sum_y \min \left\{ \mu _{x0}e^{\alpha _{xy}}, \mu _{0y}e^{\gamma _{xy}} \right\}=n_x,$
and:
b. solving in $\mu_{0y}$ the equation
$\mu_{0y}+\sum_x \min \left\{ \mu _{x0}e^{\alpha _{xy}}, \mu _{0y}e^{\gamma _{xy}} \right\}=m_y.$
until the update become below tolerance.
The solution is given by
$\mu_{xy} = \min \left\{ \mu _{x0}e^{\alpha _{xy}}, \mu _{0y}e^{\gamma _{xy}} \right\}.$
We build this into:
def solveIPFP(self,output=0,trace=False,tol = 1e-5):
self.comp_nbsteps = 0
tracemax = 100*self.nbx*self.nby
if trace:
self.traceu_x_t = np.zeros((tracemax,self.nbx))
μ_x0 = np.zeros(self.nbx)
μ_0y = self.m_y.astype(np.float)
while True :
for x in range(self.nbx):
μ_x0[x] = opt.brentq (lambda theμ : (theμ+ np.minimum(theμ*np.exp(self.α_x_y[x,0:-1]),μ_0y*np.exp(self.γ_x_y[x,:])).sum() - self.n_x[x]) ,0,1.1*self.n_x[x])
μP_x_y = np.minimum(μ_x0.reshape((-1,1))*np.exp(self.α_x_y[:,0:-1]),μ_0y*np.exp(self.γ_x_y[0:-1,:]))
self.eq_u_x = -np.log(μ_x0 / self.n_x)
for y in range(self.nby):
μ_0y[y] = opt.brentq(lambda theμ : (theμ+ np.minimum(μ_x0*np.exp(self.α_x_y[:,y]),theμ*np.exp(self.γ_x_y[0:-1,y])).sum() - self.m_y[y]) ,0,1.1*self.m_y[y])
μE_x_y = np.minimum(μ_x0.reshape((-1,1))*np.exp(self.α_x_y[:,0:-1]),μ_0y*np.exp(self.γ_x_y[0:-1,:]))
self.eq_v_y = -np.log(μ_0y / self.m_y)
rej_x_y = μP_x_y - μE_x_y # compute rejected offers
if output >= 2:
print('μP_x_y=\n',μP_x_y)
print('μE_x_y=\n',μE_x_y)
if trace and self.comp_nbsteps < tracemax:
self.traceu_x_t[self.comp_nbsteps,:] = self.eq_u_x
self.comp_nbsteps +=1
if np.max(np.abs(rej_x_y)) < tol:
if trace:
self.traceu_x_t = self.traceu_x_t[0:min(self.comp_nbsteps,tracemax),:]
break # if all offers have been accepted (within tolerance), then algorithm stops
self.eq_μ_x_y = μE_x_y
return (0)
NTU_market.solveIPFP = solveIPFP
We run the IPFP on the one-car example:
one_car_ex.solveIPFP()
print('Solution from IPFP: μ=',one_car_ex.eq_μ_x_y,', u=', one_car_ex.eq_u_x,', v=',one_car_ex.eq_v_y)
print('Predicted solution: μ=',np.exp(1)/(1+np.exp(1)),', u=', np.log(1+np.exp(1)),', v=',np.log(1+np.exp(1))-np.log(1+np.exp(1)/2))
Solution from IPFP: μ= [[0.73105858]] , u= [1.31326169] , v= [0.45496415] Predicted solution: μ= 0.7310585786300049 , u= 1.3132616875182228 , v= 0.45496415414611713
On the running_mkt example:
running_mkt.solveDARUM(het1='logit',het2='logit',trace=True)
traceD = running_mkt.traceu_x_t
running_mkt.eq_μ_x_y,running_mkt.eq_u_x,running_mkt.eq_v_y
(array([[0.21901611, 0.20862926],
[0.26548752, 0.32541622],
[0.25866091, 0.29466345]]),
array([0.55799648, 0.8938048 , 0.80592259]),
array([1.35931964, 1.76439096]))
running_mkt.solveIPFP()
traceI = running_mkt.traceu_x_t
running_mkt.eq_μ_x_y,running_mkt.eq_u_x,running_mkt.eq_v_y
(array([[0.21901611, 0.20862926],
[0.26548752, 0.32541622],
[0.25866091, 0.29466345]]),
array([0.55799648, 0.8938048 , 0.80592259]),
array([1.35931964, 1.76439096]))
Interestingly, in this example, the steps taken by both algorithms are exactly the same.
np.min(np.abs(traceD-traceI))
0.0
However, this is far from being the case in general. This is not surprising given that IPFP is a Gauss-Seidel algorithm, while DARUM is a deferred acceptance algorithm. See for example:
np.random.seed(77)
nbx,nby=3,2
mkt_tst = NTU_market(np.random.rand(nbx,nby)-0.2,np.random.rand(nbx,nby)-0.1)
mkt_tst.solveIPFP(trace =True,output=0)
t1 = mkt_tst.traceu_x_t
μ1 = mkt_tst.eq_μ_x_y
u1 = mkt_tst.eq_u_x
n1 = mkt_tst.comp_nbsteps
mkt_tst.solveDARUM(het1='logit',het2='logit',trace=True)
t2 = mkt_tst.traceu_x_t
μ2 = mkt_tst.eq_μ_x_y
u2 = mkt_tst.eq_u_x
n2 = mkt_tst.comp_nbsteps
print('Discrepancy betweeo outcomes of IPFP and DARUM = '+str(np.sum(np.abs(μ1.flatten()-μ2.flatten()))))
print('IPFP took '+str(n1)+' steps, while DARUM tool '+str(22)+' steps.')
print('trace IPFP\n',t1)
print('trace Darum\n',t2)
Discrepancy betweeo outcomes of IPFP and DARUM = 2.50716114535976e-13 IPFP took 3 steps, while DARUM tool 22 steps. trace IPFP [[1.52795096 1.30313513 1.3066622 ] [0.78064364 0.73407044 0.95041139] [0.74283247 0.70050346 0.90179972]] trace Darum [[1.52795096 1.30313513 1.3066622 ] [1.52795096 1.30313513 1.3066622 ] [1.52795096 1.30313513 1.3066622 ] [1.52795096 1.09519935 1.28411052] [1.34391957 0.95677379 1.00962647] [1.24457769 0.94029626 1.00962647] [0.9906809 0.94029626 1.00962647] [0.74283247 0.94029626 1.00962647] [0.74283247 0.94029626 1.00962647] [0.74283247 0.94029626 1.00962647] [0.74283247 0.82087269 1.00962647] [0.74283247 0.70050346 1.00962647] [0.74283247 0.70050346 1.00962647] [0.74283247 0.70050346 1.00962647] [0.74283247 0.70050346 1.00962647] [0.74283247 0.70050346 1.00962647] [0.74283247 0.70050346 1.00962647] [0.74283247 0.70050346 1.00962647] [0.74283247 0.70050346 1.00962647] [0.74283247 0.70050346 1.00962647] [0.74283247 0.70050346 1.00962647] [0.74283247 0.70050346 0.92860775] [0.74283247 0.70050346 0.90179972]]