<jupyter_start><jupyter_text>Block 6: Semi-discrete optimal transport A<NAME> (NYU) `math+econ+code' masterclass on matching models, optimal transport and applications© 2018-2019 by <NAME>. Support from NSF grant DMS-1716489 is acknowledged. <NAME> contributed.  Optimal transport I: Semidiscrete matching  References* [OTME], Ch. 5* Anderson, <NAME>, 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*.* <NAME> (2007). The Pure Characteristics Demand Model. *IER*.* <NAME> (1995). Estimating Markups and Market Conduct with Multidimensional Product Attributes. *ReStud*.* <NAME> (2017). Yoghurts Choose Consumers. Identification of Random Utility Models via Two-Sided Matching. *Mimeo*.  Motivation and setting   MotivationToday 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.  SettingConsider 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 according to a absolutely continuous distribution $P$ whose support is included on $\mathcal{X}$, normalized to $1$, so $P$ is a probability distribution.There are $M$ fountains, located at points $y_{j}\in\mathbb{R}^{d}$, $1\leq j\leq M$. Fountain $j$ is assumed to have capacity $q_{j}$, and 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 $j$ 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,...,M\right\}  }\left\{\tilde{\Phi}\left(  x,y_{j}\right)  -\tilde{v}_{j}\right\} \label{ustar}\end{align*}  A reformulationWithout 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,...,M\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 [this](PWAu) 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 forreformulating the problem as we did.  Power DiagramsThe 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$.  Voronoi tesselationsIf 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. We will generate a Voronoi tesselations where $10$ fountains are distributed uniformly on $[0,1]^2$, and $\tilde{v} = 0$.<jupyter_code>library("magick")  # For plotting animations
library("transport")
library("geometry")

SEED = 777
MAX_ITER = 1000
PREC = 0.01

set.seed(SEED)
nCells = 10

y1 = runif(nCells)  # coordinate 1
y2 = runif(nCells)  # coordinate 2
vtilde = rep(0, nCells)
q = rep(1/nCells, nCells)
demand = rep(0, nCells)
eps = 0.1  #tuning parameter

pwd = power_diagram(y1, y2, vtilde, rect = c(0, 1, 0, 1))

plot(pwd, weights = FALSE)<jupyter_output>Warning message:
"package 'magick' was built under R version 3.5.2"Linking to ImageMagick 6.9.9.14
Enabled features: cairo, freetype, fftw, ghostscript, lcms, pango, rsvg, webp
Disabled features: fontconfig, x11
Warning message:
"package 'transport' was built under R version 3.5.2"Warning message:
"package 'geometry' was built under R version 3.5.2"Loading required package: magic
Warning message:
"package 'magic' was built under R version 3.5.2"Loading required package: abind
Warning message:
"package 'abind' was built under R version 3.5.2"<jupyter_text>Demand zone of a fountainThe demand for fountain $j$ is given by $P\left(\mathcal{X}_{j}\right)  =\Pr\left(  X\in\mathcal{X}_{j}\right)  $ where $X\sim P$, which is the mass of consumers who prefer fountain $j$ over the others.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*}  Equilibrium prices: Aurenhammer's methodIntroduce the social welfare of producers and consumers as\begin{align*}S\left(  v\right) :=\sum_{j}q_{j}v_{j}+\mathbb{E}_{P}\left[  \max_{j\in\left\{  1,...,M\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}%_{P}\left[  1\left\{  \nabla u\left(  X\right)  =y_{k}\right\}  \right]=q_{k}-P\left(  \mathcal{X}_{k}^{v}\right)  .\end{align*}Thus, the excess supply for fountain $j$ is given by\begin{align*}q_{k}-P\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 welfate](defSbis) above.Hence, market clearing prices, or equilibrium prices are prices $v$ such that demand and supply clear, that is, such that $q_{k}=P\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*}  Central planner's problemThe central planner may decide arbitrarily on assigning to each inhabitant $x$ a fountain $T\left(  x\right)  \in\left\{  y_{1},...,y_{M} \right\}  $, in a such way that each fountain $j$ is used to its full capacity, that is\begin{align*}P\left(  T\left(  X\right)  =y_{j}\right)  =q_{j},~\forall j\in\left\{1,...,M\right\}  . \label{massBalance}%\end{align*}The planner seeks to maximize the total surplus subject to capacity constraints; hence\begin{align*}&  \max\mathbb{E}_{P}\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,...,M\right\}\end{align*}This is a Monge problem, whose Kantorovich relaxation is\begin{align*}\max_{\pi\in\mathcal{M}\left(  P,Q\right)  }\mathbb{E}_{\pi}\left[X^{\intercal}Y\right]  .\end{align*}  DualityBy the Monge-Kantorovich theorem, the dual problem is\begin{align*}&  \min_{u,v}\mathbb{E}_{P}\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_{M}\right\}  $, and the constraint+optimality implies $u\left(  x\right) =\max_{j\in\left\{  1,...,M\right\}  }\left\{  \Phi\left(  x,y_{j}\right) -v_{j}\right\}  $. Thus, the [dual problem](dualKantoContDiscr) rewrites as\begin{align*}\min_{v\in\mathbb{R}^{M}}\mathbb{E}_{P}\left[  \max_{j\in\left\{1,...,M\right\}  }\left\{  X^{\intercal}y_{j}-v_{j}\right\}  \right] +\sum_{j=1}^{M}q_{j}v_{j} \label{MKfiniteDim}\end{align*}which is the minimum of $S$ over $v\in\mathbb{R}^{M}.$As a result:1. There exist equilibrium prices, which are the minimizers of $S$.2. The total welfare at equilibrium coincides with the optimal welfare.  Monge vs KantorovichNote that\begin{align*}\arg\max_{j\in\left\{  1,...,M\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 Monge and the Kantorovich problems have the same solutions.  Determination of pricesWe 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 $F$ 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(  P\left(  \nabla u\left(  X\right)=y_{k}\right)  -q_{k}\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.  ComputationThe prices adjust over time so that the area of the demand cells adjust to supply. For each cell j, set `cellj=pwd$cells[[j]]`, and demand is computed using<jupyter_code>for (j in 1:nCells) {
    cellj = pwd$cells[[j]]
    demand[j] = polyarea(cellj[, 1], cellj[, 2])
}
print(demand - 0.1)<jupyter_output>[1] -0.030069573 -0.016570247 -0.038802835 -0.043641234  0.002014428
 [6] -0.024359557  0.035695538  0.025425664  0.077657409  0.012650408<jupyter_text>where `polyarea`, from the ‘geometry` package, computes the area of a cell. Finally, the prices are adjusted using:<jupyter_code>vtilde = vtilde - eps * (demand - q)<jupyter_output><empty_output><jupyter_text>where `eps` is a tuning parameter which controls the speed of adjustment of the process, and the configuration converges to equilibrium. Running this until convergence<jupyter_code># Loop it

