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.
55,426 characters · 13 sections · 70 citation commands
Tangential Wasserstein Projections
The concept of projections, that is, approximating a target quantity of interest by an optimally weighted combination of other quantities, is of fundamental relevance in statistics. Projections are generally defined between random variables in appropriately defined linear spaces van2000asymptotic. In modern statistics and machine learning applications, the objects of interest are often probability measures themselves. Examples range from object- and functional data marron2014overview to causal inference with individual heterogeneity athey2015machine.
We introduce a notion of projection between sets of probability measures supported on Euclidean spaces. The proposed definition is applicable between sets of general probability measures with different supports and possesses good computational and statistical properties. It also provides a unique solution to the projection problem under mild conditions and can replicate the geometric properties of the target measure, such as its shape and support. To achieve this, we work in the $2$-Wasserstein space, that is, the set of all probability measures with finite second moments equipped with the $2$-Wasserstein distance villani_optimal_2009.
Importantly, we focus on the multivariate setting, i.e. we consider the Wasserstein space over some Euclidean space $\mathbb{R}^d$, denoted by $\operatorname{\mathcal{W}}_2$, where the dimension $d$ can be high. The multivariate setting poses particular challenges from a mathematical, computational, and statistical perspective. In particular, $\operatorname{\mathcal{W}}_2$ is a positively curved metric space for $d>1$ ambrosio_gradient_2008, kloeckner2010geometric, which requires us to develop a definition of projection on positively curved metric spaces. Moreover, the $2$-Wasserstein distance between two probability measures is defined as the value function of the Monge-Kantorovich optimal transportation problem villani_topics_2003, which does not have a closed-form solution in multivariate settings. This is coupled with a well-known statistical curse of dimensionality for general measures ajtai1984optimal,dudley1969speed, fournier2015rate, talagrand1992matching, talagrand1994transportation, weed2019sharp.
These challenges have impeded the development of a method of projections between potentially high-dimensional probability measures. A focus so far has been on the univariate and low-dimensional setting. In particular, chen2021wasserstein, ghodrati2022distributions, and pegoraro2021fast introduced frameworks for distribution-on-distribution regressions in the univariate setting for object data. bigot_geodesic_2014, cazelles2017log developed principal component analyses on the space of univariate probability measures using geodesics on the Wasserstein space.
The most closely related works to ours are bonneel_wasserstein_2016 and werenski2022measure. The former develops a regression approach in barycentric coordinates with applications in computer graphics as well as color and shape transport problems. Their method requires solving a bilevel optimization problem, which is computationally costly and does not need to achieve a global solution. The latter works on a tangential structure like we do, but is based on “Karcher means” karcher2014riemannian, zemel2019frechet. This implies that their method works between absolutely continuous measures with densities that are bounded away from zero, with the target measure lying in the convex hull of the control measures.
The notion of projection we propose in this article circumvents these challenges in the multivariate setting by lifting the projection problem to the regular tangent space to $\operatorname{\mathcal{W}}_2$ at the target measure based on generalized geodesics. These tangent spaces exist when the target measure is regular, i.e., if it does not give mass to sets of lower Hausdorff dimension. For general measures, we construct a regular tangent space using barycentric projections ambrosio_gradient_2008. In contrast to the existing approaches, our method works for general probability measures, allows for the target measure to be outside the generalized geodesic convex hull of the control measures, and can be implemented by a constrained linear regression, which minimizes computational costs. In particular, Propositions (ref) and (ref) show that our method is a projection of the target onto the generalized geodesic convex hull of the control measures.
The method hence transforms the projection problem on the positively curved Wasserstein space into a linear optimization problem in the regular tangent space, which provides a unique solution to the projection problem under mild assumptions. This problem takes the form of a deformable template boissard2015distribution, yuille_deformable_1991, which connects our approach to this literature. The method can be implemented in three steps: (i) obtain the general tangent cone structure at the target measure, (ii) construct a regular tangent space if it does not exist, and (iii) perform a linear regression to carry out the projection in the tangent space.
The challenging part of the implementation is lifting the problem to the tangential structure: this requires computing the corresponding optimal transport plans between the target and each measure used in the projection. Many methods have been developed for this, see for instance benamou2000computational, jacobs2020fast, makkuva2020optimal, peyre2019computational, ruthotto2020machine and references therein. Other alternatives compute approximations of the optimal transport plans via regularized optimal transport problems peyre2019computational, such as entropy regularized optimal transport galichon2010matching, cuturi2013sinkhorn. The proposed projection approach is compatible with any such method. We provide results for the statistical consistency when estimating the measures via their empirical counterparts in practice.
To demonstrate the efficiency and utility of the proposed method, we extend the classical synthetic control estimator abadie2003economic,abadie2010synthetic to settings with observed individual heterogeneity in multivariate outcomes. This lets us perform the synthetic control method on the joint distribution of several outcomes and also estimate one set of optimal weights over all pre-intervention time periods. This complements the recently introduced method in gunsilius_distributional_2021, which is designed for univariate outcomes.
We apply this synthetic controls estimator to estimate the causal effect of Medicaid expansion on the population in individuals states. Specifically, we exploit the fact that the Affordable Care Act (ACA) allows individual states to decide whether to adopt such expansion. We use Montana as our target state, which adopted Medicaid expansion in 2016. Using the American Community Survey (ACS; ruggles2019ipums) data collected between 2010 and 2016, we evaluate such effects by estimating a counterfactual (“synthetic”) Montana, had the state not adopted Medicaid expansion. We conclude that the policy induces nontrivial, positive effects on Medicaid enrollment, earnings, and labor supply, while its effect on employment agrees with some previous estimates, but much less in magnitude compared to the other effects we estimated.
For probability measures $P_X, P_Y\in \mathscr{P}(\mathbb{R}^d)$ with supports $\operatorname{\mathcal{X}}, \operatorname{\mathcal{Y}} \subseteq \operatorname{\mathbb{R}}^d$, respectively, the $2$-Wasserstein distance $W_2 (P_X, P_Y)$ is defined as
Here, $|\cdot|$ denotes the Euclidean norm on $\mathbb{R}^d$ and \[\Gamma(P_X,P_Y)\coloneqq \left\{\gamma\in \mathscr{P}(\mathbb{R}^d\times\mathbb{R}^d): (\pi_1)_\#\gamma = P_X, \thickspace (\pi_2)_\#\gamma = P_Y\right\}\] is the set of all couplings of $P_X$ and $P_Y$. The maps $\pi_1$ and $\pi_2$ are the projections onto the first and second coordinate, respectively, and $T_\# P$ denotes the pushforward measure of $P$ via $T$, i.e. for any measurable $A\subseteq\mathcal{Y}$, $T_\# P(A) \equiv P(T^{-1}(A))$. An optimal coupling $\gamma\in \Gamma(P_X,P_Y)$ solving the optimal transport problem (ref) is an optimal transport plan. By Prokhorov's theorem, a solution always exists in our setting.
When $P_X$ is regular, i.e. when it does not give mass to sets of lower Hausdorff dimension in its support, then the optimal transport plan $\gamma$ solving (ref) is unique and takes the form $\gamma = (\operatorname{\mathrm{Id}} \times \nabla\varphi)_\# P_X$, where $\operatorname{\mathrm{Id}}$ is the identity map on $\mathbb{R}^d$ and $\nabla\varphi(x)$ is the gradient of some convex function. This result is known as Brenier's theorem brenier1991polar, mccann1997convexity, villani_topics_2003. By definition, all measures that possess a density with respect to Lebesgue measure are regular. In the following we will distinguish between regular and general measures.
The $2$-Wasserstein space $\operatorname{\mathcal{W}}_2\equiv\operatorname{\mathcal{W}}_2(\mathbb{R}^d)\equiv (\operatorname{\mathscr{P}}_2 (\operatorname{\mathbb{R}}^d), W_2)$ is the metric space defined on the set $\mathscr{P}_2(\mathbb{R}^d)$ of all probability measures with finite second moments supported on $\mathbb{R}^d$, with the $2$-Wasserstein distance as the metric. It is a complete and separable metric space ambrosio_gradient_2008 and also possesses a geometric structure that we exploit. In particular, it is a geodesically complete space in the sense that between any two measures $P,P'\in\operatorname{\mathcal{W}}_2$, one can define a geodesic $P_t: [0,1] \to\operatorname{\mathcal{W}}_2$ via the interpolation ambrosio_gradient_2008, mccann1997convexity $P_t\coloneqq (\pi_t)_\# \gamma$, where $\gamma$ is an optimal transport plan and $\pi_t: \mathbb{R}^d\times\mathbb{R}^d\to\mathbb{R}^d$ is defined through $\pi_t(x,y)\coloneqq (1-t)x + ty$. Using this, it can be shown that $\operatorname{\mathcal{W}}_2$ is a positively curved metric space $d>1$ ambrosio_gradient_2008 and flat for $d=1$ kloeckner2010geometric, where curvature is defined in the sense of Aleksandrov aleksandrov1951theorem. This difference in the curvature properties is the main reason for why the multivariate setting requires different approaches compared to the established results for measures on the real line.
In linear spaces, the idea of projecting an element on a set of other elements is a linear or convex combination of these elements in the set. The analogues of averages in vector spaces are Fr\'echet means or barycenters in metric spaces.
For Wasserstein spaces, this concept has been introduced in our setting in agueh_barycenters_2011 and in a more abstract sense in carlier2010matching. For any collection of probability measures $\set{P_j}_{1 \leqslant j \leqslant J}\subseteq\operatorname{\mathcal{W}}_2$, their weighted barycenter $\bar{P}(\lambda)$ for given weights $\lambda\equiv (\lambda_1,\ldots,\lambda_J)\in\Delta^J$ is a solution of the following minimization problem:
The weights $\lambda$ are defined to lie in the $J$-dimensional probability simplex $\Delta^{J}$ of nonnegative vectors in $\mathbb{R}^J$ that sum to unity. Prokhorov's theorem implies that a solution to (ref) exists agueh_barycenters_2011.
Throughout, we focus on the case where $\lambda\in \Delta^{J}$, because this provides a natural notion of interpolation between the given measures $P_j$. We could also construct a linear projection that corresponds to a “geodesic extrapolation” by relaxing the requirement that all weights $\lambda_j$ need to be non-negative. In the classical setting of random variables mapping to some Euclidean space and not random measures, this relaxation is the natural analogue to a linear regression, see abadie2015comparative for a discussion. We focus on the interpolation setting because it is in line with the notion of a projection onto the convex hull spanned by other elements. All of our results can be extended to the extrapolation setting in principle.
One way to extend the notion of projection between a target probability measure $P_0$ and a set of control measures $\{P_j\}_{j=1,\ldots,J}$ in $\operatorname{\mathcal{W}}_2$ would hence consist in finding the optimal weights $\lambda^*\in\Delta^{J}$ for which an induced barycenter $\bar{P}(\lambda^*)$ solving (ref) is as close as possible to $P_0$. This would lead to the following bi-level optimization problem, assuming that the barycenter $\bar{P}(\lambda)$ is unique for given $\lambda$:
A version of this approach is used in bonneel_wasserstein_2016 to define a notion of regression between probability measures on rectangular supports in low dimensions. The challenges with this approach are mathematical and computational. Importantly, the optimal weights $\lambda^*$ need not be unique. This is not an issue for the applications considered in bonneel_wasserstein_2016, like color transport; however, it is important in statistical settings when the weights convey information used in further procedures, like causal inference via synthetic controls, where the optimal weights are used to introduced a counterfactual outcome of a treated unit had it not been treated abadie2003economic, abadie2010synthetic, abadie2021using. Moreover, the bi-level optimization structure makes solving the problem prohibitively costly for higher-dimensional distributions. bonneel_wasserstein_2016 introduce a gradient descent approach based on an entropy-regularized analogue of $W_2$ cuturi2013sinkhorn, peyre2019computational that can be implemented in settings with low-dimensional empirical measures of rectangular support. However, in higher dimensions and with general measures, an efficient implementation of (ref) that can produce unique weights is absent.
We circumvent these difficulties by exploiting a tangential structure that can be defined on $\operatorname{\mathcal{W}}_2$ ambrosio_gradient_2008,otto2001geometry. In particular, it allows us to entirely circumvent solving a bi-level optimization problem as the one in (ref).
The tangential structure relies on the fact that geodesics $P_t$ in $\operatorname{\mathcal{W}}_2$ are linear in the transport plans $(\pi_t)_\#\gamma$. This implies a geometric tangent cone structure at each measure $P\in\operatorname{\mathcal{W}}$ that can be defined as the closure in $\operatorname{\mathscr{P}}_2(\mathbb{R}^d)$ of the set \[\mathcal{G}(P)\coloneqq \left\{\gamma\in \operatorname{\mathscr{P}}_2(\mathbb{R}^d\times\mathbb{R}^d): (\pi_1)_\#\gamma = P,\thickspace\medspace (\pi_1,\pi_1+\varepsilon\pi_2)_\#\gamma\thickspace\text{is optimal for some $\varepsilon>0$}\right\}\] with respect to the local distance
where $\gamma_{12}$ and $\gamma_{13}$ are couplings between $P$ and some other measures $P_2$ and $P_3$, respectively, and $\Gamma_1(\gamma_{12},\gamma_{13})$ is the set of all $3$-couplings $\gamma_{123}$ such that the projection of $\gamma_{123}$ onto the first two elements is $\gamma_{12}$ and the projection onto the first and third element is $\gamma_{13}$ ambrosio_gradient_2008. We can then define the exponential map at $P$ with respect to some tangent element $\gamma\in\mathcal{G}(P)$ by
This tangent cone can be constructed at every $P\in\operatorname{\mathcal{W}}$, irrespective of its support.
In the case where $P$ is absolutely continuous with respect to Lebesgue measure, the definition simplifies. The corresponding optimal transport plan between $P$ and any other measure measure $Q\in\operatorname{\mathcal{W}}_2$ is then supported on the graph of the gradient of a convex function $\nabla\varphi$ by Brenier's theorem. This allows the introduction of a regular tangent cone at $P$, defined via ambrosio_gradient_2008
where $\overline{A}^{L^2(P)}$ defines the closure of the set $A$ with respect to the distance induced by the $L^2$-norm on $P$. Interestingly, this tangent cone is actually a linear tangent space ambrosio_gradient_2008. In the following, we therefore call $\mathcal{T}_P\operatorname{\mathcal{W}}_2$ the regular tangent space (ref).
We want to define a natural analogue to the barycenter projection problem (ref) in the regular tangent space of the Wasserstein space. A starting point for this is to consider a characterization of the barycenter $\bar{P}(\lambda)$ for fixed weights of a set $\{P_j\}_{j\in\llbracket J\rrbracket}$ in regular tangent spaces. agueh_barycenters_2011 show that if at least one of the measures is absolutely continuous with respect to Lebesgue measure, then $\bar{P}(\lambda)$ can be characterized via
where $\{\varphi_j\}_{j \in \llbracket J \rrbracket}$ are the optimal transport maps from the barycenter to the respective measure $P_j$, i.e. $(\varphi_j)_\# \bar{P}(\lambda) = P_j$. Each term of the summand in (ref) is an element in $\operatorname{\mathcal{T}}_{\bar{P}(\lambda)} \operatorname{\mathcal{W}}_2 (\mathbb{R}^d)$ by construction.
More generally, the condition (ref) is a sufficient condition for $\bar{P}(\lambda)$ to be a “Karcher mean” karcher2014riemannian in $\operatorname{\mathcal{W}}_2$ zemel2019frechet. In fact, a “Karcher mean” of a set of measures $\{P_j\}_{j\in\llbracket J\rrbracket}$ is defined as the gradient of the Fr\'echet functional in $\operatorname{\mathcal{W}}_2$ and is characterized through (ref) holding $\bar{P}(\lambda)$-almost everywhere. (ref) is a stronger condition because it is assumed to hold at every point in the support of $\bar{P}(\lambda)$, not just almost every point. alvarez2016fixed use this characterization to introduce a fixed-point approach to compute Wasserstein barycenters, and werenski2022measure use this structure to introduce a replication approach for absolutely continuous measures whose densities are bounded away from zero and whose target measure lies inside the convex hull of the control measures. Related is the recent definition of weak barycenters in cazelles2021novel, where the authors replace the optimal transport maps from the classical optimal transport problem by the weak optimal transport problem introduced in gozlan2017kantorovich. Heuristically, this characterization is that of a deformable template. A measure $P$ is a deformable template if there exists a set of deformations $\set{\psi_{j}}_{j=1,\dots,J}$ such that ${\psi_{j}}_{\#} P = P_j$, in a way that their weighted average is “as close to the identity” as possible. In our setting $\psi_j\equiv\nabla\varphi_j-\operatorname{\mathrm{Id}}$ anderes_discrete_2015,boissard2015distribution,yuille_deformable_1991.
In our setting of interest, we are given a target measure $P_0$ that we want to replicate given a set $\{P_j\}_{j\in\llbracket J\rrbracket}$. The key idea for this is to adapt the characterization (ref) with respect to the target measure $P_0$. That is, we aim to find the optimal weights $\lambda^*\in\Delta^{J}$ that satisfy
where $\nabla\varphi_j$ are the optimal transport maps between the target $P_0$ and the control measures $P_j$, $j\in\llbracket J\rrbracket$.
We now show that this approach is in fact a projection of $P_0$ onto the generalized geodesic convex hull of the control measures $P_j$ with respect to $P_0$ as illustrated in Figure (ref). To define this notion of convex hull, we extend the definition of generalized geodesics ambrosio_gradient_2008, and in particular the definition of $W_P$ to the multimarginal setting, by defining, for given couplings $\gamma_{0j}\in \Gamma(P_0,P_j)$, $j\in\llbracket J\rrbracket$
where $\Gamma_1(\gamma_{01},\ldots,\gamma_{0J})\subseteq\Gamma(P_0,P_1,\ldots,P_J)$ is the set of all $(J+1)$-couplings $\bm{\gamma}$ such that the projection of $\bm{\gamma}$ onto the first- and $j$-th element is $\gamma_{0j}$. Note that this definition is similar to the multimarginal definition of the $2$-Wasserstein barycenter agueh_barycenters_2011, gangbo1998optimal, but “centered” at $P_0$. Based on this, we define the generalized geodesic convex hull of measures $\{P_j\}_{j\in\llbracket J\rrbracket}$ with respect to the measure $P_0$ as
Based on these definitions we can show that our approach is a projection of the target $P_0$ onto $\operatorname{\mathfrak{Co}}_{P_0}\left(\{P_j\}_{j=1}^J\right)$.
Proposition (ref) holds for a regular target $P_0$. In many practical settings, however, the target outcome is not a regular measure, as in our application in Section (ref). In such settings, the corresponding optimal transport problems (ref) between the target and the respective control measures $P_j$ is only achieved via optimal transport plans $\gamma_{0j}$, not maps $\nabla\varphi_j$. In contrast to the regular setting, these transport maps also do not need to be unique.
Still, the fundamental idea of the tangential projection can be extended to the more general setting, using transport plans instead of maps. The implementation for regular targets (ref) is a special case of (ref) where all plans $\gamma_{0j}$ are achieved via transport maps $\nabla\varphi_j$. This implies that the optimal weights $\lambda^*\in\Delta^{J}$ could in principle be obtained by
where $\gamma_{0j}$ are optimal transport plans between the target $P_0$ and the respective control measure $P_j$. By definition, a solution to (ref) will provide a projection of the target $P_0$ onto $\operatorname{\mathfrak{Co}}_{P_0}\left(\{P_j\}_{j=1}^J\right)$. However, solving (ref) is computationally prohibitive in practice for two reasons. First, it is again a bilevel problem, similar to the direct approach (ref). Second, $W_{P_0;\lambda}^2$ requires computing a joint coupling over $J+1$ marginal distributions, which is computationally infeasible in practice for a reasonably large $J$.
We therefore rely on barycentric projections to reduce the complex general setting to the regular tangent space and subsequently apply the projection (ref). This is computationally inexpensive, as it amounts to computing
where \[b_{\gamma_{0j}}(x_1)\coloneqq \int_{\mathbb{R}^d} x_2\dif \gamma_{0j,x_1}(x_2)\] are the barycentric projections of optimal transport plans $\gamma_{0j}$ between $P_0$ and $P_j$. Here, $\gamma_{x_1}$ denotes the disintegration of the optimal transport plan $\gamma$ with respect to $P_0$.
This approach is a natural extension of the regular setting to general probability measures for two reasons. First, if the optimal transport plans $\gamma_{0j}$ are actually induced by some optimal transport map $\nabla\gamma_j$, then $b_{\gamma_{0j}}$ reduces to this optimal transport map; in this case the general tangent cone $\mathcal{G}(P_0)$ reduces to the regular tangent cone $\mathcal{T}_{P_0}\operatorname{\mathcal{W}}_2$ ambrosio_gradient_2008. Second, by the definition of $b_{\gamma}$ and disintegrations in conjunction with Jensen's inequality it holds for all $\lambda\in\Delta^{J}$ that
This implies that for general $P_0$ we can also define a convex hull based on barycentric projections, which is of the form
Furthermore, the contraction property (ref) implies that $\operatorname{\mathfrak{Co}}_{P_0}\subseteq\widetilde{\operatorname{\mathfrak{Co}}}_{P_0}$, with equality when all transport plans are achieved via maps $\nabla\varphi_j$. Using these definitions, the following defines our notion of projection for general $P_0$ and shows that it projects onto $\widetilde{\operatorname{\mathfrak{Co}}}_{P_0}$.
Note that in contrast to the regular case in Proposition (ref), the optimal plans $\gamma_{0j}$ transporting $P_0$ to $P_j$ need not be unique, i.e., the measures $P_j$ might lie outside the cut locus of $P_0$. However, the projection for fixed $\gamma_{0j}$ is unique.
The proposed method of projections is hence a well-defined notion of a geodesic metric projection: in the case of a regular measure, our approach is a metric projection of the target measure onto the generalized geodesic convex hull made up of the control measures. In the case where the target measure is not regular, we project onto a slight extension of the generalized geodesic convex hull, which we construct by a barycentric projection. The actual projections (ref) and (ref) are simple regression problems, which are easy to compute in practice once the tangent structure has been constructed.
We now provide statistical consistency results for our method when the corresponding measures $\{P_j\}_{j\in \llbracket J\rrbracket}$ are estimated from data. We consider the case where the measures $P_j$ are replaced by their empirical counterparts \[\operatorname{\mathds{P}}_{N_j}(A)\coloneqq N_j^{-1}\sum_{n=1}^{N_j} \delta_{X_{n}}(A)\] for every measurable set $A$ in the Borel $\sigma$-algebra on $\mathbb{R}^d$, where $\delta_x(A)$ is the Dirac measure and $\left(X_{1j},\ldots, X_{N_j,j}\right)$ is an independent and identically distributed set of random variables whose distribution is $P_j$. We explicitly allow for different sample sizes $\bigcup_{j=0}^J N_j = N$ for the different measures. To save on notation we write $\widehat{\varphi}_{N_j}\equiv \widehat{\varphi}_j$, $\widehat{b}_{0j}\equiv \widehat{b}_{\gamma_{0j}, N_j}$ and $\widehat{\gamma}_{0j}\equiv \widehat{\gamma}_{N_j,N_0}$ in the following.
This consistency result directly implies consistency of the optimal weights in case the optimal transport problems between $\operatorname{\mathds{P}}_{N_0}$ and each $\operatorname{\mathds{P}}_{N_j}$ are achieved by optimal transport maps $\nabla\widehat{\varphi}_{N_j}$. Based on this we also have a consistency result for the empirical counterparts $\widetilde{\operatorname{\mathds{P}}}_{\pi,N}$ of the optimal projection $\widetilde{P}_\pi$.
Proposition (ref) and Corollary (ref) hold in all generality and without any assumptions on the corresponding measures $P_j$, except that they possess finite second moments. To get stronger results, for instance parametric rates of convergences or even asymptotic Gaussianity, one needs to make strong regularity assumptions on the measures $P_j$. Without these, the rate of convergence of optimal transport maps in terms of expected square loss is as slow as $n^{-2/d}$ hutter2021minimax. Under such additional regularity conditions, the results for the asymptotic properties are standard, because the proposed method reduces to a classical semiparametric estimation problem, as the weights $\lambda_j$ are finite-dimensional. In particular, the setting is that of a MINPIN estimator, as defined in andrews1994asymptotics.
The fact that the actual tangential projection is a simple constrained linear regression for the weights $\lambda$ implies that the regularity of $b_j$ is what drives the statistical properties of the weights $\lambda$. For instance, if the barycentric projections $b_j$ are regular in the sense that they lie in a Donsker class for their given dimension $d$ wellner2013weak and converge to their population counterpart at at least the rate of $n^{1/4}$, then the optimal weights will converge at the parametric rate, which follows directly from the arguments in andrews1989asymptotics, andrews1994asymptotics. In the regular setting, i.e. when the target $P_0$ is a regular measure so that $b_j=\nabla\varphi_j$ for some convex functions $\varphi_j$, such regularity conditions can be derived from classical regularity theory caffarelli1990interior, caffarelli1992regularity, de2013w, and have been used in deriving rates of convergences of optimal transport maps by deb2021rates, forrow2019statistical, gunsilius2021convergence, hutter2021minimax, manole2021plugin, weed2019sharp. The same holds for estimators that use barycentric projections after solving an entropy-regularized analogue to the optimal transport problem seguy2018large, pooladian_entropic_2021. Without these regularity assumptions, the curse of dimensionality outlined in these statistical results implies that the rate of convergence for the optimal weights will in general be slower than the parametric rate, especially in higher dimensions.
In this section, we provide some simulations and an application to the synthetic controls estimator to demonstrate the computational properties of the method. We use the POT package flamary2021pot to obtain the optimal transport plans, and CVXPY diamond2016cvxpy, agrawal2018rewriting to compute (ref). Additional details of our applications are contained in Appendix (ref).
We apply our estimator to Gaussian distributions in dimension $d = 10$. We draw from the following Gaussians: \[\mathbf{X}_j \sim \operatorname{\mathcal{N}} \left( \mu_j, \Sigma \right), \quad j=0,1,2, 3 ~, \] where $\mu_0 = [10, 10, \dots, 10]$, $\mu_1 = [50, 50, \dots, 50]$, $\mu_2 = [200, 200, \dots, 200]$, $\mu_3 = [-50, -50, \dots, -50]$ and $\Sigma = \operatorname{\mathrm{Id}}_{10}+ 0.8 \operatorname{\mathrm{Id}}^{-}_{10}$, with $\operatorname{\mathrm{Id}}_{10}^{-}$ the $10\times 10$ matrix with zeros on the main diagonal and ones on all off-diagonal terms. We take $\mathbf{X}_0$ as target, and $\mathbf{X}_1$, $\mathbf{X}_2$, $\mathbf{X}_3$ as controls. The optimal weights are $\lambda^* =[0.3643, 0.0943, 0.5414]$, meaning $\mathbf X_1$ and $\mathbf X_3$ receive substantial weights, while $\mathbf X_2$ only receives a small amount. This weight distribution can be understood by looking at the differences between $\mu_0$ and $\mu_1$, $\mu_2$ and $\mu_3$, separately---which, in this case, is based on the distance between the means of these distributions, since we work with Gaussian distributions of the same variance. The mean of $\mathbf X_2$ is significantly further away from the mean of the target than the other two means, so it is to be expected that $\mathbf X_2$ receives little weight. Table (ref) suggests the projection is close to the target distribution when only considering the mean, despite only having three control units.
To demonstrate the property that our method provides sparse weights in general, we provide an application on replicating a target image of an object using images of the same object taken from different angles. We use the Lego Bricks dataset available from Kaggle, which contains approximately 12,700 images of 16 different Lego bricks. These images are in the RGBA format, despite being grayscale, which provides a good resolution. For our application, the target image is contained in Figure (ref)(A), and the control images are in Figure (ref).
The result, as shown in Figure (ref), indicates the optimally-weighted projection replicates the target image well, even in a setting with only a few controls. We note two interesting findings. One, only the control images in the first row of Figure (ref) received nontrivial weights; all others received essentially zero weights. This suggests that our method use most information from controls which look sufficiently like the target and demonstrates that our method provides sparse weights. Two, even with few controls, the optimally weighted projection approximates the target well in this application.
In this section, we apply the method to extend the classical notion of synthetic controls abadie2003economic, abadie2010synthetic, abadie2021using and its generalization to univariate distributions gunsilius_distributional_2021 to general multivariate outcome distributions. Moreover, this generalization allows to estimate one set of optimal weights over all time periods jointly, while existing methods need to apply the method in every time period before averaging. This removes sparsity in the optimal weights estimated. For more details, we refer to abadie2021using and gunsilius_distributional_2021.
An application to a causal inference setting is studying the effect of health insurance coverage following state-level Medicaid expansion in the United States. A provision within the ACA allows states to decide whether to expand Medicaid for low-income households. Some states decide to adopt such expansions early on, while others did not (and still have not done so). We investigate the economic and behavioral effects of Medicaid expansion. Specifically, we consider some first-order effects (i.e. the extensive margin of Medicaid enrollment post-expansion and disemployment effects) and second-order effects (i.e. income and labor supply effects) of expanded Medicaid access.
The observational data we use is the ACS. From it, we collect the following variables:
We measure labor hours supplied and labor income in logs instead of levels. We consider Montana as the treated unit for this application. Montana adopted such expansion in 2016. For control units, we use the twelve states for which such expansion has never occurred. As of 2022, the twelve states are: Alabama, Florida, Georgia, Kansas, Mississippi, North Carolina, South Carolina, South Dakota, Tennessee, Texas, Wisconsin, Wyoming. We use data from 2010 to 2016 to estimate optimal weights for the control states in order to create the “synthetic Montana”, i.e. Montana had it not adopted Medicaid expansion. Details of sample selection and estimating “synthetic Montana” are described in Appendix (ref). As we show in the Appendix, the combination of control states replicating Montana is close to actual Montana in the pre-intervention period, implying that the observed differences in post-intervention periods can be attributed to the causal effect of the Medicaid expansion.
Consistent with findings in courtemanche2017early, mazurenko2018effects, we find significant first-order effects of Medicaid expansion, which are summarized in Figure (ref). We note that the “synthetic Montana” has much lower proportion of individuals insured under Medicaid, suggesting that expanding Medicaid eligibility directly affects the extensive margin of Medicaid enrollment. The disemployment effect is much less pronounced in comparison to the enrollment effect we estimated, but nonetheless positive and nontrivial; this is consistent with the findings in, e.g., peng2020effects, but inconsistent with those in, e.g., gooptu2016medicaid.
We also find nontrivial, positive second-order effects, summarized in Figure (ref). In both earnings and labor hours supplied, we see the “synthetic Montana” has lower averages, and narrower supports compared to the observed distributions. This suggests Medicaid expansion improves earnings, and widens the intensive margin of labor supply.
We further note two crucial findings from our application. One, the optimal weights estimated here is, again, sparse. Based on our results, summarized in Table (ref), we observed only five control states---Alabama, Mississippi, South Carolina, South Dakota, and Wisconsin---with nonzero weights. South Dakota alone constitutes over half of the total weight, suggesting it best approximates Montana compared to all other control states. Two, we estimated these weights using data from all years covered in the pre-intervention period, which provides one sparse set of weights over all time periods. This stands in contrast to the standard synthetic controls method abadie2003economic,abadie2010synthetic, where the optimal weights are obtained from taking some weighted average of weights estimated in every time unit during the pre-intervention period; this averaging over weights in each time period generates a non-sparse weight.
We have developed a projection method between sets of probability measures supported on $\mathbb{R}^d$ based on the tangential structure of the 2-Wasserstein space. Our method seeks to best approximate some target distribution that is potentially multivariate, using some chosen set of control distributions. We provide an implementation which gives unique, interpretable weights in a setting of regular probability measures. For general probability measures, we construct our projection by first creating a regular tangent space through applying barycentric projection to optimal transport plans. Our application to evaluating the first- and second-order effects of Medicaid expansion in Montana via an extension of the synthetic controls estimator demonstrates the method's efficiency and the necessity to have a method that is applicable for general proabbility measures. The approach still works without restricting optimal weights to be in the unit simplex, which would allow for extrapolation beyond the convex hull of the control units, providing a notion of tangential regression. It can also be extended to a continuum of measures, using established consistency results of barycenters le2017existence.