# A Convexified Matching Approach to Imputation and Individualized Inference

YoonHaeng Hur  University of Chicago  Tengyuan Liang Liang acknowledges the
generous support from the NSF Career Grant (DMS-2042473), and the William
Ladany Faculty Fellowship from the University of Chicago Booth School of
Business. University of Chicago

###### Abstract

We introduce a new convexified matching method for missing value imputation
and individualized inference inspired by computational optimal transport. Our
method integrates favorable features from mainstream imputation approaches:
optimal matching, regression imputation, and synthetic control. We impute
counterfactual outcomes based on convex combinations of observed outcomes,
defined based on an optimal coupling between the treated and control data
sets. The optimal coupling problem is considered a convex relaxation to the
combinatorial optimal matching problem. We estimate granular-level individual
treatment effects while maintaining a desirable aggregate-level summary by
properly constraining the coupling. We construct transparent, individual
confidence intervals for the estimated counterfactual outcomes. We devise fast
iterative entropic-regularized algorithms to solve the optimal coupling
problem that scales favorably when the number of units to match is large.
Entropic regularization plays a crucial role in both inference and
computation; it helps control the width of the individual confidence intervals
and design fast optimization algorithms.

Keywords— Missing value imputation, optimal coupling, matching, propensity
score, individualization, entropic regularization.

##  1 Introduction

One central topic in applied econometric research is to assess the effects of
policy interventions reliably. On the one hand, to analyze nonexperimental
data, advanced econometric estimates are devised to provide granular
counterfactual answers by leveraging specific structural models. In an
influential paper in 1986 [LaL86], LaLonde questioned whether such
sophisticated econometric estimates are credible. He compared them with the
experimental benchmarks on some coarse summary—for example, the average
treatment effect (ATE)—and concluded unfavorably. On the other hand, when
experimental data, as in randomized controlled trials, are available,
elementary statistical estimates can determine the ATE and eliminate
confounding explanations. The ATE estimate, however, does not answer whether
the treatment works for an individual. In modern applications such as
personalized medicine and online marketing, the treatment effects vary across
individuals; the treatment might be beneficial for some individuals but
ineffective for others [LM16]. One potentially costly approach is to conduct
individualized experiments in different time windows, known as N-of-1 trials
[Hil55, LR23]. Therefore, it is desirable to develop nonexperimental methods
that conform to the coarse estimates—for example, those for ATE or the average
treatment effect on the treated (ATT)—at the aggregate level, while delivering
more granular, individualized inference.

Individualized inference is about imputing missing counterfactual outcomes—the
outcomes that would have been observed had the subjects received an opposite
treatment—and their uncertainty quantifications. Three popular approaches are
available in the literature: matching, regression imputation, and synthetic
control. This paper proposes a new convexified matching method by integrating
these three approaches. The central quantity we play is a coupling matrix
between the treated and control data clouds, which resembles matching and is
used to synthesize the counterfactual outcomes as in classic nonparametric
regression.

Matching [Rub06, Ros20] is a widely adopted way to impute counterfactual
outcomes. For each subject, the counterfactual outcome is estimated by
identifying units in the opposite treatment group similar to the subject in
covariates and then averaging their outcomes. Due to its simplicity, nearest
neighbor matching, which finds subjects closest under a suitable distance
between the covariates, is favored in practice. Exact or approximate matching
with high-dimensional covariates can be difficult, and thus propensity score
matching was proposed as a remedy [RR83]. There are more sophisticated
matching methods to improve balance or forbid certain matches, often referred
to as the optimal matching [Ros20, Zub+23]; they require solving combinatorial
optimization problems by network optimization techniques or mixed integer
programming, which can be computationally expensive. However, even analyzing
the aggregate ATE estimates resulting from combinatorial optimization problems
is highly nontrivial, and these estimates may differ from simple ATE estimates
based on propensity score reweighting [Ros87, HI01, HIR03]. As in matching, we
solve for a coupling matrix—a convex relaxation to the combinatorial
constraints—without access to outcomes, separating the design and analysis
phases in an observational study. For a large-scale matching problem with $N$
units, we devise entropic-regularized algorithms from optimal transport
[Vil03, PC19] to solve convexified matching, calling a matrix scaling
subroutine with a number of times nearly independent of $N$; see Section 4. In
contrast, in a conventional matching problem, solvability within polynomial
time of $N$ may not be feasible.

Regression imputation is another approach to estimating counterfactual
outcomes. The idea is simple: under the ignorability or unconfoundedness
assumption [RR83], the underlying functions that map the covariates to the
responses under the treatment and the control are identifiable and estimable
based on observational data. It translates a causal inference problem into a
regression problem. Parametric or nonparametric regression techniques estimate
the conditional response function under the treatment and the control, given
the covariates, and in turn find the conditional average treatment effect
(CATE), see [Hah98, HIT97, HV05, MCH11]. At the aggregate level,
semiparametric efficient estimation of ATE leveraging regression imputation
techniques is studied in [Hah98]. Admittedly, if the quantity of interest is
at the aggregate-level—a one-dimensional parameter such as ATE or ATT—the
specific nonparametric regression technique matters less, as long as it
estimates the function reasonably well. There has been a fruitful line of
literature on the use of flexible machine learning methods in estimating CATE
[BCH14, AI16, Che+17, Far15, FLM21, WA18], in place of classic Nadaraya-Watson
(NW) kernel nonparametric regression. Albeit naive, let us use the NW
estimator as an example to illustrate a shared feature by regression
imputation and matching: to compute a counterfactual outcome for a treated
unit, NW finds local neighbors in the control group and uses convex weights
based on the nonparametric kernel to synthesize the outcome. We will adopt
this aggregation by convex weights feature into our convexified matching,
whereas our coupling weights are optimized globally. Our coupling weights are
determined to optimize imputation quality and individual uncertainty
quantification; see Section 3.

Hybrid methods were proposed to combine matching/propensity score weighting
[Ros87, HI01, HIR03, LMZ18] with regression imputation techniques to
efficiently estimate ATE, notably the doubly robust estimators [VR03, Cat10,
Far15, FLM21, Che+22]. The estimator compares favorably as it is valid when
either the propensity score function or the regression function is consistent,
robust to bias due to misspecification. Similar in spirit, bias correction
using regression to improve the nearest neighbor matching was studied in
[AI11, AI06].

From a different vein, synthetic control [AG03, ADH10, Aba21] provides a fresh
look at counterfactual imputation: a convex combination of control units may
synthesize a subject under treatment better than any one of the control unit.
Mathematically, a convex combination of data may simultaneously reduce the
approximation error or bias, and, at the same time, reduce the variance driven
by the idiosyncratic errors. This convex relaxation idea also eases the
computation for the optimization problem in matching. We adopt this convex
combination idea to approximate the overall covariate information under
treatment by convex combinations of the covariate information under control,
defined based on a data set-to-data set coupling; see Section 2.

We develop a convexified matching method that solves an optimal coupling,
which in turn defines convex combinations for missing value imputation. Unlike
combinatorial optimization problems, where the optimization variable must be
an extreme point of a certain polytope to represent an assignment, we optimize
over the whole polytope as the variable represents the weights of the convex
combination. Moreover, we add an entropic regularization to the optimization,
where the strength of regularization controls the bias and variance tradeoff.
The resulting formulation is a smooth convex optimization problem confined to
the polytope that can be solved efficiently. Lastly, by properly specifying
the constraints, we can pair with propensity score weighting [Ros87, HI01,
HIR03] estimators. We can guarantee that the aggregate summary of individual
treatment effects coincides with the desired ATE or ATT estimates.

We introduce an inference procedure to construct individual confidence
intervals for the estimated counterfactual outcomes based on the convexified
matching method. Providing a credible confidence interval around the
individual treatment effect can help decide whether to adopt the treatment for
each individual. The width of the confidence interval is determined by the
approximation error in covariate balancing and the entropy of the convex
weights, which comprise the objective function of the optimization. Therefore,
how we formulate the convexified matching is directly targeted at
individualized inference. Entropic regularization plays a crucial role in
inference, namely, controlling the width of the individual confidence
intervals, and in computation, namely, designing fast iterative algorithms to
solve the optimization.

In summary, we integrate important features from matching, regression
imputation, and synthetic control. We impute counterfactual outcomes by convex
combinations defined based on an optimal coupling. The coupling is a convex
relaxation to optimal matching and can be solved efficiently using iterative
matrix scaling subroutines called the Sinkhorn algorithm [Sin67, Cut13]. We
estimate granular individual treatment effects while maintaining a desirable
aggregate-level summary by properly constraining the coupling. We construct
transparent, individual confidence intervals for the estimated counterfactual
outcomes, where the optimization objective controls the width of the
confidence intervals.

#### Notation

We denote by $\|\cdot\|_{2}$ and $\langle\cdot,\cdot\rangle$ the Euclidean
norm and inner product, respectively. For a square matrix $A$, let
$\mathrm{tr}(A)$ denote its trace. For matrices $A,B$ of the same dimension,
$\langle A,B\rangle:=\mathrm{tr}(A^{\top}B)$ is the Frobenius inner product.
For any integer $n\in\mathbb{N}$, let $1_{n}=(1,\ldots,1)\in\mathbb{R}^{n}$
denote the vector whose entries are all $1$. Let
$\Delta_{n}:=\\{a\in\mathbb{R}_{+}^{n}:\sum_{i=1}^{n}a_{i}=1\\}$ denote the
simplex of probability vectors and let
$\Delta_{n}^{+}:=\\{a\in\Delta_{n}:a_{i}>0~{}~{}\forall i\\}$ denote the
interior of $\Delta_{n}$.

##  2 Imputation Method: Synthetic Coupling

We first introduce Kernel Synthetic Coupling (KSC), a convexified matching
method for missing value imputation. Later in Section 3, we build individual
confidence intervals around the imputed values. To streamline the exposition,
we cast the KSC method following the notations in the potential outcome
framework [Ney23, Rub74] for a binary treatment—a leading example for missing
value imputation.

We consider $N$ units indexed by $1,\ldots,N$, where each unit is associated
with two potential outcome variables $Y_{i}(1),Y_{i}(0)\in\mathbb{R}$ under
treatment and control, respectively, and a covariate vector
$x_{i}\in\mathbb{R}^{d}$. We can observe only one of $Y_{i}(1),Y_{i}(0)$
depending on the treatment assignment $Z_{i}=1$ or $Z_{i}=0$, denoting that
unit $i$ is treated or not. Accordingly, the observed outcome $Y_{i}$ of unit
$i$ is given as $Y_{i}=Z_{i}Y_{i}(1)+(1-Z_{i})Y_{i}(0)$. We let
$\mathcal{T}=\\{i\in[N]:Z_{i}=1\\}$ and $\mathcal{C}=\\{i\in[N]:Z_{i}=0\\}$
denote the sets of indices of the treated units and control units,
respectively, so that $\mathcal{T}\cup\mathcal{C}=\\{1,\ldots,N\\}$. Also, let
$N_{t}:=|\mathcal{T}|$ and $N_{c}:=|\mathcal{C}|$ denote the numbers of the
treated and control units. The individual treatment effect of unit $i$ is
defined as $\tau_{i}=Y_{i}(1)-Y_{i}(0)$, which is unobservable as only one of
$Y_{i}(1),Y_{i}(0)$ can be observed.

###  2.1 A Simple Formulation

Without loss of generality, we focus on imputing the missing potential
outcomes of the treated units; counterfactual outcomes of the control units
can be similarly obtained. We estimate $Y_{j}(0)$ by some $\widehat{Y}_{j}(0)$
for each treated unit $j$. This, in turn, allows us to estimate the individual
treatment effect $\tau_{j}$ by
$Y_{j}-\widehat{Y}_{j}(0)=Y_{j}(1)-\widehat{Y}_{j}(0)$. The proposed
convexified matching combines two ideas. First, we want to find a match
between the data cloud in the treated group and that in the control using the
covariate information. Second, as in the synthetic control method, the
matching quality is determined by how well it approximates the covariate of
each treated unit by a convex combination of the covariates of control units.
As a result, this procedure involves finding a matrix, which we call a
coupling, indexed by a pair $(i,j)$ for $i\in\mathcal{C}$ and
$j\in\mathcal{T}$ denoting the weight assigned to the covariate of control
unit $i$ to approximate the covariate of treated unit $j$. We solve the
following optimization to obtain the optimal coupling:

$\displaystyle\min_{\pi\in\mathbb{R}_{+}^{N_{c}\times N_{t}}}\quad$
$\displaystyle\frac{1}{2N_{t}}\sum_{j\in\mathcal{T}}\|x_{j}-\sum_{i\in\mathcal{C}}\tfrac{\pi_{ij}}{1/N_{t}}x_{i}\|_{2}^{2}+\lambda\sum_{i\in\mathcal{C}}\sum_{j\in\mathcal{T}}\pi_{ij}\log\tfrac{\pi_{ij}}{e},$
(1) $\displaystyle\mathrm{subject~{}to}\quad$
$\displaystyle\sum_{i\in\mathcal{C}}\pi_{ij}=1/N_{t}\quad\forall
j\in\mathcal{T},$ (2)
$\displaystyle\sum_{j\in\mathcal{T}}\pi_{ij}=1/N_{c}\quad\forall
i\in\mathcal{C}.$ (3)

The first term in the objective function (1) is the average of the squared
approximation errors for covariate balancing. In the first constraint (2), for
each treated unit $j$, the synthetic weights
$(\tfrac{\pi_{ij}}{1/N_{t}})_{i\in\mathcal{C}}$ sum to $1$ and thus the term
$\sum_{i\in\mathcal{C}}\tfrac{\pi_{ij}}{1/N_{t}}x_{i}$ in (1) is a convex
combination to approximate $x_{j}$.

The second term in (1) is an entropic regularization with regularization
parameter $\lambda>0$. Increasing the parameter $\lambda$ encourages the
solution $\pi$ to be uniform, namely, $\pi_{ij}=\tfrac{1}{N_{c}\cdot N_{t}}$
for all $(i,j)$’s, which essentially corresponds to increasing the number of
neighbors in the nearest neighbor matching.