# create canvas
frames = image_graph(width = 800, height = 800, res = 150)

t = 1
cont = TRUE
while ((cont == TRUE) && (t < MAX_ITER)) {
    if (t%%10 == 0) {
        # Print on the screen some message
        cat(paste0("iteration: ", t, "\n"))
    }
    for (j in 1:nCells) {
        cellj = pwd$cells[[j]]
        demand[j] = polyarea(cellj[, 1], cellj[, 2])
    }
    if (max(abs(demand - q)) < PREC/nCells) {
        cont = FALSE
    } else {
        t = t + 1
        vtilde = vtilde - 0.1 * (demand - q)
        pwd = power_diagram(y1, y2, vtilde, rect = c(0, 1, 0, 1))
        plot(pwd, weights = FALSE)
    }
}

# animate
animation = image_animate(frames, 10)
image_write(animation, "powerdiagram.gif")
animation<jupyter_output>iteration: 10
iteration: 20
iteration: 30
iteration: 40
iteration: 50
{'fraction_non_alphanumeric': 0.1660135375846099, 'fraction_numerical': 0.02144638403990025, 'mean_word_length': 4.50647312671636, 'pattern_counts': {'":': 0, '<': 23, '<?xml version=': 0, '>': 21, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '411695', 'license': 'AFL-1.1', 'n_tokens_mistral': 4994, 'n_tokens_neox': 4639, 'n_words': 1517, 'path': 'slides_ipynb_mec_optim/B06_semidiscretetransport.ipynb', 'repo_name': 'math-econ-code/optim_2019-06', 'size': '14095'}