Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
212,354 characters · 67 sections · 338 citation commands
An econometrician's guide to optimal transport
\address{NYU and Penn State}
\nocite{*}
Optimal transport traces its origins to 18th Century civil engineering. It draws its central ideas from several diverse areas of mathematics. It finds application in fields as diverse as non-Euclidean geometry, cosmology, fluid dynamics, traffic management, machine learning, economics and computational biology. There is a rapidly growing community of researchers who work {\em on} optimal transport as well as researchers who work {\em with} optimal transport. There is also a growing body of reference books, textbooks and surveys about various aspects of optimal transport theory and its applications. Among the latter, in economics alone, optimal transport has wide ranging applications, from family economics to mechanism design, with recent books and surveys to account for them. The objective of this new review is to give an account of optimal transport theory and methods specifically designed for applications in econometrics. The target audience includes all producers and consumers of econometric methods. Econometric theorists will find here new proof techniques, new computational tools, new ways to model phenomena and new ways to perform inference. Applied economists will find here many new econometric techniques applicable to a wide range of fields with step-by-step guides to implementation.
We start with a brief overview of optimal transport designed mostly for readers who were previously unacquainted with the subject. We provide a gentle general introduction, which includes a short history of the subject, the basic formulations of the problem, a short description of the main aspects of the theory that are relevant to econometrics as well as its classical computational tools. We then devote the second section to describing a collection of algorithms that are the main building blocks in most of the procedures presented in this review. We propose a simplex-based algorithm to derive optimal matchings, which are optimal transport plan between two discrete distributions. We derive the optimal transport map between two Gaussian distributions. We present the solution to optimal transport in $\mathbb R$ and their relation with rearrangement inequalities and quantiles. We describe algorithms from computational geometry to compute optimal transport maps from a discrete distribution to a uniform in $[0,1]^d$, and the iterated proportional fitting procedure to compute entropic relaxations of the optimal transport problem. Finally, we discuss methods to compute the Wasserstein distance between two distributions, including recent work on Wasserstein generative adversarial networks.
We then proceed to the review of econometric papers, where the contribution relies in a significant way on optimal transport theory. The applications of optimal transport we review here range across most areas of econometrics, including causal inference, discrete choice, structural estimation, data combination, program evaluation and treatment assignment, risk and inequality measurement, nonparametric identification, mis-specification and robustness, distributional regression, and machine learning. We organize them not by area but according to the aspect of optimal transport theory that underlies the main innovation in each paper. We identify three main aspects of the theory as well as several extensions. The main aspects of the theory we identify are: First, the existence of optimal transport plans and Kantorovich duality theory; Second, the uniqueness of optimal transport maps and generalized monotonicity; Third, optimal transport distance between distributions. They correspond broadly, but not exactly with: First data combination applications; Second, multivariate quantiles and identification; Third, distributional robustness. For each aspect of the theory, we explain how it is useful in econometrics, we review econometric innovations based on this aspect, and finally describe in full detail a few contributions that we find particularly representative or salient with a guide to implementation, including algorithms and code.
Optimal transport has natural connections with the economics of matching. General applications to economic theory are beyond the scope of this review. However, an aspect of the economics of matching is directly related to econometric theory and practice. The problem of recovering transport costs from realized transport plans, called inverse optimal transport problem, is directly relevant to the identification and estimation of matching surpluses and complementarities in empirical matching theory. In section (ref), we present the inverse optimal transport problem; we review its applications in empirical matching and its connections with parametric estimation of matching, Poisson regressions, gravity equations and metric learning.
For those readers who wish to work with optimal transport theory to develop new econometric methods, this review is an introduction, not a substitute for more advanced texts. Beyond the pioneering books of rachev1998mass [rachev1998mass,rachev1998mass-app], which are treasure troves of useful results, but somewhat difficult to navigate, we recommend villani2003topics [villani2003topics,villani2008optimal], santambrogio2015optimal, galichon2016optimal, peyre2019computational and chewi2024statistical. villani2003topics is the book that ushered us both into the world of optimal transport. It's hard to overstate its elegance and influence. It provides a concise overview of the whole theory up to that point, suitable for any audience with a modicum of mathematical sophistication. It is, and is likely to remain, the best introduction to optimal transport theory. villani2008optimal expands considerably on villani2003topics. For the purposes of applications to econometrics, it is particularly valuable for the deep insights it provides into solutions to the Monge problem, i.e., existence, uniqueness and cyclical monotonicity of optimal transport maps. santambrogio2015optimal is also a very valuable resource. Unlike villani2008optimal which is primarily designed as a reference for researchers who wish to work on optimal transport and advance the theory, santambrogio2015optimal targets researchers who wish to work with optimal transport, such as applied mathematicians and econometricians among others. galichon2016optimal introduces optimal transport to economists as a unifying tool for modeling, identification, and computation, and sketches its core applications across matching, equilibrium analysis, and empirical economics; it is addressed to a general economics, rather than econometrics, audience, and it predates many of the examples we review. peyre2019computational is a great resource on all computational aspects of optimal transport. Finally chewi2024statistical contains an account of recent developments on the estimation of optimal transport plans and maps and optimal transport distances, which are particularly relevant to applications in econometrics. More recently, Galichon2026DCM focuses on discrete choice models in connection with multiple aspects of optimal transport.
In all this review, unless specified otherwise, the results are valid when the bases spaces $\mathcal X$ and $\mathcal Y$ are Polish spaces. However, in most applications, they are subsets of $\mathbb R^d$.
The origin of optimal transport theory is always traced back to monge1781memoire, with the original stylized formulation of the problem of moving specific piles of dirt to build specific structures at minimum cost. Its author was Gaspard Monge, influential scientist, engineer and politician during the French Revolution and Empire, and one of the three founding fathers of the \'Ecole polytechnique in Paris. The Russian mathematician Leonid Kantorovich rediscovered\footnote{Kantorovich made the connection with the Monge problem in 1948.} the problem in 1938 and proposed a probabilistic formulation with a duality theorem in kantorovich1940 [kantorovich1940,kantorovich1942translocation], which laid the foundations of linear programming methods, and for which he was awarded the 1975 Nobel Prize for economics, jointly with Tjalling Koopmans, “for their contributions to the theory of optimum allocation of resources”.\footnote{The American mathematician Frank Hitchcock independently formulated a special case of the optimal transport problem in hitchcock1941distribution, together with an early version of the simplex algorithm to solve it.} He also used optimal transport to define a distance between probability measures, which, through the vagaries of misattribution, is now commonly known as Wasserstein distance. The solution to the original Monge problem, when the transportation cost is quadratic, was given by brenier1987decomposition [brenier1987decomposition,brenier1991polar] and ruschendorf1990characterization based on earlier work by sudakov1979geometric and knott1984optimal in particular. They show that solutions to the optimal transport problem with quadratic cost exist, are unique and are cyclically monotone (gradients of convex functions). mccann1995existence extended this result by showing that two probability measures can always be mapped into each other with a cyclically monotone map, under minimal regularity and without the finite variance condition necessary to define quadratic optimal transport. Chapter 10 of villani2008optimal gives comprehensive results on existence and uniqueness of solutions to optimal transport with more general costs based on work by mccann2001polar and bernard2007optimal. More recent developments that are relevant to econometric applications are displacement convexity in mccann1997convexity, multi-marginal optimal transport in pass2015multi, and weak optimal transport in gozlan2017kantorovich. Optimal transport theory and its applications are now in an explosive phase of development, one of the latest twists being the massive interest in Wasserstein barycenters and Wasserstein Generative Adversarial Networks in machine learning and artificial intelligence.
Optimal Transport is a coupling problem. It is the problem of finding a coupling of two probability measures to minimize an aggregate cost. Let $\mu$ and $\nu$ be probability measures on $\mathcal X$ and $\mathcal Y$ respectively. A coupling of $(\mu,\nu)$ is the construction of a pair of random elements $(X,Y)$ with probability distribution $\pi$ on $\mathcal X\times\mathcal Y$ such that $X$ has distribution $\mu$ and $Y$ has distribution $\nu$. The probability distribution $\pi$ is then said to have marginals $\mu$ and $\nu$, which is equivalent to
for all integrable functions $\varphi$ and $\psi$ on $\mathcal X$ and $\mathcal Y$ respectively. By extension, a distribution $\pi$ with marginals $\mu$ and $\nu$ is also called a coupling of $(\mu,\nu)$ and the collection of all such couplings is denoted $\mathcal M(\mu,\nu)$.
The coupling is called deterministic if there is a measurable function $T:\mathcal X\rightarrow\mathcal Y$ such that $Y=T(X)$. The function $T$ is then called a change of variables, which is equivalent to
for all integrable function $\psi$ on $\mathcal Y$.
The optimal transport problem minimizes an aggregate cost. Call $c: \mathcal X\times\mathcal Y\rightarrow \mathbb R$ the cost function. Optimal transport has two main formulations.
The Kantorovich formulation of the optimal transport problem is that of finding a coupling $(X,Y)$ of $(\mu,\nu)$ with distribution $\pi$ that minimizes the aggregate cost
Such a distribution $\pi$ with marginals $\mu$ and $\nu$ is called an optimal transport plan from distribution $\mu$ to distribution $\nu$.
The Monge formulation of the optimal transport problem is that of finding a deterministic coupling $(X,T_{\mu\rightarrow\nu}(X))$ of $(\mu,\nu)$ that minimizes the aggregate cost
The map $T_{\mu\rightarrow\nu}$ that defines such a deterministic coupling of $(\mu,\nu)$ is called an optimal transport map from distribution $\mu$ to distribution $\nu$.
The main aspects of optimal transport theory that are relevant to econometrics are related to the following basic questions about optimal transport plans, optimal transport maps and the value of the optimal transport objective.
The answers to these questions may naturally depend on the cost function $c$ and the marginal distributions $\mu$ and $\nu$. Conversely, we may also ask:
The answer to question (1) is straightforward for optimal transport plans, which exist in great generality and are generically non unique. Optimal transport plans exist under the following condition on the cost function, which we will assume in the rest of this review.
Existence of an optimal transport plan holds because $\mathcal M(\mu,\nu)$ is compact, and, as C\'edric Villani puts it, “has probably been known since time immemorial” villani2008optimal. The answer to question (1) for optimal transport maps is much more involved and is one of the most significant aspects of optimal transport theory. Existence of optimal transport maps is much trickier because the set of deterministic couplings is not compact. Proofs of existence and uniqueness of optimal transport maps for certain classes of transport cost functions and under certain regularity conditions on the probability distributions $\mu$ and $\nu$ rely on some of the answers to question (2), so we will return to them later in this section.
A first characterization of optimal transport plans and answer to question (2) is the equivalence between cyclical monotonicity and optimality of a transport plan. Equivalence between cyclical monotonicity and optimality can be understood in the following way. Recall that a transport plan $\pi$ is a joint distribution on $\mathcal X\times\mathcal Y$. If a pair $(x,y)$ is in the support of $\pi$, then some mass is being transferred from $x$ to $y$. If $\pi$ is an optimal transport plan, then the mass transferred from $x$ to $y$ cannot be more cheaply transferred via $x\rightarrow y^\prime\rightarrow x^\prime\rightarrow y$, say. More generally, optimality implies the property of cyclical monotonicity of the support of $\pi$.
Conversely, it can be shown that cyclical monotonicity of the support implies optimality of the transport plan under assumption (ref). Equivalence of cyclical monotonicity and optimality, true for real-values and lower semicontinuous cost functions is proved in schachermayer2009characterization. See equivalence of (a) and (b) in Theorem 5.10(ii) of villani2008optimal.
Cyclical monotonicity also relates to the central answer to question (2) and a major aspect of optimal transport theory, namely Kantorovich duality, which also holds in great generality. The dual Kantorovich problem is that of maximizing loading and unloading costs
under the constraint that loading cost $\varphi(x)$ at $x$ and unloading cost $\psi(y)$ at $y$ are continuous and bounded and do not exceed the total transportation cost $c(x,y)$. Hence the supremum in the dual Kantorovich problem is taken over the set
Since all transport plans are couplings of $(\mu,\nu)$, integrating the constraint yields
Hence the Kantorovich dual $\tilde{\mathcal C}(\mu,\nu)$ is no larger than the primal $\mathcal C(\mu,\nu)$ (weak duality). It turns out they are equal (strong duality) under very general conditions, which is the content of the Kantorovich duality theorem.
The dual constraint can be rewritten $\psi(y)\leq c(x,y)-\varphi(x)$. Hence, function $\psi$ can be taken equal to
which is called the $c$-conjugate of $\varphi$. Similarly, the dual objective can be further improved if we replace $\varphi$ with
which is the double conjugate of $\varphi$. If we iterate the procedure, the pair $(\varphi^{cc},\varphi^c)$ stays unchanged. Hence, for a pair $(\varphi,\psi)$ to achieve the dual optimal value, $\psi=\varphi^c$ and $\varphi$ must be $c$-concave in the sense of the following definition, due to ruschendorf1990characterization.
Then, if $\varphi$ is indeed $c$-concave, and $(\varphi,\varphi^c)$ achieves the dual optimal value, the set where the constraints bind
is the cyclically monotonic support of an optimal transport plan, i.e., of a solution to the primal Kantorovich optimal transport problem. It is called the $c$-subdifferential of $\varphi$, which explains the notation. The terminology will become clear when we consider the special case of quadratic cost later in this subsection. With definition (ref) and ((ref)), we can state the second part of the Kantorovich duality, which gives necessary and sufficient conditions for optimality.
This brings us back to the second part of question (1): when do optimal transport maps $T$ exist? In other words, when is there an optimal coupling $\pi$ that is deterministic? Heuristically, for this to happen, the support of an optimal transport plan must be concentrated on a curve. Since the support of an optimal transport plan is $\partial_c\varphi$ for some $c$-concave function $\varphi$, for it to be concentrated on a curve, the set $\partial_c\varphi(x):=\{ y\in\mathcal Y: \varphi(x)+\varphi^c(y)=c(x,y)\}$ must be reduced to a single $y$ for almost all $x$. Since the pair $(\varphi,\varphi^c)$ satisfies the dual constraint, for any $x^\prime$, we have $\varphi(x^\prime)+\varphi^c(y)\leq c(x^\prime,y)$. Eliminating $\varphi^c(y)$ yields $\varphi(x^\prime) - \varphi(x) \leq c(x^\prime,y) - c(x,y).$ Since this is true for all $x^\prime$ approaching $x$, assuming suitable differentiability, we obtain
The latter has a unique solution $y$, as desired, if $\nabla_xc(x,y)$ in injective:
Assumption (ref) rules out, for instance, additive costs of the form $c(x,y)=a(x)+b(y)$, which would lead to multiple optimal couplings. In order for equation ((ref)) to make sense, we need differentiability of $\varphi$. All we know about $\varphi$ is that it is $c$-concave, as part of dual optimal pair. Hence, a sufficient condition is the ($\mu$-almost sure) differentiability of all $c$-concave functions, which brings us to formulate:
The conditions of theorem (ref) are satisfied\footnote{Note that theorem (ref) does not apply to distance costs $c(x,y)=d(x,y)$.} in the special case, where the cost function takes the form $c(x,y)=h(x-y)$ with $h$ strictly convex. From ((ref)), the optimal transport map then takes the special form $y=T(x)=x-\left(\nabla h\right)^{-1}\left(\nabla\varphi(x)\right)$. The even more special case of quadratic cost $c(x,y)=\|x-y\|^2/2$ is particularly important both historically and in terms of applications. In that case, the optimal transport map takes the form $y=T(x)=x-\nabla\varphi(x)$, which is the gradient of the convex function $\vartheta (x) := x^\top x/2-\varphi(x)$.
The first part of theorem (ref) says that for any two probability distributions $\mu$ and $\nu$ on $\mathbb R^d$, if $\mu$ is absolutely continuous, there is a unique gradient of a convex function that pushes $\mu$ to $\nu$ in the sense of the change of variables formula ((ref)). When the optimal transport problem from $\mu$ to $\nu$ is well defined, the second part of theorem (ref) states that this gradient of a convex function is the unique optimal transport map from $\mu$ to $\nu$ from theorem (ref). This relates to a classical result in convex analysis, which states that gradients of convex functions are characterized by cyclical monotonicity of their graph\footnote{With strictly convex costs in $\mathbb R^n$, $\nabla T$ is diagonalizable with nonnegative eigenvalues. See the remark on the bottom of page 278 of villani2008optimal.}. See for instance section 24 of rockafellar1970convex.
Turning now to question (3), we examine the aspects of optimal transport theory relating to measuring discrepancy between probability distributions. The optimal transport cost between two probability distributions can be used to measure the discrepancy between the two. More precisely, optimal transport can be used to define a family of distances on the space of probability distributions called Wasserstein distances.
Definition (ref) is a theorem-definition, in that it implies that the optimal transport functional $W_p$ satisfies the properties of a distance, i.e., it is non negative, it equals zero if and only if the two probability distributions are identical, and it satisfies the triangle inequality. In case $p=1$, which corresponds to cost function $c(x,y)=\|x-y\|$, the $c$-concave functions (definition (ref)) are the 1-Lipschitz functions. Hence, the Kantorovich dual takes the special form
which is called Kantorovich-Rubinstein dual. Wasserstein distance $W_1$ was historically named Kantorovich-Rubinstein distance, and is now commonly called Earth Mover's Distance in computer science. The latter name refers to the fact that Wasserstein distances, by definition of optimal transport, define the distance with the amount of mass that needs to be shifted to move from one distribution to another. This is in sharp contrast with traditional distances between absolutely continuous probability distributions that aggregate the difference between densities at each point.
Under the conditions of theorem (ref), there is a unique optimal transport map $T$ and a deterministic coupling $(X,T(X))$ of $(\mu,\nu)$ that solves the optimal transport problem of $\mu$ to $\nu$. The conditions of theorem (ref) are satisfied in particular if $\mu$ is absolutely continuous and the cost function is $c(x,y)=\|x-y\|^p$, $p>1$. Hence, the Wasserstein distance has the following expression:
Wasserstein distances $W_p$ induce a notion of convergence of probability distributions.
A first observation from theorem (ref) is that convergence in $W_p$ implies convergence in $W_q$ when $p>q$, and convergence in $W_1$ is the weakest. The second observation is that convergence in $W_p$ is stronger then traditional convergence in distribution, since it also requires convergence of moments of order $p$. Another feature of Wasserstein distances is that they satisfy the Kantorovich duality theorem, which provides many computational and technical advantages. The space of probability distributions also inherits geometric features from the base space $\mathcal X$. A first way to see this is to compare Wasserstein distance with the total variation distance $\mbox{TV}(\mu,\nu)=\sup_{A\subset\mathcal X}|\mu(A)-\nu(A)|$. It can be shown by Kantorovich duality that $\mbox{TV}(\mu,\nu)$ is the optimal cost of transporting $\mu$ to $\nu$ with cost function $c(x,y)=1\{x\ne y\}$. Hence, none of the geometry of $\mathcal X$ is preserved, since the indicator $1\{x\ne y\}$ says nothing about how far apart $x$ and $y$ are. In contrast, the Wasserstein distance $W_p(\delta_x,\delta_y)$ between the two Dirac masses at $x$ and $y$ is equal to the distance between $x$ and $y$ themselves. More generally, the Wasserstein distances endow the space of probability distributions with a rich geometry, that involves notions of geodesics, interpolation, convexity, and barycenters.
Among geometric concepts, interpolation and barycenters are the most directly relevant to econometric applications. Most of the relevant theory on barycenters is with respect to the $2$-Wasserstein distance $W_2$ on $\mathcal W_2:=\{\mu:{\textstyle\int_{\mathcal X}}\|x\|^2\;d\mu(x)<\infty\}$. A barycenter of $K$ probability distributions $\mu_1,\ldots,\mu_K$ on $\mathcal X_1\subset\mathbb R^d,\ldots,\mathcal X_K\subset\mathbb R^d$ respectively, is defined, in analogy with the Euclidean case, as the minimizer
given a collection $\lambda_1,\ldots,\lambda_K$ of non negative weights that sum to $1$. Let $x$ be a point in $\mathcal X_1\times\ldots\times\mathcal X_K$. In analogy with the fact that the Euclidean barycenter $\Sigma_{k=1}^K\lambda_kx_k$ of $(x_1,\ldots,x_K)$ satisfies
it can be shown that the Wasserstein barycenter is obtained from the following optimization problem:
Problem ((ref)) is called a multi-marginal optimal transport problem with cost function $c(x_1,\ldots,x_K)=\Sigma_{k=1}^K\|x_k-\Sigma_{j=1}^K\lambda_jx_j\|^2$. In the particular case of problem ((ref)), by theorem (ref), there is a unique solution, which yields a useful characterization of the Wasserstein barycenter.
This section briefly describes the main variants of the duality of optimal transport. First, entropic optimal transport is the case, where the primal Kantorovich formulation contains an entropic regularization that penalizes departures of the optimal transport plan from the independent coupling. Second, unbalanced and weak optimal transport relax the sharp constraint on the marginals. Finally, multi-marginal optimal transport extends optimal transport to the problem of finding a coupling (or {\em teaming}) of more than two marginals, while minimizing an aggregate cost.
Entropic optimal transport is a regularization of the Kantorovich formulation of the optimal transport problem that was originally proposed for computational purposes. We discuss this in section (ref). However, some theoretical features of entropic optimal transport are also useful in econometric applications. We describe the basic formulation and duality in what follows.
Let $H(\pi\vert \pi^\prime)$ be the relative entropy between probability distributions $\pi$ and $\pi^\prime$, defined as follows.
Entropic optimal transport solves the following problem:
Entropic optimal transport penalizes joint distributions that are far from the independent coupling $\mu\otimes\nu$ of $(\mu,\nu)$. It's a convex problem, so it provides a computationally convenient approximation of $\mathcal C(\mu,\nu)$ when $\varepsilon$ is small. The dual version of EOT is also useful in econometric applications.
Optimal pairs satisfy the first order optimality conditions
Moreover, if $\pi$ is an entropic optimal transport plan, i.e., solves the primal ((ref)) and $(\varphi,\psi)$ is a pair of entropic optimal transport potentials, i.e., solve the dual in theorem (ref), then they are related by
for all integrable functions $\zeta$. This can be written $d\pi=e^{\frac{\varphi+\psi-c}{\varepsilon}}d(\mu\otimes\nu)$ in short.
Partial and unbalanced optimal transport are variants of optimal transport motivated by situations in which the origin and destination marginals have different masses, or the particular application prescribes only part of the mass be transported. It can also accommodate applications where there is some flexibility in the specification of the marginals. In the most common version, unbalanced optimal transport solves the following problem:
where $H$ denotes relative entropy, defined in ((ref)). Unbalanced optimal transport penalizes marginals that are far from $\mu$ and $\nu$. The dual version of UOT is also useful in econometric applications.
The basic Kantorovich formulation of the optimal transport problem from section (ref) can be rewritten:
The term in brackets in the previous display is a function of both $x$ and the conditional distribution $\pi_x$ of $Y$ conditional on $X=x$ when the joint distribution of $(X,Y)$ is $\pi$. Weak optimal transport generalizes this to more general functions $\mathcal F(x,\pi_x)$ of $x$ and $\pi_x$.
In addition to traditional optimal transport, many extensions are special cases of weak optimal transport.
Duality theory for weak optimal transport closely follows Kantorovich duality.
The optimal transport problem has a natural extension with more than two prescribed marginals $\mu_1,\ldots,\mu_K$ on $\mathcal X_1,\ldots,\mathcal X_K$. The Kantorovich formulation is the following.
Considered as a problem of shifting mass $\mu_1$ on $\mathcal X_1$ to $\mathcal X_2\times\ldots\times\mathcal X_K$, it differs from a standard optimal transport problem in two ways. First, the origin and destination dimensions are different. Second, the probability distribution on $\mathcal X_2\times\ldots\times\mathcal X_K$ is not fully specified, only its marginals $\mu_2,\ldots,\mu_K$. The Kantorovich duality extends straightforwardly, with dual
where the supremum is taken over continuous and bounded functions satisfying the dual constraints. The theory on existence and uniqueness of optimal transport maps, on the other hand, is more involved. Up to relabeling of the variables, the Monge version of the problem takes the following form.
where the infimum is taken over functions $F_k$ that push $\mu_1$ forward to $\mu_k$ in the sense of ((ref)). Two cases are well understood: the case with quadratic cost, and the case with marginals on $\mathbb R$. In the first case, the cost function is taken to be:
The following theorem shows uniqueness of the optimal transport maps $(F_2,\ldots,F_K)$.
The second case, which is well understood is the case with univariate marginals with submodular transport cost. We first define submodularity.
Theorem (ref) tells us that the Monge multi-marginal optimal transport problem admits the deterministic coupling $(X_1,F_2(X_1),\ldots,F_K(X_1))$ as solution. This coupling is obtained by monotone rearrangements of each component of the random vector. As a corollary, we have the following rearrangement inequality:
for any joint distribution $\pi$ with marginals $(\mu_1,\ldots,\mu_K)$, where $Q_\mu$ is the quantile function of probability distribution $\mu$.
Important special cases include optimal transport on the real line, between Gaussians and between discrete distributions.
The case $\mathcal X\subseteq\mathcal R$ and $\mathcal Y\subseteq\mathbb R$ yields rich connections between optimal transport theory and the theory of monotone rearrangements of probability distributions. As a result, it yields closed form solutions for optimal transport maps and Wasserstein distances.
The class of cost functions relevant to econometric applications is the class of submodular costs. For convenience, we give a useful equivalent definition of submodularity in the case of $\mathbb R^2$:
The class of (strictly) submodular functions includes functions of the form $c(x,y)=h(x-y)$, with $h$ (strictly) convex, or $h(x+y)$ with $h$ (strictly) concave. This includes the case $c(x,y)=|x-y|^p,$ for $p\geq1$, which is relevant to Wasserstein distances in particular. The general theory simplifies greatly in case of submodular costs in $\mathbb R^2$ because $c$-cyclical monotonicity of a subset $\Gamma$ of $\mathbb R^2$ simplifies to monotonicity. A subset $\Gamma$ of $\mathbb R^2$ is said to be monotone if $(x,y)\in\Gamma$ and $(x^\prime,y^\prime)\in\Gamma$ imply $(x-x^\prime)(y-y^\prime)\geq0$. Therefore, optimal transport plans for submodular costs in $\mathbb R^2$ have monotone supports. Heuristically, this means that mass is transported from the lowest $x$ to the lowest $y$, then the second lowest $x$ to the second lowest $y$, and so on from left to right. Hence, an optimal transport plan is a monotone rearrangement of mass $\mu$ into mass $\nu$. This means that an optimal coupling $(X,Y)$ of $(\mu,\nu)$ is such that $X$ and $Y$ are comonotone random variables.
By the probability integral transform, it is always possible to write $X=Q_X(U)$ and $Y=Q_Y(V)$, for a pair $(U,V)$ of uniform random variables. The random variables $X$ and $Y$ are comonotone, when $U=V$, so that they are ordered in the same way. Optimality of the comonotone coupling can be formalized with the following result.
The Cambanis-Simons-Stout monotone rearrangement inequalities are a direct corollary of theorem (ref). So is the closed form solution for Wasserstein distances in $\mathbb R$:
When $\mu$ is absolutely continuous with respect to Lebesgue measure, the comonotonic coupling $(X,Y):=(Q_\mu(U),Q_\nu(U))$ is also equal to $(X,Q_\nu(F_\mu(X)))$, where $F_\mu$ is the cumulative distribution function of $\mu$. Therefore, the monotone non decreasing map $T:=Q_\nu\circ F_\mu$ is the unique optimal transport map. Notice that the optimal transport plan and the optimal transport map are independent of the submodular cost function in the definition of the optimal transport problem.
Beyond the scalar case described in the previous section, closed form solutions for optimal transport maps are rare. One such closed form solution is the optimal transport map between Gaussian distributions with quadratic transportation cost $c(x,y)=\|x-y\|^2$. Let $\mu$ and $\nu$ be Gaussian probability distributions on $\mathbb R^d$, with means $m_1$ and $m_2$ and covariance matrices $\Sigma_1$ and $\Sigma_2$ respectively. If we can find a deterministic coupling $(X,T(X))$ of $(\mu,\nu)$, such that $T$ is the gradient of a convex function, we know from theorem (ref) that such a coupling is unique and that $T$ is the optimal transport map.
Let $b$ be a vector in $\mathbb R^d$ and $A$ a $d\times d$ matrix. Since affine maps of the form $x\mapsto b+Ax$, preserve Gaussianity, it is natural to look for $T$ among affine maps. If $A$ is symmetric and positive semidefinite, then $T:x\mapsto b+Ax$ is the gradient of the convex function $x\mapsto x^\top b+x^\top Ax/2$. There remains to find $A$ and $b$ such that $A$ is symmetric positive semidefinite and $(X,b+AX)$ is a coupling of $(\mu,\nu)$.
If $A$ is symmetric positive semidefinite and $X$ has Gaussian distribution $\mu$ with mean $m_1$ and covariance $\Sigma_1$, then $b+AX$ has Gaussian distribution with mean $b+Am_1$ and covariance $A\Sigma_1 A$. Hence, $(X,b+AX)$ is a coupling of $(\mu,\nu)$ if $m_2=b+Am_1$ and $\Sigma_2=A\Sigma_1A$. A symmetric positive semidefinite solution of $\Sigma_2=A\Sigma_1A$ can be found by noticing that
From the latter, we obtain $A$ and therefore the optimal transport map $T$ is
Using expression ((ref)) for the Wasserstein distance, we have the closed form
The case, where the cost function $c(x,y):=1\{(x,y)\in \Gamma\}$ is binary is very useful in probability theory and its applications. We will show later that it has several uses in econometrics, particularly in data combination and partial identification of structural models.
With cost function $c(x,y)=1\{(x,y)\in \Gamma\}$, the Kantorovitch dual of the optimal transport problem
takes the form
It can be shown that the value of the dual is unchanged if the potential functions $\varphi$ and $\psi$ are restricted to be indicator functions themselves, specifically if $\varphi=1\{x\in A\}$ and $\psi(y)=-1\{y\in B\}$, and the supremum is taken over pairs of sets $(A,B)$. The dual constraint then becomes
which implies
The right-hand side of the latter expression is equal to $1$ if and only if $y$ is in the set
Therefore, the value of the dual is unchanged if the potential functions $\varphi$ and $\psi$ are restricted to $\varphi=1\{x\in A\}$ and $\psi(y)=-1\{y\in A^\Gamma\}$, and the supremum is taken over sets $A$.
We therefore have the following result (theorem 2.27 in villani2003topics):
Taking $\Gamma:=\{(x,y):d(x,y)\leq\varepsilon\}$ for some distance $d$ and some $\varepsilon\geq0$, a direct application of theorem (ref) is one of Volker Strassen's many celebrated theorems, namely the characterization of the L\'evy-Prokhorov distance between probability distributions (see corollary 1.28 and remark 1.29 in villani2003topics for details).
We now investigate the formulations of the optimal transport problem and aspects of the theory when the transported distributions have finite support. In that case, the optimal transport problem is also called the optimal assignment problem. We now have finite spaces $\mathcal X=\{1,\ldots,N\}$ and $\mathcal Y=\{1,\ldots,M\}$. The distributions $\mu$ and $\nu$ on $\mathcal X$ and $\mathcal Y$ are now characterized by probability mass vectors $(\mu_1,\ldots,\mu_N)$ and $(\nu_1,\ldots,\nu_M)$. A coupling of $(\mu,\nu)$ is a probability mass vector $\pi=(\pi_{xy})$ on $\mathcal X\times\mathcal Y$ with marginals $\mu$, and $\nu$, which is now equivalent to the adding up constraints
In the finite case, a coupling $\pi$ of $(\mu,\nu)$ is deterministic if $M=N$, $\mu_x=\nu_y=1/N$ and $\pi$ is the matrix of a permutation of $\{1,\ldots,N\}$. Each $x\in\mathcal X$ is assigned to exactly one $y\in\mathcal Y$. The discrete version of the Kantorovich formulation of the optimal transport problem is that of finding a coupling $\pi$ of $(\mu,\nu)$, i.e., a distribution on $\mathcal X\times\mathcal Y$ that satisfies ((ref)), that achieves
The discrete version of the Monge formulation of the optimal transport problem is that of finding a deterministic coupling, i.e., a permutation $\sigma$ of $\{1,\ldots,N\}$, that achieves
The latter is also called the pure assignment problem. Kantorovich duality in case of discrete distributions is a special instance of the duality of linear programming. Once specialized to distributions with finite support, theorem (ref) states equality of the primal problem ((ref)) under the adding up constraints ((ref)) with the dual
Both the primal and the dual have solutions $\pi$ and $(\varphi,\psi)$ respectively. Finally, in the language of linear programming, the concentration of a solution $\pi$ to the primal on the cyclically monotone set $\partial_x\varphi$ of ((ref)), where dual constraints are binding, i.e.,
is called complementary slackness. By theorem (ref), if $\pi\in\mathcal M(\mu,\nu)$ and $\varphi$ is $c$-concave, then the plan $\pi$ and the pair of potentials $(\varphi,\varphi^c)$ are optimal transport solutions if and only if they are complementary in the sense of ((ref)).
When $M=N$ and $\mu_x=\nu_y=1/N$, the Kantorovich form of the discrete optimal transport problem has a deterministic coupling as a solution. This coupling is concentrated on the graph of the map $y=\sigma(x)$, where $\sigma$ is a permutation of $\{1,\ldots,N\}$. The Wasserstein metric also specializes easily to discrete distributions in case $N=M$. Letting $d(x,y)$ be a distance on $\{1,\ldots,N\}$, and letting the cost function be $c(x,y)=d(x,y)^p$ for some $p\geq1$, then $W_p(\mu,\nu)=(\mathcal C(\mu,\nu))^{1/p}$ defines a distance on the set of probability measures on $\{1,\ldots,N\}$.
The discrete optimal transport problem in its Monge formulation is also called the optimal assignment problem. The classical algorithms to solve the assignment problem are the Hungarian algorithm of kuhn1955hungarian [kuhn1955hungarian,kuhn1956variants], munkres1957algorithms, edmonds1972theoretical and the auction algorithm of bertsekas1979distributed. All standard computing packages contain efficient implementations of these algorithms. Their generic complexity is $O(n^3)$.
However, econometric applications rely mostly on the Kantorovich formulation of the optimal transport problem. Extensions of the Hungarian and the auction algorithms exist for the Kantorovich formulation, but variants of the network simplex algorithm hitchcock1941distribution, and shortest path algorithms, tomizawa1971some, jonker1987shortest, are often preferable. An even more favored alternative is entropy regularization. The discrete Kantorovich optimal transport problem is replaced with the discrete version of the entropic optimal transport with a suitably small $\varepsilon$. The problem is solved using the iterated proportional fitting (IPFP), or Sinkhorn algorithm of deming1940least, sinkhorn1964relationship, and ruschendorf1995convergence. The latter iterates between $\varphi$ and $\psi$ to solve ((ref)).
In what follows, we describe a specific form of the network simplex algorithm for the Kantorovich formulation of the optimal transport problem. We also propose a Python implementation of the algorithm for illustration purposes. However, for best results, we recommend using optimized code. For instance, the Python optimal transport library POT, see flamary2021pot, implements a variant of the network simplex algorithm, while the Julia optimal transport library OptimalTransport.jl implements a variant of the shortest path algorithm of jonker1987shortest. Both also implement the IPFP algorithm for entropy regularized optimal transport.
Recall the notation from section (ref). The origin set is $\mathcal X:=\{x_1,\ldots,x_M\}$ and the destination set is $\mathcal Y:=\{y_1,\ldots,y_N\}$. The marginal probability mass vectors are $\mu:=(\mu_i:=\mu_{x_i})_{i=1}^M$ and $\nu:=(\nu_i:=\nu_{x_i})_{i=1}^N$ respectively. We look for a joint probability $\pi:=(\pi_{ij}:=\pi_{x_iy_j})_{ij}$ with marginals $\mu$ and $\nu$ that minimizes transport cost $C$ with elements $C_{ij}:=c(x_i,y_j)$ . Consider the bipartite graph $G$ with $M+N$ nodes in $\mathcal X\cup\mathcal Y$ and $MN$ edges in $\mathcal X\times\mathcal Y$. The network simplex algorithm finds an optimal $\pi$ by finding optimal weights on the edges of graph $G$.
The computation of optimal transport maps often involves a discrete marginal, typically an empirical distribution, and a continuous marginal, often the uniform on $[0,1]^d$. The computation of the optimal transport map in such contexts is called the semi-discrete optimal transport problem. Let random vector $X$ have absolutely continuous distribution $\mu$ on a closed and compact $\mathcal X$, and $Y$ have discrete distribution $\nu$ with probability mass $((y_1,q_1),\ldots,(y_K,q_K))$. The mass points $y_1,\ldots,y_K$ are in $\mathbb R^d$, and the weights $q_1,\ldots,q_K$ are in $[0,1]$ and sum to $1$. The Kantorovich dual of the optimal transport problem from $\mu$ to $\nu$ is
The $c$-conjugate of the vector $\psi\in\mathbb R^K$ is the function $\psi^c(x)=\min_{1\leq j\leq K}\{c(x,y_j)-\psi_j\}.$ The Kantorovich semi-dual given by
is therefore a finite dimensional optimization problem. For each $j\leq K$, on the set $\mathcal L_j(\psi)$, called {\em Laguerre cell} and defined by
the $c$-conjugate of $\psi$ is equal to $c(x,y_j)-\psi_j$. Hence, the Kantorovich semi-dual can be rewritten
which is a convex program with optimizer $\hat\psi$. The first order conditions are $q_j=\mu(\mathcal L_j(\hat\psi))$ for all $j$ and the optimal transport map is the piecewise constant map, which takes values $y_i$ on each Laguerre cell $\mathcal L_j(\hat\psi)$. The partition of $\mathcal X$ into Laguerre cells is called a power diagram, and the original algorithm is due to aurenhammer1998minkowski. A variant of this algorithm is implemented in the power diagram component of the Julia DelauneyTriangulation.jl library, see VandenHeuvel2024DelaunayTriangulation.
In many applications in economics, such as matching markets and trade gravity equations, the computational problem is reversed: instead of computing the optimal transport plan given marginals and cost, one is concerned with recovering transport costs from the optimal transport plan. Here we address the issue of computing the optimal transport cost given perfect knowledge of the optimal transport plan. Here we give details of the algorithm proposed in carlier2023sista.
Assume that the transport plan solves the discrete version of entropic optimal transport with transport cost $c(x,y)$:
with dual (where we have slightly redefined the dual variables in order to drop the $-1$ inside the exponential)
Optimal $\pi$ and $(\varphi,\psi)$ are related by
Assume that the true transport cost belongs to a parametric family $c(x,y;\beta)$ with finite dimensional parameter $\beta$. The object of the procedure is to recover $\beta$ from the knowledge of $\pi$ and the optimality conditions. Hence we are looking for $\beta$ such that
satisfies the marginal constraints
and moment constraints
Conditions ((ref)-(ref)) can be seen as the first order conditions of a Poisson maximum likelihood, i.e., the maximization of
The algorithm builds on the Sinkhorn algorithm in the following way. Choose step size $\rho$ and initial values $\beta^0,\varphi^0$ and $\psi^0$. The algorithm alternates Sinkhorn updates on the potentials $\varphi$ and $\psi$ (to enforce the marginals):
with a gradient descent to update $\beta$:
carlier2023sista add a LASSO regularization and show convergence of the proximal version of the algorithm for a class of sparse parametric cost functions.
The $1$-Wasserstein distance is increasingly used in machine learning algorithms. Conversely, deep neural nets can be used to compute the $1$-Wasserstein distance between two probability distributions. The starting point is always the Kantorovich-Rubinstein dual ((ref)) of $W_1(\mu,\nu)$.
The idea of generative adversarial networks is the competition between a {\em generator} who tries to mimic samples from the true data generating process, and a {\em critic}, who tries to discriminate between the generator's data and data from the true data generating process. The generator uses data $(Z_1,\ldots,Z_n)$ sampled from a prior $P$, and transforms them through a parametric function $g_\theta$ to mimic true data. The critic tries to discriminate between the generator data $(g_\theta(Z_1),\ldots,g_\theta(Z_n))$ drawn from $\nu_\theta$ (the push-forward of $P$ by $g_\theta$) and true data $(X_1,\ldots,X_n)$ drawn from the true data generating process $\mu$. In Wasserstein generative adversarial networks (hereafter WGAN) of arjovsky2017wasserstein, the critic uses a parametric family of Lipschitz functions $\varphi_\eta$ and the criterion
The critic chooses $\eta$ to maximize ((ref)) and therefore maximize discrimination. The generator then chooses parameter $\theta$ to make discrimination as difficult as possible, hence solving
As indicated in the display above, the inner maximization is an approximation of the $1$-Wasserstein distance between two fixed probability distributions $\mu$ and $\nu_\theta$. Parameter values $\eta^\ast$ that maximize $\mathcal D(\eta,\theta)$ in ((ref)) and hence approximate $W_1(\mu,\nu_\theta)$ are typically computed using root mean square propagator (RMSProp). Efficient implementations of RMSProp can be found in the machine learning libraries such as PyTorch in Python and Optimizers.jl in Julia. One aspect of the algorithm that remains unsatisfactory is the way to impose the Lipschitz constraint on the family of neural nets $\varphi_\eta$ to conform with the Kantorovitch-Rubinstein dual ((ref)). The current preferred solution is to add a gradient penalty to the critic's objective ((ref)). This term directly penalizes the norm of the derivative of the potential $\varphi_\eta$:
where the $\hat X_i=\varepsilon_i X_i+(1-\varepsilon_i)g_\theta(Z_i)$ are random convex combinations of real and generated observations.
The primitives in optimal transport problems are the marginal distributions to be coupled and the transport cost. In many applications of optimal transport, at least one of the marginals must be estimated from data. In this section, we survey the theory that addresses the question of how the estimation of (one of) the marginals affects:
We address each question in turn. In all this section, unless otherwise specified, $\mu_n$ denotes the empirical distribution based on a size $n$ i.i.d. sample of observations $(X_1,\ldots,X_n)$ with distribution $\mu$. Most of the results in this section extend to more general estimators of the marginals, denoted $\hat\mu_n$, as long as they converge in distribution with a suitable rate to the true marginals.
Under very general conditions, optimal transport plans between estimated marginals converge to an optimal transport plan between the true marginals.
In order to quantify the rate of convergence of optimal transport, we need to specify particular classes of transport costs. Most of the results in the literature relate to Wasserstein distances $W_p$.
We consider the rate of convergence of Wasserstein distance $W_p(\mu_n,\mu)$ to zero. Rates of convergence of $W_p(\mu_n,\nu_m)$ to $W_p(\mu,\nu)$ follow from the triangle inequality. First consider $W_1$. Let $X$ be a random vector with distribution $\mu$ and $\mu_n$ be the empirical distribution relative to the i.i.d. sample $(X_1,\ldots,X_n)$. By the Kantorovich-Rubinstein duality ((ref)),
The right-hand side is the supremum over the class of Lipschitz functions of the empirical process $\varphi\mapsto\Sigma_{i=1}^n (\varphi(X_i) - \mathbb E \varphi(X_i))/n$. Hence, asymptotic distributions and rates of convergence can be obtained with empirical process methods. Asymptotic bounds on the rate of convergence of $\mathbb EW_1(\mu_n,\mu)$ obtained using empirical process theory are the following:
The rates above are the best attainable rates of convergence, except when $d=2$, where they are off by a $\sqrt{\ln n}$ factor. The rate $n^{-1/d}$ for $d\geq3$ is an instance of the curse of dimensionality. The rate of convergence slows exponentially with the dimension of the space $d$. This is related to the fact that the supremum of the empirical process in ((ref)) is taken over a very large class of functions.
Some ways to circumvent the curse of dimensionality are the following. In each case, the curse of dimensionality is circumvented in the sense that the rate of convergence is $n^{-1/2}$ irrespective of dimension $d$.
Existing theory on the estimation of optimal transport maps concerns mostly the case of optimal transport with at least one absolutely continuous marginal and quadratic cost. We will concentrate on this case here. As we have seen in theorem (ref), the optimal transport map from $\mu$ to $\nu$ with transport cost $c(x,y)=\|x-y\|^2$ is $T=x-\nabla\varphi(x)$, where $\varphi$ is a $c$-concave function that achieves the dual. Equivalently, the optimal transport map is $T=\nabla\vartheta(x)$, where $\vartheta$ is a convex function that achieves
where $\tilde\vartheta^\ast$ is the convex conjugate of $\tilde\vartheta$. This form is called semi-dual, because the dual constraints are incorporated in the dual objective using the fact that any optimal dual pair must be of the form $(\vartheta,\vartheta^\ast)$. The optimal transport map can be estimated as the gradient $\hat T:=\nabla\hat\vartheta$ of a solution $\hat\vartheta$ to
with the usual normalization $\vartheta(0)=0$. Computational issues are taken up in section (ref). To derive statistical properties of this estimator of the optimal transport map, we need to specify the class $\mathcal F$ of convex functions over which the maximum in ((ref)) is taken. We also need to specify the smoothness of marginals $\mu$ and $\nu$. There exists classes $\mathcal F$ of convex functions such that the solution of ((ref)) achieves the following rates
In the display above, $s$ is the number of bounded derivatives of $\nabla\vartheta$. Up to logarithm factors, the rates above are minimax rates for estimators of $s$-smooth optimal transport maps between marginals $\mu$ and $\nu$ with densities bounded away from $0$ and $\infty$ on a compact support. When $\nabla\hat\vartheta$ maximizes ((ref)) where $\mu_n$ and $\nu_n$ are replaced with Gaussian kernel estimators with bandwidth $h_n$, the following type of pointwise central limit theorem
can be shown on the Torus $\mathbb R^d/ \mathbb Z^d$ (to avoid boundary issues) with $d\geq3$.
The most complete account of the history of research on optimal transport theory up to that point is given in chapter 3 of villani2008optimal. The history of some more recent developments can be found in the preface of santambrogio2015optimal. Each chapter of villani2008optimal also contains very detailed bibliographical notes that precisely trace the history of ideas and contributions in each aspect of the theory. The basic formulations and main mathematical questions relating to optimal transport are elegantly presented in the introduction of villani2003topics. Theorem (ref) (the principal Kantorovich duality theorem) is taken from Theorem 1.3 page 19 of villani2003topics and theorem 5.10 page 58 of villani2008optimal.
Chapter 4 of villani2008optimal gives a proof of existence of optimal transport plans, Theorem 4.1. Chapter 5 of villani2008optimal gives a very thorough introduction to $c$-cyclical monotonicity, $c$-concavity, and their relation to optimality and the dual Kantorovich problem. The relation between $c$-cyclical monotonicity and duality is also detailed in section 1.6 of santambrogio2015optimal. Chapter 5 of villani2008optimal also presents as complete a picture of Kantorovich duality as needed for applications to econometrics. Theorem (ref) is taken from theorem 5.10 of villani2008optimal. Theorem (ref) is taken from theorem 5.20 page 77 of villani2008optimal.
Chapter 10 of villani2008optimal gives the best account of the theory underlying solutions to the Monge problem, namely existence of optimal transport maps, and sufficient conditions on the cost function and the marginal distributions. Theorem (ref) is implied by theorem 10.28 page 243 of villani2008optimal. Section 1.3 of santambrogio2015optimal is a good resource to understand the results in the special case, where the cost function has the form $c(x,y)=h(x-y)$, with $h$ strictly convex. Theorem (ref) can be found as theorem 2.9 page 63 of santambrogio2015optimal for instance. This includes the quadratic cost, which is particularly relevant to econometric applications. Theorem (ref) is a combination of theorem 9.4 page 209 and theorem 10.42, corollary 10.44 and particular case 10.45 on page 256 of villani2008optimal. See also chapter 3 of villani2003topics.
Sections 5.1 and 5.2 of santambrogio2015optimal give a thorough treatment of Wasserstein distances between probability distributions on $\mathbb R^d$. Theorem (ref) is theorem 5.9 page 184 of santambrogio2015optimal. Chapter 6 of villani2008optimal considers the more general case of distributions on a Polish space, which can be useful in applications with distributions over the space of distributions, for example. Section 1.2 of chewi2024statistical is also a very good introduction. Section 7.1 of chewi2024statistical is a very accessible introduction to geodesics and interpolation in the Wasserstein metric. Chapter 5 of villani2003topics introduces displacement convexity. Chapters 16 and 17 of villani2008optimal are much more thorough and include more recent material. The original paper agueh2011barycenters is still the best reference for the theory of Wasserstein barycenters.
nutz2021introduction gives a very complete account of the mathematics of entropic optimal transport. Theorem (ref) is derived from theorem 7 page 35 of nutz2021introduction. Theorem (ref) is a special case of theorem 1.2.19 page 28 of chizat2017unbalanced. pass2015multi is an in-depth survey of the theory of multi-marginal optimal transport up to that point. Theorem (ref) is taken from theorem 2.1 page 27 of gangbo1998optimal and theorem (ref) is taken from theorem 4.1 page 527 of carlier2003class. gozlan2017kantorovich introduce weak optimal transport\footnote{See also chone2023weak.}, and veraguas2018existence prove the duality theorem presented in section (ref). More precisely, theorem (ref) is taken from theorem 3.1 page 203 of veraguas2018existence.
Chapter 3 of galichon2016optimal gives a concise presentation of discrete optimal transport duality and a simple constructive proof of the existence of a deterministic solution to the discrete Kantorovich optimal transport problem. Section 2.4 of peyre2019computational gives a proof of the triangle inequality for the discrete Wasserstein distance which helps understand the proof for continuous distributions in theorem 7.3 of villani2003topics.
The classical algorithms to compute optimal assignments are presented and reviewed in chapter 3 of peyre2019computational: The network simplex algorithm in section 3.5, the Hungarian algorithm in section 3.6 and the auction algorithm in section 3.7.
chewi2024statistical is an excellent recent resource on the estimation of optimal transport. Chapter 2 of chewi2024statistical covers the estimation of Wasserstein distances (mostly $W_1$), rates of convergence and ways of circumventing the curse of dimensionality. The semi-dual estimator of optimal transport maps was introduced in chernozhukov2017monge, with a proof of uniform convergence to the population map. Chapter 3 of chewi2024statistical derives rates of convergence for the semi-dual estimator of optimal transport maps. The minimax rates of section (ref) are derived in hutter2021minimax. The central limit theorem for optimal transport maps can be found in manole2023central.
In this second part of our review, we examine applications of optimal transport methods in econometrics around three major aspects of the theory: first, the duality of optimal transport and its entropic, unbalanced, weak and multi-marginal extensions; second, the uniqueness and cyclical monotonicity of optimal transport maps; third, the metric on the space of probability distributions based on optimal transport. Broadly, the corresponding categories of econometric applications are the following: for the first aspect, partial identification and data combination problems; for the second one, multivariate quantiles and their applications; and for the third one, distributional robustness.
This section reviews econometric applications that involve direct use of optimal transport formulation and the duality of optimal transport for computational advantage in a variety of problems or for the characterization of identified sets and sharp bounds in partially identified problems, specifically in incomplete models and broadly defined data combination problems. These applications mostly rely on the convenience of the simplex algorithm, which relies on theorem (ref), on rearrangement inequalities of theorem (ref), and on the Kantorovich duality theorem (ref) and its extensions to entropic (ref), unbalanced (ref), weak (ref), and multi-marginal (ref) optimal transport. Econometric applications often involve conditioning on exogenous covariates, which, as we shall see, can be handled in a variety of ways.
In a variety of different applications, optimal transport is used to formulate the problem and unlock computational advantages. This includes applications to discrete games with multiple equilibria, in galichon2011set, henry2015combinatorial, gu2023dual, to discrete choice models in galichon2022cupid, chiong2016duality, shi2018estimating, bonnet2022yogurts, to treatment assignment under budget constraints in sunada2025optimal. Optimal transport duality is used to transform an infinite dimensional optimization problem to a finite dimensional one, as in the partial identification of incomplete models in galichon2006inference (galichon2006inference,galichon2009test,galichon2011set), to transform optimization of expectations with respect to the joint distribution to expectations with respect to the marginals only, in the references above, as in the estimation of distributional treatment effects in ji2023model. Monotone rearrangement inequalities of theorem (ref) are applied to the characterization of treatment effect bounds in fan2009partial, fan2010sharp, fan2014identifying, and in data combination problems as in gechter2024generalizing,fan2024multidimensional,d2024linear,meango2025combining.
In the rest of this section, we review several classes of applications: treatment effects, data combination, incomplete models, and discrete choice and matching estimation. We explain in detail how optimal transport is applied in each and we provide algorithms.
Optimal transport methods have been used in solving optimal assignment problems, but there is no established literature yet. kitagawa2025who consider the problem of optimally matching one population with another to minimize a matching cost (the opposite of a matching surplus). They estimate the entropy regularized optimal assignment based on estimated cost and marginals. sunada2025optimal consider a utilitarian planner seeking to maximize utility $w(x,t)$ from giving treatment $t\in\{0,1\}$ to individuals with characteristics $x$. Only a proportion $p$ of the total population can be treated. What is the propensity score $\mu(t\vert x)$ that maximizes this constrained optimization problem
subject to $\int \mu(1\vert x) dF_X(x)=p$? Since the treatment assignment must be Bernoulli random variable with probability of success $p$ (distribution $F_T$), the maximum above solves the optimal transport problem
Here the optimal transport formulation is used for computational convenience.
A large class of applications of optimal transport duality concerns treatment effects. The treatment effects model we consider has the following main ingredients. A sample of variables $(Y,D,X)$ is observed. The variable $D$ is a binary treatment indicator. Unobserved potential outcomes $(Y_0,Y_1)$ determine the observed outcome through the relation $Y=Y_0(1-D)+Y_1D$. Outcomes $Y$ may be scalar or multivariate, depending on the application. Finally, $X$ denotes a vector of exogenous covariates. Assume selection on observables $(Y_0,Y_1)\perp\!\!\!\!\perp D\vert X$ or any other condition that guarantees identification of the joint distributions of $(Y_0,X)$ and $(Y_1,X)$. In all cases reviewed below, $(Y_i,D_i,X_i)_{i=1}^n$ is an i.i.d. sample of observations. We denote $P_{Y_0,X}$ and $P_{Y_1,X}$ the joint distributions of $(Y_0,X)$ and $(Y_1,X)$ respectively, and $P_{Y_0,X;n}$ and $P_{Y_1,X;n}$ their empirical counterparts based on the sample. Similarly, we denote $P_{Y_0\vert X}$ and $P_{Y_1\vert X}$ the conditional distributions of $Y_0\vert X$ and $Y_1\vert X$ respectively, and $P_{Y_0\vert X;n}$ and $P_{Y_1\vert X;n}$ their empirical counterparts based on the sample.
Optimal transport is applied to this framework in two very different ways. First, optimal transport is applied to covariate matching procedures in conditional average treatment effects estimation in gunsilius2021matching, dunipace2021optimal, and charpentier2023optimal, and optimal experimental design, such as site selection for external validity in bouyamourn2025where . Second, fan2009partial, fan2010sharp, fan2014identifying, fan2017partial, ji2023model, ober2023estimating, kaji2023assessing, lin2025estimation and fan2024multidimensional apply optimal transport methods to the (partial) identification and estimation of functionals of the joint distribution of $(Y_0,Y_1,X)$. This is an instance of data combination problem, since $Y_0$ and $Y_1$ are never observed for the same individual. A variety of applications to other types of data combination problems will be reviewed in the next subsection.
gunsilius2021matching propose an alternative to propensity score matching for the estimation of average treatment effects in a setting with selection on observables. The idea is to match each treated observation with an average of control observations weighted by the similarity in their covariates to the treated individual. We have
The idea of the covariate matching procedure is to replace the term $\mathbb E[Y_1\vert D=0,X]$ by a weighted average of controls. Let $\mu_0$ and $\mu_1$ be the covariate distributions of control and treated units respectively, and let $\mu_{0;n}$ and $\mu_{1;n}$ be their empirical counterparts based on the sample. Let $\hat\pi$ be the coupling of $\mu_{0;n}$ and $\mu_{1;n}$ that solves the unbalanced optimal transport problem ((ref)). As an optimal transport solution, this matching of covariates does not require the support of covariates to be the same for treated and control units. In addition, the unbalanced optimal transport problem discards both control and treated units that don't have sufficiently close covariate matches in the opposite group. gunsilius2021matching show that it reduces estimation bias. With the unbalanced optimal transport solution $\hat\pi$, we can construct the weights as the conditional probability distribution $\hat\pi_{1\vert 0}(x_1\vert x_0)$ of treated group covariates given the value of the control group covariate. Each treated unit $i$ with covariate $X_i$ is matched with control unit outcomes $Y_k$ with a weight $\hat\pi_{1\vert 0}(X_k\vert X_i)$. Finally, $\mathbb EY_1$ is estimated with the following empirical version of ((ref)):
We turn now to the second major application of optimal transport to treatment effects, which is partial identification and estimation of functionals of the joint distribution of potential outcomes. When outcomes are scalar, and the function $h$ is submodular or supermodular, fan2014identifying propose to apply monotone rearrangement inequalities of theorem (ref) to derive closed form sharp bounds for $\mathbb E[h(Y_0,Y_1)]$. These closed form solutions can be extended to conditional sharp bounds on $\mathbb E[h(Y_0,Y_1)\vert X=x]$ based on conditional versions of the monotone rearrangement inequalities of theorem (ref). kaji2023assessing extend monotone rearrangement inequalities to derive sharp bounds on subgroup treatment effects. They also derive sharp bounds on the subgroup proportion of winners (those who benefit from the treatment). The latter can be derived using optimal transport duality with zero-one costs functions as we show below. Their bounds analysis is motivated by the policy relevance of the average effect of a treatment or the sign of the treatment effect for the section of the population with the lowest outcomes before treatment. The subgroup treatment effects are defined as $\mathbb E[Y_1-Y_0\vert a<U_0<b]$. In the latter expression, $U_0$ is the rank of an individual in the untreated outcome distribution. The latter is defined by $Y_0=Q_0(U_0)$, where $Q_0$ is the quantile function of $Y_0$. The sharp upper bound on $\mathbb E[Y_1-Y_0\vert a<U_0<b]$ is obtained in the rearrangement where individuals ranked from rank $a$ to rank $b$ in the untreated distribution are ranked from rank $1-b+a$ to rank $1$ in the treated distribution. Similarly, the sharp lower bound is obtained in the rearrangement where individuals ranked from rank $a$ to rank $b$ in the untreated distribution are ranked from rank $0$ to rank $b-a$ in the treated distribution. The sharp bounds are therefore given by
The subgroup proportion of winners from theorem 3 in kaji2023assessing, can be obtained with an application of theorem (ref) on Kantorovich duality with binary transport cost functions. The subgroup proportion of winners is defined as $\mathbb P(Y_1>Y_0\vert a< U_0 <b)$ and can be obtained by dividing $\mathbb P(Y_1>Y_0, a< U_0 <b)$ by $b-a$. By construction, $U_0$ is the uniform random variable on $[0,1]$ such that $Y_0=Q_0(U_0)$. Define $U_1$ similarly, so that $Y_1=Q_1(U_1)$ and define $V:=F_0(Q_1(U_1))$, so that by definition $Y_1>Y_0$ if and only if $V_1>U_0$. The lower bound on $\mathbb P(Y_1>Y_0, a< U_0 <b)$ is equal to the solution of the binary cost optimal transport problem
Since $\mathbb P(Y_1\leq Y_0, a< U_0 <b)=b-a-\mathbb P(Y_1>Y_0, a< U_0 <b)$, the upper bound can be obtained analogously. By theorem (ref), the lower bound $\inf_\pi \pi((U_0,V_1)\in\Gamma)$ is equal to its dual $\sup\left[P_{U_0}(A)-P_{V_1}\left(A^\Gamma\right)\right]$, where $A^\Gamma:=\{v:\exists u\in A,(u,v)\notin\Gamma\}$. Hence, when $A\not\subseteq(a,b)$, $v$ is always in $A^\Gamma$, and $A$ therefore does not contribute to the dual. Hence, take the supremum over $A\subseteq(a,b)$. Note that there exists $u\in A, u\geq v$ if and only if $v=\bar a:=\sup A$. Hence $A^\Gamma=[0,\bar a]$. The dual value is non smaller if $A$ is restricted to take the form $(a,\bar a]$ (in the terminology of galichon2011set, the sets $A=(a,\bar a]$ form a {\em core determining class}). Finally, since $V_1=F_0(Q_1(U_1))$, with $U_1$ uniform on $[0,1]$, its cdf is $[F_0\circ Q_1]^{-1}=F_1\circ Q_0$. Hence, the dual is equal to
Dividing by $b-a$ and adding the non-negativity constraint, we obtain the lower bound in theorem 3 of kaji2023assessing.
Monotone rearrangement inequalities rely on scalar potential outcomes. More generally, one of the major advantages of optimal transport methods is the ability to handle multivariate outcomes. ji2023model, fan2024multidimensional and lin2025estimation apply optimal transport methods to derive sharp bounds for functionals of the joint distribution of possibly multivariate potential outcomes. fan2024multidimensional and lin2025estimation rely on the primal formulation of the bounds and the ruschendorf1991bounds method of conditioning: The lower bound on $\mathbb E[h(Y_0,Y_1,X)]$ is obtained by taking the minimum over joint distributions for $(Y_0,Y_1,X)$ with given multivariate marginals for $(Y_0,X)$ and $(Y_1,X)$. This situation is different from the standard coupling problem because of the overlapping marginals: $X$ is common to both multivariate marginals. The method of conditioning in ruschendorf1991bounds consists in writing
to remove the overlapping marginals. lin2025tightening propose an alternative approach to remove the overlapping marginals. They relax the problem to
under the constraint $\pi_{Y_0,X_0}\sim P_{Y_0,X}$ and $\pi_{Y_1,X_1}\sim P_{Y_1,X}$. This relaxation is similar to the approach to conditioning in li2022finite. ji2023model apply a conditional version of Kantorovich duality to derive and estimate sharp bounds on $\mathbb E[h(Y_0,Y_1,X)$. For instance, the lower bound is equal to
where the last inequality holds for any $\varphi_{0,x},\varphi_{1,x}$ such that
for all $x,y_0,y_1$. For each $x$, the dual functions solve
The dual formulation ((ref)) has two main benefits. First, the minimization problem over the joint distribution is turned into a problem involving only the marginals. Second, the inequality in ((ref)) is valid for any choice of function pair $(\varphi_{0,x},\varphi_{1,x})$ that satisfies the dual constraint ((ref)). Therefore the lower bound $\mathbb E[\varphi_{0,X}(Y_0)]+\mathbb E[\varphi_{1,X}(Y_1)]$ is valid (if potentially conservative) even if $(\varphi_{0,x},\varphi_{1,x})$ is not a pair of dual solutions, i.e., it doesn't solve ((ref)). Using this dual representation, ji2023model propose the following procedure to conduct inference on the lower bound $\theta_L$ (and symmetrically for the upper bound).
The treatment effect problems described in the previous subsection is a particular instance of data combination, since both potential outcomes are never observed for the same individual. A variety of other data combination problems can also be addressed with optimal transport methods. Various instances of the short and long regression problem of cross2002regressions can be solved with the monotone rearrangement inequalities of theorem ((ref)), as in fan2014identifying, gechter2024generalizing, d2024linear and meango2025combining, and with Kantorovich duality, as in meango2025combining. fan2024multidimensional provide a unified conditional optimal transport framework to address these and more general data combination problems. d2025partially derive sharp bounds in partially linear models with data combination using weak optimal transport.
Two recent applications of the short and long regression framework of cross2002regressions directly use optimal transport methods. The first, gechter2024generalizing, considers external validity of the results of randomized experiments. Population $e$ is treated at random, so $F^e_{Y_0}$ and $F^e_{Y_1}$ are identified. In population $a$ (for alternative), no one is treated, so $F^a_{Y_0}$ is identified. Call $f^e_0$ and $f^a_0$ the identified densities of $Y_0$ in the experimental and alternative populations respectively. What remains to be identified for the treatment effect in the alternate population is $\mathbb E^a[Y_1]$. Under the population similarity assumption $\mathbb E^a[Y_1\vert Y_0]=\mathbb E^e[Y_1\vert Y_0]$, we have
The monotone rearrangement theorem ((ref)) then provides sharp bounds on $\mathbb E^a[Y_1]$ as desired. In the second application of optimal transport to short and long regressions, meango2025combining propose to use stated preferences from surveys to identify revealed preferences from actual choices. Actual binary choices $D\in\{0,1\}$ are observed together with an endogenous driver of choice $X$ in the revealed preference data set. A different data set contains $X$ and the variable $P$, which is the stated expected probability of choosing $D=1$. The quantity of interest is the counterfactual (or potential) choice $D(x)$ when the driver of choice is externally set to $x$. The main identifying assumption is that stated preferences reveal the unobserved heterogeneity relevant to choices, i.e., $D(x)\perp\!\!\!\!\perp X\vert P$. Under this assumption,
Using the same strategy as in the external validity example, we can write the previous expression
and use the monotone rearrangement theorem to derive sharp bounds. However, the latter involve the quantile function of a ratio of density functions, which is not conducive to inference. meango2025combining therefore provide an alternative characterization of the bounds based on a direct application of Kantorovich duality.
fan2025partial develop a general framework that includes many data combination problems. They analyze the (partial) identification of finite dimensional $\theta$ in the moment equality model $\mathbb Em(Y_0,Y_1,X;\theta)=0$, where $m$ is a vector of $d_m$ moment functions. The distributions of $(Y_0,X)$ and $(Y_1,X)$ are identified in different data sets, but the joint distribution of $(Y_0,Y_1,X)$ is not. While they use the potential outcomes notation $(Y_0,Y_1)$ of treatment effects, their framework is more general. The identified set $\Theta_I$ is defined as the set of parameter values $\theta$ such that $\mathbb E_\pi m(Y_0,Y_1,X;\theta)=0$ for some joint probability distribution $\pi$, which is compatible with the overlapping marginals $P_{Y_0,X}$ and $P_{Y_1,X}$. By the ruschendorf1991bounds method of conditioning, the identified set is characterized by the existence of a joint probability distribution $\pi$ with marginals $P_{Y_0\vert X=x}$ and $P_{Y_1\vert X=x}$ such that
The latter is a problem of existence of a joint probability distribution. It can be transformed into an optimal transport problem, as in galichon2006inference and ekeland2010optimal (see section (ref) below). See the related result in theorem 1 of franguridi2025inference. Heuristically, the existence of a joint probability $\pi$ that makes ${\textstyle\int} m\;d\pi$ equal to zero is equivalent to $0$ being larger than $\inf_\pi{\textstyle\int} \lambda^\top m\;d\pi$ and smaller than $\sup_\pi{\textstyle\int} \lambda^\top m\;d\pi$ for any direction $\lambda$ in the unit sphere $\mathcal S^{d_m}$ ($d_m$ is the dimension of the vector $m$ of moment functions). Hence, the identified set $\Theta_I$ is also characterized by
for all $\lambda\in\mathcal S^{d_m}$. The special case of the framework in fan2024multidimensional, where the function $m(y_0,y_1,x;\theta)$ is linear in $\theta$ yields a useful characterization of the identified set $\Theta_I$ in a variety of empirically relevant settings. We illustrate the procedure and the characterization of the identified set $\Theta_I$ in the special case of the linear projection model. Let $X=(X_p,X_{np})$, where covariates in $X_p$ appear in the projection and $X_{np}$ collects additional covariates that do not appear in the projection. Consider the model with scalar $Y_1$ satisfying
The parameter of interest is $\vartheta$. The model can be rewritten in the general framework of fan2024multidimensional with
In this case, characterization ((ref)) can immediately be rewritten
for all $\lambda\in\mathcal S^{d_m}$ (here $d_m$ is equal to the dimension of $\theta$, which is also the dimension of $Y_0$). When $Y_0$ is scalar, the monotone rearrangement theorem ((ref)) applies directly to the term in brackets in ((ref)), which is equal to
When $Y_0$ is multivariate, it can be shown that the knowledge of the distribution of $Y_0$ gives no additional information relative to the knowledge of the distribution of $\lambda^\top Y_0$, so that bound ((ref)) is still sharp. See the discussion below theorem 1 in d2024linear. Bound ((ref)) leads directly to the characterization of the identified set for the parameter of interest $\vartheta$ in theorems 1 and 2 of d2024linear and proposition 4.2 in fan2024multidimensional. The monotone rearrangement theorem ((ref)) applies in ((ref)) because $Y_1$ is scalar. When $Y_1=(Y_{1l},Y_{1r}^\top)^\top$ is multivariate, as in the case of the linear projection model
the optimal transport characterization of the identified set still provides computational advantages. An important feature of all these bounds in linear projection models is that conditioning on the full vector $(X_p,X_{np})$ of covariates yields tighter bounds than conditioning only on $X_p$. Hence, covariates that are irrelevant in the point identified linear projection are no longer irrelevant in case of data combination.
d2025partially consider the partially linear regression model
In the display above, $Y$ is the scalar outcome, and $X_c$, $X_{nc}$ are vectors of covariates. The unknown function $f$ and finite dimensional parameter $\beta$ are the objects of inference. Both $Y$ and $X_c$ are observed in one data set, and $X_{nc}$ and $X_c$ are observed in another data set. Hence, the distributions of $(Y,X_c)$ and $(X_{nc},X_c)$ are identified, but the joint distribution of $(Y,X_{nc},X_c)$ is not. The identification procedure in d2025partially follows the steps:
The crucial step is the characterization of the identified set for $\beta$ in point 3 of the previous list. Start with the case $X$ is scalar. Assume $Y$ has mean zero without loss of generality. The random variables $Y$ and $X$ can be coupled in any way that preserves the marginal distribution and the martingale constraint $\mathbb E[Y\vert X]=\beta X$. The first observation is that we cannot expect much informativeness on $\beta$ with only marginal information. In particular, $X\perp\!\!\!\!\perp Y$ rationalizes the data. However, the model is not devoid of empirical content, and the bounds can be informative, when tightened with additional restrictions. With the notation $X^\beta:= \beta X$, the identified set is the set of $\beta$'s such that $\mathbb E[\widetilde Y\vert \widetilde{X^\beta}]=\widetilde{X^\beta}$ for some $\widetilde Y$ distributed like $Y$ and some $\widetilde{X^\beta}$ distributed like $X^\beta$. The marginal distributions are known for $Y$ and $X^\beta$. But nothing is known of their joint distribution, except the martingale constraint. This is classic a problem of coupling under constraints, which is equivalent to Lorenz dominance (see for instance chapter 17.C of marshall1979inequalities):
Of course, things are not that simple when $X$ is not scalar, because $\mathbb E[Y\vert X]$ is no longer the same as $\mathbb E[Y\vert \beta^\top X]$. There remains to show that the set of $\beta$'s such that $\mathbb E[\widetilde Y\vert \widetilde{X^\beta}]=\widetilde{X^\beta}$ for some $\widetilde Y$ distributed like $Y$ and some $\widetilde{X^\beta}$ distributed like $X^\beta$ is the same as the set of $\beta$'s such that $\mathbb E[\widetilde Y\vert \widetilde{X}]=\beta^\top\widetilde{X}$ for some $\widetilde Y$ distributed like $Y$ and some $\widetilde{X}$ distributed like $X$. d2025partially show this using the duality of weak optimal transport.\footnote{Beyond the scope of this review, weak optimal transport is also used to characterize labor market equilibrium in chone2021matching and paty2022algorithms.}\footnote{In related work, d2021rationalizing characterize rational expectations with the existence of a martingale coupling.}
The first class of applications of optimal transport methods concerns inference in parametric incomplete structural models considered in galichon2006inference (galichon2006inference,galichon2009test,galichon2011set). The vector of variables of interest $(Y,X,U)$ satisfies support constraint $(Y,X,U) \in \Gamma(\theta)$ almost surely, and $U$ has fixed and known distribution $Q_U$.\footnote{The distribution $Q_U$ may also dependent on an unknown finite dimensional parameter.} Both vectors or variables $Y$ and $X$ are observed, in the sense that available data consists in a sample $((Y_1,X_1),\ldots,(Y_n,X_n))$. Variables in vector $U$ are unobserved. Variables in vector $X$ are exogenous\footnote{The distribution of $U$ may also depend on $X$ as long as it is the conditional distribution $Q_{U\vert X}$ is known up to a finite dimensional parameter vector.} (in the sense that $U\perp X$), and there are no restrictions on the process generating $(X_1,\ldots,X_n)$. All endogenous variables are subsumed in vector $Y$. The model is incomplete in that multiple values of endogenous variables may be consistent with a single value of exogenous and unobserved variables. This can be seen in the fact that the set $\{ y: (y,x,u)\in\Gamma(\theta) \}$ may not be a singleton for all $(u,x)$. This corresponds to the fact that the model fails to produce a unique prediction. Recent empirical examples of such models include discrete choice with unobserved heterogeneity in consideration sets in barseghyan2021heterogeneous and market structure and competition in airline markets in ciliberto2021market.
The object of inference is the finite dimensional parameter $\theta\in\Theta$. The identified set for $\theta$ (also known as sharp identified region) is the set $\Theta_I$ of all $\theta$ such that the joint distribution $P_{Y,X}$ of the data is one of the data generating processes predicted by the model under $\theta$. As galichon2006inference (galichon2006inference) and ekeland2010optimal point out, this is equivalent to the existence of joint distribution with marginals $P_{Y,X}$ and $Q_U$ with support contained in $\Gamma(\theta)$, which in turn is equivalent to $0$ being the value of the optimal transport problem
This can be equivalently written
Since the cost function is an indicator function, we can directly apply theorem (ref) to conclude that ((ref)) is equivalent to
where $A^\Gamma=\{u\in\mathcal U:\,\exists y\in A,\,(y,x,u)\in\Gamma(\theta)\}$, and where the supremum is over all measurable subsets of the set of outcomes $y$. This gives a characterization of the identified set with a collection of conditional moment inequalities, called {\em Strassen-Artstein inequalities}. The same characterization of the identified set was obtained by beresteanu2011sharp using random set theory and the Artstein theorem (corollary 1.4.11 page 83 of molchanov2005theory). When the set $\mathcal Y$ of outcomes is finite, such as the set of strategy profiles in finite action finite player games, this duality is an instance of transformation of an infinite dimensional into a finite dimensional problem. ekeland2010optimal propose an extension of the Kantorovich duality to characterize the identified set in incomplete models, where the probability distribution of the latent variable is not specified, but is known to satisfy a set of moment constraints. schennach2014entropic and li2018identification offer alternative approaches to the same general question based on different optimization problems.
galichon2011set show that the identified set for $\theta$ can be characterized by the solution of a finite optimal transport problem when the set of outcomes is discrete. Let $\mathcal Y:=\{y_1,\ldots,y_L\}$ be the finite set of outcomes. Then the set of predicted outcomes (set of equilibria in a finite game) when $u$ ranges over its (usually continuous) domain is also finite. Call it $\mathcal U^\ast(x;\theta):=\{u^\ast_1(x;\theta),\ldots,u^\ast_K(x;\theta)\}$. Let $P_{Y\vert X=x}:=(p(y_1\vert x),\ldots,p(y_L\vert x))$ be the probability mass function of $Y$, and let $Q_{U^\ast\vert X=x;\theta}:=(q(u^\ast_1\vert x;\theta),\ldots,q(u^\ast_K\vert x;\theta))$, where for each $k$,
The characterization ((ref)) of the identified set for $\theta$ can now be rewritten as the finite optimal transport problem
where the minimization is over $\pi$ satisfying the margin constraints
for each $k\leq K$ and $l\leq L$.
li2022finite also apply optimal transport method to derive a finite sample valid inference method for parameter vector $\theta$ or a lower dimensional transformation of $\theta$ such as a subvector. They therefore consider the problem of testing $H_0(S):\Theta_I\cap S\ne\varnothing$, for some region $S$ of the parameter space. In the most common case, where a component, say $\theta_1$, of the vector of structural parameters, with true value $\theta_{01}$, is of interest, $S:=\{\theta\in\Theta: \theta_1=\theta_{10}\}$, and a confidence region for $\theta_1$ is obtained by inverting the test, i.e., including all values $\theta_{10}$ such the test fails to reject at the chosen significance level.
The test statistic in li2022finite is inspired by the characterization ((ref)) of the identified set. However, the conditioning over covariates $X$ is handled by replacing the indicator transport cost function $1_{\{(y,x,u)\notin\Gamma(\theta)\}}$ with a discrepancy $\delta$ between $(u,x)$ and $(y,x^\prime)$ which is non negative, lower semi-continuous and equal to zero if and only if $x=x^\prime$ and $(y,x,u)\in\Gamma(\theta)$. The idea of the discrepancy is to simultaneously penalize departures from the model structure and bad covariate matches. li2022finite recommend to construct the discrepancy as follows:
In the expression above, $\Sigma_U$ is the (known) covariance matrix of the random vector $U$ with distribution $Q_U$ and $\hat\Sigma_X$ the empirical covariance matrix of the sample $(X_1,\ldots,X_n)$. This relaxation approach to conditional optimal transport is similar to the one later proposed in lin2025tightening.
The test statistic in li2022finite is based on a discrete optimal transport problem based on ((ref)) with the indicator replaced with the discrepancy $\delta$. Then the value of the optimal transport problem, which is a function of $\theta$, is profiled by taking the infimum over all $\theta\in S$. The inner optimal transport problem is discrete. When the support of outcomes $y$ is finite, this can be achieved with the galichon2011set strategy described above. More generally, it is achieved with a low discrepancy sequence $(\tilde u_1,\ldots,\tilde u_n)$ that approximates the distribution $Q_U$. Call $\Pi_n$ the set of non negative matrices such that $\Sigma_i\pi_{ij}=\Sigma_j\pi_{ij}=1/n$, for all $i,j \leq n.$ The test statistic is defined as
Replace in ((ref)) the sample $(Y_1,\ldots,Y_n)$ with an arbitrary sequence $\tilde y:=(\tilde y_1,\ldots,\tilde y_n)$ of elements of the outcome space and call the resulting statistic $T_n(\tilde y;S)$. Call $\tilde U:=(\tilde U_1,\ldots,\tilde U_n)$ a sample of independent draws from $Q_U$. Critical values are obtained as quantiles of the statistic
In the construction of ((ref)), the null hypothesis is enforced because $(\tilde y_j,X_j,\tilde U_j)\in\Gamma(\theta')$ for some $\theta^\prime$ in $S$. The test has exact validity, because the supremum over $\theta^\prime$ and $\tilde y$ means that $\tilde T_n(S)$ achieves the worst case (largest in first order stochastic dominance) distribution for $T_n(S)$ under the null. In practice, the critical values are obtained numerically using a large number of replications of $\tilde T_n(S)$ each with an independent sample $\tilde U:=(\tilde U_1,\ldots,\tilde U_n)$.
chiong2016duality and bonnet2022yogurts formulate the estimation of discrete choice models as an optimal transport problem, which provides a computationally attractive solution to the demand inversion problem. Let agent $i$ choose option $y$ among a finite set $\mathcal Y$. Agent $i$ chooses $y$ to maximize utility
In the display above, the mean utility $\delta_y$ of option $y$ depends on observable characteristics $X=x$ of the collection of choices $\mathcal Y$, and $\varepsilon_{iy}$ is the random utility component. As in berry1993automobile, the conditional distribution of the random utility shocks $\varepsilon$ given $X=x$ is known and equal to $P_{\varepsilon\vert X=x}$. Market shares are known and represented by a probability mass function $q_{Y\vert X=x}$ over $\mathcal Y$ conditional on $X=x$. The demand inversion problem is the problem of finding the set of mean utility vectors $\delta(x)=(\delta_y(x))_{y\in\mathcal Y}$ compatible with market shares $q_{Y\vert X=x}$. Hence, the object of interest is the set
We let $c(y,\varepsilon_i\vert x_i)=\varepsilon_{iy}$, and define $G(\delta(x)):=\mathbb {\textstyle\int}[\max_y(\delta_y(x)+c(y,\varepsilon\vert x))]dP_{\varepsilon\vert X=x}(\varepsilon\vert X=x)$. By the envelope theorem $\delta(x)\in D(q_{Y\vert X=x})$ if and only if $q_{Y\vert X=x}$ is in the subdifferential $\partial G(\delta(x))$. By Theorem 23.5 page 218 of rockafellar1970convex, the latter is equivalent to $\delta(x)\in\partial G^\ast(q_{Y\vert X=x})$. Hence, by definition of the convex conjugate $G^\ast$,
Therefore $\delta(x)=-\psi(x)\in D(q_{Y\vert X=x})$ if and only if it minimizes
which is the dual of the optimal transport problem with cost $-c(y,\varepsilon\vert x)$. In this optimal transport problem, the transport potential $\psi(x)$ is equal to minus the mean utility $\delta(x)$. Hence, this optimal transport is a parametric conditional optimal transport problem, where the parametric model is imposed on the potential. galichon2022cupid estimate matching surplus in matching markets using a parameterization of the potential in an entropic optimal transport formulation of the problem.
This section reviews econometric applications of optimal transport maps. The aspects of optimal transport theory involved here relate to cyclical monotonicity of optimal transport plans and maps (definition (ref) and theorem (ref)), the uniqueness of optimal transport maps (theorem (ref)), and optimal transport theory with quadratic costs (theorem (ref)). Cyclical monotonicity provides a multivariate notion of monotonicity. This allows to define multivariate monotone rearrangements in ekeland2012comonotonic and multivariate comonotonicity in ekeland2012comonotonic and puccetti2010multivariate. Multivariate rearrangemenmts are used to define vector quantile functions in ekeland2012comonotonic, galichon2012dual, and chernozhukov2017monge (see also hallin2021distribution) and vector quantile regressions in carlier2016vector. In turn, vector quantile functions are used in econometrics to identify nonlinear simultaneous equations in chernozhukov2021identification and gunsilius2023condition, define robust risk measures in ekeland2012comonotonic and ruschendorf2012law, copulas between multivariate marginals in fan2023vector, multidimensional inequality measures in fan2022lorenz and hallin2025multiple, distributional difference-in-differences in torous2024optimal, principal component analysis in gunsilius2023independent and rank based distribution-free nonparametric inference in deb2023multivariate.
As explained in section (ref), the theory of optimal transport in $\mathbb R$ is closely related to the theory of monotone rearrangements. Consider the problem of transporting the uniform distribution $\mu$ on $[0,1]$ to a given distribution $\nu$ with quadratic cost $c(t,y)=|t-y|^2$ for $t\in[0,1]$ and $y\in\mathcal Y\subseteq\mathbb R$. Let $t\mapsto Q_\nu(t)$ be the quantile function of $\nu$. The cost function is submodular, hence by theorem (ref), the value of the optimal transport problem is
and the optimal transport map from the uniform $\mu$ to $\nu$ is simply $t\mapsto Q_\nu(t)$, i.e., the quantile function of $\nu$. The optimal coupling is the comonotone coupling $(T,Q_\nu(T))$, where $T$ is uniformly distributed and $Q_\nu(T)$ is the monotone rearrangement of $Y\sim\nu$ that minimizes the transport cost from $\mu$ to $\nu$.
Consider now the case of a distribution $\nu$ in $\mathbb R^d$. By theorem (ref), we know that there is a unique map $Q_\nu$ on $[0,1]^d$ that minimizes the optimal transport problem
which is the cost of transporting the uniform distribution $\mu$ on $[0,1]^d$ to $\nu$ with quadratic cost function $c(u,y):=\|u-y\|^2$. As an optimal transport map, $Q_\nu$ satisfies a multivariate notion of monotonicity, i.e., $c$-cyclical monotonicity of definition (ref). With quadratic cost $c(u,y)=\|u-y\|^2$, $c$-cyclical monotonicity boils down to the traditional cyclical monotonicity of rockafellar1966characterization, which requires
for any $K$, and any collection of vectors $(u_1,\ldots,u_K)$, setting $u_{K+1}=u_1$. Cyclical monotonicity characterizes maps that are the gradient of a convex function, which is another way to define multivariate monotonicity (since the derivatives of convex functions are the non-decreasing functions in $\mathbb R$). Hence, optimal transport maps can be interpreted as multivariate monotone rearrangements. Optimal transport maps from the uniform on $[0,1]^d$ can therefore define vector quantiles. If $\nu$ does not have finite second moments, then the unique gradient of a convex function $Q_\nu$ that transports mass from $\mu$ to $\nu$ can still define the vector quantile, even though the optimal transport problem with quadratic cost is no longer well defined. chernozhukov2017monge therefore define vector quantiles in the following way.
When $d=1$, $Q_\nu$ is the traditional quantile function (automatically satisfied when (1) and (2) hold), which explains why this notion is proposed as an extension of the traditional notion of quantile. When $d\geq 1$, and if, in addition, $\nu$ is absolutely continuous, $Q_\nu$ is invertible and $Q_\nu^{-1}(Y)$ is uniformly distributed on $[0,1]^d$, whenever $Y$ has distribution $\nu$. Thus, when $\nu$ is absolutely continuous, the map $Q_\nu^{-1}$ is a cyclically monotone map that transports $\nu$ into the uniform $\mu$ on $[0,1]^d$. Hence, it can be construed as a cardinal to ordinal transformation: It transforms a random vector $Y$ with distribution $\nu$ into a uniformly distributed random vector $U=Q_\nu^{-1}(Y)$ while minimizing distortion in the transformation (because of cyclical monotonicity). This motivates the definition of multivariate ranks in chernozhukov2017monge: The multivariate rank of vector $Y$ is $U=R_\nu(Y):=Q_\nu^{-1}(Y)$. This also motivates the extension of definition (ref) to random vectors:
Definition (ref) was orginally proposed by ekeland2012comonotonic. puccetti2010multivariate discuss alternative notions of multivariate comonotonicity, and torous2024optimal call it {\em cyclical comonotonicity} and use it to define an analogue of parallel trends for distributional difference-in-differences.
The vector quantile of a random vector with distribution $\nu$ can be estimated using the semi-discrete algorithm of section (ref) to compute the optimal transport map between the uniform distribution on $[0,1]^d$ and the empirical distribution $\nu_n$ of a sample $(Y_1,\ldots,Y_n)$. The vector rank map can be estimated with the optimal transport map from the empirical distribution $\nu_n$ to the empirical distribution $\mu_n$ of a Halton low discrepancy sequence of size $n$ on $[0,1]^d$. See section (ref) for the definition of the discrete optimal transport map and the simplex algorithm (ref) for its computation.
In addition to distributional difference-in-difference already mentioned, vector quantiles and ranks are applied to quantile regression in carlier2016vector, multivariate stochastic orders in ekeland2012comonotonic, galichon2012dual, charpentier2016local and fan2022lorenz, the measurement of risk and inequality in ekeland2012comonotonic, fan2022lorenz and hallin2025multiple, the characterization of dependence between random vectors in fan2023vector, distribution-free multivariate inference in hallin2021distribution and deb2023multivariate, the identification of nonlinear simultaneous equations in chernozhukov2021identification and nonlinear principal component analysis in gunsilius2023independent. We examine each in turn.
As noted before, in econometric applications, we are often interested in the conditional version of the optimal transport problem. As vector quantiles are defined based on theorem (ref), conditional vector quantiles can be defined based on a conditional version of theorem (ref). We call conditional vector quantile of random vector $Y$ on $\mathbb R^d$ conditional on $X=x$ the almost surely unique cyclically monotone map $Q_{Y\vert X=x}$ such that $Q_{Y\vert X=x}(U)$ has the same distribution as $Y\vert X=x$ for any uniformly distributed $U$ on $[0,1]^d$. carlier2016vector propose a parametric version of this conditional vector quantile, which extends traditional quantile regression to the multivariate outcome case. The parametric vector quantile takes the form
where $x$ is a $p$-valued vector of covariates (including an intercept), and $\beta(u)$ is the $p\times d$ valued matrix of regression coefficients and $\beta(u)^\top x$ is constrained to be cyclically monotone in $u$ for each $x$. When $d=1$, ((ref)) reduces to the traditional quantile regression. We now explain how the optimal transport formulation allows efficient computation of the regression weights $\beta(u)$. By definition of the vector quantile, we know from theorem (ref) that $Q_{Y\vert X=x}(u)=\nabla\varphi(u,x)$, where $\varphi$ solves the semi-dual optimization program
where the supremum is taken over all convex functions, and where $\varphi^\ast(y,x)$ is understood as the convex conjugate at $y$ of $u \mapsto \varphi(u,x)$ for fixed $x$. Under the assumption that $Q_{Y\vert X=x}(u)=\beta(u)^\top x$, we have $\nabla\varphi(u,x)=\beta(u)^\top x$, hence the solution to the semi-dual program ((ref)) is $\varphi(u,x):=B(u)^\top x$, where $B$ solves the following program among functions $B$ such that $B(u)^\top x$ is a convex function of $u$.
After integration over $X$, the function $B$ solves
Finally, the solution is $\beta(u)=D B(u)$, the Jacobian of $B$.
Multivariate quantiles and ranks are natural building blocks of multivariate risk and inequality measures. Consider a random vector $Y$ on $\mathbb R^d$ with probability distribution $\nu$ and vector quantile function $Q_\nu$. In risk analysis, vector $Y$ is interpreted as a vector of uncertain prospects or financial exposures. In inequality analysis, the vector $Y$ is interpreted as the vector of attributes (wealth, health, education level) of an individual randomly drawn from the population of interest. A risk or inequality measure is a functional $\varrho$ that takes random vector $Y$ and returns a real number. We consider only {\em law invariant} functionals, i.e., functionals that only depend on the distribution of $Y$ or, equivalently, on its vector quantile function. galichon2012dual show that functionals $\varrho$ are comonotonic additive, i.e., $\varrho(Y+Y^\prime)=\varrho(Y)+\varrho(Y^\prime)$ for comonotonic $Y$ and $Y^\prime$ if and only if they are of the form $\varrho(Y)={\textstyle\int_{[0,1]^d}}\left(\phi(u)^\top Q_Y(u)\right)du$, for some vector function $\phi:\mathbb R^d\rightarrow\mathbb R^d$. Such a functional $\varrho$ is called a rank dependent functional, since it is a weighted average of vector ranks only. Rank dependent risk and inequality measures are of this form.
charpentier2016local extend the Bickel-Lehmann dispersive order and its characterization in landsberger1994co. A distribution $\nu$ is more dispersed than a distribution $\nu^\prime$ if there are random vectors $Y\sim\nu$, $Y^\prime\sim\nu^\prime$ and $Z$ vector comonotonic with $Y^\prime$ such that $Y=Y^\prime+Z$. They show that it is equivalent to cyclical monotonicity of the difference in quantile maps $Q_\nu-Q_{\nu^\prime}$. fan2022lorenz extend the Lorenz order of increasing inequality to multi-attribute inequality. A population with resource distribution $\nu$ is more unequal than a population with resource distribution $\nu^\prime$ if
for all $r=(r_1,\ldots,r_d)\in[0,1]^d$. Similarly to the traditional scalar case, this has the following simple interpretation. For any rank in the population, the cumulative resource in each attribute held by all individuals below that rank is larger in the less unequal population. fan2022lorenz show that the increasing inequality ordering shares the traditional interpretation in terms of inequality increasing transfers and preference by any rank dependent inequality measure.
fan2023vector propose an extension of the Sklar theorem and a definition of copula to characterize dependence between two or more random vectors. Let $Y=(Y_1^\top,\ldots,Y_K^\top)^\top$ be a random vector with distribution $\nu$ on $\mathbb R^{d_1\times \ldots\times d_K}$. Each component $Y_k$ is a random vector on $\mathbb R^{d_k}$ with vector rank map $R_k$. Then there is a unique probability distribution $\nu_C$ such that for each collection of measurable subsets $(A_1,\ldots,A_K)$
The probability distribution $\nu_C$ on $[0,1]^{d_1\times \ldots\times d_K}$ characterizes the dependence between random vectors $Y_1,\ldots,Y_K$ independently of each of their multivariate marginals.
deb2023multivariate propose rank based inference for mutual independence and for equality of two distributions. Tests are exactly distribution-free, which means that the null distribution of the test statistics are free of the unknown data generating process at all sample sizes. Let $(Y_1,\ldots,Y_n)$ be a sample of independently and identically distributed random vectors with distribution $\nu$ on $\mathbb R^d$. The vector rank map $R_\nu$ is the optimal transport map from $\nu$ to the uniform distribution $\mu$ on $[0,1]^d$. Empirical ranks are defined in the following way. Let $\mathcal V_n:=(v_1,\ldots,v_n)$ be a low discrepancy sequence, i.e., a deterministic sequence of points in $[0,1]^d$ such that the empirical distribution $\mu_n:=\Sigma_i\delta_{v_i}/n$ approximates $\mu$ and converges in distribution faster than the empirical distribution of a random sample. A popular choice is a Halton sequence. The empirical rank map $R_{\nu,n}$ is the optimal transport map from $\nu_n$ to $\mu_n$, i.e., it solves ((ref)) in section (ref). Finally, define order statistics $(Y_{(1)},\ldots,Y_{(n)})$ as any arbitrary ordering of the sample (for instance ordering with respect to the first coordinate). Then, the following theorem (proposition 2.2 in deb2023multivariate and proposition 1.6.1 of hallin2021distribution) shows that multivariate ranks inherit the distribution-free property of univariate ranks.
Using theorem (ref) and the combinatorial central limit theorem of hoeffding1951combinatorial, deb2023multivariate propose a test of independence of two random vectors, and a test of the hypothesis that two random vectors have the same distribution. In both cases, the test statistic is a function of empirical vector ranks only, and the limiting distribution is of the form $\Sigma_{j=1}^\infty\eta_jZ_j^2$, where the $Z_j$'s are standard normals and the weights $\eta_j$ do not depend on the underlying distributions, or the choice of low discrepancy sequence $\mathcal V_n$.
In traditional systems of $d$ linear simultaneous equations, the outcome variable $Y$ is a random vector in $\mathbb R^{d_y}$ that satisfies $Y=AX+\varepsilon$, where $A$ is a matrix of parameters, $X$ is a vector of covariates and $\varepsilon$ is a random vector in $\mathbb R^d_{\varepsilon}$ satisfying $\mathbb E[\varepsilon\vert X]=0$. We consider identification of the nonparametric non additively separable extension $Y=F(X,\varepsilon)$. Without further assumption, the distribution of $\varepsilon$ and the function $F$ cannot be separately identified. More precisely, for any distribution $P_{Y\vert X=x}$, and any pair of absolutely continuous distributions $(P_{\varepsilon\vert X=x},\tilde P_{\tilde\varepsilon\vert X=x})$, we know by mccann1995existence that there exists an invertible map $T$ such that $F(x,\varepsilon)$ and $\tilde F(x,\tilde\varepsilon):=F(x,T(\tilde\varepsilon))$ have the same distribution $P_{Y\vert X=x}$. Therefore, $(F,P_{\varepsilon\vert X=x})$ and $(\tilde F,\tilde P_{\tilde\varepsilon\vert X=x})$ are observationally equivalent. Hence, we normalize the conditional distribution of $\varepsilon$ to be a fixed known distribution $P_{\varepsilon\vert X=x}$. In case $d_z=d_\varepsilon$, function $F$ is identified as the unique cyclically monotone map that transports probability distribution $P_{\varepsilon\vert X=x}$ to $P_{Y\vert X=x}$. This extends the nonparametric identification result in matzkin2003nonparametric to $d_y=d_\varepsilon>1$. Note that the distribution of $Y$ need not be absolutely continuous. More generally, chernozhukov2021identification also show nonparametric identification of the model $\varepsilon=H(Y,X)$ in case $d_Y\geq d_\varepsilon\geq1$. The condition for identification of $H(Y,X)$ is that $Y$ is obtained as the hedonic demand of a consumer with type $(x,\varepsilon)$ and the identification is based on theorem (ref).
gunsilius2023independent propose an alternative to principal component analysis based on optimal transport. Consider a random vector $Y$ in $\mathbb R^d$ with absolutely continuous distribution $\nu$ and a predetermined desired number $k<d$ of principal factors. The proposed procedure is as follows:
The basis vectors $(v^1,\ldots,v^k)$ are the principal directions. The vector $Y^{*k}$ is the $k$-factor approximation of $Y$. The linear transformation in PCA is replaced with the optimal transport map $T_\nu$ (if $Y$ is normal, they coincide) and entropy replaces variance as a way to quantify the relevance of components.
The Wasserstein distance induces a geometry on the space of probability distributions. This geometry inherits features from the geometry of the base space, typically $\mathbb R^d$. Heuristically, this is because it is more costly to transport mass far away than to transport it close by. As a result, the Wasserstein distance is useful to model robustness concerns to mis-specification of baseline distributions, in blanchet2019quantifying, adjaho2022externally, kido2022distributionally, gu2023dual, and fan2025quantifying, and measurement error in schennach2022optimally and forneron2024fitting. The Wasserstein distance is based on shifts of mass around the baseline space, so it naturally allows the comparison of probability distributions with different supports. This is particularly useful in applications to treatment effect problems with limited overlap between the supports of covariate distributions for the treated and control populations. The Wasserstein barycenter is used in gunsilius2023distributional to define distributional synthetic controls and applied to distributional regression by zhu2023autoregressive. The Wasserstein distance has computational advantages that make it useful to simulate from distributions, as shown in arellano2023recovering and athey2024using, and to estimate parametric models in kaji2023adversarial and fan2024minimum.
The prototypical distributional robust optimization framework features a reference probability distribution $\mu$, which is known (possibly only up to a finite dimensional parameter vector), and an adversarial view of nature, which is supposed to choose the true data generating process $\nu$ within a Wasserstein neighborhood. In this section, we will call Wasserstein neighborhood of a probability distribution $\mu$ a set of probability distributions $\nu$ such that
for some $\rho>0$ and some cost function $\delta$, in most cases a metric on $\mathcal Y\subseteq\mathbb R^d$.
In blanchet2019quantifying, the parameter of interest is the expectation of a function $f$, so that the distributional robust optimization problem is specified as
Letting $\lambda$ denote the Lagrange multiplier associated with the Wasserstein constraint, the Lagrangian formulation of the constrained optimization problem is
The inner optimal transport problem has Kantorovich dual
and the semi-dual is achieved by taking $\varphi(y)=0$ and $\psi(\tilde y)=\sup_y(f(y)-\lambda\delta(y,\tilde y))$. Hence, the original problem is equal to
blanchet2019quantifying apply this duality result to the computation of worst-case probabilities, i.e., probability distributions $\nu_{\tiny\mbox{max}}$ and $\nu_{\tiny\mbox{min}}$ in the Wasserstein neighborhood $W_\delta(\nu,\mu)\leq\rho$ that achieve the maximum and the minimum of ${\textstyle\int} fd\nu$ respectively.
In econometric applications of distributionally robust optimization, the problem of conditioning on exogenous covariates is addressed in a variety of ways. adjaho2022externally consider the problem of choosing treatment assignment $\tau(x)\in\{0,1\}$ to maximize the utilitarian welfare under adversarial shifts in the joint probability distribution $P_{Y_0,Y_1,X}$ for the scalar potential outcomes $Y_0$, $Y_1$, and the vector of covariates $X$. They assume the planner considers a potentially mis-specified reference joint probability distribution $P_{Y_0,Y_1,X}$. This can be obtained for instance with access to a pilot treatment data set with selection on observables and an assumption of either perfect correlation of $Y_0$ and $Y_1$ conditional on $X$, or independence of $Y_0$ and $Y_1$ conditional on $X$. The planner then seeks robustness to the misspecification of the joint probability (in particular robustness to the assumption of comonotone or independent coupling) and adjaho2022externally define the robust welfare criterion as the maximin response to an adversarial nature. Nature is assumed to choose the worse case data generating process $\nu$ in the Wasserstein neighborhood $W_\delta(\nu,P)$ with
The robust welfare criterion is therefore defined as
Using the equality between ((ref)) and ((ref)) above, and the fact that $y+\lambda|y-\tilde y|$ is minimized by setting $y=-\infty$ if $\lambda\in[0,1)$ and $y=\tilde y$ if $\lambda\geq1$, the robust welfare criterion is equal to
fan2025quantifying consider a similar robustness concern in the allocation of treatment. However, they do not assume that the planner has a reference joint probability distribution fo $(Y_0,Y_1,X)$. Instead, there is a reference probability distribution for $(Y_0,X)$ and $(Y_1,X)$, justified for instance with access to a pilot treatment data set with selection on observables, but no additional information on the joint probability distribution. fan2025quantifying derive bounds for the expectation of a function of $Y_0,Y_1,$ and $X$. They combine the method of conditioning of ruschendorf1991bounds to deal with the overlapping marginals with the distributional robust optimization approach of blanchet2019quantifying to account for the uncertainty about the marginals. The method applies to a large variety of settings, including ones with multivariate potential outcomes. However, we choose to illustrate it on a special case with scalar potential outcomes. We compare the result in fan2025quantifying to adjaho2022externally in the special case of robust treatment assignment. For $l\in\{0,1\}$, define the following quantities. Let $\mu_l$ be the reference distribution for $(Y_l,X)$. Let the Wasserstein neighborhood $W_\delta(\cdot,\mu_l)\leq \rho_l$ of $\mu_l$ be defined with radius $\rho_l>0$ and
The robust treatment criterion is
where the infimum is now taken with respect to probability distribution $\pi\in\mathcal M(\nu_0,\nu_1)$ for $(Y_0,Y_1,X)$ with overlapping marginals $\nu_0$ for $(Y_0,X)$ and $\nu_1$ for $(Y_1,X)$, and each $\nu_l$ lies in the Wasserstein neighborhood $W_\delta(\nu_l,\mu_l)$. Using the same technique, the robust welfare criterion is shown to equal
with $\phi_\lambda^l:=\inf_{\tau(x)=l}\{\lambda_0\|x_0-x\|+\lambda_1\|x_1-x\|\}$. The difference with the adjaho2022externally case is the data combination problem and the duplication of covariates to deal with overlapping marginals.
gu2023dual add distributional robustness concerns to the class of incomplete models described in section (ref).\footnote{The sensitivity analysis in gu2023dual is related to prior work by christensen2023counterfactual, where the distributional robustness neighborhood are characterized by $f$-divergences instead of the Wasserstein distance.} Observed variables are collected in $(Y,X)$. In the identification analysis below, we assume the probability distribution $P_{Y,X}$ generating them is known. There is also a vector $U$ of latent variables with distribution $Q_U$. The variables in vector $X$ are exogenous in the sense that $U$ is independent of $X$. The model structure is characterized by the support restriction $(Y,U,X)\in\Gamma(\theta)$. We are interested in a functional $\mathbb E[\varphi(Y^\ast,U,X;\theta)]$ of counterfactual outcomes $Y^\ast$, which satisfy a counterfactual variant of the structural model, namely $(Y^\ast,U,X)\in\Gamma^\ast(\theta)$. Both actual structure $\Gamma(\theta)$ and counterfactual structure $\Gamma^\ast(\theta)$ are known up to the finite dimensional unknown parameter vector $\theta$. Because of distributional robustness concerns, the distribution $Q_U$ of the latent vector $U$ is only assumed to lie in a Wasserstein neighborhood of a reference distribution $Q_V$: $W_1(Q_U,Q_V)\leq\rho$. By definition of the Wasserstein distance $W_1$, there is a coupling $(U,V)$ of the vectors $U$ and $V$ with respective distributions $Q_U$ and $Q_V$ such that $\|U-V\|\leq\rho.$ Hence, the lower bound for the parameter of interest is given by
where the inner infimum is over all probability distributions $\pi$ for $(Y^\ast,Y,U,V,X)$ satisfying the following constraints:
The inner optimization program in ((ref)) has Lagrangian formulation:
subject to (1), (2) and (3) above. Call
the collection of values of the latent variables $(y^\ast,u)$ compatible with observations $(y,x)$. Define the random set
as the set of values that the function $\varphi(y^\ast,u,x;\theta)+\lambda(\|u-v\|-\rho)$ can take when latent variables vary over $\mathcal G(y,x;\theta)$ and $v$ is free. The inner optimization problem in ((ref)) is equal to the infimum over the Aumann expectation of the random set $\Phi$, i.e., over the collection of integrals ${\textstyle\int_\pi}\left[\varphi(y^\ast,u,x;\theta)+\lambda(\|u-v\|-\rho)\right]\,d\pi(y,v,x)$ subject to (2) and (3) when latent variables vary over $\mathcal G(y,x;\theta)$ and $v$ is free. Under the assumptions of theorem 2.1.20 page 236 of molchanov2005theory, the infimum of the Aumann expectation is the expectation of the infimum over $\mathcal G(y,x;\theta)$. Hence, calling
the inner optimization problem in ((ref)) is equal to the optimal transport problem with cost $c((y,x),v;\theta,\lambda)$:
Finally, the sharp lower bound for the parameter of interest $\mathbb E[] \varphi(Y^\ast,U,X;\theta)]$ is therefore equal to $\inf_\theta \sup_{\lambda\geq0}\mathcal C(\theta,\lambda)$.
Xu2025distributionally study treatment effect prediction under population shift by defining the target object through a distributionally robust optimization problem. Let $P$ denote the source distribution of potential outcomes and let $\mathcal P_\delta = \{ Q : W_p(Q,P) \le \delta \}$ be a Wasserstein neighborhood. The robust treatment effect is defined as the minimax solution \[ \tau_\delta^\ast = \mbox{arg}\min_{\tau} \sup_{Q \in \mathcal P_\delta} \mathbb E_Q \bigl[(Y(1)-Y(0)-\tau)^2\bigr]. \] Because the joint distribution of $(Y(1),Y(0))$ is not identified, the problem is coupled with partial identification over the set of copulas consistent with the marginal potential outcome distributions, yielding sharp optimistic and pessimistic minimax solutions. The resulting robust predictor coincides with the source average treatment effect when $\delta=0$ and exhibits shrinkage toward zero as $\delta$ increases.
schennach2022optimally propose to model robustness to measurement error using the Wasserstein distance. Suppose the data $(Y_1,\ldots,Y_n)$ is a sample of possibly mismeasured version of a latent variable $Z$, which satisfies moment conditions $\mathbb Eg(Z,\theta)=0$. The quantity of interest is the finite dimensional parameter $\theta$. schennach2022optimally propose the optimally transported generalized method of moments (OTGMM) estimator for $\theta$ which minimizes the solution of the following program:
For each value of $\theta$, the program finds a set of points $(\tilde Y_1,\ldots,\tilde Y_n)$ that satisfy the empirical moments, while minimizing the 2-Wasserstein distance between the empirical distribution $\Sigma_{i=1}^n\delta_{\tilde Y_i}/n$ relative to $(\tilde Y_1,\ldots,\tilde Y_n)$ and the empirical distribution $\Sigma_{i=1}^n\delta_{Y_i}/n$ relative to the actual sample. Where generalized empirical likelihood estimators re-weight observations to minimize KL discrepancy to account for sampling bias, OTGMM adjusts sample points to minimize 2-Wasserstein distance to account for measurement error. forneron2024fitting apply related ideas to robustness in dynamic state space models. daljord2021black also use an optimal transport approach to the challenge of measuring illicit trade volume.
The geometry induced by the Wasserstein distance in the space of probability measures makes it possible to work with averages of distributions for distributional regression and synthetic controls in particular. gunsilius2023distributional extends the synthetic control methodology to estimate distributional causal effects. In the traditional setting, the analyst observes an aggregate outcome $Y_{jt}$ for each unit $j=1,\ldots,J$ and time period $t=1,2$. Unit $J$ receives a treatment at period $2$ only. All other units are never treated. This is easily generalizable to more treatment units and more periods. In case we have access to disaggregated data within each unit, we can identify and estimate the probability distribution $P_{jt}$ that generates $Y_{jt}$ for each unit and time period. The synthetic control distribution $P_{J2}^N$ is a weighted average of never treated units. The averaging is with respect to Wasserstein distance: It is a Wasserstein barycenter as defined in ((ref)). The weights are chosen so that the synthetic control distribution in period $1$ is closest in Wasserstein distance to the distribution for unit $J$ in period $1$.
zhu2023autoregressive use the Wasserstein geometry to define autoregressive processes for distribution regression.
The Wasserstein distance can be used to simulate data. One example is the filtering of the latent variable $U_t$ in forneron2024fitting described above. In arellano2023recovering, observed outcomes $Y$ are assumed to be the sum $Y=\eta+\eta^\prime$ of two latent factors of interest $\eta$ and $\eta^\prime$. The latent factors are the objects of interest and we wish to simulate samples of the latter. A sample of observable outcomes $(Y_1,\ldots,Y_n)$ is available. Let $\sigma_\eta$ and $\sigma_\eta^\prime$ be independent permutations. We want to generate samples $(\hat \eta_1,\ldots,\hat\eta_n)$ and $(\hat \eta_1^\prime,\ldots,\hat\eta_n^\prime)$ to minimize the 2-Wasserstein distance
where the regularization (nonparametric) set $\mathcal E$ is defined by $|\eta| \leq\bar C_n$ and
athey2024using propose to use Wasserstein GANs (described in section (ref)) to produce synthetic data that closely mimic a given real data set.
Wassertsein GANs can also be used in estimation. kaji2023adversarial propose adversarial estimation\footnote{See also kaji2021adversarial.}. They use traditional GANs, but their method can be trivially adapted to Wasserstein GANs. Let $P_\theta$ be the parametric model entertained, such that data can be simulated via $g_\theta(Z_i)$, where $Z_i$ follows a given reference distribution (as explained in section (ref)). Let $P_{\theta,m}$ be the empirical distribution relative to a simulated sample $(g_\theta(Z_1),\ldots,g_\theta(Z_m))$. The adversarial estimator for $\theta$ is the minimizer $\hat\theta$ of the 1-Wasserstein distance between $P_{\theta,m}$ and the empirical distribution $P_n$ of the data sample $(X_1,\ldots,X_n)$:
In kaji2023adversarial, $n/m\rightarrow0$, to avoid needlessly suffering from simulation noise when $n$ is small. fan2024minimum also use a minimum Wasserstein distance estimation principal. However, to circumvent the curse of dimensionality, they propose to minimize the sliced Wasserstein distance defined in (ref). With the notation of section (ref), the sliced Wasserstein distance between the empirical distribution $\mu_n$ of a sample of observations and the parametric distribution (or the empirical distribution of a simulation sample thereof) is
The MSD (Minimum Sliced Wasserstein Distance) estimator $\hat\theta$ is defined by
One notable feature of the minimum sliced Wasserstein estimator is that it is asymptotically normal even when the support of the data generating process depends on the parameter. In a similar spirit, qu2024distributionally propose a distributionally robust instrumental variable estimation, which solves
where $\tilde P_n$ is the empirical distribution of $(\Pi_ZX_i,\Pi_ZY_i)_i$.
A final application of the Wasserstein distance to inference problems can be adapted from galichon2013dilation. The latter use a different metric on the space of probability measures, but the problem can be adapted in the following way. In the framework of section (ref), we now consider inference on the basis of an empirical distribution $P_{Y,X;n}$. Hence the true distribution $P_{Y,X}$ is unknown. A confidence region for the parameter vector $\theta$ can be characterized by
In the expression above, $\rho_n$ is calibrated to the desired confidence level.
While this article is dedicated to the use of optimal transport {in} econometrics, it is particularly relevant to discuss the {econometrics of optimal transport} itself, or, equivalently, the {econometrics of matching models with transferable utility}. In the transferable-utility (TU) framework, the total utility generated by a match between types $x$ and $y$ can be summarized by a {joint surplus} $\Phi_{xy}$, which is split between the two agents through transfers. Empirically, the econometrician often observes the {equilibrium matching pattern} between two finite populations---firms and workers, men and women, origins and destinations in migration flows---but the underlying systematic surplus $\Phi$ (or the corresponding transportation cost $c(x,y)=-\Phi_{xy}$) that rationalizes this allocation is not directly observed. Recovering this economic primitive from the observed allocation is the object of {Inverse Optimal Transport} (IOT). Because structural matching models with transferable utility are fundamentally optimal-transport problems, and equilibrium allocations in these models are the solutions to such problems, IOT provides the natural framework for structural identification and estimation.
This perspective has proved empirically powerful in a wide range of applications. In the marriage market, choo2006marries use a TU matching model with unobserved heterogeneity to measure the gains to marriage over time and to assess how the U.S.\ Supreme Court’s ruling in {Roe v.\ Wade} affected the “value of marriage” implicit in observed sorting patterns. dupuy2014personality construct {indices of mutual attractiveness} that summarize, in a low-dimensional way, how a rich set of characteristics (education, height, BMI, health, risk attitudes, personality traits) contributes to sorting, by estimating the surplus function in a TU matching model. In a related vein, chiappori2012fatter investigate how physical appearance and body shape affect the structure of matching gains, showing that the valuation of body size across genders interact with that of education and income to shape systematic sorting in marriage markets.
CiscatoGalichonGousse2020 structurally compare homogamy across same-sex and different-sex couples, quantifying how assortative mating on education, age, race, and wages differs across these groups. ChiapporiSalanieWeiss2017 study the {marital college premium}, showing how changes in education levels and educational assortative matching translate into changes in the returns to college in the marriage market. More recently, ChiapporiFiorioGalichonVerzillo2025 exploit an extremely rich administrative income dataset to measure assortative matching on income, using a flexible TU matching framework that can accommodate complex sorting patterns and their implications for household-income inequality. Beyond household and labor applications, TU matching methods have also been applied to corporate finance: GuadalupeRappoportSalanieThomas2023 use a structural assortative-matching framework to study mergers and acquisitions, showing how the matching between targets and acquirers reflects systematic complementarities in firm characteristics. Related TU matching models have also been applied to industrial organization contexts, for instance in the estimation of matching games with transfers between firms and their trading partners or suppliers Fox2018MatchingTransfers, underscoring the ubiquity of IOT-type questions whenever economic agents sort into relationships on the basis of underlying complementarities.
Throughout this section we place ourselves in the {discrete} setting with the notation of section (ref). In economic applications, it is often more natural to work not with the cost but with the {systematic surplus} $\Phi_{xy} = -c(x,y),$ since \(\Phi_{xy}\) represents the deterministic component of the joint utility that two agents of types \(x\) and \(y\) derive from matching. Given a surplus matrix \(\Phi=(\Phi_{xy})\), the classical (unregularized) optimal transport problem ((ref)) is repeated here for convenience:
where \(\Pi(\mu,\nu)\) denotes the transportation polytope, which is the set of of feasible couplings between \(\mu\) and \(\nu\), that is, the set of $\pi_{xy} \geq 0$ such that $\sum_y \pi_{xy}=\mu_x$ and $\sum_x \pi_{xy} = \nu_y$. In the inverse optimal transport, the econometrician observes the matching pattern \(\widehat{\pi}\) and seeks to infer the surplus matrix \(\Phi\) that rationalizes it.
IOT is a notoriously hard problem in its unregularized form. A detailed account of the difficulties associated with unregularized assignment and transportation problems can be found in section 6.7 of the monograph of burkard2012assignment, which documents the combinatorial structure and degeneracies of assignment polytopes.
The fundamental difficulty is simple to state: {unregularized optimal transport typically produces extremely sparse optimal solutions}. In the fully balanced assignment case ($\mu_x=\nu_y=\tfrac{1}{N}$), the optimal solution is always a {vertex} of the Birkhoff polytope of doubly stochastic matrices and therefore corresponds to a {permutation matrix}. More generally, the optimal solution is an {extreme point} of the transportation polytope (see section (ref)), and extreme points have support size at most $M+N-1$. This means the mass of $\pi$ concentrates on a sparse set of pairs.
In many important cases in the continuous setting too --- for instance, when the cost function $c$ is quadratic --- the optimizer is a {Monge map}: a deterministic matching function $y=T(x)$. Thus the mapping from $c$ to $\pi$ is highly discontinuous and set-valued. Inverting it is extremely challenging, and in most cases impossible without strong structural restrictions.
A classical example from economics makes the identification problem vivid. In Becker’s model of marriage markets, suppose types are ordered scalars and the matching is positively assortative. Then any {supermodular} surplus function $\Phi$ generates positive assortative matching as the unique stable (and optimal transport) outcome. Because the class of supermodular functions is extremely large, the observed matching pattern carries very little identifying information: the same allocation can be rationalized by an infinite-dimensional family of surplus functions. Thus, in the absence of regularization, the IOT problem is fundamentally {under-identified}. Many surplus functions $\Phi$ (or cost functions $c=-\Phi$) produce the same optimal transport plan. This motivates the need for additional structure to restore identifiability and to regularize the inverse problem. In practice, this structure may come from parametric restrictions on the surplus $\Phi$ (for example, bilinear or low-rank specifications), from sparsity-promoting priors on its parameters, or---most powerfully—from {entropic regularization}. Entropic regularization plays a dual role: it introduces a smooth, strictly convex component into the forward optimal transport problem, which makes the mapping $c\mapsto\pi$ single-valued and differentiable; and it has the natural economic interpretation of capturing the effect of unobserved heterogeneity in preferences. In the next subsection we formalize this regularized formulation of IOT, and we show how it yields a well-posed econometric problem whose structure is closely connected to Poisson pseudo–maximum likelihood (hereafter PPML) and to gravity models of bilateral trade flows.
We begin by formally defining the regularized IOT problem in the discrete setting. Let $\Phi_{xy}=c(x,y)$ be the surplus matrix, and let $\pi^\sigma(\Phi)$ denote the unique solution of the entropically regularized optimal transport problem
with marginals $(\mu,\nu)$.
The entropic term makes the forward map $\Phi \mapsto\pi^\sigma(\Phi)$ smooth, strictly convex, and single-valued, turning the inverse problem into a well-posed estimation problem. Crucially, entropic OT always yields a coupling with {full support}, i.e., $\pi^\sigma_{xy}>0$ for all $(x,y)$, which eliminates the sparsity and non-invertibility issues inherent in classical (unregularized) IOT.
The model of choo2006marries provides a prototypical empirical setting in which entropic regularization arises naturally. There are finite type sets $\mathcal{X}$ and $\mathcal{Y}$, and agents on both sides may remain unmatched. The systematic surplus shared by types $(x,y)$ is \[ \Phi_{xy}=U_{xy}+V_{xy}, \] where $U_{xy}$ and $V_{xy}$ are deterministic components of utility for each side. The systematic utilities of unmatched individuals are normalized to zero, so $U_{x0}=V_{0y}=0$. Additive i.i.d. Gumbel shocks generate a multinomial-logit structure.
Let $m_x$ and $m_y$ denote the masses of women and men of each type. The equilibrium matching flows $\pi_{xy}$ (i.e., mass of $(x,y)$ pairs, including $x=0$ and $y=0$) satisfy
where the utilities of unmatched individuals are normalized to zero. Eliminating type-specific utilities yields the well-known choo2006marries identification formula:
which allows to identify the entire $\Phi$ directly from matching flows. Moreover, as shown by galichonSalanie2009,galichon2022cupid, the equilibrium is the solution of the {unbalanced} entropically regularized optimal transport problem, which is a convex optimization problem:
decker2013unique study the uniqueness of the equilibrium and substitution effects in this model.
dupuy2014personality develop a continuous-logit matching model without singlehood, where every agent is matched. Let $\mathcal X$ and $\mathcal Y$ be continuous type spaces with densities $\mu(x)$ and $\nu(y)$. With i.i.d. Gumbel shocks, the equilibrium matching density is:
where $a(x)$ and $b(y)$ are endogenous potentials ensuring the correct marginals. Taking logs yields the identification formula
so $\Phi$ is identified up to the additive fixed effects $\varphi(x)$ and $\psi(y)$.
A central contribution of galichon2022cupid is to show that the convex-analytic formulation above extends far beyond the Gumbel (logit) specification. Let $\varepsilon=(\varepsilon_{y})_{y\in\mathcal{Y}\cup\{ 0\} } $ denote the vector of taste shocks for a type-$x$ individual, with some arbitrary joint distribution, and let $z_x=(z_{xy})_{y\in\mathcal{Y}\cup\{ 0\} }$ be a vector of {systematic utilities} associated with the various matching options (including the option of being unmatched). We classically define the social surplus function as the expected indirected utility of each type
which is finite and convex in $z_x$. Galichon and Salanié define the entropy of choice as the convex conjugate of the former, that is
is finite when $p_x=(p_{xy})_{y}$ lies in the simplex (interpreted as the vector of choice probabilities of a type-$x$ woman across all options). An analogous pair of functions $H_y(\cdot)$ and $H_y^\star(\cdot)$ can be defined for type-$y$ men.
At the aggregate level, if $\pi_{xy}$ is the matching measure and $m_x$ (resp.\ $m_y$) the mass of type-$x$ women (resp. type-$y$ men), galichon2022cupid define the entropy of matching as the weighted sums of the entropy of choice, that is
where $\pi_{x\cdot}=(\pi_{xy})_{y}$ and $\pi_{\cdot y}=(\pi_{xy})_{x}$ denote the row and column profiles of $\pi$. Equilibrium matching is then characterized as the solution of the convex program, which characterizes the social welfare as the convex conjugate of the entropy of matching, which is given by
In the Gumbel (logit) case, $G_x$ takes the log-sum-exp form, and $G_x^\star$ reduces to an entropy term so that $\mathcal{E}(\pi)$ is (a multiple of) the Kullback–Leibler divergence; this recovers the entropic OT formulation (ref). For general random-utility models, the distribution of the shocks is encoded in the convex functionals $G_x^\star$ and $H_y^\star$, yielding a large class of entropy-like regularizers. In particular, galichon2022cupid show that for any fixed heterogeneity distribution satisfying mild convexity conditions, the mean utilities and the matching measure are identified from observed matching patterns via the convex program (ref). This firmly embeds regularized IOT as a general econometric framework, not restricted to the logit case.
The nonparametric model of choo2006marries identifies the entire surplus matrix $\Phi$ up to a normalization when matching flows are observed, but in many empirical settings a {parametric} structure on $\Phi$ is both substantively meaningful and statistically advantageous. Parametric restrictions impose economic structure, reduce dimensionality, and enable the use of rich covariates on both sides of the market. They also allow the econometrician to examine how specific characteristics contribute to the joint surplus, much as parametric discrete-choice models allow one to interpret the role of covariates in choice probabilities.
galichon2022cupid suggest to parameterize the matching surplus using
where $\lambda \in \mathbb R^K$ is the parameter to be estimated, and $(\phi_{xyk})_{xy,k}$ are basis function that parameterize the surplus. For instance, the various elements $(\phi_{xyk})_{xy,k}$ associated with various indices $k$ may distances in various socio-demographic dimensions, such as age, education, or socioeconomic status, and so on. By an application of the envelope theorem in expression (ref), the derivative of the social welfare with respect to $\lambda_k$ is seen to be
where $\pi^\lambda $ is the maximizer of (ref). galichon2022cupid define the moment matching estimator $\hat{\lambda}$ as the value of the parameter vector $\lambda$ for which the predicted moments $\sum_{xy} \pi^\lambda _{xy} \phi_{xyk}$ coincide with the observed moments $\sum_{xy} \hat{\pi} _{xy} \phi_{xyk}$, and, using equation (ref), they note that the moment matching estimator is obtained as the unique solution to the convex optimization problem
dupuy2014personality consider an extension of the previous framework in the continuous case, with the surplus being specified as
where $W_x$ and $W_y$ are vectors of observable characteristics of types $x$ and $y$ respectively, and the parameter $\lambda$ is $A$, a matrix of parameters capturing complementarities, so that the vector of entries of $A$ is the parameter vector $\lambda$. Such models allow the econometrician to estimate how education, income, health, age, personality traits, or other observables contribute to sorting patterns through the entries of $A$. Imposing a bilinear structure improves statistical efficiency, supports interpretation, and, when supplemented with low-rank or sparsity restrictions, promotes parsimony.
In the case when the random utilities are Gumbel, one can show that the dual expression of $\mathcal{E}^\star (\Phi)$ in (ref) boils down to
galichon2022cupid propose several estimation methods for parametric IOT. Taking the normalization $\sigma=1/2$, the entropically regularized optimal matching induced by $\Phi_{xy}(\lambda)=\sum_{k} \phi_{xyk} \lambda_k $ satisfies
and the moment matching estimation problem (ref) becomes
Letting $\theta = (\lambda^\top,a^\top,b^\top)^\top$, this suggests the following Poisson regression with two-way fixed effects for the observed matches $\hat\pi_{xy}$:
If the observations $\hat\pi_{xy}$ were drawn from a Poisson distribution with intensity equal to these expectations, and if one assigns a weight of 1 to matched pairs, and a weight 1/2 to unmatched elements, then the weighted log-likelihood would be equal to
Therefore $\theta$ is estimated by a weighted Poisson regression. Yet since Poisson sampling is not assumed, inference for the estimator of $\theta $ that maximizes ((ref)) is a PPML.
Without unmatched agents, the entropy regularized optimal transport solution ((ref)) can be reinterpreted as a structural gravity model of bilateral trade flows, where $\pi_{xy}$ are the flows from exporter $x$ to importer $y$, and $a_x$ and $b_y$ are exporter and importer fixed effects respectively. As is well-known in the trade literature since the paper by silva2006log, the gravity equation can be obtained by an unweighted Poisson regression, where the uniform weights stems from the fact that there are no unmatched agents to account for.
While entropic regularization smooths the matching pattern and guarantees full support, empirical researchers are often interested in {sparse} or {structured} specifications of the surplus function, especially when $\Phi_{xy}$ is parameterized by a high-dimensional matrix $A$. Sparsity can be promoted by adding an $\ell_1$-penalty to the parametric IOT likelihood: \[ \max_{A} \left\{ \sum_{x,y} \widehat{\pi}_{xy} \log \pi_{xy}(A) - \sum_{x,y} \pi_{xy}(A) - \lambda \|A\|_{1} \right\}, \] which drives many entries of $A$ exactly to zero. This yields a composite optimization problem combining a smooth (entropic-OT) component and a nonsmooth ($\ell_1$) penalty.
A computationally efficient solution to this problem is the SISTA algorithm ({Sinkhorn Iterative Soft-Thresholding Algorithm}) of carlier2023sista. SISTA alternates between two steps: (i) a {Sinkhorn/Bregman iteration} that updates the dual potentials of the entropic OT problem and computes the gradient of the smooth part with respect to $A$, and (ii) a {soft-thresholding proximal step} \[ A \leftarrow \operatorname{Soft}_{\lambda \tau} \big( A + \tau\, \nabla_A \log \pi(A) \big), \] where $\tau$ is a step size. The first step exploits the structure of entropic OT to compute $\nabla_A \pi(A)$ efficiently via the dual potentials, and the second step implements the proximal operator for the $\ell_1$ penalty. SISTA is monotone, scalable, and tailored to situations where the parameter matrix is large but the OT structure provides efficient low-dimensional updates.
Conceptually, LASSO regularization plays a role opposite to that of entropy: the entropic term {diffuses} mass and generates dense matching patterns, while the $\ell_1$ penalty {concentrates} the parameters by shrinking many coefficients to zero. Together, the two regularizers offer complementary control of smoothness, sparsity, interpretability, and generalization in parametric IOT.
The IOT problem is closely related to---and in some ways a special case of---the broader field of {metric learning}. In metric learning, one seeks to infer a distance or similarity function from observed pairwise relationships; in IOT, one seeks to infer a {transport cost} (or surplus) function from observed matching flows. Both problems aim to recover the geometry of a space from behavioral or relational data. However, the constraints imposed by transferable-utility equilibrium and the convexity induced by entropic regularization give IOT several structural advantages over classical metric learning, which is typically nonconvex and statistically more delicate.
In its conventional form, metric learning seeks a matrix $M\succeq 0$ such that the Mahalanobis distance \[ d_M(x,y) = \|x - y\|_M^2 = (x-y)^\top M (x-y) \] reflects observed similarities or dissimilarities between data points. To ensure $M\succeq 0$, one writes $M = L^\top L$ and optimizes over $L$, which makes the problem inherently nonconvex. Triplet-loss objectives, hinge losses, and neighborhood constraints, as in weinberger2009distance, hoffer2015deep, or bellet2013supervised, all involve nonconvex formulations. As a result, solutions may be sensitive to initialization, prone to local minima, and difficult to interpret. In addition, classical metric learning does not enforce consistency with marginal distributions or equilibrium conditions: similarity constraints are typically local (e.g. triplets) rather than arising from a global balance condition such as feasibility of a transport plan.
A natural bridge between metric learning and IOT emerged with the use of differentiable entropy-regularized OT as a loss function, as in cuturi2014fast and courty2016optimal. By learning the cost matrix $c_\theta$ parameterized by $\theta$ so as to match observed transport plans (or distributions, in domain adaptation), these approaches cast metric learning as the minimization of \[ \mathcal{L}(\theta) = D\big( \pi^\sigma(c_\theta), \widehat{\pi} \big), \] where $\pi^\sigma(c_\theta)$ is the solution to the entropy regularized optimal transport with cost $c_\theta$, $\hat\pi$ is the observed matching pattern, and $D$ is a divergence between probability distributions. Differentiability of $\pi^\sigma$ (thanks to entropy) gives access to stochastic-gradient methods, moving metric learning closer to convex-analytic optimal transport.
However, these methods still (i) lack equilibrium structure, (ii) do not incorporate marginal constraints as equilibrium conditions, and (iii) often rely on heuristic or local triplet-loss-like objectives.
The IOT framework imposes the full set of constraints from transferable-utility equilibrium. Instead of comparing distances for isolated triplets or neighborhoods, IOT uses the {entire joint distribution} of matches to recover the structural surplus. Metric learning becomes: \[ \text{find } \theta \text{ such that } \pi^\sigma\!\left( c_\theta \right) \approx \widehat{\pi}, \] where the matching pattern is an equilibrium object satisfying marginal constraints and dual optimality conditions. This transforms metric learning from a nonconvex factorization problem into a convex estimation problem. In particular, if the parameter $\theta$ enters linearly in the cost function (as in the bilinear surplus case studied above), then, the entropically regularized OT mapping is {smooth and globally convex}, convex duality provides explicit gradients and fixed-point characterizations, and the equilibrium constraints ensure global rather than local consistency. In economic terms, IOT learns the affinity between types in a way that respects a TU equilibrium: the recovered geometry is the one generating the observed sorting pattern under optimality and feasibility.
Classical metric learning is notoriously difficult because the Mahalanobis parametrization $M=L^\top L$ renders the optimization problem nonconvex, because similarity constraints are typically local (for example, triplets) rather than global equilibrium conditions, because distances are identifiable only up to arbitrary additive and multiplicative normalizations, and because gradient-based methods can be numerically unstable in the absence of entropy or other smoothing devices. In contrast, IOT with entropic regularization is considerably more tractable: when $\theta$ enters linearly in the cost function $c_\theta$, the objective in $\theta$ is globally concave (as in log-likelihood or KL-fitting formulations), the transport plan $\pi^\sigma(c_\theta)$ is smooth, unique, and has full support, the dual potentials yield stable gradients through convex duality, and the equilibrium conditions impose a global structure that sharpens identification. Seen from this angle, IOT can be viewed as {metric learning under convexity and equilibrium constraints}, which explains why entropically regularized IOT enjoys strong computational and statistical advantages over classical metric-learning formulations.
Viewing IOT as metric learning reveals why bilinear or low-rank surplus models are so powerful in empirical matching: they learn a latent metric or affinity space in which assortative patterns become linear (or low-dimensional). Structured sparsity (group Lasso, nuclear norm, SISTA) further enhances interpretability by identifying which characteristics of agents matter for sorting and which complementarities drive equilibrium behavior. This interpretation also emphasizes the conceptual contrast with classical discrete-choice estimation: discrete-choice models learn individual preferences, whereas IOT learns {pairwise complementarities} governing equilibrium matches. Both are metric-learning problems, but only IOT uses equilibrium OT geometry.
Outside the economics literature, the paper by StuartWolfram2020 formulates IOT as a fully fledged Bayesian inverse problem, studying well-posedness and posterior contraction for cost learning from noisy observations of optimal transport plans. More recently, AndradePeyrePoon2024 analyze $\ell_1$-regularized IOT in an entropic setting and establish sparsistency of the resulting estimators, that is, consistent recovery of the support of the true cost function. They show that, as the entropic penalty varies, IOT interpolates between classical Lasso and graphical Lasso, thereby connecting sparse IOT to sparse precision-matrix and graph estimation problems that are central in modern statistics and machine learning.
Optimal transport has become an exceptionally powerful tool in econometrics because it unifies, extends, and operationalizes several fundamental ideas in a way that is both conceptually elegant and empirically tractable. First, OT generalizes the notion of distance from points to probability distributions, yielding metrics---such as Wasserstein distances---that metrize weak convergence and underlie a wide range of applications in distributional comparisons, semiparametric inference, and distributionally robust methods. Second, OT extends the classical univariate theory of monotone rearrangements to the multivariate setting, thereby providing economically meaningful notions of multivariate quantiles, ranks, copula-based dependence, and multidimensional measures of inequality and risk, all of which are central in modern empirical work with rich heterogeneity. Third, entropically regularized OT reveals a deep algebraic connection with generalized linear models featuring two-way fixed effects, placing OT squarely within the econometric panel-data tradition and allowing researchers to draw on a mature body of identification, estimation, and inference tools. Finally, OT enjoys outstanding computational properties: exact solutions reduce to linear programming in finite dimensions, while regularized variants admit scalable Sinkhorn-type algorithms in high dimensions, enabling applications to large datasets and machine-learning environments. These combined features---a rich geometric structure, powerful generalizations of classical econometric concepts, tight connections to familiar statistical models, and remarkable computational efficiency---explain why optimal transport has rapidly become a central instrument in modern econometric analysis.