The second constraint (3) enforces that all control units contribute equally
and thus seeks a matching between treated and control data sets; this
constraint can be generalized to incorporate inverse propensity scores, to be
shown in the next section. If this constraint is dropped, the optimization can
be decoupled into $N_{t}$ independent optimization of the same type, each with
one treated unit $j\in\mathcal{T}$, also solvable by our optimization
algorithm detailed in Section 4. We shall show in a second that this
constraint will be crucial. It enables the convex program to synthesize
granular information while enforcing, at the coarse level, that the answers
agree with typical estimators for the ATE or ATT.

Finally, the optimization variable $\pi\in\mathbb{R}_{+}^{N_{c}\times N_{t}}$
is a matrix whose entries are indexed by $(i,j)$ for $i\in\mathcal{C}$ and
$j\in\mathcal{T}$, instead of the usual indexing by $i\in\\{1,\ldots,N_{c}\\}$
and $j\in\\{1,\ldots,N_{t}\\}$, to denote that $\pi_{ij}$ is the weight
between control unit $i$ and treated unit $j$. We call
$\pi\in\mathbb{R}_{+}^{N_{c}\times N_{t}}$ satisfying the constraints (2) and
(3) a coupling, which generalizes doubly stochastic matrices to non-square
matrices with prescribed row and column sums.

Upon obtaining the solution, denoted by $\widehat{\pi}$, we impute the
counterfactual outcome $Y_{j}(0)$ of treated unit $j$ by the convex
combination of the outcomes of control units with the weights
$(N_{t}\widehat{\pi}_{ij})_{i\in\mathcal{C}}$:

$\widehat{Y}_{j}(0):=\sum_{i\in\mathcal{C}}\tfrac{\widehat{\pi}_{ij}}{1/N_{t}}Y_{i},$
(4)

which in turn leads to the individual treatment effect estimate
$\widehat{\tau}_{j}:=Y_{j}-\widehat{Y}_{j}(0)$. By the constraint (3), the
average of the estimated individual treatment effects coincides with the
difference in the mean outcomes between the treatment and control groups
denoted as $\widehat{\tau}^{\mathsf{DiM}}$:

$\frac{1}{N_{t}}\sum_{j\in\mathcal{T}}\widehat{\tau}_{j}=\frac{1}{N_{t}}\sum_{j\in\mathcal{T}}(Y_{j}-\sum_{i\in\mathcal{C}}\tfrac{\widehat{\pi}_{ij}}{1/N_{t}}Y_{i})=\frac{1}{N_{t}}\sum_{j\in\mathcal{T}}Y_{j}-\frac{1}{N_{c}}\sum_{i\in\mathcal{C}}Y_{i}=:\widehat{\tau}^{\mathsf{DiM}}.$
(5)

Therefore, the proposed convexified matching allows for estimating individual
treatment effects while maintaining the desired ATT estimate
$\widehat{\tau}^{\mathsf{DiM}}$.

The resulting optimization is a smooth convex optimization problem. To see
this, notice that the first term of the objective function (1) is the
following quadratic function of $\pi$:

$\frac{N_{t}}{2}\langle\pi,K_{cc}\pi\rangle-\langle\pi,K_{ct}\rangle+\frac{\mathrm{tr}(K_{tt})}{2N_{t}},$
(6)

where $K_{cc}\in\mathbb{R}^{N_{c}\times N_{c}}$,
$K_{ct}\in\mathbb{R}^{N_{c}\times N_{t}}$, and
$K_{tt}\in\mathbb{R}^{N_{t}\times N_{t}}$ are the Gram matrices whose entries
are inner products of the covariates, namely, $K_{cc}=(\langle
x_{i},x_{i^{\prime}}\rangle)_{i,i^{\prime}\in\mathcal{C}}$, $K_{ct}=(\langle
x_{i},x_{j}\rangle)_{(i,j)\in\mathcal{C}\times\mathcal{T}}$, and
$K_{c,t}=(\langle x_{j},x_{j^{\prime}}\rangle)_{j,j^{\prime}\in\mathcal{T}}$.
This quadratic function is convex as the Gram matrix $K_{cc}$ is positive
semidefinite. The second term—the entropy term—is convex. The constraint set
resulting from (2) and (3) is a convex polytope consisting of couplings. In
Section 4, we introduce and analyze efficient algorithms to solve this
optimization problem based on a simple iterative matrix scaling procedure
called the Sinkhorn algorithm [Sin67].

Lastly, we discuss several existing methods and concepts related to the
proposed method. [AL21] considers multiple treated units ($N_{t}>1$) and
proposes running synthetic controls separately for treated units with a
different regularization function. The proposed method differs from running
synthetic controls separately because of the constraint (3). If there is only
one treated unit ($N_{t}=1$) and no constraint (3), the proposed method is
equivalent to the standard synthetic control with entropic regularization as
in [Hai12] which proposes entropic regularization for covariate balancing.
Also, notice that the first term of (1) is upper bounded under the constraint
(2) as follows:

$\frac{1}{2N_{t}}\sum_{j\in\mathcal{T}}\|x_{j}-\sum_{i\in\mathcal{C}}\tfrac{\pi_{ij}}{1/N_{t}}x_{i}\|_{2}^{2}=\frac{1}{2N_{t}}\sum_{j\in\mathcal{T}}\|\sum_{i\in\mathcal{C}}\tfrac{\pi_{ij}}{1/N_{t}}(x_{j}-x_{i})\|_{2}^{2}\leq\frac{1}{2N_{t}}\sum_{j\in\mathcal{T}}\sum_{i\in\mathcal{C}}\tfrac{\pi_{ij}}{1/N_{t}}\|x_{j}-x_{i}\|_{2}^{2},$

where the inequality uses Jensen’s inequality. This upper bound analytically
shows that the average approximation error based on convex combinations is
smaller than that based on pairwise distances, which aligns with the
motivation mentioned in the previous section. Replacing the first term of the
objective function (1) by the above upper bound—while maintaining the
constraints (2), (3)—leads to the optimal transport problem [Vil03] with
entropic regularization [Cut13].

###  2.2 Extensions

We present two extensions of the simple formulation introduced in Section 2.1.
The first is extending the main insights to cover the nonlinear case when the
potential outcomes can be nonlinear functions of covariates, using a kernel
trick [SC04, SC08]. The second is incorporating notions of propensity scores
[RR83] into our synthetic coupling formulation.

#### Kernel Synthetic Coupling

We can kernelize the presented simple formulation to extend to nonlinear
cases. The kernel trick is widely used in machine learning to modify methods
that originated in the linear setting to accommodate nonlinearity, such as
ridge regression, support vector machines, and principal component analysis.
Recall that our objective (1) relies on the Gram matrices as shown in (6).
Therefore, to kernelize the simple formulation, we can replace the Gram
matrices with the kernel matrices, namely, letting
$K_{cc}=(k(x_{i},x_{i^{\prime}}))_{i,i^{\prime}\in\mathcal{C}}$,
$K_{ct}=(k(x_{i},x_{j}))_{(i,j)\in\mathcal{C}\times\mathcal{T}}$, and
$K_{tt}=(k(x_{j},x_{j^{\prime}}))_{j,j^{\prime}\in\mathcal{T}}$ for some
kernel function $k\colon\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}$. The
kernelized formulation is essentially the same convex optimization problem. It
can be tackled by the algorithms we introduce in Section 4.

Applying the above kernel trick is equivalent to replacing the first term of
(1) with the approximation error in the reproducing kernel Hilbert space
(RKHS) $\mathcal{H}$ associated with the kernel $k$, namely,
$\frac{1}{2N_{t}}\sum_{j\in\mathcal{T}}\|\phi_{x_{j}}-\sum_{i\in\mathcal{C}}\tfrac{\pi_{ij}}{1/N_{t}}\phi_{x_{i}}\|_{\mathcal{H}}^{2}$,
where $\|\cdot\|_{\mathcal{H}}$ is the norm in the RKHS $\mathcal{H}$ and
$\phi_{x}=k(x,\cdot)$ denotes the canonical feature of $x\in\mathbb{R}^{d}$.
In other words, we generalize the matching of covariates in the Euclidean
space to the matching of embedded canonical features in the RKHS
$\mathcal{H}$, which leads to the following kernelized formulation:

$\displaystyle\mathrm{KSC}^{v,w}_{\lambda}:=\min_{\pi\in\mathbb{R}_{+}^{N_{c}\times
N_{t}}}$
$\displaystyle\frac{1}{2}\sum_{j\in\mathcal{T}}v_{j}\|\phi_{x_{j}}-\sum_{i\in\mathcal{C}}\tfrac{\pi_{ij}}{v_{j}}\phi_{x_{i}}\|_{\mathcal{H}}^{2}+\lambda\sum_{i\in\mathcal{C}}\sum_{j\in\mathcal{T}}\pi_{ij}\log\tfrac{\pi_{ij}}{e},$
(7) $\displaystyle\mathrm{subject~{}to}$
$\displaystyle\sum_{i\in\mathcal{C}}\pi_{ij}=v_{j}\quad\forall
j\in\mathcal{T}\quad\text{and}\quad\sum_{j\in\mathcal{T}}\pi_{ij}=w_{i}\quad\forall
i\in\mathcal{C},$

where $v=(v_{j})_{j\in\mathcal{T}}\in\Delta_{N_{t}}^{+}$ and
$w=(w_{i})_{i\in\mathcal{C}}\in\Delta_{N_{c}}^{+}$ are strictly positive
probability vectors that generalize the constraints (2) and (3) to allow
treated and control units to have different weights.111Whenever
$v=\frac{1}{N_{t}}1_{N_{t}}$, we will drop the dependency on $v$ and
abbreviate as $\mathrm{KSC}^{w}_{\lambda}$. With the optimal coupling
$\widehat{\pi}$, we impute using convex combination same as in (4), a classic
idea in nonparametric regression.

#### Propensity Score

Using arbitrary weights through
$v=(v_{j})_{j\in\mathcal{T}}\in\Delta_{N_{t}}^{+}$ and
$w=(w_{i})_{i\in\mathcal{C}}\in\Delta_{N_{c}}^{+}$ allows for various average
treatment effect estimators, particularly including those incorporating
propensity scores. To estimate the ATT, we can set
$v=\frac{1}{N_{t}}1_{N_{t}}$ as before, and the average of the estimated
individual treatment effects of the treated leads to

$\frac{1}{N_{t}}\sum_{j\in\mathcal{T}}\widehat{\tau}_{j}=\frac{1}{N_{t}}\sum_{j\in\mathcal{T}}(Y_{j}-\sum_{i\in\mathcal{C}}\tfrac{\widehat{\pi}_{ij}}{1/N_{t}}Y_{i})=\frac{1}{N_{t}}\sum_{j\in\mathcal{T}}Y_{j}-\sum_{i\in\mathcal{C}}w_{i}Y_{i}.$

Now, let $\widehat{p}\colon\mathbb{R}^{d}\to(0,1)$ be a suitable
estimator—often estimated by logistic regression—of the propensity score,
$x\mapsto\mathbb{P}(W=1\,|\,X=x)$. If the weights $w$ are chosen based on the
propensity scores as

$w_{i}=\tfrac{\widehat{p}(x_{i})}{1-\widehat{p}(x_{i})}/\sum_{i^{\prime}\in\mathcal{C}}\tfrac{\widehat{p}(x_{i^{\prime}})}{1-\widehat{p}(x_{i^{\prime}})}\quad\forall
i\in\mathcal{C},$ (8)

the aggregate effects match the normalized Inverse Probability Weighting (IPW)
estimator [HI01], that is,

$\frac{1}{N_{t}}\sum_{j\in\mathcal{T}}\widehat{\tau}_{j}=\frac{\sum_{j\in\mathcal{T}}Y_{j}}{N_{t}}-\frac{\sum_{i\in\mathcal{C}}\tfrac{\widehat{p}(x_{i})}{1-\widehat{p}(x_{i})}Y_{i}}{\sum_{i\in\mathcal{C}}\tfrac{\widehat{p}(x_{i})}{1-\widehat{p}(x_{i})}}=:\widehat{\tau}^{\mathsf{IPW}}_{\mathsf{ATT}}.$
(9)

Namely, the proposed KSC method is guaranteed to yield individual treatment
effect estimates that conform to (9) at the aggregate level.

Similarly, solving $\mathrm{KSC}^{v,w}_{\lambda}$ in (7) by setting

$\displaystyle
v_{j}=\tfrac{1}{\widehat{p}(x_{j})}/\sum_{j^{\prime}\in\mathcal{T}}\tfrac{1}{\widehat{p}(x_{j^{\prime}})}\quad\forall
j\in\mathcal{T}\quad\text{and}\quad
w_{i}=\tfrac{1}{1-\widehat{p}(x_{i})}/\sum_{i^{\prime}\in\mathcal{C}}\tfrac{1}{1-\widehat{p}(x_{i^{\prime}})}\quad\forall
i\in\mathcal{C},$

we obtain an coupling $\widehat{\pi}$ that matches the IPW estimator for ATE

$\displaystyle\sum_{j\in\mathcal{T}}v_{j}\widehat{\tau}_{j}=\frac{\sum_{j\in\mathcal{T}}\tfrac{1}{\widehat{p}(x_{j})}Y_{j}}{\sum_{j\in\mathcal{T}}\tfrac{1}{\widehat{p}(x_{j})}}-\frac{\sum_{i\in\mathcal{C}}\tfrac{1}{1-\widehat{p}(x_{i})}Y_{i}}{\sum_{i\in\mathcal{C}}\tfrac{1}{1-\widehat{p}(x_{i})}}=:\widehat{\tau}^{\mathsf{IPW}}_{\mathsf{ATE}}.$

###  2.3 Numeric Illustration

Now, we demonstrate our imputation method, KSC, using a real-world data set.
This section aims to concisely present the features of KSC and compare it with
other imputation methods; further details and robust evaluations are provided
in Section 5 together with the inference procedure, to be presented in the
next section.

