© 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
Galichon (2016). Optimal Transport Methods in Economics. Chapter 5. Princeton University Press.
Anderson, de Palma, and Thisse (1992). Discrete Choice Theory of Product Differentiation. MIT.
Aurenhammer (1987). Power Diagrams: Properties, Algorithms and Applications. SIAM J Computing.
Lancaster (1966). A New Approach to Consumer Theory. JPE.
Berry, Pakes (2007). The Pure Characteristics Demand Model. IER.
Feenstra, Levinsohn (1995). Estimating Markups and Market Conduct with Multidimensional Product Attributes. ReStud.
Bonnet, Galichon, Shum (2017). Yoghurts Choose Consumers. Identification of Random Utility Models via Two-Sided Matching. Mimeo.
Leclerc, Merigot. pysdot library. https://github.com/sd-ot/pysdot
Today we'll consider a version of the transportation problem where we seek to match a continuous distribution on $\mathbb{R}^{d}$ with a discrete distribution. This problem is called a semi-discrete transportation problem.
Actually, we will introduce this problem not as a matching problem, but as a demand problem. We'll model the demand for facilities (such as schools, stores) in the physical space. The same approach applies to the demand for products (e.g. cars) in the characteristics space, see e.g. Lancaster (1966), Feenstra and Levinsohn (1995), and Berry and Pakes (2007).
We'll simulate fountain locations on a city represented by the two dimensional square.
We shall now load the libraries that we need. They are a bit specific, as they require combinatorial geometry routines. The pysdot library by Hugo Leclerc and Quentin Mérigot is still at an early stage of development, but is quite promising and easy to pip install, so we will adopt it for this course.
# !pip install pysdot
from pysdot import PowerDiagram
from pysdot.radial_funcs import RadialFuncUnit
from pysdot import OptimalTransport
from pysdot.domain_types import ConvexPolyhedraAssembly
import numpy as np
import random as rd
Consider inhabitants of a city whose geographic coordinates are $x\in\mathcal{X}=\left[0,1\right]^{2}$. More generally, $\mathcal{X}$ will be a convex subset of $\mathbb{R}^{d}$ ($d=2$ is only to fix ideas). The location of inhabitants is distributed with a density of mass $n(x)$ which is positive on $\mathcal{X}$. $n$ is assumed to have unit total mass: $\int_\mathcal{X} n(x)dx=1$, so it is a probability density function.
There are $J$ fountains, located at points $y_{j}\in\mathbb{R}^{d}$, $1\leq j\leq J$. Fountain $j$ is assumed to have capacity $q_{j}$, which means it can serve a mass $q_j$ of inhabitants. It is assumed that $\sum_{j}q_{j}=1$, which means that total supply equals the total demand.
An inhabitant at $x$ has a transportation cost associated with using fountain located at $y$ which is proportional to the square distance to the fountain
\begin{align*} \tilde{\Phi}\left( x,y\right) :=-\left\vert x-y\right\vert ^{2}/2. \label{Phistar} \end{align*}Let $\tilde{v}_{j}$ be the price charged by fountain $j$. The utility of the consumer at location $x$ is therefore $\tilde{\Phi}\left( x,y_{j}\right) -\tilde{v}_{j}$, and the indirect surplus of the consumer at $x$ is given by
\begin{align*} \tilde{u}\left( x\right) =\max_{j\in\left\{ 1,...,J\right\} }\left\{ \tilde{\Phi}\left( x,y_{j}\right) -\tilde{v}_{j}\right\} \label{ustar} \end{align*}Without loss of generality, one can replace the quadratic surplus $\tilde{\Phi}\left( x,y\right) =-\left\vert x-y\right\vert ^{2}/2$ by the scalar product surplus
\begin{align*} \Phi\left( x,y\right) :=x^{\intercal}y. \label{PhiScalProd} \end{align*}Indeed, note that $\tilde{\Phi}\left( x,y\right) =\Phi\left( x,y\right) - \left\vert x\right\vert ^{2}/2 - \left\vert y\right\vert ^{2}/2$, and introduce the reduced indirect surplus $u\left( x\right)$ and the $v_{j}$'s the reduced prices as
\begin{align*} u\left( x\right) =\tilde{u}\left( x\right) +\left\vert x\right\vert ^{2}/2\text{, and }v_{j}=\tilde{v}_{j}+\left\vert y_{j}\right\vert ^{2}/2, \label{uandv} \end{align*}One immediately sees that $\tilde{u}\left( x\right) +\tilde{v}_{j}\geq \tilde{\Phi}\left( x,y_{j}\right) $ if and only if $u\left( x\right) +v_{j}\geq\Phi\left( x,y_{j}\right) $. It follows that the consumer at location $x$ chooses fountain $j$ that maximizes
\begin{align*} u\left( x\right) =\max_{j\in\left\{ 1,...,J\right\} }\left\{ \Phi\left(x,y_{j}\right) -v_{j}\right\} . \label{PWAu} \end{align*}
Hence the problem can be reexpressed so that the surplus of consumer $x$ at fountain $j$ is simply $x^{\intercal}y_{j}-v_{j}$. It is clear from inspection that (unlike $\tilde{u}$), the reduced surplus $u$ is a piecewise affine and convex function from $\mathbb{R}^{d}$ to $\mathbb{R}$. The connection with convex and piecewise affine functions is the reason for reformulating the problem as we did.
The demand set of fountain $j$ is
\begin{align*} \mathcal{X}_{j}^{v}:=\left\{ x\in\mathcal{X}:\tilde{\Phi}\left(x,y_{j}\right) -\tilde{v}_{j}\geq\tilde{\Phi}\left( x,y_{k}\right) -\tilde{v}_{k},~\forall k\right\} \end{align*}which is equivalent to
\begin{align*} \mathcal{X}_{j}^{v}=\left\{ x\in\mathcal{X}:x^{\intercal}\left( y_{j} -y_{k}\right) \geq v_{j}-v_{k},~\forall k\right\}. \end{align*}Basic properties:
$\mathcal{X}_{j}$ is a convex polyhedron;
The intersection of $\mathcal{X}_{j}$ and $\mathcal{X}_{k}$'s lies in the hyperplane of equation $\{x:x^{\intercal}\left( y_{j}-y_{k}\right) +v_{k}-v_{j}=0\}$;
The set $\mathcal{X}_{j}$ weakly increases when $v_{k}$ ($k\neq j$) increases, and strictly decreases when $v_{j}$ decreases.
The system of sets $\left( \mathcal{X}_{j}^{v}\right) _{j}$ is called the power diagram associated to the price system $v$.
If fountains do not charge any fee, that is, if $\tilde{v}_{j}=0$, or equivalently if $v_{j}=\left\vert y_{j}\right\vert ^{2}/2$, then $\mathcal{X}_{j}^{0}$ is the set of consumers who are closer to fountain $j$ than to any other fountain. The cells $\mathcal{X}_{j}^{0}$ form a partition of $\mathcal{X}$ called Voronoi tesselation, which is a very particular case of a power diagram. The Voronoi diagrams have the property that fountain $j$ belongs to cell $\mathcal{X}_{j}^{0}$; when $\tilde{v}\neq0$, this property may no longer hold for more general power diagrams.
Example. We will generate a Voronoi tesselations where $10$ fountains are distributed uniformly on $[0,1]^2$, and $\tilde{v} = 0$.
rd.seed(777)
nCells = 10
Ys = np.random.uniform(0,1,2*nCells).reshape((-1,2))
vor_dia = PowerDiagram(Ys)
vor_dia.display_jupyter()
2
The demand for fountain $j$ is given by $\mathbb{P}_n\left(\mathcal{X}_{j}\right) = \int_\mathcal{X} 1\{x\in\mathcal{X}_{j}\} n(x) dx $ where $\mathbb{P}_n$ is the probability distribution of consumer locations.
Note that in general $x^{\intercal}y_{j}-u\left( x\right) \leq v_{j}$; yet if consumer $x$ chooses fountain $j$, then this inequality holds as an equality. Hence, the set of consumer who prefer fountain $j$ is given by
\begin{align*} \mathcal{X}_{j}=\arg\max_{x\in\mathcal{X}}\left\{ x^{\intercal}y_{j}-u\left( x\right) \right\} \label{defXj} \end{align*}By first order conditions $x\in\mathcal{X}_{j}$ if and only if $\nabla u\left(x\right) = y_{j}$ (assuming $u$ is differentiable at $x$). Therefore
\begin{align*} \mathcal{X}_{j}:=\nabla u^{-1}\left( \left\{ y_{j}\right\} \right) . \label{Demand} \end{align*}Fountain example. We see in the picture above that cells have different areas. The areas of the cells are given by:
vor_dia.integrals()
array([0.06924676, 0.14662215, 0.0255011 , 0.22922231, 0.04548061,
0.04697387, 0.06799278, 0.26857535, 0.07354108, 0.02684399])
Introduce the social welfare of producers and consumers as
\begin{align*} S\left( v\right) :=\sum_{j}q_{j}v_{j}+\mathbb{E}_{\mathbb{P}_n}\left[ \max_{j\in\left\{ 1,...,J\right\} }\left\{ X^{\intercal}y_{j}-v_{j}\right\} \right] . \label{defSbis} \end{align*}
We have
\begin{align*} \frac{\partial S\left( v\right) }{\partial v_{k}}=q_{k}-\mathbb{E}% _{\mathbb{P}_n}\left[ 1\left\{ \nabla u\left( X\right) =y_{k}\right\} \right] =q_{k}-\mathbb{P}_n\left( \mathcal{X}_{k}^{v}\right) . \end{align*}Thus, the excess supply for fountain $j$ is given by
\begin{align*} q_{k}-\mathbb{P}_n\left( \mathcal{X}_{k}^{v}\right) =\frac{\partial S\left( v\right) }{\partial v_{k}} \label{exprDemand}% \end{align*}where $S$ is defined by the social welfare above.
Hence, market clearing prices, or equilibrium prices are prices $v$ such that demand and supply clear, that is, such that $q_{k}=\mathbb{P}_n\left(\mathcal{X}_{k}^{v}\right)$ for each $k$; in other words
\begin{align*} \frac{\partial S\left( v\right) }{\partial v_{k}}=0. \end{align*}The central planner may decide arbitrarily on assigning to each inhabitant $x$ a fountain $T\left( x\right) \in\left\{ y_{1},...,y_{J} \right\} $, in a such way that each fountain $j$ is used to its full capacity, that is
\begin{align*} \mathbb{P}_n\left( T\left( X\right) =y_{j}\right) =q_{j},~\forall j\in\left\{ 1,...,J\right\} . \label{massBalance}% \end{align*}The planner seeks to maximize the total surplus subject to capacity constraints; hence
\begin{align*} & \max\mathbb{E}_{\mathbb{P}_n}\left[ X^{\intercal}T\left( X\right) \right] \label{welfare}\\ & s.t. P\left( T\left( X\right) =y_{j}\right) =q_{j},~\forall j\in\left\{ 1,...,J\right\} \end{align*}This is a Monge problem, whose Kantorovich relaxation is
\begin{align*} \max_{\mu\in\mathcal{M}\left( \mathbb{P}_n,q\right) }\mathbb{E}_{\mu}\left[X^{\intercal}Y\right] . \end{align*}By the Monge-Kantorovich theorem, the dual problem is
\begin{align*} & \min_{u,v}\mathbb{E}_{\mathbb{P}_n}\left[ u\left( X\right) \right] +\mathbb{E} _{q}\left[ v\left( Y\right) \right] \label{dualKantoContDiscr}\\ & s.t. u\left( x\right) +v\left( y\right) \geq x^{\intercal}y, \end{align*}
where the constraint should hold almost surely with respect to $P$ and $Q$.
The constraint should be verified for $y\in\left\{ y_{1},...,y_{J}\right\} $, and the constraint+optimality implies $u\left( x\right) =\max_{j\in\left\{ 1,...,J\right\} }\left\{ \Phi\left( x,y_{j}\right) -v_{j}\right\} $. Thus, the dual problem rewrites as
\begin{align*} \min_{v\in\mathbb{R}^{J}}\mathbb{E}_{\mathbb{P}_n}\left[ \max_{j\in\left\{1,...,J\right\} }\left\{ X^{\intercal}y_{j}-v_{j}\right\} \right] +\sum_{j=1}^{J}q_{j}v_{j} \label{MKfiniteDim} \end{align*}which is the minimum of $S$ over $v\in\mathbb{R}^{J}.$
As a result:
There exist equilibrium prices, which are the minimizers of $S$.
The total welfare at equilibrium coincides with the optimal welfare.
Note that
\begin{align*} \arg\max_{j\in\left\{ 1,...,J\right\} }\left\{ \Phi\left( x,y_{j}\right) -v_{j}\right\} \end{align*}is a singleton for almost every $x$ (it is not a singleton when $x$ is at the boundary between two cells). The assumption that $P$ is absolutely continuous is crucial here.
Hence the map
\begin{align*} T\left( x\right) =\nabla u\left( x\right) \end{align*}is defined almost everywhere and coincides with $\arg\max$ whenever it is defined. Thus the solution does not involve to split mass.
We turn to a discussion on the numerical determination of the prices (we discuss the determination of the $v$'s, as the expression for the $w$'s immediately follows). The function $S$ to minimize being convex, we can use a standard gradient descent algorithm in which the increase in prices is given by
\begin{align*} v_{j}^{t+1}-v_{j}^{t}=\varepsilon\left( \mathbb{P}_n\left( \nabla u\left( X\right) =y_{j}\right) -q_{j}\right) , \label{tatonnement}% \end{align*}which has immediately an economic interpretation: the fountains that are over-demanded raise their prices, while the fountains that are under-demanded lower their prices. This a tâtonnement process.
Algorithm Take an initial guess of $v^{0}$. At step $t$, define $v^{t+1}$ by
\begin{align*} v_{j}^{t+1}=v_{j}^{t}-\varepsilon_{t}\frac{\partial S}{\partial v_{j}}\left(v^{t}\right), \end{align*}for $\varepsilon_{t}$ small enough. Stop when $\frac{\partial S}{\partial v_{j}}\left( v^{t+1}\right) $ is sufficiently close to zero.
The gradient descent method is implemented as follows.
rel_tol = 1e-4
q_j = np.ones(nCells)/nCells
vtilde_j = np.zeros(nCells)
cont = True
pow_dia = PowerDiagram(Ys,vtilde_j)
while cont:
demand_j = pow_dia.integrals()
if ((demand_j - q_j)/q_j).max()<rel_tol:
cont=False
else:
vtilde_j = vtilde_j - 0.1 * (demand_j - q_j)
pow_dia.set_weights(vtilde_j)
pow_dia.display_jupyter()
2
We can display the prices:
print(vtilde_j - vtilde_j[0])
[ 0. 0.01314871 0.04433732 -0.06718019 0.11014237 0.11008925 0.08778759 -0.13640295 0.04872085 0.12719081]
Alternatively, one could have used the following routine:
def make_square(box=[0, 0, 1, 1]):
density = ConvexPolyhedraAssembly()
density.add_box([box[0], box[1]], [box[2], box[3]])
return density
dens = make_square()
def optimal_transport(density, Y, masses, psi0 = None, err=1e-8):
center = (density.min_position() + density.max_position())/2
halfsides = (density.max_position() - density.min_position())/2
ratio = 1/np.max(np.abs((Y-center)/halfsides))
psi = (1-ratio)*np.sum((Y-center)**2, axis=-1)
ot = OptimalTransport(Y, psi, density, radial_func=RadialFuncUnit(), obj_max_dw=err)
ot.set_masses(masses)
ot.adjust_weights()
return ot.pd
pow_dia2 = optimal_transport(dens, Ys, q_j)
vtilde_j2 = pow_dia2.weights
pow_dia2.display_jupyter()
print(vtilde_j2 - vtilde_j2[0])
print((vtilde_j2 -vtilde_j + vtilde_j[0]- vtilde_j2[0]).max())
[ 0. 0.01316282 0.04434139 -0.0671193 0.11013903 0.11015548 0.08778734 -0.13640282 0.04878224 0.12725911] 6.829379990777595e-05
Exercise. Implement a coordinate descent version of the algorithm.