We apply the proposed method to evaluate the National Supported Work (NSW)
demonstration program, first analyzed by [LaL86]. This data set has been
widely studied in the program evaluation literature. The NSW program was
conducted from 1975 to 1978 in the United States, which aimed at providing a
job training program to disadvantaged workers, where the treatment, namely,
the training, was randomly assigned to some of them. The data we consider here
is a specific subset of the experimental data used in [DW99], which consists
of $N=445$ subjects with $N_{t}=185$ treated units and $N_{c}=260$ control
units. The outcome of interest $Y_{i}$ is the post-treatment earnings recorded
in 1978. Each individual is associated with six variables: age, years of
education, and four indicator variables denoting whether the individual is
black, Hispanic, married, and a high school dropout. Following [DW99], the
pretreatment outcomes, the earnings in 1974 and 1975, and two indicator
variables denoting whether these pretreatment earnings are zero, are included
as covariates. Accordingly, we have $x_{i}\in\mathbb{R}^{d}$ with $d=10$.

(a) NSW $\lambda=0.001$ (b) NSW $\lambda=0.01$

(c) PSID $\lambda=0.001$ (d) PSID $\lambda=0.01$

Figure 1: Scatter plots of the imputed counterfactual outcomes of the treated
along with the histograms of the marginal distributions. The $x$-coordinate is
$\widehat{Y}_{j}(0)$ imputed by the proposed method with
$\lambda\in\\{0.001,0.01\\}$, while the $y$-coordinate is the imputed
counterfactual outcomes by the nearest neighbor matching with replacement or
the regression method. (a) and (b) are based on the NSW experimental data,
while (c) and (d) are based on the PSID data after trimming. The proposed
method uses the linear kernel for both data. For the NSW data, we use uniform
weights $v,w$ in (7), while for the PSID data, we set $w$ to be the propensity
score-based weights following (8).

Figure 1(a) and Figure 1(b) show the scatter plots of the imputed
counterfactual outcomes of the treated units, where the $x$-coordinate is
$\widehat{Y}_{j}(0)$ imputed by the proposed method, while the $y$-coordinate
is the imputed counterfactual outcomes by the nearest neighbor matching with
replacement or the regression imputation method. The proposed method relies on
the linear kernel following the simple formulation in Section 2.1 with
$\lambda\in\\{0.001,0.01\\}$, the nearest neighbor matching uses $5$ and $10$
nearest neighbors, and the regression method fits a linear regression model to
the observed outcome using all the covariates and the treatment indicator. At
first glance, the scatter plots visually suggest that all three imputation
methods provide comparable imputed counterfactual outcomes. One crucial
difference is that the average of the imputed counterfactual outcomes and,
thus, the estimated individual treatment effects vary across the methods. The
ATT estimate computed by the proposed method always matches—regardless of the
choice of $\lambda$—the mean difference between the treatment and control
groups, which is roughly 1794.3, as shown in (5)—the unbiased estimator of the
ATT in the experimental setup. In contrast, the $k$-nearest neighbor matching
with $k=5$, with $k=10$, and the regression imputation method produce the ATT
estimates 2030.5, 1776.6, and 1706.2, respectively. While the average of the
imputed counterfactual outcomes does not depend on the choice of $\lambda$, we
can see that it determines the level of individualization. The histograms
along the $x$-axix show that the imputed counterfactual outcomes based on the
proposed method with $\lambda=0.001$ are more dispersed than those based on
$\lambda=0.01$; the latter is based on larger regularization leading to more
uniform weights, which is similar to the fact that $k=5$ is more
dispersed—along the $y$-axis—than $k=10$ in the nearest neighbor matching.

Since [LaL86], there has been a debate on the credibility of the average
treatment effect estimation using observational control groups. Since such
control groups can be significantly different from the control group of the
experimental data, the resulting ATE estimates can deviate much from the
experimental benchmark, say, the mean difference of 1794.3 computed earlier,
which has been the central argument of [LaL86]. Later, several methods [DW99,
DW02] have been proposed to recover the experimental benchmark estimate using
observational control groups. Here, we focus on one such method based on the
normalized inverse probability weighting (IPW) estimator (9) and apply the KSC
method with the weights based on the propensity scores. To this end, we
consider a control group from the Panel Study of Income Dynamics (PSID) data,
consisting of 2490 units, with which we estimate the propensity scores by
logistic regression. The resulting ATT estimate by
$\widehat{\tau}^{\mathsf{IPW}}_{\mathsf{ATT}}$ is 2579.7, which is far from
the benchmark of 1794.3, rooted in the stark difference between the
experimental and observational control groups. To remedy this, we trim the
PSID data by only taking the units whose propensity scores are in
$[0.05,0.95]$, yielding 214 units. The resulting ATT estimate by
$\widehat{\tau}^{\mathsf{IPW}}_{\mathsf{ATT}}$ is then 1748.0, which resembles
the experimental benchmark. Taking this trimmed PSID data as a control group,
we apply the proposed method with the weights $w_{i}$ based on the propensity
scores as in (8) and compare with the matching and regression methods. The
results are shown in Figure 1(c) and Figure 1(d). Again, though visually
similar in distributions, the ATT estimates differ across the methods, where
the proposed method with the propensity score-based weights provides the ATT
estimate 1748.0, regardless of $\lambda$. In contrast, the $k$-nearest
neighbor matching with $k=5$, with $k=10$, and the regression method produce
the ATT estimates 1567.3, 1003.0, and 1065.5, respectively. This example
illustrates that the proposed imputation method allows for producing granular
information, namely, the individual treatment effects, from a coarse level
estimate by $\widehat{\tau}^{\mathsf{IPW}}_{\mathsf{ATT}}$. The coarse level
estimate is closer to the experimental benchmark estimate than inherently
imputation-based methods such as the nearest neighbor matching or regression.

##  3 Inference: Individual Confidence Intervals

Constructing confidence intervals for the estimated individual treatment
effects is conducive to reliable decision-making. This section introduces and
analyzes the method to construct individual confidence intervals based on the
proposed convexified matching. To this end, we build a confidence interval for
the missing counterfactual outcome $Y_{j}(0)$ for each treated unit
$j\in\mathcal{T}$. We first introduce appropriate model assumptions. Then, we
elucidate the construction of individual confidence intervals, followed by a
numerical example, and discuss the coverage and efficiency of the proposed
method.

###  3.1 Setup

#### The Model

Let $\mathcal{H}$ be the reproducing kernel Hilbert space (RKHS) associated
with a kernel function
$k\colon\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}$, where
$\langle\cdot,\cdot\rangle_{\mathcal{H}}$ denotes the RKHS inner product,
$\|\cdot\|_{\mathcal{H}}$ denotes the RKHS norm, and
$\phi_{x}:=k(x,\cdot)\in\mathcal{H}$ is the corresponding reproducing kernel
function. We postulate the following model for the potential outcomes: for
each $i=1,\ldots,N$, letting $x_{i}\in\mathbb{R}^{d}$ be the covariate of unit
$i$, we have

$\displaystyle Y_{i}(0)$
$\displaystyle=f_{0}(x_{i})+\varepsilon_{0,i},~{}\varepsilon_{0,i}\sim\mathcal{N}(0,\sigma_{0}^{2}),$
$\displaystyle Y_{i}(1)$
$\displaystyle=f_{1}(x_{i})+\varepsilon_{1,i},~{}\varepsilon_{1,i}\sim\mathcal{N}(0,\sigma_{1}^{2}),$

where $f_{0},f_{1}\in\mathcal{H}$ are some functions in the RKHS and
$\sigma_{0},\sigma_{1}\geq 0$. Let $\mathbf{X}:=[x_{1},\ldots,x_{N}]$ denote
the fixed design matrix and let $\mathbf{Z}:=[Z_{1},\ldots,Z_{N}]$ denote the
possibly random treatment assignment vector.

######  Assumption 1.

$\\{\varepsilon_{0,i},\varepsilon_{1,i}\\}_{i=1}^{N}$ are independent and

$\mathbf{Z}\mathrel{\perp\\!\\!\\!\perp}\\{\varepsilon_{0,i},\varepsilon_{1,i}\\}_{i=1}^{N}~{}|~{}\mathbf{X},$
(10)

that is, the treatment assignment and the errors are independent given the
covariates.

######  Remark 1.

Typical assumptions in causal inference involve the following: (i)
$(Y_{i}(0),Y_{i}(1),X_{i})_{i=1}^{N}$ are independently drawn from some
population distribution defined on
$\mathbb{R}\times\mathbb{R}\times\mathbb{R}^{d}$ and (ii) the treatment
assignment $Z_{i}$ is independent of $(Y_{i}(0),Y_{i}(1))$ conditional on
$X_{i}$. Unlike (i), our model postulates a signal-plus-noise structure for
the potential outcomes, where the covariates are from a fixed design matrix.
All our results are conditioned on the fixed design $\mathbf{X}$, which we
highlight by using the lowercase $\\{x_{i}\\}$ notation instead of
$\\{X_{i}\\}$. On the one hand, Assumption 1 serves the same purpose as (ii)
in spirit but is stronger in the sense that the errors are assumed to be
Gaussian. On the other hand, we do not require $\mathbf{X}$ nor $\mathbf{Z}$
to be i.i.d. drawn, allowing $Z_{i}$’s to have arbitrary dependence across
units or even to be adversarially chosen conditioned on the fixed design
$\mathbf{X}$.

Under this setup, for each treated unit $j\in\mathcal{T}$, the goal is to
build an individual confidence interval for its missing counterfactual outcome
$Y_{j}(0)$. More specifically, we aim to construct a confidence interval
$[\widehat{L}_{j},\widehat{U}_{j}]$ that contains the conditional expectation
$f_{0}(x_{j})=\operatorname*{\mathbb{E}}[Y_{j}(0)\,|\,x_{j}]=\operatorname*{\mathbb{E}}[Y_{j}(0)\,|\,\mathbf{X},\mathbf{Z}]$
with a desired level of confidence as follows:

$\mathbb{P}\left(f_{0}(x_{j})\in[\widehat{L}_{j},\widehat{U}_{j}]~{}|~{}\mathbf{X},\mathbf{Z}\right)\geq
1-\alpha,$ (11)

namely, the probability that the interval $[\widehat{L}_{j},\widehat{U}_{j}]$
contains the conditional expectation $f_{0}(x_{j})$ is at least $1-\alpha$,
where the probability is conditional on the design matrix $\mathbf{X}$ and the
treatment assignments $\mathbf{Z}$.

###  3.2 Individual Confidence Intervals

#### Construction

We consider the following Kernerlized Synthetic Coupling (KSC) discussed in
Section 2.2. To construct an individual interval for the counterfactual
control outcome of each treated unit, we solve (7) with uniform weights
$v=\frac{1}{N_{t}}$ on the treated, which we rewrite in the matrix form as
follows: let $(w_{i})_{i\in\mathcal{C}}$ be positive numbers such that
$\sum_{i\in\mathcal{C}}w_{i}=1$,

$\displaystyle\mathrm{KSC}^{w}_{\lambda}=\min_{\pi\in\mathbb{R}_{+}^{N_{c}\times
N_{t}}}$
$\displaystyle\tfrac{N_{t}}{2}\langle\pi,K_{cc}\pi\rangle-\langle\pi,K_{ct}\rangle+\tfrac{\mathrm{tr}(K_{tt})}{2N_{t}}+\lambda\langle\pi,\log\tfrac{\pi}{e}\rangle,$
(12) $\displaystyle\mathrm{subject~{}to}$
$\displaystyle\sum_{i\in\mathcal{C}}\pi_{ij}=\tfrac{1}{N_{t}}\quad\forall
j\in\mathcal{T}\quad\text{and}\quad\sum_{j\in\mathcal{T}}\pi_{ij}=w_{i}\quad\forall
i\in\mathcal{C}.$

Let $\widehat{\pi}$ be the solution of the above optimization problem and
define a probability transition matrix
$\widehat{P}=(\widehat{p}_{ij})_{(i,j)\in\mathcal{C}\times\mathcal{T}}$, where
$\widehat{p}_{ij}=N_{t}\widehat{\pi}_{ij}$ so that
$\sum_{i\in\mathcal{C}}\widehat{p}_{ij}=1$ for any $j\in\mathcal{T}$. For each
treated unit $j$, recall the imputed value $\widehat{Y}_{j}(0)$ defined in
(4), for which we construct a confidence interval
$[\widehat{L}_{j},\widehat{U}_{j}]$, where $\widehat{L}_{j}$ and
$\widehat{U}_{j}$ are defined as follows:

$\displaystyle\widehat{L}_{j}$
$\displaystyle:=\widehat{Y}_{j}(0)-\widehat{\theta}\cdot\sqrt{\left(K_{tt}+\widehat{P}^{\top}K_{cc}\widehat{P}-2K_{ct}^{\top}\widehat{P}\right)_{jj}}-z_{1-\frac{\alpha}{2}}\cdot\widehat{\sigma}_{0}\sqrt{\left(\widehat{P}^{\top}\widehat{P}\right)_{jj}},$
(13) $\displaystyle\widehat{U}_{j}$
$\displaystyle:=\widehat{Y}_{j}(0)+\widehat{\theta}\cdot\sqrt{\left(K_{tt}+\widehat{P}^{\top}K_{cc}\widehat{P}-2K_{ct}^{\top}\widehat{P}\right)_{jj}}+z_{1-\frac{\alpha}{2}}\cdot\widehat{\sigma}_{0}\sqrt{\left(\widehat{P}^{\top}\widehat{P}\right)_{jj}}.$
(14)

Here, $z_{1-\frac{\alpha}{2}}$ is the $(1-\frac{\alpha}{2})$-quantile of the
standard normal distribution and $\widehat{\theta},\widehat{\sigma}_{0}$ are
suitable estimates of the norm $\|f_{0}\|_{\mathcal{H}}$ and the standard
deviation $\sigma_{0}$, respectively.

We apply kernel ridge regression to the control group data to estimate
$\|f_{0}\|_{\mathcal{H}}$ and $\sigma_{0}$. We can estimate
$f_{0}\in\mathcal{H}$ by
$\widehat{f}_{0}(x)=\sum_{i\in\mathcal{C}}\widehat{\beta}_{i}k(x_{i},x)$,
where $\widehat{\beta}=(K_{cc}+\rho I_{N_{c}})^{-1}Y_{c}\in\mathbb{R}^{N_{c}}$
with a kernel matrix $K_{cc}$ defined in Section 2.2, a vector of control
outcomes $Y_{c}=(Y_{i})_{i\in\mathcal{C}}\in\mathbb{R}^{N_{c}}$, and a ridge
regularization parameter $\rho>0$. Then, we estimate $\|f_{0}\|_{\mathcal{H}}$
by

$\widehat{\theta}:=\|\widehat{f}_{0}\|_{\mathcal{H}}=\sqrt{\langle\widehat{\beta},K_{cc}\widehat{\beta}\rangle},$

and the variance $\sigma_{0}^{2}$ by

$\widehat{\sigma}_{0}^{2}=\frac{1}{N_{c}}\sum_{i\in\mathcal{C}}(Y_{i}-\widehat{f}_{0}(x_{i}))^{2}=\frac{1}{N_{c}}\langle
Y_{c}-K_{cc}\widehat{\beta},Y_{c}-K_{cc}\widehat{\beta}\rangle.$

Before diving into each component of the construction, we first look at a
simple numerical example to illustrate the performance of the proposed method.

#### Numerical Example

We consider a set of fixed one-dimensional covariates
$x_{1},\ldots,x_{N}\in[0,1]$ with $N=500$ and we randomly choose $N_{t}=200$
treatment units. Let $\mathcal{H}$ be the RKHS associated with the Gaussian
kernel $k(x,x^{\prime})=\exp(-\gamma\cdot|x-x^{\prime}|^{2})$, where
$\gamma=2.5$, and suppose $f_{0}(x)=k(0.5,x)$, which yields
$\|f_{0}\|_{\mathcal{H}}=(k(0.5,0.5))^{1/2}=1$. Then, we generate the
potential outcomes by $Y_{i}(0)=f_{0}(x_{i})+\varepsilon_{0,i}$, where
$\varepsilon_{0,1},\ldots,\varepsilon_{0,N}$ are independently drawn from
$N(0,\sigma_{0}^{2})$. As $x_{1},\ldots,x_{N}$ and $Z_{1},\ldots,Z_{N}$ are
fixed, the randomness is solely from the residuals $\varepsilon_{i}$’s, which
we sample 1000 times.

(a) $\sigma_{0}=0.1$ (b) $\sigma_{0}=1$

(c) $\sigma_{0}=3$

Figure 2: Confidence intervals. The potential outcomes are generated by
$Y_{i}(0)=f_{0}(x_{i})+\varepsilon_{0,i}$, where
$\varepsilon_{0,1},\ldots,\varepsilon_{0,N}$ are independently drawn from
$N(0,\sigma_{0}^{2})$. The black dashed curve is the ground truth conditional
expectation function $f_{0}(x)=\exp(-2.5\times|x-0.5|^{2})$. The blue solid
curve shows the estimated counterfactual outcome
$\widehat{Y}_{j}(0)=\sum_{i\in\mathcal{C}}\widehat{p}_{ij}Y_{i}$ for each
treated unit $j$. The blue shaded region represents
$[\widehat{L}_{j},\widehat{U}_{j}]$ defined in (13) and (14), where
$\widehat{\theta}$ and $\widehat{\sigma}_{0}$ are estimated by the kernel
ridge regression as explained in Section 3.2. The oracle interval
$[L_{j}^{\ast},U_{j}^{\ast}]$ defined in (19) and (20) is shown using the red
error bars.

Figure 2 shows the constructed confidence intervals
$[\widehat{L}_{j},\widehat{U}_{j}]$ for each treated unit $j$ for different
values of $\sigma_{0}$. For each $\sigma_{0}\in\\{0.1,1,3\\}$, we solve (7)
for $\lambda\in\\{0.1,0.01,0.001\\}$ with uniform weights, run the kernel
ridge regression to obtain $\widehat{\theta},\widehat{\sigma}_{0}$, and
produce $[\widehat{L}_{j},\widehat{U}_{j}]$ which is shown as blue shaded
regions in the plots. We compare this interval with an “oracle” interval
$[L_{j}^{\ast},U_{j}^{\ast}]$, which—as we will explain in the next section—
is guaranteed to have an exact $1-\alpha$ coverage, namely,

$\mathbb{P}\left(f_{0}(x_{j})\in[L_{j}^{\ast},U_{j}^{\ast}]~{}|~{}\mathbf{X},\mathbf{Z}\right)=1-\alpha.$
(15)

By construction, if the estimates $\widehat{\theta},\widehat{\sigma}_{0}$ are
accurate enough, the interval $[\widehat{L}_{j},\widehat{U}_{j}]$ is supposed
to contain the oracle interval $[L_{j}^{\ast},U_{j}^{\ast}]$, which can be
verified in Figure 2. For each $\sigma_{0}$, decreasing $\lambda$ results in
smaller approximation errors, namely, the bias correction reflected in
$\left(K_{tt}+\widehat{P}^{\top}K_{cc}\widehat{P}-2K_{ct}^{\top}\widehat{P}\right)_{jj}^{1/2}$,
and the estimated counterfactual outcomes (solid blue curves) that are more
wiggly as the obtained weights have more variations. As we will see,
$[\widehat{L}_{j},\widehat{U}_{j}]$ is wider than
$[L_{j}^{\ast},U_{j}^{\ast}]$ by the twice of the bias correction, namely,

$(\widehat{U}_{j}-\widehat{L}_{j})-(U_{j}^{\ast}-L_{j}^{\ast})=2\cdot\widehat{\theta}\cdot\left(K_{tt}+\widehat{P}^{\top}K_{cc}\widehat{P}-2K_{ct}^{\top}\widehat{P}\right)_{jj}^{1/2}.$

We can verify from Figure 2 that this gap between the two intervals gets
smaller as $\lambda$ decreases, consistent with the above observation that the
bias decreases.

Figure 3 plots the coverage of $[\widehat{L}_{j},\widehat{U}_{j}]$ and
$[L_{j}^{\ast},U_{j}^{\ast}]$, namely, the probability that they contain
$f_{0}(x_{j})$ estimated using the 1,000 samples. When $\lambda=0.1$, we can
see that the coverage of $[\widehat{L}_{j},\widehat{U}_{j}]$ is almost $1$,
meaning that the interval is too conservative, which results from the fact
that the bias term is not small enough compared to the variance of the
residuals. As $\lambda$ decreases, we can see that the coverage of
$[\widehat{L}_{j},\widehat{U}_{j}]$ gets closer to $1-\alpha=0.95$, which is
consistent with the fact that the bias decreases. Particularly, when
$\sigma_{0}=3$ and $\lambda=0.001$, the coverage of
$[\widehat{L}_{j},\widehat{U}_{j}]$ is closest to the coverage of the oracle
interval $[L_{j}^{\ast},U_{j}^{\ast}]$, suggesting that the bias is dominated
by the variance. Meanwhile, for $(\sigma_{0},\lambda)=(1,0.01)$ and
$(\sigma_{0},\lambda)=(3,0.1)$, the coverage of
$[\widehat{L}_{j},\widehat{U}_{j}]$ is larger than the coverage of the oracle
interval $[L_{j}^{\ast},U_{j}^{\ast}]$ by approximately a constant, which
suggests that the bias term and the variance term are nearly balanced. In the
next section, we provide theoretical guidance to choose the regularization
parameter $\lambda$ to balance the bias and variance tradeoff.

(a) $\sigma_{0}=0.1$ (b) $\sigma_{0}=1$

(c) $\sigma_{0}=3$

Figure 3: Coverage of the confidence intervals
$[\widehat{L}_{j},\widehat{U}_{j}]$ and $[L_{j}^{\ast},U_{j}^{\ast}]$
estimated using the 1,000 samples. The red dashed line shows the nominal
coverage $1-\alpha=0.95$.

###  3.3 Coverage and Efficiency

Now, we dive into each component in constructing individual confidence
intervals to dissect the coverage and efficiency. We discuss how the entropic
regularization parameter $\lambda$ drives the bias and variance tradeoff, as
seen in the numerical simulations. Furthermore, we will show that optimizing
the entropic-regularized convexified matching program (7) with the desired
$\lambda$ chosen above is closely connected to optimizing the average of the
squared lengths of the individual confidence intervals
$\frac{1}{N_{t}}\sum_{j\in\mathcal{T}}\mathrm{Len}_{j}^{2}$.

#### Point Estimate

For each treated unit $j\in\mathcal{T}$, the center of the confidence interval
$[\widehat{L}_{j},\widehat{U}_{j}]$ is the imputed counterfactual outcome
$\widehat{Y}_{j}(0)$:

$\displaystyle\widehat{Y}_{j}(0)=\sum_{i\in\mathcal{C}}\widehat{p}_{ij}Y_{i}(0).$
(16)

This weighted average estimator leverages all the information in $\mathbf{X}$
to anchor an optimal coupling (or matching) between the treated and control
units for synthesizing the controls. It uses convex combinations of the
control units to synthesize the granular counterfactual outcomes for the
treated units, while balancing the aggregate level statistics, the ATT. We
remark that the optimal coupling $\widehat{\pi}$—and thus
$\widehat{P}=N_{t}\widehat{\pi}$ —is calculated based solely on $\mathbf{X}$
and $\mathbf{Z}$, not the outcomes $Y_{i}$’s. Recall the following bias-
variance decomposition of the point estimate $\widehat{Y}_{j}(0)$:

$\widehat{Y}_{j}(0)-\operatorname*{\mathbb{E}}[Y_{j}(0)\,|\,\mathbf{X},\mathbf{Z}]=\underbrace{\widehat{Y}_{j}(0)-\operatorname*{\mathbb{E}}[\widehat{Y}_{j}(0)\,|\,\mathbf{X},\mathbf{Z}]}_{:=\mathcal{V}_{j}}+\underbrace{\operatorname*{\mathbb{E}}[\widehat{Y}_{j}(0)\,|\,\mathbf{X},\mathbf{Z}]-\operatorname*{\mathbb{E}}[Y_{j}(0)\,|\,\mathbf{X},\mathbf{Z}]}_{:=\mathcal{B}_{j}}.$

#### Bias-Awareness

The point estimate could be biased as

$\mathcal{B}_{j}=\operatorname*{\mathbb{E}}[\widehat{Y}_{j}(0)\,|\,\mathbf{X},\mathbf{Z}]-\operatorname*{\mathbb{E}}[Y_{j}(0)\,|\,\mathbf{X},\mathbf{Z}]=\sum_{i\in\mathcal{C}}\widehat{p}_{ij}f_{0}(x_{i})-f_{0}(x_{j})=\langle
f_{0},\sum_{i\in\mathcal{C}}\widehat{p}_{ij}\phi_{x_{i}}-\phi_{x_{j}}\rangle_{\mathcal{H}},$

where the last equality uses the reproducing property of the RKHS. Since
$f_{0}$ is unknown and estimating such a function uniformly on the support of
the covariate vector can be difficult, we rely on a bias-awareness correction
by invoking the Cauchy-Schwarz inequality:

$|\mathcal{B}_{j}|\leq\|f_{0}\|_{\mathcal{H}}\cdot\|\phi_{x_{j}}-\sum_{i\in\mathcal{C}}\widehat{p}_{ij}\phi_{x_{i}}\|_{\mathcal{H}}=\|f_{0}\|_{\mathcal{H}}\cdot\left(K_{tt}+\widehat{P}^{\top}K_{cc}\widehat{P}-2K_{ct}^{\top}\widehat{P}\right)_{jj}^{1/2}.$
(17)

This explains subtracting and adding the term
$\widehat{\theta}\cdot\left(K_{tt}+\widehat{P}^{\top}K_{cc}\widehat{P}-2K_{ct}^{\top}\widehat{P}\right)_{jj}^{1/2}$
in (13) and (14), respectively, where $\widehat{\theta}$ estimates the unknown
quantity $\|f_{0}\|_{\mathcal{H}}$.

#### Variance

Under the model,

$\mathcal{V}_{j}=\widehat{Y}_{j}(0)-\sum_{i\in\mathcal{C}}\widehat{p}_{ij}f_{0}(x_{i})=\sum_{i\in\mathcal{C}}\widehat{p}_{ij}\varepsilon_{0,i}\sim
N(0,\sigma_{0}^{2}\sum_{i\in\mathcal{C}}\widehat{p}_{ij}^{2}).$ (18)

If we knew the bias $\mathcal{B}_{j}$, we could construct the following oracle
confidence interval $[L_{j}^{\ast},U_{j}^{\ast}]$:

$\displaystyle L_{j}^{\ast}$
$\displaystyle:=\widehat{Y}_{j}(0)-\mathcal{B}_{j}-z_{1-\frac{\alpha}{2}}\cdot\sigma_{0}\sqrt{\sum_{i\in\mathcal{C}}\widehat{p}_{ij}^{2}},$
(19) $\displaystyle U_{j}^{\ast}$
$\displaystyle:=\widehat{Y}_{j}(0)-\mathcal{B}_{j}+z_{1-\frac{\alpha}{2}}\cdot\sigma_{0}\sqrt{\sum_{i\in\mathcal{C}}\widehat{p}_{ij}^{2}},$
(20)

Due to (18), the oracle interval $[L_{j}^{\ast},U_{j}^{\ast}]$ is guaranteed
to have the exact $1-\alpha$ coverage, namely, (15) holds. However, since
$\mathcal{B}_{j}$ is unknown, we rely on the aforementioned bias-awareness
correction to construct the confidence interval
$[\widehat{L}_{j},\widehat{U}_{j}]$, obtained by inflating the oracle interval
using (17).

#### Efficiency: Interval Length

We now explain how the width of the confidence interval is related to our
convexified matching objective. As a result, the explanation will shed light
on why the entropic regularization serves as a tuning parameter to trade off
bias and variance for overall interval length efficiency. First, we need the
following simple facts on the Kullback-Leibler divergence, Hellinger distance,
and $\chi^{2}$ distance.

######  Proposition 1.

For any probability vector $p\in\Delta_{n}$,

$\frac{1}{2(n\|p\|_{\infty}+1)}\leq\frac{\sum_{i=1}^{n}p_{i}\log(p_{i})+\log(n)}{n\sum_{i=1}^{n}p_{i}^{2}-1}\leq
1.$

The above fact will help us relate the entropy term and the standard Euclidean
norm of the weights. For a fixed $j\in\mathcal{T}$, let us consider that
$(\widehat{p}_{ij})_{i\in\mathcal{C}}\in\Delta_{N_{c}}$, the $j$-th column of
$\widehat{P}$, is delocalized in the sense
$\max_{i\in\mathcal{C}}\widehat{p}_{ij}\leq(M/2-1)/N_{c}$ with some constant
$M>4$. The interval length is as follows:

$\displaystyle\mathrm{Len}_{j}:=2\cdot\|f_{0}\|_{\mathcal{H}}\cdot\|\phi_{x_{j}}-\sum_{i\in\mathcal{C}}\widehat{p}_{ij}\phi_{x_{i}}\|_{\mathcal{H}}+2\cdot
z_{1-\frac{\alpha}{2}}\cdot\sigma_{0}\sqrt{\sum_{i\in\mathcal{C}}\widehat{p}_{ij}^{2}}.$

From $\widehat{\pi}_{ij}=\widehat{p}_{ij}/N_{t}$ and Proposition 1, we deduce
the following bounds:

$\displaystyle\mathrm{Len}_{j}^{2}$ $\displaystyle\leq
8\left\\{\|f_{0}\|_{\mathcal{H}}^{2}\|\phi_{x_{j}}-\sum_{i\in\mathcal{C}}\tfrac{\pi_{ij}}{1/N_{t}}\phi_{x_{i}}\|_{\mathcal{H}}^{2}+M\cdot
z^{2}_{1-\frac{\alpha}{2}}\sigma_{0}^{2}\frac{N_{t}}{N_{c}}\left[\sum_{i\in\mathcal{C}}\widehat{\pi}_{ij}\log\tfrac{\widehat{\pi}_{ij}}{e}+\frac{\tfrac{1}{M}+1+\log(N_{c}N_{t})}{N_{t}}\right]\right\\},$
$\displaystyle\mathrm{Len}_{j}^{2}$ $\displaystyle\geq
4\left\\{\|f_{0}\|_{\mathcal{H}}^{2}\|\phi_{x_{j}}-\sum_{i\in\mathcal{C}}\tfrac{\pi_{ij}}{1/N_{t}}\phi_{x_{i}}\|_{\mathcal{H}}^{2}+z^{2}_{1-\frac{\alpha}{2}}\sigma_{0}^{2}\frac{N_{t}}{N_{c}}\left[\sum_{i\in\mathcal{C}}\widehat{\pi}_{ij}\log\tfrac{\widehat{\pi}_{ij}}{e}+\frac{2+\log(N_{c}N_{t})}{N_{t}}\right]\right\\}.$

Averaging over $j\in\mathcal{T}$, and barring the constant $M$, the upper and
lower bounds on the right-hand sides remind us of the entropic regularized
convexified matching objective in (7), with the regularization parameter

$\lambda\asymp\frac{z^{2}_{1-\frac{\alpha}{2}}}{N_{c}}\frac{\sigma_{0}^{2}}{\|f_{0}\|_{\mathcal{H}}^{2}}.$

The above tradeoff confirms the empirical findings: for a high signal-to-noise
problem, to balance bias and variance, we need a smaller $\lambda$ for better
approximation error; for a low signal-to-noise problem, a larger $\lambda$ is
preferred to reduce variance. Such a phenomenon confirms the observations in
Figure 2. Searching for a coupling minimizing the objective of convexified
matching is closely related to minimizing the interval length, provided the
regularization parameter $\lambda$ is appropriately chosen. Therefore,
optimizing the entropic regularized convexified matching program (7) with the
desired $\lambda$ chosen above is closely connected to optimizing the average
of the squared lengths of the individual confidence intervals
$\frac{1}{N_{t}}\sum_{j\in\mathcal{T}}\mathrm{Len}_{j}^{2}$.

##  4 Optimization: Algorithms and Analysis

This section introduces and analyzes optimization algorithms to solve the
convexified matching problem, KSC defined in (7) with uniform weights
$v=\frac{1}{N_{t}}$ on the treated, as shown in (12). Notice that it is an
instance of the constrained nonlinear optimization problem: for
$n,m\in\mathbb{N}$, $a\in\Delta_{n}^{+}$, $b\in\Delta_{m}^{+}$, $\lambda>0$,
and $g\colon\mathbb{R}^{n\times m}\to\mathbb{R}$,

$\min_{\pi\in\Pi_{a,b}}~{}g(\pi)+\lambda h(\pi),$ (21)

where $\Pi_{a,b}=\\{\pi\in\mathbb{R}_{+}^{n\times m}:\pi
1_{m}=a~{}\text{and}~{}\pi^{\top}1_{n}=b\\}$ and
$h\colon\mathbb{R}_{+}^{n\times m}\to\mathbb{R}$ is the entropy function
defined by
$h(\pi)=\sum_{i=1}^{n}\sum_{j=1}^{m}\pi_{ij}\log\frac{\pi_{ij}}{e}$. The
convexified matching amounts to the case where $g$ is a convex quadratic
function, in which (21) is a convex program taking the form (12). It can be
solved using the interior point method in $\mathrm{poly}(nm)$ time, possibly
slow in practice for coupling matrices with $nm$ large.

In this section, we devise a faster algorithm that requires an almost
dimension-free number of oracle calls to a subroutine—a simple iterative
matrix scaling procedure called the Sinkhorn algorithm [Sin67]. When $\lambda$
is sufficiently large, the algorithm compares favorably to the interior point
method, requiring only a $\log(1/\varepsilon)$ number of Sinkhorn subroutines.
Later, we modify and extend the algorithm to cover all
$\lambda\in\mathbb{R}_{+}$ by connecting to the steepest descent method under
the Kullback-Leibler divergence. The convergence is admittedly slower for
small $\lambda$: $\log(nm)/\varepsilon$ number of Sinkhorn subroutines is
required, but it still compares favorably to the interior point method.

#### Additional notation

For a matrix $A\in\mathbb{R}^{n\times m}$, let $\|A\|_{\mathrm{F}}$ denote its
Frobenius norm, $\|A\|_{1}:=\sum_{i=1}^{n}\sum_{j=1}^{m}|A_{ij}|$ denote its
$L^{1}$ norm, and $\|A\|_{\infty}:=\max_{1\leq i\leq n}\max_{1\leq j\leq
m}|A_{ij}|$ denote its $L^{\infty}$ norm. Let $h$ keep its definition as the
entropy function, and let $\nabla h(A)$ denote the gradient of $h$, which is
$\log(A)$, the entrywise logarithm of $A$, provided all entries of $A$ are
positive. For any $a\in\Delta_{n}$ and $b\in\Delta_{m}$, let
$\Pi_{a,b}^{+}:=\\{A\in\Pi_{a,b}:A_{ij}>0~{}~{}\forall i,j\\}$. Let
$\Delta_{n,m}=\\{A\in\mathbb{R}_{+}^{n\times
m}:\sum_{i=1}^{n}\sum_{j=1}^{m}A_{ij}=1\\}$. For vectors
$u,v\in\mathbb{R}^{n}$, let $\frac{u}{v}$ denote the entrywise division
provided all the entries of $v$ are nonzero.

###  4.1 Optimality Conditions and the Sinkhorn Algorithm

When $g$ is convex and $\lambda>0$, the strong convexity of $h$ on $\Pi_{a,b}$
implies that (21) admits a unique minimizer. Moreover, as $h$ prevents the
minimizer from having a zero entry, (21) admits a unique minimizer on
$\Pi_{a,b}^{+}$, which is the interior of $\Pi_{a,b}$. It turns out that this
unique minimizer of (21) is the fixed point of an operator related to the
entropic regularized optimal transport problem [Cut13], as we shall show in
Proposition 2. Later, in Theorem 1, we derive convergence to the minimizer via
the iterative matrix scaling subroutine. This subroutine is referred to as the
Sinkhorn algorithm [Sin67], which we introduce below.

######  Definition 1 (Sinkhorn [Sin67, Cut13]).

Fix $\lambda>0$. The following operator
$\Phi_{\lambda}\colon\mathbb{R}^{n\times m}\to\Pi_{a,b}^{+}$ is well-defined:

$\Phi_{\lambda}(C):=\operatorname*{arg\,min}_{\pi\in\Pi_{a,b}}\left(\langle
C,\pi\rangle+\lambda h(\pi)\right)\quad\forall C\in\mathbb{R}^{n\times m}.$

The right-hand side is often called the entropic regularized optimal transport
problem, given a cost matrix $C$. For any $C\in\mathbb{R}^{n\times m}$,
suppose there are $\mu\in\mathbb{R}^{n}$ and $\nu\in\mathbb{R}^{m}$ such that

$\exp\left(-\frac{\mu
1_{m}^{\top}+1_{n}\nu^{\top}+C}{\lambda}\right)\in\Pi_{a,b},$

where $\exp$ is applied entrywise. Then, the following must hold:

$\exp\left(-\frac{\mu
1_{m}^{\top}+1_{n}\nu^{\top}+C}{\lambda}\right)=\Phi_{\lambda}(C).$

Moreover, for any $C\in\mathbb{R}^{n\times m}$, one can obtain
$\Phi_{\lambda}(C)$ by iteratively scaling the rows and columns of the matrix
$\exp(-C/\lambda)$ using the Sinkhorn algorithm summarized in Algorithm 1,
namely, $\Phi_{\lambda}(C)$ is the limit of the sequence produced from
$\mathrm{Sinkhorn}(e^{-C/\lambda},a,b)$.

Algorithm 1 $\mathrm{Sinkhorn}(P,a,b)$

0: A matrix $P\in\mathbb{R}^{n\times m}$ with positive entries,
$a\in\Delta_{n}^{+}$, and $b\in\Delta_{m}^{+}$.

1: Initialize $k\leftarrow 0$, $x^{(k)}\leftarrow 1_{n}$, and
$y^{(k)}\leftarrow\frac{b}{P^{\top}x^{(k)}}$.

2: repeat

3: $x^{(k+1)}\leftarrow\frac{a}{Py^{(k)}}\quad\text{and}\quad
y^{(k+1)}\leftarrow\frac{b}{P^{\top}x^{(k+1)}}.$

4: $P^{(k+1)}=\mathrm{diag}(x^{(k+1)})P\mathrm{diag}(y^{(k+1)})$.

5: $k\leftarrow k+1$

6: until Discrepancy between $P^{(k)}1_{m}$ and $a$ reaches a desired level.

7: return $P^{(k)}$.

For simplicity, Algorithm 1 states the Sinkhorn algorithm such that the column
sum of $P^{(k)}$ matches with $b$, namely, $(P^{(k)})^{\top}1_{n}=b$ for any
$k\geq 0$. Then, to assess the accuracy of the scaled matrix $P^{(k)}$, we may
only need to check a suitable discrepancy between $a$ and the row sum of
$P^{(k)}$, that is, $P^{(k)}1_{m}$.

The following proposition shows that when $g$ is convex, the unique minimizer
of (21) is the fixed point of an operator given by the composition of
$\Phi_{\lambda}$ and $\nabla g$.

######  Proposition 2.

Fix $\lambda>0$. For a function $g\colon\mathbb{R}^{n\times m}\to\mathbb{R}$
such that $\nabla g$ exists, define an operator
$T_{\lambda}\colon\mathbb{R}^{n\times m}\to\Pi_{a,b}^{+}$ defined by
$T_{\lambda}=\Phi_{\lambda}\circ\nabla g$, that is,

$T_{\lambda}(A)=\operatorname*{arg\,min}_{\pi\in\Pi_{a,b}}(\langle\nabla
g(A),\pi\rangle+\lambda h(\pi))\quad\forall A\in\mathbb{R}^{n\times m}.$ (22)

If $g$ is convex, (21) admits a unique minimizer contained in $\Pi_{a,b}^{+}$,
say, $\pi^{\star}\in\Pi_{a,b}^{+}$, and $\pi^{\star}$ is the unique fixed
point of $T_{\lambda}$, namely, $\pi^{\star}=T_{\lambda}(\pi^{\star})$.

Based on the established connection between the optimality condition of (21)
and the Sinkhorn algorithm, we propose two algorithms to solve (21) in the
following sections.

###  4.2 Fixed-Point Algorithm

By Proposition 2, finding the fixed point of $T_{\lambda}$ is equivalent to
solving (21). Therefore, we propose solving (21) by the fixed-point iterations
based on the operator $T_{\lambda}\colon\mathbb{R}^{n\times
m}\to\Pi_{a,b}^{+}$, which is simply choosing the initial point
$\pi^{(0)}\in\Pi_{a,b}$ properly and iterating $\pi^{(k+1)}\leftarrow
T_{\lambda}(\pi^{(k)})$ for $k\geq 0$, where
$T_{\lambda}(\pi^{(k)})=\Phi_{\lambda}(\nabla g(\pi^{(k)}))$ is approximated
by $\mathrm{Sinkhorn}(e^{-\nabla g(\pi^{(k)})/\lambda},a,b)$ as explained in
Definition 1. Algorithm 2 summarizes this procedure.

Algorithm 2 $\mathrm{FixedPoint}(g,\lambda,a,b)$

0: A map $g\colon\mathbb{R}^{n\times m}\to\mathbb{R}$, a parameter
$\lambda>0$, $a\in\Delta_{n}^{+}$, and $b\in\Delta_{m}^{+}$.

1: Pick any $\pi^{(0)}\in\Pi_{a,b}$ and set $k\leftarrow 0$.

2: repeat

3: $\pi^{(k+1)}\leftarrow\mathrm{Sinkhorn}(e^{-\nabla
g(\pi^{(k)})/\lambda},a,b)$.

4: $k\leftarrow k+1$.

5: until Discrepancy between $\pi^{(k)}$ and $\pi^{(k-1)}$ reaches a desired
level.

6: return $\pi^{(k)}$.

We show that $T_{\lambda}$ is a contraction if the gradient of $g$ is
Lipschitz and $\lambda$ is larger than the Lipschitz constant. In such a case,
Algorithm 2 converges to the fixed point quickly. When $g$ is a quadratic
function, we can write the Lipschitz constant in terms of the $L^{\infty}$
norm of the Hessian matrix, which is independent of the input dimension $nm$.
We do not require $g$ to be convex for this result; $T_{\lambda}$ admits a
unique fixed point due to the contraction argument, which is independent of
Proposition 2. The role of Proposition 2 is to translate this convergence
result into the convergence to the minimizer of (21) for convex $g$ by
equating the fixed point to the minimizer.

######  Theorem 1.

Let $g\colon\mathbb{R}^{n\times m}\to\mathbb{R}$ be a quadratic function, with
some $H\in\mathbb{R}^{n\times n}$ that is symmetric and
$C\in\mathbb{R}^{n\times m}$,

$g(\pi)=\frac{1}{2}\langle\pi,H\pi\rangle+\langle
C,\pi\rangle\quad\forall\pi\in\mathbb{R}^{n\times m}.$

Assume $\lambda>\|H\|_{\infty}$. Then, the operator $T_{\lambda}$ defined in
(22) is a contraction under the distance defined by $\|\cdot\|_{1}$ and has a
unique fixed-point, say, $\pi^{\star}\in\Pi_{a,b}^{+}$. Moreover, Algorithm 2,
assuming the inner loop $\mathrm{Sinkhorn}$ is always exact, outputs a
sequence $(\pi^{(k)})_{k\geq 0}$ such that for any $T\in\mathbb{N}$,

$\|\pi^{(T)}-\pi^{\star}\|_{1}\leq\frac{(\|H\|_{\infty}/\lambda)^{T}}{1-(\|H\|_{\infty}/\lambda)}\|\pi^{(1)}-\pi^{(0)}\|_{1}.$
(23)

######  Remark 2.

Theorem 1 shows that when $\lambda>\|H\|_{\infty}$, to obtain an
$\varepsilon$-approximate fixed point
$\|\pi^{(T)}-\pi^{\star}\|_{1}^{2}\leq\varepsilon$ , we need
$\log(1/\varepsilon)/2\log(\lambda/\|H\|_{\infty})$ number of Sinkhorn
routines—a dimension-free quantity.

###  4.3 Local Algorithm with Kullback-Leibler Geometry

The fixed-point algorithm may not converge for sufficiently small $\lambda$.
To fill in the gap when $\lambda\in[0,\|H\|_{\infty})$ left unanswered by
Theorem 1, we employ the steepest descent method under Bregman divergence to
solve (21) as follows: let $f=g+\lambda h$,

$\pi^{(k+1)}=\operatorname*{arg\,min}_{\pi\in\Pi_{a,b}}\left(\langle\nabla
f(\pi^{(k)}),\pi\rangle+\frac{D_{h}(\pi,\pi^{(k)})}{\tau_{k}}\right),$ (24)

where $D_{h}(\pi,\pi^{(k)})=h(\pi)-h(\pi^{(k)})-\langle\nabla
h(\pi^{(k)}),\pi-\pi^{(k)}\rangle$ is the Bregman divergence, which equals the
Kullback-Leibler divergence on $\Pi_{a,b}$, and $\tau_{k}>0$ is a suitable
step size. With the entropic regularization $h$, we implement the steepest
descent over the polytope $\Pi_{a,b}$ under the Kullback-Leibler divergence.
We emphasize that this can be solved using the Sinkhorn algorithm as well
since (24) is equivalent to

$\begin{split}\pi^{(k+1)}&=\operatorname*{arg\,min}_{\pi\in\Pi_{a,b}}\left(\langle\nabla
g(\pi^{(k)})+(\lambda-\tau_{k}^{-1})\nabla
h(\pi^{(k)}),\pi\rangle+\tau_{k}^{-1}h(\pi)\right)\\\
&=\Phi_{\tau_{k}^{-1}}\left(\nabla g(\pi^{(k)})+(\lambda-\tau_{k}^{-1})\nabla
h(\pi^{(k)})\right).\end{split}$

The above can be approximated by applying the Sinkhorn algorithm to $\nabla
g(\pi^{(k)})+(\lambda-\tau_{k}^{-1})\nabla h(\pi^{(k)})$. Algorithm 3
summarizes this procedure. Note that letting $\tau_{k}=\lambda^{-1}$ for all
$k\geq 0$ in Algorithm 3 recovers Algorithm 2.

Algorithm 3 $\mathrm{SteepestDescentKL}(g,\lambda,\\{\tau_{k}\\}_{k\geq
0},a,b)$

0: A map $g\colon\mathbb{R}^{n\times m}\to\mathbb{R}$, a parameter
$\lambda>0$, $a\in\Delta_{n}^{+}$, and $b\in\Delta_{m}^{+}$.

0: Step sizes $\\{\tau_{k}\\}_{k\geq 0}$.

1: Pick any $\pi^{(0)}\in\Pi_{a,b}$ and set $k\leftarrow 0$.

2: repeat

3: $\pi^{(k+1)}\leftarrow\mathrm{Sinkhorn}(e^{-(\nabla
g(\pi^{(k)})+(\lambda-\tau_{k}^{-1})\nabla h(\pi^{(k)}))/\tau_{k}^{-1}},a,b)$.

4: $k\leftarrow k+1$.

5: until Discrepancy between $\pi^{(k)}$ and $\pi^{(k-1)}$ reaches a desired
level.

6: return $\pi^{(k)}$.

We show that Algorithm 3 outputs a sequence that converges to the minimum of
(21) provided the step size is below a certain threshold.

######  Theorem 2.

Let $g\colon\mathbb{R}^{n\times m}\to\mathbb{R}$ be a convex quadratic
function, for some $H\in\mathbb{R}^{n\times n}$ that is symmetric and positive
semidefinite and $C\in\mathbb{R}^{n\times m}$,

$g(\pi)=\frac{1}{2}\langle\pi,H\pi\rangle+\langle
C,\pi\rangle\quad\forall\pi\in\mathbb{R}^{n\times m}.$

Then, for any $\lambda\geq 0$, if Algorithm 3, assuming the inner loop
$\mathrm{Sinkhorn}$ is always exact, with constant stepsize $\tau_{k}=\tau$
for all $k\geq 0$, where $\tau^{-1}\geq\|H\|_{\infty}+\lambda$, it outputs a
sequence $(\pi^{(k)})_{k\geq 0}$ such that for any $T\in\mathbb{N}$,

$g(\pi^{(T)})+\lambda h(\pi^{(T)})-\min_{\pi\in\Pi_{a,b}}(g(\pi)+\lambda
h(\pi))\leq\frac{1}{T}\frac{\log(nm)}{\tau}.$ (25)

For $\lambda>0$, the following holds as well:

$\|\pi^{(T)}-\pi^{\star}\|_{1}^{2}\leq\frac{1}{T}\frac{2\log(nm)}{\lambda\tau},$
(26)

where $\pi^{\star}\in\Pi_{a,b}^{+}$ is the unique minimizer of (21).

######  Remark 3.

The $\lambda\in[\|H\|_{\infty},\infty)$ case has been studied in Theorem 1.
Theorem 2 shows that for any $\lambda\in[0,\|H\|_{\infty})$, with the choice
$\tau^{-1}=2\|H\|_{\infty}$, we need $1/\varepsilon\cdot
2\|H\|_{\infty}\log(nm)$ number of Sinkhorn subroutines to get
$\varepsilon$-close to the objective value. Further more, when $\lambda$ is
strictly positive, we can assure that after $1/\varepsilon\cdot
4\|H\|_{\infty}\log(nm)/\lambda$ number of Sinkhorn subroutines,
$\|\pi^{(T)}-\pi^{\star}\|_{1}^{2}\leq\varepsilon$, with a logarithmic
dependence on the problem dimension.

##  5 Applications: Revisit the NSW Data

In this section, we numerically investigate our imputation and individualized
inference method using the NSW data and compare it with other methods,
continued from Section 2.3.

###  5.1 Imputation

We conduct a robustness check of the imputation results in Section 2.3 on the
NSW data by specifying different nonlinear kernels. We consider the RBF kernel
and the polynomial kernel of degree two. The results are shown in Figure 4 and
Figure 5. Though the overall patterns are similar to those using the linear
kernel, we can observe some differences. First, for the experimental data, the
imputed counterfactual outcomes via the RBF kernel with $\lambda=0.001$ show
the largest variability. In Figure 4(a), there are two treated units whose
imputed counterfactual outcomes are above 15000, which do not occur in the
results via the linear or polynomial kernel; as a result, the right tail of
the histogram along the $x$-axis is longer than the other cases. For the PSID
data, we can see that using the polynomial kernel leads to the most spread out
imputed outcomes, visualized by the right tails of the histograms along the
$x$-axis in (a) and (b) of Figure 5 which are longer than the other cases. We
point out that unlike the NN-matching method and our KSC method, whose
counterfactual distributions are supported on $\mathbb{R}_{+}$ even for the
nonexperimental PSID data, the regression imputation method presents a wider
support for imputed values, even with negative outcomes as seen across Figures
1, 4 and 5.

(a) NSW $\lambda=0.001$ (b) NSW $\lambda=0.01$

(c) PSID $\lambda=0.001$ (d) PSID $\lambda=0.01$

Figure 4: Scatter plots of the imputed counterfactual outcomes of the treated
along with the histograms of the marginal distributions. The setup is the same
as in Figure 1, but the results are based on the RBF kernel for both NSW and
PSID data.

(a) NSW $\lambda=0.001$ (b) NSW $\lambda=0.01$

(c) PSID $\lambda=0.001$ (d) PSID $\lambda=0.01$

Figure 5: Scatter plots of the imputed counterfactual outcomes of the treated
along with the histograms of the marginal distributions. The setup is the same
as in Figure 1, but the results are based on the polynomial kernel of degree
two for both NSW and PSID data.

###  5.2 Inference

We apply the method presented in Section 3 to construct confidence intervals
around the imputed counterfactual outcomes. Figure 6 (a) shows the results for
the experimental data using the linear kernel with $\lambda\in\\{0.001,1\\}$,
which plots the estimated individual treatment effects (ITEs),
$\widehat{\tau}_{j}=Y_{j}-\widehat{Y}_{j}(0)$, against the logarithm of
earnings in 1975, where the confidence intervals are shown as error bars
around the estimated ITEs. Here, only 74 treated units (out of 185 treated
units) with positive earnings in 1975 are shown. We pick this index merely for
visualization of individualization. For comparison, we also show the blue
points representing the ITEs based on imputing all the treated units with the
mean of the control group outcomes, along with the confidence intervals based
on the standard error of the mean; the case where there is no
individualization in the imputed control outcomes, and thus no meaningful
individual confidence intervals addressing the bias. Let us first focus on the
point estimates of the ITEs. We can see that most of the estimated ITEs fall
in the interval $(-10000,10000)$, seemingly uncorrelated with the earnings in
1975, while there are seven individuals whose ITEs are above 10000. These
seven units with the largest ITEs are likely to be the ones who benefit most
from the treatment, standing out from the rest of the treatment group. The
confidence intervals help us decide how sure we are about these seemingly
significant ITEs. We can confidently say that the top two ITEs are extremely
large compared to the rest of the ITEs since the lower ends of their
individual confidence intervals are larger than the upper ends of the rest of
the ITEs, regardless of the choice of $\lambda$. Now, looking at the
confidence intervals, rather than only looking at the point estimates, the
four units whose ITEs are between 10000 and 30000 are less distinguishable
from the rest. The error bars do help to suggest which units benefit the more
here. For the third largest ITE, around 20000, the conclusion may vary
depending on $\lambda$. If we read the confidence intervals corresponding to
$\lambda=1$ (red, dashed), this unit seems to be more significant than the
rest, while this conclusion is less assertive for $\lambda=0.001$ as the
confidence interval of this unit overlaps with the interval $(-15000,20000)$
which contains the rest of the confidence intervals. In other words, using
$\lambda=0.001$ leads to the most conservative conclusion regarding who
benefits the most from the treatment as the confidence intervals are wider.
Unlike the confidence intervals obtained by the proposed method, the intervals
based on the mean imputation (shown in blue) are extremely narrow, which does
not take into account individual heterogeneity, overlooking a potentially
large bias; the confidence interval here is useless for individual decision
making and may result in overly optimistic conclusions, for instance,
essentially the majority benefits significantly from the program. (b) and (c)
of Figure 6 repeat the same analysis using the RBF kernel and the polynomial
kernel of degree two, respectively. The overall patterns are similar to those
using the linear kernel, but the confidence intervals tend to be wider.
Accordingly, even the second largest ITE above 30000 is not as distinguishable
from the rest as in the linear kernel case.

The results based on the PSID data are shown in Figure 7. One notable
difference is that the intervals are generally much wider than those of the
NSW experimental data. One possible explanation is that the PSID data, even
after trimmed, differs significantly from the experimental data, leading to
the matching with more conservative error estimates. Accordingly, the
confidence intervals are wider, which makes it harder to make a decisive
conclusion with confidence.

Lastly, we briefly comment on the choice of $\rho$, the regularization
parameter of the kernel ridge regression mentioned in Section 3.2. Increasing
$\rho$ leads to a smaller estimate of $\|f_{0}\|_{\mathcal{H}}$, leading to a
smaller bias correction. Increasing $\rho$ will typically result in a larger
estimate of $\sigma_{0}$, implying a large variability. The correct choice of
$\rho$ depends on the underlying signal-to-noise ratio for the data-generating
process in the control outcomes, namely, the complexity of the function vs.
the amount of noise. In all examples of this section, we choose $\rho$ using
5-fold cross-validation.

(a) Linear (b) RBF

(c) Polynomial

Figure 6: NSW experimental data: estimated individual treatment effects (ITEs)
of the treated together with the confidence intervals along the earnings
before the treatment. Confidence intervals are shown as error bars around the
estimated ITEs. 74 treated units (out of 185 treated units) with positive
earnings in 1975 are shown. The $x$-axis is the logarithm of earnings in 1975,
while the $y$-axis is the estimated ITEs. The results are based on the
experimental data using uniform weights $v,w$ in (7) with three different
kernels: linear, RBF, and polynomial of degree two. The blue points correspond
to imputing all the treated units with the mean of the control group outcomes,
where the error bars show the confidence interval based on the standard error
of the mean.

(a) Linear (b) RBF

(c) Polynomial

Figure 7: PSID nonexperimental data: estimated individual treatment effects
(ITEs) of the treated along the earnings before the treatment. The setup is
the same as in Figure 6, but the results are based on the trimmed PSID data.
For the weights $v,w$ in (7), we let $v$ be uniform, while $w$ is the
propensity score-based weights as in (8).

##  6 Conclusion

We proposed a convexified matching method for missing value imputation and
individualized inference, integrating favorable features from optimal
matching, regression imputation, and synthetic control. We impute
counterfactual outcomes based on convex combinations of observed outcomes,
defined by an optimal coupling between the treated and control data sets.
Finding an optimal coupling is a convex relaxation to the combinatorial
optimal matching problem, for which we propose efficient algorithms based on
matrix scaling. Unlike existing imputation methods, we begin with a desirable
aggregate-level summary and estimate granular-level individual treatment
effects by properly constraining the coupling so that the estimated individual
effects are consistent with the aggregate summary. We provided a method to
construct individual confidence intervals for the estimated counterfactual
outcomes, along with a simulation study to demonstrate the effectiveness of
our method. The simulation confirms our theoretical insight on the level of
entropic regularization needed. We establish that entropic regularization
plays a crucial role in both efficiency in inference and computation: first,
the convexified matching objective is gauged for minimizing the width of the
individual confidence intervals, trading off bias and variance; second, the
entropic regularization enables us to design fast algorithms. We demonstrated
the empirical performance of our method on the NSW data, using both
experimental and nonexperimental data sets, with a comparison to existing
methods and robustness checks.

## References

  * [Aba21] Alberto Abadie “Using Synthetic Controls: Feasibility, Data requirements, and Methodological Aspects” In _Journal of Economic Literature_ 59.2, 2021, pp. 391–425
  * [ADH10] Alberto Abadie, Alexis Diamond and Jens Hainmueller “Synthetic control methods for comparative case studies: Estimating the effect of California’s tobacco control program” In _Journal of the American statistical Association_ 105.490, 2010, pp. 493–505
  * [AG03] Alberto Abadie and Javier Gardeazabal “The economic costs of conflict: A case study of the Basque Country” In _American Economic Review_ 93.1, 2003, pp. 113–132
  * [AI06] Alberto Abadie and Guido W. Imbens “Large Sample Properties of Matching Estimators for Average Treatment Effects” In _Econometrica_ 74.1, 2006, pp. 235–267
  * [AI11] Alberto Abadie and Guido W. Imbens “Bias-Corrected Matching Estimators for Average Treatment Effects” In _Journal of Business & Economic Statistics_ 29.1, 2011, pp. 1–11
  * [AL21] Alberto Abadie and Jérémy L’hour “A penalized synthetic control estimator for disaggregated data” In _Journal of the American Statistical Association_ 116.536, 2021, pp. 1817–1834
  * [AI16] Susan Athey and Guido Imbens “Recursive partitioning for heterogeneous causal effects” In _Proceedings of the National Academy of Sciences_ 113.27, 2016, pp. 7353–7360
  * [BCH14] Alexandre Belloni, Victor Chernozhukov and Christian Hansen “Inference on treatment effects after selection among high-dimensional controls” In _Review of Economic Studies_ 81.2 Oxford University Press, 2014, pp. 608–650
  * [Cat10] Matias D Cattaneo “Efficient semiparametric estimation of multi-valued treatment effects under ignorability” In _Journal of Econometrics_ 155.2 Elsevier, 2010, pp. 138–154
  * [Che+17] Victor Chernozhukov et al. “Double/debiased/neyman machine learning of treatment effects” In _American Economic Review_ 107.5 American Economic Association 2014 Broadway, Suite 305, Nashville, TN 37203, 2017, pp. 261–265
  * [Che+22] Victor Chernozhukov et al. “Locally robust semiparametric estimation” In _Econometrica_ 90.4 Wiley Online Library, 2022, pp. 1501–1535
  * [Cut13] Marco Cuturi “Sinkhorn distances: lightspeed computation of optimal transport” In _Advances in Neural Information Processing Systems_ 26, 2013
  * [DW99] Rajeev H Dehejia and Sadek Wahba “Causal effects in nonexperimental studies: Reevaluating the evaluation of training programs” In _Journal of the American statistical Association_ 94.448 Taylor & Francis, 1999, pp. 1053–1062
  * [DW02] Rajeev H Dehejia and Sadek Wahba “Propensity score-matching methods for nonexperimental causal studies” In _Review of Economics and Statistics_ 84.1, 2002, pp. 151–161
  * [Far15] Max H Farrell “Robust inference on average treatment effects with possibly more covariates than observations” In _Journal of Econometrics_ 189.1 Elsevier, 2015, pp. 1–23
  * [FLM21] Max H Farrell, Tengyuan Liang and Sanjog Misra “Deep neural networks for estimation and inference” In _Econometrica_ 89.1 Wiley Online Library, 2021, pp. 181–213
  * [Hah98] Jinyong Hahn “On the role of the propensity score in efficient semiparametric estimation of average treatment effects” In _Econometrica_ JSTOR, 1998, pp. 315–331
  * [Hai12] Jens Hainmueller “Entropy balancing for causal effects: A multivariate reweighting method to produce balanced samples in observational studies” In _Political Analysis_ 20.1 Cambridge University Press, 2012, pp. 25–46
  * [HIT97] James J Heckman, Hidehiko Ichimura and Petra E Todd “Matching as an econometric evaluation estimator: Evidence from evaluating a job training programme” In _Review of Economic Studies_ 64.4 Wiley-Blackwell, 1997, pp. 605–654
  * [HV05] James J Heckman and Edward Vytlacil “Structural equations, treatment effects, and econometric policy evaluation” In _Econometrica_ 73.3 Wiley Online Library, 2005, pp. 669–738
  * [Hil55] Austin Bradford Hill “Principles of Medical Statistics” Oxford University Press, 1955
  * [HI01] Keisuke Hirano and Guido W Imbens “Estimation of causal effects using propensity score weighting: An application to data on right heart catheterization” In _Health Services and Outcomes Research Methodology_ 2, 2001, pp. 259–278
  * [HIR03] Keisuke Hirano, Guido W Imbens and Geert Ridder “Efficient estimation of average treatment effects using the estimated propensity score” In _Econometrica_ 71.4, 2003, pp. 1161–1189
  * [LaL86] Robert J LaLonde “Evaluating the econometric evaluations of training programs with experimental data” In _American Economic Review_ JSTOR, 1986, pp. 604–620
  * [LMZ18] Fan Li, Kari Lock Morgan and Alan M Zaslavsky “Balancing covariates via propensity score weighting” In _Journal of the American Statistical Association_ 113.521, 2018, pp. 390–400
  * [LR23] Tengyuan Liang and Benjamin Recht “Randomization inference when n equals one” In _arXiv preprint arXiv:2310.16989_ , 2023
  * [LM16] Keli Liu and Xiao-Li Meng “There is individualized treatment. Why not individualized inference?” In _Annual Review of Statistics and Its Application_ 3, 2016, pp. 79–111
  * [MCH11] Thomas MaCurdy, Xiaohong Chen and Han Hong “Flexible estimation of treatment effect parameters” In _American Economic Review_ 101.3 American Economic Association, 2011, pp. 544–551
  * [Ney23] Jerzy Neyman “On the application of probability theory to agricultural experiments: Essay on principles, Section 9”, In Polish, but translated in English with discussion in Statistical Science, 5(4):465–472, 1990, 1923
  * [PC19] Gabriel Peyré and Marco Cuturi “Computational optimal transport: With applications to data science” In _Foundations and Trends® in Machine Learning_ 11.5-6, 2019, pp. 355–607
  * [Ros87] Paul R Rosenbaum “Model-based direct adjustment” In _Journal of the American statistical Association_ 82.398, 1987, pp. 387–394
  * [RR83] Paul R Rosenbaum and Donald B Rubin “The central role of the propensity score in observational studies for causal effects” In _Biometrika_ 70.1 Oxford University Press, 1983, pp. 41–55
  * [Ros20] Paul R. Rosenbaum “Modern Algorithms for Matching in Observational Studies” In _Annual Review of Statistics and Its Application_ 7, 2020, pp. 143–176
  * [Rub74] Donald B Rubin “Estimating Causal Effects of Treatments in Randomized and Nonrandomized Studies” In _Journal of Educational Psychology_ 66.5, 1974, pp. 688
  * [Rub06] Donald B. Rubin “Matched Sampling for Causal Effects” Cambridge University Press, 2006
  * [SC04] John Shawe-Taylor and Nello Cristianini “Kernel Methods for Pattern Analysis” Cambridge University Press, 2004
  * [Sin67] Richard Sinkhorn “Diagonal equivalence to matrices with prescribed row and column sums” In _The American Mathematical Monthly_ 74.4, 1967, pp. 402–405
  * [SC08] Ingo Steinwart and Andreas Christmann “Support Vector Machines” Springer, 2008
  * [Tsy09] Alexandre B. Tsybakov “Introduction to Nonparametric Estimation” Springer, 2009
  * [VR03] Mark J Van der Laan and James M Robins “Unified Methods for Censored Longitudinal Data and Causality” Springer, 2003
  * [Vil03] Cédric Villani “Topics in Optimal Transportation” American Mathematical Society, 2003
  * [WA18] Stefan Wager and Susan Athey “Estimation and inference of heterogeneous treatment effects using random forests” In _Journal of the American Statistical Association_ 113.523, 2018, pp. 1228–1242
  * [Zub+23] José R Zubizarreta, Elizabeth A Stuart, Dylan S Small and Paul R Rosenbaum “Handbook of Matching and Weighting Adjustments for Causal Inference” CRC Press, 2023

##  Appendix A Proofs

###  A.1 Proof of Proposition 1

###### Proof of Proposition 1.

For any $p=(p_{1},\ldots,p_{n})\in\Delta_{n}$, we have

$\sum_{i=1}^{n}p_{i}\log(p_{i})+\log(n)\leq n\sum_{i=1}^{n}p_{i}^{2}-1,$

where the left-hand side is the Kullback-Leibler (KL) divergence of $p$ from
$\frac{1}{n}1_{n}$ which is known to be bounded by the $\chi^{2}$ distance of
$p$ from $\frac{1}{n}1_{n}$ on the right-hand side; see Lemma 2.7 of [Tsy09].
Meanwhile, the KL divergence is bounded below by the Hellinger distance; see
Lemma 2.4 of [Tsy09]. Therefore, we have

$\sum_{i=1}^{n}p_{i}\log(p_{i})+\log(n)\geq\sum_{i=1}^{n}(\sqrt{p_{i}}-\sqrt{1/n})^{2},$

where the right-hand side is the Hellinger distance between $p$ and
$\frac{1}{n}1_{n}$. The right-hand side can be further lower bounded,

$\sum_{i=1}^{n}(\sqrt{p_{i}}-\sqrt{1/n})^{2}=\sum_{i=1}^{n}\frac{(p_{i}-1/n)^{2}}{(\sqrt{p_{i}}+\sqrt{1/n})^{2}}\geq\frac{n\sum_{i=1}^{n}p_{i}^{2}-1}{2(n\|p\|_{\infty}+1)}.$

Combining these inequalities, we have

$\frac{1}{2(n\|p\|_{\infty}+1)}\leq\frac{\sum_{i=1}^{n}p_{i}\log(p_{i})+\log(n)}{n\sum_{i=1}^{n}p_{i}^{2}-1}\leq
1.$

∎

###  A.2 Proof of Proposition 2

###### Proof of Proposition 2.

We have already discussed that (21) must admit a unique minimizer
$\pi^{\star}\in\Pi_{a,b}^{+}$. We show that $\pi^{\star}$ is a unique fixed
point of $T_{\lambda}$, which follows from the Karush-Kuhn-Tucker (KKT)
conditions. To see this, rewrite (21) as

$\displaystyle\min_{\pi\in\mathbb{R}_{+}^{n\times m}}$
$\displaystyle(g(\pi)+\lambda h(\pi))$ $\displaystyle\mathrm{subject~{}to}$
$\displaystyle\pi 1_{m}=a\quad\text{and}\quad\pi^{\top}1_{n}=b.$

Then, define the Lagrangian $L\colon\mathbb{R}_{+}^{n\times
m}\times\mathbb{R}^{n}\times\mathbb{R}^{m}\to\mathbb{R}$ by

$L(\pi,\mu,\nu)=g(\pi)+\lambda h(\pi)+\langle\mu,\pi
1_{m}-a\rangle+\langle\nu,\pi^{\top}1_{n}-b\rangle.$

The KKT conditions for $(\pi,\mu,\nu)\in\mathbb{R}_{+}^{n\times
m}\times\mathbb{R}^{n}\times\mathbb{R}^{m}$ are $\pi\in\Pi_{a,b}$ and $\pi$ is
a minimizer of $L(\cdot,\mu,\nu)$ over $\mathbb{R}_{+}^{n\times m}$. Due to
$h$, we can see that $L(\cdot,\mu,\nu)$ must admit a unique minimizer at the
interior of $\mathbb{R}_{+}^{n\times m}$, namely, all the entries of the
unique minimizer are strictly positive. Accordingly, the minimizer must
satisfy the first-order condition:

$\nabla_{\pi}L(\pi,\mu,\nu)=\nabla g(\pi)+\lambda\nabla h(\pi)+\mu
1_{m}^{\top}+1_{n}\nu^{\top}=0,$

which leads to

$\pi=\exp\left(-\frac{\mu 1_{m}^{\top}+1_{n}\nu^{\top}+\nabla
g(\pi)}{\lambda}\right).$ (27)

In summary, the KKT conditions are $\pi\in\Pi_{a,b}$ and (27), which—by
Definition 1—leads to the following:

$\pi=\Phi_{\lambda}(\nabla g(\pi))=T_{\lambda}(\pi).$

Hence, we have shown that $\pi$ is a minimizer of (21) if and only if $\pi$ is
a fixed point of $T_{\lambda}$. As we have already proved that (21) admits a
unique minimizer $\pi^{\star}$, we conclude that $T_{\lambda}$ has a unique
fixed point $\pi^{\star}$. ∎

###  A.3 Proof of Theorem 1

Theorem 1 is based on the following lemma that analyzes the Lipschitz property
of the operator $\Phi_{\lambda}$—defined in Proposition 1—w.r.t. the
$L^{\infty}$ metric.

######  Lemma 1.

Fix $\lambda>0$. Define a map $\Phi_{\lambda}\colon\mathbb{R}^{n\times
m}\to\Pi_{a,b}^{+}$ by

$\Phi_{\lambda}(C)=\operatorname*{arg\,min}_{\pi\in\Pi_{a,b}}\left(\langle
C,\pi\rangle+\lambda h(\pi)\right).$

Then, for any $C_{1},C_{2}\in\mathbb{R}^{n\times m}$, we have

$\|\Phi_{\lambda}(C_{1})-\Phi_{\lambda}(C_{2})\|_{1}\leq\frac{\|C_{1}-C_{2}\|_{\infty}}{\lambda}.$
(28)

###### Proof of Lemma 1.

We first claim that

$\langle C_{1}+\lambda\nabla
h(\Phi_{\lambda}(C_{1})),\pi-\Phi_{\lambda}(C_{1})\rangle\geq
0\quad\forall\pi\in\Pi_{a,b}.$ (29)

Suppose not, namely, there exists $\pi\in\Pi_{a,b}$ such that

$\langle C_{1}+\lambda\nabla
h(\Phi_{\lambda}(C_{1})),\pi-\Phi_{\lambda}(C_{1})\rangle<0.$

In other words, letting $q(\pi):=\langle C_{1},\pi\rangle+\lambda h(\pi)$, we
have $\pi\in\Pi_{a,b}$ such that

$\langle\nabla q(\Phi_{\lambda}(C_{1})),\pi-\Phi_{\lambda}(C_{1})\rangle<0.$

As $q$ is differentiable at $\Phi_{\lambda}(C_{1})$ and $q$ is convex on
$\Pi_{a,b}$, this means that we can find a point $\pi^{\prime}$ on the line
segment between $\Phi_{\lambda}(C_{1})$ and $\pi$ such that
$q(\pi^{\prime})<q(\Phi_{\lambda}(C_{1}))$, which contradicts that
$\Phi_{\lambda}(C_{1})$ is a minimizer of $q$ over $\Pi_{a,b}$. Hence, (29)
must hold. Letting $\pi=\Phi_{\lambda}(C_{2})$ in (29), we have

$\langle C_{1}+\lambda\nabla
h(\Phi_{\lambda}(C_{1})),\Phi_{\lambda}(C_{2})-\Phi_{\lambda}(C_{1})\rangle\geq
0.$

Changing the role of $C_{1}$ and $C_{2}$,

$\langle C_{2}+\lambda\nabla
h(\Phi_{\lambda}(C_{2})),\Phi_{\lambda}(C_{1})-\Phi_{\lambda}(C_{2})\rangle\geq
0.$

Combining the above two inequalities, we have

$\langle\nabla h(\Phi_{\lambda}(C_{1}))-\nabla
h(\Phi_{\lambda}(C_{2})),\Phi_{\lambda}(C_{1})-\Phi_{\lambda}(C_{2})\rangle\leq-\frac{1}{\lambda}\langle
C_{1}-C_{2},\Phi_{\lambda}(C_{1})-\Phi_{\lambda}(C_{2})\rangle.$ (30)

By applying Hölder’s inequality to the right-hand side of (30), we have

$-\frac{1}{\lambda}\langle
C_{1}-C_{2},\Phi_{\lambda}(C_{1})-\Phi_{\lambda}(C_{2})\rangle\leq\frac{\|C_{1}-C_{2}\|_{\infty}\cdot\|\Phi_{\lambda}(C_{1})-\Phi_{\lambda}(C_{2})\|_{1}}{\lambda}.$
(31)

Also, as $h$ is $1$-strongly convex with respect to $\|\cdot\|_{1}$ on
$\Delta_{n,m}$, we have

$\langle\nabla h(P_{1})-\nabla
h(P_{2}),P_{1}-P_{2}\rangle\geq\|P_{1}-P_{2}\|_{1}^{2}\quad\forall
P_{1},P_{2}\in\Delta_{n,m},$

which allows us to lower bound the left-hand side of (30) as follows:

$\langle\nabla h(\Phi_{\lambda}(C_{1}))-\nabla
h(\Phi_{\lambda}(C_{2})),\Phi_{\lambda}(C_{1})-\Phi_{\lambda}(C_{2})\rangle\geq\|\Phi_{\lambda}(C_{1})-\Phi_{\lambda}(C_{2})\|_{1}^{2}.$
(32)

Combining (30), (31), and (32), we have (28). ∎

###### Proof of Theorem 1.

By (28) of Lemma 1, for any $A_{1},A_{2}\in\mathbb{R}^{n\times m}$.

$\begin{split}\|T_{\lambda}(A_{1})-T_{\lambda}(A_{2})\|_{1}&=\|\Phi_{\lambda}(\nabla
g(A_{1}))-\Phi_{\lambda}(\nabla g(A_{2}))\|_{1}\\\ &\leq\frac{\|\nabla
g(A_{1})-\nabla g(A_{2})\|_{\infty}}{\lambda}\\\
&=\frac{\|H(A_{1}-A_{2})\|_{\infty}}{\lambda}\\\
&\leq\frac{\|H(A_{1}-A_{2})\|_{\infty}}{\lambda}\\\
&\leq\frac{\|H\|_{\infty}}{\lambda}\|A_{1}-A_{2}\|_{1},\end{split}$

where the last inequality follows from
$\|H(A_{1}-A_{2})\|_{\infty}\leq\|H\|_{\infty}\|A_{1}-A_{2}\|_{1}$. Therefore,
$T_{\lambda}$ is a contraction on $\Pi_{a,b}$ equipped with $\|\cdot\|_{1}$ if
$\lambda>L=\|H\|_{\infty}$. By the Banach–Caccioppoli theorem, $T_{\lambda}$
has a unique fixed point, and we have (23). ∎

###  A.4 Proof of Theorem 2

###### Proof of Theorem 2.

We first show the following from (24): for any $k\geq 0$ and
$\pi\in\Pi_{a,b}$,

$\langle\nabla f(\pi^{(k)}),\pi^{(k+1)}-\pi\rangle\leq\frac{\langle\nabla
h(\pi^{(k)})-\nabla h(\pi^{(k+1)}),\pi^{(k+1)}-\pi\rangle}{\tau_{k}}.$ (33)

To see this, pick $\pi\in\Pi_{a,b}$ and define
$q(\varepsilon)=f((1-\varepsilon)\pi^{(k+1)}+\varepsilon\pi)+\tau_{k}^{-1}D_{h}((1-\varepsilon)\pi^{(k+1)}+\varepsilon\pi,\pi^{(k)})$
for $\varepsilon\in[0,1]$. As $q(\varepsilon)\geq q(0)$ for
$\varepsilon\in[0,1]$, the right derivative of $q$ at $\varepsilon=0$ must be
nonnegative, which leads to (33). Next, letting $D_{f}$ be the Bregman
divergence with respect to $f$, we have

$\begin{split}f(\pi^{(k+1)})-f(\pi)&=f(\pi^{(k+1)})-f(\pi^{(k)})+f(\pi^{(k)})-f(\pi)\\\
&=D_{f}(\pi^{(k+1)},\pi^{(k)})+\langle\nabla
f(\pi^{(k)}),\pi^{(k+1)}-\pi^{(k)}\rangle+f(\pi^{(k)})-f(\pi)\\\ &\leq
D_{f}(\pi^{(k+1)},\pi^{(k)})+\langle\nabla
f(\pi^{(k)}),\pi^{(k+1)}-\pi^{(k)}\rangle+\langle\nabla
f(\pi^{(k)}),\pi^{(k)}-\pi\rangle\\\ &\leq
D_{f}(\pi^{(k+1)},\pi^{(k)})+\frac{\langle\nabla h(\pi^{(k)})-\nabla
h(\pi^{(k+1)}),\pi^{(k+1)}-\pi\rangle}{\tau_{k}},\end{split}$

where the first inequality uses the convexity of $f$ and the second inequality
follows from (33). For the last term on the right-hand side, we have

$\begin{split}\langle\nabla h(\pi^{(k)})-\nabla
h(\pi^{(k+1)}),\pi^{(k+1)}-\pi\rangle&=\langle\log\frac{\pi^{(k)}}{\pi^{(k+1)}},\pi^{(k+1)}-\pi\rangle\\\
&=\langle\log\frac{\pi}{\pi^{(k)}},\pi\rangle-\langle\log\frac{\pi}{\pi^{(k+1)}},\pi\rangle-\langle\log\frac{\pi^{(k+1)}}{\pi^{(k)}},\pi^{(k+1)}\rangle\\\
&=D_{h}(\pi,\pi^{(k)})-D_{h}(\pi,\pi^{(k+1)})-D_{h}(\pi^{(k+1)},\pi^{(k)}).\end{split}$

Hence,

$f(\pi^{(k+1)})-f(\pi)\leq
D_{f}(\pi^{(k+1)},\pi^{(k)})+\frac{D_{h}(\pi,\pi^{(k)})-D_{h}(\pi,\pi^{(k+1)})-D_{h}(\pi^{(k+1)},\pi^{(k)})}{\tau_{k}}.$
(34)

Meanwhile,

$\begin{split}D_{f}(\pi^{(k+1)},\pi^{(k)})&=\frac{1}{2}\langle\pi^{(k+1)}-\pi^{(k)},H(\pi^{(k+1)}-\pi^{(k)})\rangle+\lambda
D_{h}(\pi^{(k+1)},\pi^{(k)})\\\
&\leq\frac{\|H\|_{\infty}}{2}\|\pi^{(k+1)}-\pi^{(k)}\|_{1}^{2}+\lambda
D_{h}(\pi^{(k+1)},\pi^{(k)})\\\
&\leq(\|H\|_{\infty}+\lambda)D_{h}(\pi^{(k+1)},\pi^{(k)}),\end{split}$

where the first inequality follows from Hölder’s inequality and the second
inequality uses Pinsker’s inequality $\|\pi^{(k+1)}-\pi^{(k)}\|_{1}^{2}\leq
2D_{h}(\pi^{(k+1)},\pi^{(k)})$; see, Lemma 2.5 of [Tsy09]. Now, as
$\tau_{k}^{-1}=\tau^{-1}\geq\|H\|_{\infty}+\lambda$, we have

$D_{f}(\pi^{(k+1)},\pi^{(k)})\leq\frac{D_{h}(\pi^{(k+1)},\pi^{(k)})}{\tau},$

which, together with (34), leads to

$f(\pi^{(k+1)})-f(\pi)\leq\frac{D_{h}(\pi,\pi^{(k)})-D_{h}(\pi,\pi^{(k+1)})}{\tau}\quad\forall\pi\in\Pi_{a,b}.$

By letting $\pi=\pi^{(k)}$, we have $f(\pi^{(k+1)})\leq f(\pi^{(k)})$ for all
$k\geq 0$. Hence,

$f(\pi^{(T)})-f(\pi)\leq\frac{1}{T}\sum_{k=0}^{T-1}(f(\pi^{(k+1)})-f(\pi))\leq\frac{1}{T}\frac{D_{h}(\pi,\pi^{(0)})-D_{h}(\pi,\pi^{(T)})}{\tau}\leq\frac{1}{T}\frac{D_{h}(\pi,\pi^{(0)})}{\tau}.$

As $f$ always admits a minimizer on $\Pi_{a,b}$ (even when $\lambda=0$),
plugging a minimizer to $\pi$ and using the fact that
$D_{h}(\pi,\pi^{(0)})\leq\log(nm)$, we obtain (25). Lastly,

$\begin{split}f(\pi^{(T)})-f(\pi^{\star})&=\langle\nabla
f(\pi^{\star}),\pi^{(T)}-\pi^{\star}\rangle+D_{f}(\pi^{(T)},\pi^{\star})\\\
&=\langle\nabla
f(\pi^{\star}),\pi^{(T)}-\pi^{\star}\rangle+\frac{1}{2}\langle\pi^{(T)}-\pi^{(k)},H(\pi^{(T)}-\pi^{(k)})\rangle+\lambda
D_{h}(\pi^{(T)},\pi^{\star})\\\
&\geq\frac{\lambda}{2}\|\pi^{(T)}-\pi^{\star}\|_{1}^{2},\end{split}$

where the last inequality uses Pinsker’s inequality and $\langle\nabla
f(\pi^{\star}),\pi^{(T)}-\pi^{\star}\rangle\geq 0$, which holds as
$\pi^{\star}$ is a minimizer of $f$. Hence, we have (26). ∎