The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
63,723 characters
Kernel Instrumental Variable Regression
\maketitle
\begin{abstract}
Instrumental variable (IV) regression is a strategy for learning causal relationships in
observational data. If measurements of input $X$ and output $Y$ are confounded,
the causal relationship can nonetheless be identified if an instrumental variable $Z$
is available that influences $X$ directly, but is conditionally
independent of $Y$ given $X$ and the unmeasured confounder. The classic two-stage least squares algorithm (2SLS) simplifies the estimation problem by modeling all relationships as linear functions. We propose kernel instrumental variable regression (KIV), a nonparametric generalization of 2SLS, modeling relations among $X$, $Y$, and $Z$ as nonlinear functions in reproducing kernel Hilbert spaces (RKHSs). We prove the consistency of KIV under mild assumptions, and derive conditions under which convergence occurs at the minimax optimal rate for unconfounded, single-stage RKHS regression. In doing so, we obtain an efficient ratio between training sample sizes used in the algorithm's first and second stages. In experiments, KIV outperforms state of the art alternatives for nonparametric IV regression.
\end{abstract}
\section{Introduction}
Instrumental variable regression is a method in causal statistics
for estimating the counterfactual effect of input $X$ on output $Y$
using observational data \cite{stock2003retrospectives}. If measurements of $(X,Y)$ are confounded,
the causal relationship--also called the structural relationship--can nonetheless be identified if an instrumental variable $Z$ is
available, which is independent of $Y$ conditional on $X$ and the unmeasured confounder. Intuitively, $Z$ only influences $Y$ via $X$, identifying the counterfactual relationship of interest.
Economists and epidemiologists use instrumental variables to overcome issues of strategic interaction, imperfect compliance, and selection bias.
The original application is demand estimation: supply cost shifters ($Z$) only influence sales ($Y$) via price ($X$), thereby identifying counterfactual demand even though prices reflect both supply and demand market forces \cite{wright1928tariff,blundell2012measuring}. Randomized assignment of a drug ($Z$) only influences patient health ($Y$) via actual consumption of the drug ($X$), identifying the counterfactual effect of the drug even in the scenario of imperfect compliance \cite{angrist1996identification}. Draft lottery number ($Z$) only influences lifetime
earnings ($Y$) via military service ($X$), identifying the
counterfactual effect of military service on earnings despite selection bias
in enlistment \cite{angrist1990lifetime}.
The two-stage least squares algorithm (2SLS), widely used in economics,
simplifies the IV estimation problem by assuming linear relationships: in \textit{stage 1}, perform linear regression to obtain the conditional
means $\bar{x}(z):=\mathbb{E}_{X|Z=z}(X)$; in \textit{stage 2},
linearly regress outputs $Y$ on these conditional means.
2SLS works well when the underlying assumptions hold. In practice, the relation between $Y$ and $X$ may not be linear, nor may be the relation between $X$
and $Z$.
In the present work, we introduce kernel instrumental variable regression
(KIV), an easily implemented nonlinear generalization of 2SLS (Sections~\ref{sec:framework} and~\ref{sec:problemAndAlgo}).\footnote{Code: \url{https://github.com/r4hu1-5in9h/KIV}} In
\textit{stage 1} we learn a conditional mean embedding,
which is the conditional expectation $\mu(z):=\mathbb{E}_{X|Z=z}\psi(X)$
of features $\psi$ which map $X$ to a reproducing kernel Hilbert
space (RKHS) \cite{song2009hilbert}. For a sufficiently rich RKHS, called a characteristic RKHS, the mean embedding of a random variable
is injective \cite{sriperumbudur2010relation}. It follows that the conditional mean embedding characterizes the full distribution of $X$ conditioned on $Z$, and not just the conditional mean. We then implement \textit{stage 2} via kernel ridge regression of outputs $Y$ on these conditional mean embeddings, following the two-stage distribution regression approach described by \cite{szabo2015two,szabo2016learning}. As in our work, the inputs for \cite{szabo2015two,szabo2016learning} are distribution embeddings. Unlike our case, the earlier work uses unconditional embeddings computed from independent samples.
As a key contribution of our work, we provide consistency guarantees
for the KIV algorithm for an increasing number of training samples
in stages 1 and 2 (Section~\ref{sec:consistency}). To establish stage 1 convergence, we note that
the conditional mean embedding \cite{song2009hilbert} is the solution to a regression problem \cite{grunewalder2012conditional,grunewalder2012modelling,grunewalder2013smooth}, and thus equivalent to kernel dependency estimation \cite{ciliberto2016consistent,cortes2005general}. We prove that the kernel estimator of the conditional mean embedding (equivalently, the conditional expectation operator) converges in RKHS-norm, generalizing classic results by \cite{smale2005shannon,smale2007learning}.
We allow the conditional mean embedding RKHS to be infinite-dimensional, which presents specific challenges that we carefully address in our analysis. We also discuss previous approaches to establishing consistency in both finite-dimensional \cite{grunewalder2012modelling} and infinite-dimensional \cite{song2009hilbert,song2010nonparametric,fukumizu2013kernel,hefny2015supervised,ciliberto2016consistent} settings.
We embed the stage 1 rates into stage 2 to get end-to-end guarantees for the two-stage procedure, adapting
\cite{caponnetto2007optimal,szabo2015two,szabo2016learning}. In particular, we provide
a ratio of stage 1 to stage 2 samples required for minimax optimal rates in
the second stage, where the ratio depends on the difficulty of
each stage. We anticipate that these proof strategies will apply generally
in two-stage regression settings.
\section{Related work}
Several approaches have been proposed to generalize 2SLS to the
nonlinear setting, which we will compare in our experiments (Section~\ref{sec:experiments}). A first generalization is via basis function approximation \cite{newey2003instrumental}, an approach called sieve IV, with uniform convergence rates in \cite{chen2018optimal}. The challenge in \cite{chen2018optimal} is how to define an appropriate finite dictionary of basis functions. In a second approach, \cite{carrasco2007linear,darolles2011nonparametric} implement stage 1 by computing the conditional distribution of the input $X$
given the instrument $Z$ using a ratio of Nadaraya-Watson density
estimates. Stage 2 is then ridge regression in the space
of square integrable functions. The overall algorithm has a finite sample consistency guarantee, assuming smoothness of the $(X,Z)$ joint density in stage 1 and the regression in stage 2 \cite{darolles2011nonparametric}. Unlike our bound, \cite{darolles2011nonparametric} make no claim about the optimality of the result. Importantly, stage 1 requires the
solution of a statistically challenging problem: conditional density
estimation. Moreover, analysis assumes the same number of training samples used in both stages. We will discuss this bound in more
detail in Appendix~\ref{sec:comparisonDarollesNonparametric}
(we suggest that the reader first cover Section~\ref{sec:consistency}).
Our work also relates to kernel and IV approaches to learning dynamical systems, known in machine learning as predictive state representation models (PSRs) \cite{boots2013hilbert,hefny2015supervised,downey2017predictive} and in econometrics as panel data models \cite{anderson1981estimation,arellano1991some}. In this setting, predictive states (expected future features given
history) are updated in light of new observations. The calculation of
the predictive states corresponds to stage 1 regression, and the
states are updated via stage 2 regression. In the kernel case, the
predictive states are expressed as conditional mean embeddings
\cite{boots2013hilbert}, as in our setting. Performance of the kernel PSR method is guaranteed by a finite sample bound \cite[Theorem 2]{hefny2015supervised}, however this bound is not minimax optimal. Whereas \cite{hefny2015supervised} assume an equal number of training samples in stages 1 and 2, we find that unequal numbers of training samples matter for minimax optimality. More importantly, the bound makes strong
smoothness assumptions on the inputs to the stage 1 and stage 2 regression
functions, rather than assuming smoothness of the regression
functions as we do. We show that the smoothness assumptions on
the inputs made in \cite{hefny2015supervised} do not hold in our
setting, and we obtain stronger end-to-end bounds under more realistic
conditions. We discuss the PSR bound in more detail
in Appendix~\ref{sec:comparisonHefnyNonparametric}.
Yet another recent approach is deep IV, which uses neural networks in both stages and permits learning even for complex high-dimensional data such as images \cite{hartford2017deep}. Like \cite{darolles2011nonparametric}, \cite{hartford2017deep} implement stage 1 by estimating a conditional density. Unlike \cite{darolles2011nonparametric}, \cite{hartford2017deep} use a mixture density network \cite[Section 5.6]{Bishop06}, i.e. a mixture model parametrized by
a neural network on the instrument $Z$. Stage 2 is neural network regression, trained using stochastic gradient descent (SGD). This presents a challenge: each step of SGD requires expectations using the stage 1 model, which are computed by drawing samples and averaging. An unbiased gradient estimate requires two independent sets of samples from the stage 1 model \cite[eq. 10]{hartford2017deep}, though a single set of samples may be used if an upper bound on the loss is optimized \cite[eq. 11]{hartford2017deep}. By contrast, our stage 1 outputs--conditional mean embeddings--have a closed form solution and exhibit
lower variance than sample averaging from a conditional density model. No theoretical guarantee on the consistency of the neural network
approach has been provided.
In the econometrics literature, a few key assumptions make learning a nonparametric IV model tractable. These include the completeness condition \cite{newey2003instrumental}: the structural relationship between $X$ and $Y$ can be identified only if the stage 1 conditional expectation is injective. Subsequent works impose additional stability and link assumptions \cite{blundell2007semi,chen2012estimation,chen2018optimal}: the conditional expectation of a function of $X$ given $Z$ is a smooth function of $Z$. We adapt these assumptions to our setting, replacing the completeness condition with the characteristic property \cite{sriperumbudur2010relation}, and replacing the stability and link assumptions with the concept of prior \cite{smale2007learning,caponnetto2007optimal}. We describe the characteristic and prior assumptions in more detail below.
Extensive use of IV estimation in applied economic research has revealed a common pitfall: weak instrumental variables. A weak instrument satisfies Hypothesis 1 below, but the relationship between a weak instrument $Z$ and input $X$ is negligible; $Z$ is essentially irrelevant. In this case, IV estimation becomes highly erratic \cite{bound1995problems}. In \cite{staiger1997instrumental}, the authors formalize this phenomenon with local analysis. See \cite{murray2006avoiding,stock2002survey} for practical and theoretical overviews, respectively. We recommend that practitioners resist the temptation to use many weak instruments, and instead use few strong instruments such as those described in the introduction.
Finally, our analysis connects early work on the RKHS with recent developments in the RKHS literature. In \cite{nashed1974generalized}, the authors introduce the RKHS to solve known, ill-posed functional equations. In the present work, we introduce the RKHS to estimate the solution to an uncertain, ill-posed functional equation. In this sense, casting the IV problem in an RKHS framework is not only natural; it is in the original spirit of RKHS methods.
For a comprehensive review of existing work and recent advances in kernel mean embedding research, we recommend \cite{muandet2017kernel,gretton2018notes}.
\section{Problem setting and definitions}\label{sec:framework}
\textbf{Instrumental variable:}
We begin by introducing our causal assumption about the instrument. This prior knowledge, described informally in the introduction, allows us to recover the counterfactual effect of $X$ on $Y$. Let $(\mathcal{X},\mathcal{B}_{\mathcal{X}})$, $(\mathcal{Y},\mathcal{B}_{\mathcal{Y}})$, and $(\mathcal{Z},\mathcal{B}_{\mathcal{Z}})$ be measurable spaces. Let $(X,Y,Z)$ be a random variable on $\mathcal{X}\times \mathcal{Y}\times \mathcal{Z}$ with distribution $\rho$.
\begin{hypothesis}\label{iv}
Assume
\begin{enumerate}
\item $Y=h(X)+e$ and $\mathbb{E}[e|Z]=0$
\item $\rho(x|z)$ is not constant in $z$
\end{enumerate}
\end{hypothesis}
We call $h$ the \textit{structural function} of interest. The error term $e$ is unmeasured, confounding noise. Hypothesis~\ref{iv}.1, known as the exclusion restriction, was introduced by \cite{newey2003instrumental} to the nonparametric IV literature for its tractability.
Other hypotheses are possible, although a very different approach is then needed \cite{imbens2009identification}. Hypothesis~\ref{iv}.2, known as the relevance condition, ensures that $Z$ is actually informative. In Appendix~\ref{sec:simpleIVadditiveNoise}, we compare Hypothesis~\ref{iv} with alternative formulations of the IV assumption.
We make three observations. First, if $X=Z$ then Hypothesis~\ref{iv} reduces to the standard regression assumption of unconfounded inputs, and $h(X)=\mathbb{E}[Y|X]$; if $X=Z$ then prediction and counterfactual prediction coincide. The IV model is a framework that allows for causal inference in a more general variety of contexts, namely when $h(X)\neq \mathbb{E}[Y|X]$ so that prediction and counterfactual prediction are different learning problems. Second, Hypothesis~\ref{iv} will permit identification of $h$ even if inputs are confounded, i.e. $X\cancel\raisebox{0.05em}{\rotatebox[origin=c]{90}{$\models$}} e$. Third, this model includes the scenario in which the analyst has a combination of confounded and unconfounded inputs. For example, in demand estimation there may be confounded price $P$, unconfounded characteristics $W$, and supply cost shifter $C$ that instruments for price. Then $X=(P,W)$, $Z=(C,W)$, and the analysis remains the same.
\begin{wrapfigure}{R}{0.25\textwidth}
\vspace{-15pt}
\begin{center}
\begin{adjustbox}{width=.25\textwidth}
\begin{tikzpicture}[->,>=stealth',shorten >=1pt,auto,node distance=2.8cm,
semithick]
\tikzstyle{every state}=[draw=black,text=black]
\node[state] (z) [fill=gray] {$\mathcal{Z}$};
\node[state] (hz) [right of=z] {$\mathcal{H}_{\mathcal{Z}}$};
\node[state] (x) [below of=z, fill=gray] {$\mathcal{X}$};
\node[state] (hx) [right of=x] {$\mathcal{H}_{\mathcal{X}}$};
\node[state] (y) [below of=x, fill=gray] {$\mathcal{Y}$};
\path (z) edge node {$ $} (x)
edge node {$\mu\in \mathcal{H}_{\Xi}$} (hx)
edge node {$\phi$} (hz)
(x) edge node {$h\in\mathcal{H}_{\mathcal{X}}$} (y)
edge node {$\psi$} (hx)
(hz) edge node {$E^*\in\mathcal{H}_{\Gamma ^*}$} (hx)
(hx) edge node {$H\in\mathcal{H}_{\Omega}$} (y);;
\end{tikzpicture}
\end{adjustbox}
\vspace{-10pt}
\caption{The RKHSs}
\label{rkhs_s}
\end{center}
\vspace{-15pt}
\end{wrapfigure}
Hypothesis~\ref{iv} provides the operator equation $
\mathbb{E}[Y|Z]=\mathbb{E}_{X|Z}h(X)
$ \cite{newey2003instrumental}. In the language of 2SLS, the LHS is the \textit{reduced form}, while the RHS is a composition of \textit{stage 1} linear compact operator $\mathbb{E}_{X|Z}$ and \textit{stage 2} structural function $h$. In the language of functional analysis, the operator equation is a Fredholm integral equation of the first kind \cite{nashed1974generalized,kress1989linear,newey2003instrumental,florens2003inverse}. Solving this operator equation for $h$ involves inverting a linear compact operator with infinite-dimensional domain; it is an ill-posed problem \cite{kress1989linear}. To recover a well-posed problem, we impose smoothness and Tikhonov regularization.
\textbf{RKHS model:}
We next introduce our RKHS model. Let $k_{\mathcal{X}}:\mathcal{X}\times \mathcal{X}\rightarrow\mathbb{R}$ and $k_{\mathcal{Z}}:\mathcal{Z}\times \mathcal{Z}\rightarrow\mathbb{R}$ be measurable positive definite kernels corresponding to scalar-valued RKHSs $\mathcal{H}_{\mathcal{X}}$ and $\mathcal{H}_{\mathcal{Z}}$. Denote the feature maps
\begin{align*}
\psi&:\mathcal{X}\rightarrow\mathcal{H}_{\mathcal{X}} ,\enskip x\mapsto k_{\mathcal{X}}(x,\cdot)\quad\quad \phi:\mathcal{Z}\rightarrow\mathcal{H}_{\mathcal{Z}} ,\enskip z\mapsto k_{\mathcal{Z}}(z,\cdot)
\end{align*}
Define the \textit{conditional expectation operator} $E:\mathcal{H}_{\mathcal{X}}\rightarrow \mathcal{H}_{\mathcal{Z}}$ such that
$
[Eh](z)=\mathbb{E}_{X|Z=z}h(X)
$. $E$ is the natural object of interest for stage 1. We define and analyze an estimator for $E$ directly. The conditional expectation operator $E$ conveys exactly the same information as another object popular in the kernel methods literature, the \textit{conditional mean embedding} $\mu:\mathcal{Z}\rightarrow \mathcal{H}_{\mathcal{X}}$ defined by
$
\mu(z)=\mathbb{E}_{X|Z=z}\psi(X)
$ \cite{song2009hilbert}. Indeed, $\mu(z)=E^*\phi(z)$ where $E^*:\mathcal{H}_{\mathcal{Z}}\rightarrow\mathcal{H}_{\mathcal{X}}$ is the adjoint of $E$. Analogously, in 2SLS $\bar{x}(z)=\pi' z$ for stage 1 linear regression parameter $\pi$.
The structural function $h:\mathcal{X}\rightarrow \mathcal{Y}$ in Hypothesis \ref{iv} is the natural object of interest for stage 2. For theoretical purposes, it is convenient to estimate $h$ indirectly. The structural function $h$ conveys exactly the same information as an object we call the \textit{structural operator} $H:\mathcal{H}_{\mathcal{X}}\rightarrow\mathcal{Y}$. Indeed, $h(x)=H\psi(x)$. Analogously, in 2SLS $h(x)=\beta' x$ for structural parameter $\beta$. We define and analyze an estimator for $H$, which in turn implies an estimator for $h$. Figure~\ref{rkhs_s} summarizes the relationships among equivalent stage 1 objects $(E,\mu)$ and equivalent stage 2 objects $(H,h)$.
Our RKHS model for the IV problem is of the same form as the model in \cite{nashed1974convergence,nashed1974generalized, nashed1974regularization} for general operator equations.
We begin by choosing RKHSs for the structural function $h$ and the reduced form $\mathbb{E}[Y|Z]$, then
construct a tensor-product RKHS for the conditional expectation operator $E$.
Our model differs from the RKHS model proposed by \cite{carrasco2007linear,darolles2011nonparametric}, which directly learns the conditional expectation operator $E$ via Nadaraya-Watson density estimation. The RKHSs of \cite{engl1996regularization,carrasco2007linear,darolles2011nonparametric} for the structural function $h$ and the reduced form $\mathbb{E}[Y|Z]$
are defined from the right and left singular functions of $E$, respectively. They appear in the consistency argument, but not in the ridge penalty.
\section{Learning problem and algorithm}\label{sec:problemAndAlgo}
2SLS consists of two stages that can be estimated separately. Sample splitting in this context means estimating stage 1 with $n$ randomly chosen observations and estimating stage 2 with the remaining $m$ observations. Sample splitting alleviates the finite sample bias of 2SLS when instrument $Z$ weakly influences input $X$ \cite{angrist1995split}. It is the natural approach when an analyst does not have access to a single data set with $n+m$ observations of $(X,Y,Z)$ but rather two data sets: $n$ observations of $(X,Z)$, and $m$ observations of $(Y,Z)$. We employ sample splitting in KIV, with an efficient ratio of $(n,m)$ given in Theorem~\ref{rate}. In our presentation of the general two-stage learning problem, we denote stage 1 observations by $(x_i,z_i)$ and stage 2 observations by $(\tilde{y}_i,\tilde{z}_i)$.
\subsection{Stage 1}
We transform the problem of learning $E$ into a vector-valued kernel ridge regression following \cite{grunewalder2012conditional,grunewalder2013smooth,ciliberto2016consistent}, where the hypothesis space is the vector-valued RKHS $\mathcal{H}_{\Gamma}$ of operators mapping $\mathcal{H}_{\mathcal{X}}$ to $\mathcal{H}_{\mathcal{Z}}$.
In Appendix \ref{sec:vectorValuedRKHS}, we review the theory of vector-valued RKHSs as it relates to scalar-valued RKHSs and tensor product spaces. The key result is that the tensor product space of $\mathcal{H}_{\mathcal{X}}$ and $\mathcal{H}_{\mathcal{Z}}$ is isomorphic to $\mathcal{L}_2(\mathcal{H}_{\mathcal{X}},\mathcal{H}_{\mathcal{Z}})$, the space of Hilbert-Schmidt operators from $\mathcal{H}_{\mathcal{X}}$ to $\mathcal{H}_{\mathcal{Z}}$.
If we choose the vector-valued kernel $\Gamma$ with feature map $(x,z)\mapsto [\phi(z)\otimes \psi(x)](\cdot)=\phi(z)\langle \psi(x),\cdot \rangle_{\mathcal{H}_{\mathcal{X}}}$, then $\mathcal{H}_{\Gamma}= \mathcal{L}_2(\mathcal{H}_{\mathcal{X}},\mathcal{H}_{\mathcal{Z}})$ and it shares the same norm.
We now state the objective for optimizing $E\in \mathcal{H}_{\Gamma}$.
The optimal $E$ minimizes the expected discrepancy
\begin{align*}
E_{\rho}&=\operatorname*{\arg\!\min} \mathcal{E}_1(E),\quad \mathcal{E}_1(E)= \mathbb{E}_{(X,Z)} \|\psi(X)-E^*\phi(Z)\|^2_{\mathcal{H}_{\mathcal{X}}}
\end{align*}
Both \cite{grunewalder2013smooth} and \cite{ciliberto2016consistent} refer to $\mathcal{E}_1$ as the surrogate risk.
As shown in \cite[Section 3.1]{grunewalder2012conditional} and \cite{grunewalder2013smooth}, the surrogate risk upper bounds the natural risk for the conditional expectation,
where the bound becomes tight when $\mathbb{E}_{X|Z=(\cdot)}f(X)\in \mathcal{H}_{\mathcal{Z}},\:\forall f\in\mathcal{H}_{\mathcal{X}}$.
Formally, the target operator is the constrained solution $
E_{\mathcal{H}_{\Gamma}}=\operatorname*{\arg\!\min}_{E\in\mathcal{H}_{\Gamma}}\mathcal{E}_1(E)
$.
We will assume $E_{\rho}\in \mathcal{H}_{\Gamma}$ so that $E_{\rho}=E_{\mathcal{H}_{\Gamma}}$.
Next we impose Tikhonov regularization. The regularized target operator and its empirical analogue are given by
\begin{align*}
E_{\lambda}&=\operatorname*{\arg\!\min}_{E\in \mathcal{H}_{\Gamma}} \mathcal{E}_{\lambda}(E),\quad
\mathcal{E}_{\lambda}(E)=\mathcal{E}_1(E)+\lambda\|E\|^2_{\mathcal{L}_2(\mathcal{H}_{\mathcal{X}},\mathcal{H}_{\mathcal{Z}})} \\
E^n_{\lambda}&=\operatorname*{\arg\!\min}_{E\in \mathcal{H}_{\Gamma}} \mathcal{E}_{\lambda}^n(E),\quad
\mathcal{E}_{\lambda}^n(E)=\dfrac{1}{n}\sum_{i=1}^n \|\psi(x_i)-E^*\phi(z_i)\|^2_{\mathcal{H}_{\mathcal{X}}}+\lambda\|E\|^2_{\mathcal{L}_2(\mathcal{H}_{\mathcal{X}},\mathcal{H}_{\mathcal{Z}})}
\end{align*}
Our construction of a vector-valued RKHS $\mathcal{H}_{\Gamma}$ for the conditional expectation operator $E$ permits us to estimate stage 1 by kernel ridge regression. The stage 1 estimator of KIV is at once novel in the nonparametric IV literature and fundamentally similar to 2SLS. Basis function approximation \cite{newey2003instrumental,chen2018optimal} is perhaps the closest prior IV approach, but we use infinite dictionaries of basis functions $\psi$ and $\phi$. Compared to density estimation \cite{carrasco2007linear, darolles2011nonparametric,hartford2017deep}, kernel ridge regression is an easier problem.
Alternative stage 1 estimators in the literature estimate the singular system of $E$ to ensure that the adjoint of the estimator equals the estimator of the adjoint. These estimators differ in how they estimate the singular system: empirical distribution \cite{darolles2011nonparametric}, Nadaraya-Watson density \cite{darolles2004kernel}, or B-spline wavelets \cite{chen1997shape}. The KIV stage 1 estimator has the desired property by construction; $(E_{\lambda}^n)^*=(E^*)_{\lambda}^n$. See Appendix~\ref{sec:vectorValuedRKHS} for details.
\subsection{Stage 2}
Next, we transform the problem of learning $h$ into a scalar-valued kernel ridge regression that respects the IV problem structure. In Proposition \ref{prop:existenceOfConditionalMeanElement} of Appendix~\ref{sec:vectorValuedRKHS}, we show that under Hypothesis \ref{hyp:measurableBoundedFeatures} below,
$$
\mathbb{E}_{X|Z=z}h(X)=[Eh](z)=\langle h,\mu(z)\rangle_{\mathcal{H}_{\mathcal{X}}}=H\mu(z)
$$
where $h\in\mathcal{H}_{\mathcal{X}}$, a scalar-valued RKHS; $E\in \mathcal{H}_{\Gamma}$, the vector-valued RKHS described above; $\mu\in\mathcal{H}_{\Xi}$, a vector-valued RKHS isometrically isomorphic to $\mathcal{H}_{\Gamma}$; and $H\in \mathcal{H}_{\Omega}$, a scalar-valued RKHS isometrically isomorphic to $\mathcal{H}_{\mathcal{X}}$.
It is helpful to think of $\mu(z)$ as the embedding into $\mathcal{H}_{\mathcal{X}}$ of a distribution on $\mathcal{X}$ indexed by the conditioned value $z$. When $k_\mathcal{X}$ is characteristic, $\mu(z)$ uniquely embeds the conditional distribution, and $H$ is identified. The kernel $\Omega$ satisfies $k_{\mathcal{X}}(x,x')=\Omega(\psi(x),\psi(x'))$. This expression establishes the formal connection between our model and \cite{szabo2015two,szabo2016learning}. The choice of $\Omega$ may be more general; for nonlinear examples see \cite[Table 1]{szabo2016learning}.
We now state the objective for optimizing $H\in \mathcal{H}_{\Omega}$. Hypothesis~\ref{iv} provides the operator equation, which may be rewritten as the regression equation
$$
Y=\mathbb{E}_{X|Z}h(X)+e_Z=H\mu(Z)+e_Z,\quad \mathbb{E}[e_Z|Z]=0
$$
The unconstrained solution is
\begin{align*}
H_{\rho}&=\operatorname*{\arg\!\min} \mathcal{E}(H),\quad
\mathcal{E}(H) = \mathbb{E}_{(Y,Z)}\|Y-H\mu(Z)\|_{\mathcal{Y}}^2
\end{align*}
The target operator is the constrained solution $
H_{\mathcal{H}_{\Omega}}=\operatorname*{\arg\!\min}_{H\in\mathcal{H}_{\Omega}}\mathcal{E}(H)
$. We will assume $H_{\rho}\in \mathcal{H}_{\Omega}$ so that $H_{\rho}=H_{\mathcal{H}_{\Omega}}$. With regularization,
\begin{align*}
H_{\xi}&=\operatorname*{\arg\!\min}_{H\in \mathcal{H}_{\Omega}}\mathcal{E}_{\xi}(H),\quad
\mathcal{E}_{\xi}(H)=\mathcal{E}(H)+\xi\|H\|^2_{\mathcal{H}_{\Omega}} \\
H^{m}_{\xi}&=\operatorname*{\arg\!\min}_{H\in \mathcal{H}_{\Omega}}\mathcal{E}^{m}_{\xi}(H),\quad
\mathcal{E}^{m}_{\xi}(H)=\dfrac{1}{m}\sum_{i=1}^{m}\|\tilde{y}_i-H\mu(\tilde{z}_i)\|_{\mathcal{Y}}^2+\xi\|H\|^2_{\mathcal{H}_{\Omega}}
\end{align*}
The essence of the IV problem is this: we do not directly observe the conditional expectation operator $E$ (or equivalently the conditional mean embedding $\mu$) that appears in the stage 2 objective. Rather, we approximate it using the estimate from stage 1. Thus our KIV estimator is $\hat{h}^{m}_{\xi}=\hat{H}^{m}_{\xi}\psi$ where
\begin{align*}
\hat{H}^{m}_{\xi}&=\operatorname*{\arg\!\min}_{H\in \mathcal{H}_{\Omega}}\hat{\mathcal{E}}^{m}_{\xi}(H),\quad
\hat{\mathcal{E}}^{m}_{\xi}(H)=\dfrac{1}{m}\sum_{i=1}^{m}\|\tilde{y}_i-H\mu^n_{\lambda}(\tilde{z}_i) \|_{\mathcal{Y}}^2+\xi\|H\|^2_{\mathcal{H}_{\Omega}}
\end{align*}
and $\mu^n_{\lambda}=(E_{\lambda}^n)^*\phi$. The transition from $H_{\rho}$ to $H^{m}_{\xi}$ represents the fact that we only have $m$ samples. The transition from $H^{m}_{\xi}$ to $\hat{H}^{m}_{\xi}$ represents the fact that we must learn not only the structural operator $H$ but also the conditional expectation operator $E$. In this sense, the IV problem is more complex than the estimation problem considered by \cite{nashed1974convergence,nashed1974regularization} in which $E$ is known.
\subsection{Algorithm}
We obtain a closed form expression for the KIV estimator. The apparatus introduced above is required for analysis of consistency and convergence rate. More subtly, our RKHS construction allows us to write kernel ridge regression estimators for both stage 1 and stage 2, unlike previous work. Because KIV consists of repeated kernel ridge regressions, it benefits from repeated applications of the representer theorem \cite{wahba1990spline,scholkopf2001generalized}. Consequently, we have a shortcut for obtaining KIV's closed form; see Appendix~\ref{sec:alg_deriv} for the full derivation.
\begin{algorithm}\label{alg} Let $X$ and $Z$ be matrices of $n$ observations. Let $\tilde{y}$ and $\tilde{Z}$ be a vector and matrix of $m$ observations.
\begin{align*}
W&=K_{XX}(K_{ZZ}+n\lambda I)^{-1}K_{Z\tilde{Z}},\quad
\hat{\alpha}= (WW'+m\xi K_{XX})^{-1}W\tilde{y},\quad
\hat{h}_{\xi}^m(x)=(\hat{\alpha})'K_{Xx}
\end{align*}
where $K_{XX}$ and $K_{ZZ}$ are the empirical kernel matrices.
\end{algorithm}
Theorems~\ref{stage1} and~\ref{rate} below theoretically determine efficient rates for the stage 1 regularization parameter $\lambda$ and stage 2 regularization parameter $\xi$, respectively. In Appendix~\ref{sec:validation}, we provide a validation procedure to empirically determine values for $(\lambda,\xi)$.
\section{Consistency}\label{sec:consistency}
\subsection{Stage 1}
\textbf{Integral operators:} We use integral operator notation from the kernel methods literature, adapted to the conditional expectation operator learning problem. We denote by $L^2(\mathcal{Z},\rho_{\mathcal{Z}})$ the space of square integrable functions from $\mathcal{Z}$ to $\mathcal{Y}$
with respect to measure $\rho_{\mathcal{Z}}$, where $\rho_{\mathcal{Z}}$ is the restriction of $\rho$ to $\mathcal{Z}$.
\begin{definition}
The stage 1 (population) operators are
\begin{align*}
S_1^*&:\mathcal{H}_{\mathcal{Z}}\hookrightarrow L^2(\mathcal{Z},\rho_{\mathcal{Z}}) ,\enskip \ell\mapsto \langle \ell,\phi(\cdot) \rangle_{\mathcal{H}_{\mathcal{Z}}}\quad
S_1:L^2(\mathcal{Z},\rho_{\mathcal{Z}})\rightarrow\mathcal{H}_{\mathcal{Z}},\enskip \tilde{\ell}\mapsto \int\phi(z)\tilde{\ell}(z)d\rho_{\mathcal{Z}}(z)
\end{align*}
\end{definition}
$T_1=S_1\circ S_1^*$ is the uncentered covariance operator of \cite[Theorem 1]{fukumizu2004dimensionality}. In Appendix~\ref{sec:cov_technical}, we prove that $T_1$ exists and has finite trace even when $\mathcal{H}_{\mathcal{X}}$ and $\mathcal{H}_{\mathcal{Z}}$ are infinite-dimensional. In Appendix~\ref{sec:cov_review}, we compare $T_1$ with other covariance operators in the kernel methods literature.
\textbf{Assumptions:} We place assumptions on the original spaces $\mathcal{X}$ and $\mathcal{Z}$, the scalar-valued RKHSs $\mathcal{H}_{\mathcal{X}}$ and $\mathcal{H}_{\mathcal{Z}}$, and the probability distribution $\rho(x,z)$. We maintain these assumptions throughout the paper. Importantly, we assume that the vector-valued RKHS regression is correctly specified: the true conditional expectation operator $E_{\rho}$ lives in the vector-valued RKHS $\mathcal{H}_{\Gamma}$. In further research, we will relax this assumption.
\begin{hypothesis}Suppose that $\mathcal{X}$ and $\mathcal{Z}$ are Polish spaces, i.e. separable and completely metrizable topological spaces
\end{hypothesis}
\begin{hypothesis}\label{hyp:measurableBoundedFeatures} Suppose that
\begin{enumerate}
\item $k_{\mathcal{X}}$ and $k_{\mathcal{Z}}$ are continuous and bounded:
$
\sup_{x\in\mathcal{X}}\|\psi(x)\|_{\mathcal{H}_{\mathcal{X}}}\leq Q$, $
\sup_{z\in\mathcal{Z}}\|\phi(z)\|_{\mathcal{H}_{\mathcal{Z}}}\leq \kappa
$
\item $\psi$ and $\phi$ are measurable
\item $k_{\mathcal{X}}$ is characteristic \textup{\cite{sriperumbudur2010relation}}
\end{enumerate}
\end{hypothesis}
\begin{hypothesis}Suppose that $E_{\rho}\in \mathcal{H}_{\Gamma}$. Then $
\mathcal{E}_1(E_{\rho})=\inf_{E\in\mathcal{H}_{\Gamma}}\mathcal{E}_1(E)
$
\end{hypothesis}
Hypothesis 3.3 specializes the completeness condition of \cite{newey2003instrumental}. Hypotheses 2-4 are sufficient to bound the sampling error of the regularized estimator $E_{\lambda}^n$. Bounding the approximation error requires a further assumption on the smoothness of the distribution $\rho(x,z)$.
We assume $\rho(x,z)$ belongs to a class of distributions parametrized by $(\zeta_1,c_1)$, as generalized from \cite[Theorem 2]{smale2007learning} to the space $\mathcal{H}_\Gamma$.
\begin{hypothesis}\label{hyp:smaleSmoothness}
Fix $\zeta_1<\infty$. For given $c_1\in(1,2]$, define the prior $\mathcal{P}(\zeta_1,c_1)$ as the set of probability distributions $\rho$ on $\mathcal{X}\times \mathcal{Z}$ such that a range space assumption is satisfied: $\exists G_1\in\mathcal{H}_{\Gamma}$ s.t.
$
E_{\rho}=T_1^{\frac{c_1-1}{2}}\circ G_1$ and $
\|G_1\|^2_{\mathcal{H}_{\Gamma}}\leq \zeta_1
$
\end{hypothesis}
We use composition symbol $\circ$ to emphasize that $G_1:\mathcal{H}_{\mathcal{X}}\rightarrow\mathcal{H}_{\mathcal{Z}}$ and $T_1:\mathcal{H}_{\mathcal{Z}}\rightarrow \mathcal{H}_{\mathcal{Z}}$. We define the power of operator $T_1$ with respect to its eigendecomposition; see Appendix~\ref{sec:cov_technical} for formal justification. Larger $c_1$ corresponds to a smoother conditional expectation operator $E_{\rho}$. Proposition~\ref{cef} in Appendix~\ref{sec:stage1_reg} shows $E_{\rho}^*\phi(z)=\mu(z)$, so Hypothesis 5 is an indirect smoothness condition on the conditional mean embedding $\mu$.
\textbf{Estimation and convergence:} The estimator has a closed form solution, as noted in \cite[Section 3.1]{grunewalder2012conditional} and \cite[Appendix D]{grunewalder2012modelling}; \cite{ciliberto2016consistent} use it in the first stage of the structured prediction problem. We present the closed form solution in notation similar to \cite{caponnetto2007optimal} in order to elucidate how the estimator simply generalizes linear regression. This connection foreshadows our proof technique.
\begin{theorem}\label{sol_1}
$\forall \lambda>0$, the solution $E^n_{\lambda}$ of the regularized empirical objective $\mathcal{E}^n_{\lambda}$ exists, is unique, and
$$
E^n_{\lambda}= (\mathbf{T}_1+\lambda)^{-1}\circ\mathbf{g}_1,\quad \mathbf{T}_1=\dfrac{1}{n}\sum_{i=1}^n \phi(z_i)\otimes\phi(z_i),\quad
\mathbf{g}_1=\dfrac{1}{n}\sum_{i=1}^n \phi(z_i)\otimes \psi(x_i)
$$
\end{theorem}
We prove an original, finite sample bound on the RKHS-norm distance of the estimator $E^n_{\lambda}$ from its target $E_{\rho}$.
The proof is in Appendix~\ref{sec:stage1smaleconvergenceproofs}.
\begin{theorem}\label{stage1}
Assume Hypotheses 2-5. $\forall \delta\in(0,1)$, the following holds w.p. $1-\delta$:
\begin{align*}
\|E^n_{\lambda}-E_{\rho}\|_{\mathcal{H}_{\Gamma}}&\leq r_E(\delta,n,c_1):= \dfrac{ \sqrt{\zeta_1}(c_1+1)}{4^{\frac{1}{c_1+1}}} \bigg(\dfrac{4\kappa(Q+\kappa \|E_{\rho}\|_{\mathcal{H}_{\Gamma}}) \ln(2/\delta)}{ \sqrt{n\zeta_1}(c_1-1)}\bigg)^{\frac{c_1-1}{c_1+1}} \\
\lambda&=\bigg(\dfrac{8\kappa(Q+\kappa \|E_{\rho}\|_{\mathcal{H}_{\Gamma}}) \ln(2/\delta)}{ \sqrt{n\zeta_1}(c_1-1)}\bigg)^{\frac{2}{c_1+1}}
\end{align*}
\end{theorem}
The efficient rate of $\lambda$ is $n^{\frac{-1}{c_1+1}}$. Note that the convergence rate of $E^n_{\lambda}$ is calibrated by $c_1$, which measures the smoothness of the conditional expectation operator $E_{\rho}$.
\subsection{Stage 2}
\textbf{Integral operators:} We use integral operator notation from the kernel methods literature, adapted to the structural operator learning problem. We denote by $L^2(\mathcal{H}_{\mathcal{X}},\rho_{\mathcal{H}_{\mathcal{X}}})$ the space of square integrable functions from $\mathcal{H}_{\mathcal{X}}$ to $\mathcal{Y}$
with respect to measure $\rho_{\mathcal{H}_{\mathcal{X}}}$, where $\rho_{\mathcal{H}_{\mathcal{X}}}$ is the extension of $\rho$ to $\mathcal{H}_{\mathcal{X}}$ \cite[Lemma A.3.16]{steinwart2008support}. Note that we present stage 2 analysis for general output space $\mathcal{Y}$ as in \cite{szabo2015two,szabo2016learning}, though in practice we only consider $\mathcal{Y}\subset \mathbb{R}$ to simplify our two-stage RKHS model.
\begin{definition}\label{def:stage2operators}
The stage 2 (population) operators are
\begin{align*}
S^*&:\mathcal{H}_{\Omega}\hookrightarrow L^2(\mathcal{H}_{\mathcal{X}}, \rho_{\mathcal{H}_{\mathcal{X}}}),\enskip H\mapsto \Omega^*_{(\cdot)}H \\
S&:L^2(\mathcal{H}_{\mathcal{X}}, \rho_{\mathcal{H}_{\mathcal{X}}})\rightarrow\mathcal{H}_{\Omega},\enskip \tilde{H}\mapsto \int \Omega_{\mu(z)} \circ\tilde{H}\mu(z) d\rho_{\mathcal{H}_{\mathcal{X}}}(\mu(z))
\end{align*}
\end{definition}
where $\Omega_{\mu(z)}: \mathcal{Y}\rightarrow\mathcal{H}_{\Omega}$ defined by $y\mapsto \Omega(\cdot,\mu(z))y$ is the point evaluator of \cite{micchelli2005learning,carmeli2006vector}. Finally define $T_{\mu(z)}=\Omega_{\mu(z)}\circ \Omega^*_{\mu(z)}$ and covariance operator $T=S\circ S^*$.
\textbf{Assumptions:} We place assumptions on the original space $\mathcal{Y}$, the scalar-valued RKHS $\mathcal{H}_{\Omega}$, and the probability distribution $\rho$. Importantly, we assume that the scalar-valued
RKHS regression is correctly specified: the true structural operator $H_{\rho}$ lives in the scalar-valued RKHS $\mathcal{H}_{\Omega}$.
\begin{hypothesis}\label{hyp:polish}Suppose that $\mathcal{Y}$ is a Polish space
\end{hypothesis}
\begin{hypothesis}\label{hyp:hilbertSchmidt} Suppose that
\begin{enumerate}
\item The $\{\Omega_{\mu(z)}\}$ operator family is uniformly bounded in Hilbert-Schmidt norm: $\exists B$ s.t. $\forall \mu(z)$,
$
\|\Omega_{\mu(z)}\|^2_{\mathcal{L}_2(\mathcal{Y},\mathcal{H}_{\Omega})}=Tr(\Omega_{\mu(z)}^*\circ \Omega_{\mu(z)})\leq B
$
\item The $\{\Omega_{\mu(z)}\}$ operator family is H\"older continuous in operator norm: $\exists L>0$, $\iota \in(0,1]$ s.t. $\forall \mu(z),\mu(z')$,
$
\|\Omega_{\mu(z)}-\Omega_{\mu(z')}\|_{\mathcal{L}(\mathcal{Y},\mathcal{H}_{\Omega})} \leq L \|\mu(z)-\mu(z')\|^{\iota}_{\mathcal{H}_{\mathcal{X}}}
$
\end{enumerate}
\end{hypothesis}
Larger $\iota$ is interpretable as smoother kernel $\Omega$.
\begin{hypothesis}\label{hyp:boundedY} Suppose that
\begin{enumerate}
\item $H_{\rho}\in \mathcal{H}_{\Omega}$. Then $
\mathcal{E}(H_{\rho})=\inf_{H\in\mathcal{H}_{\Omega}}\mathcal{E}(H)
$
\item $Y$ is bounded, i.e. $\exists C<\infty$ s.t. $\|Y\|_{\mathcal{Y}}\leq C$ almost surely
\end{enumerate}
\end{hypothesis}
The convergence rate from stage 1 together with Hypotheses 6-8 are sufficient to bound the excess error of the regularized estimator $\hat{H}_{\xi}^{m}$ in terms of familiar objects in the kernel methods literature, namely the residual, reconstruction error, and effective dimension. We further assume $\rho$ belongs to a stage 2 prior to simplify these bounds. In particular, we assume $\rho$ belongs to a class of distributions parametrized by $(\zeta,b,c)$ as defined originally in \cite[Definition 1]{caponnetto2007optimal}, restated below.
\begin{hypothesis}\label{hyp:caponnettoSmoothness}
Fix $\zeta<\infty$. For given $b\in(1,\infty]$ and $c\in(1,2]$, define the prior $\mathcal{P}(\zeta,b,c)$ as the set of probability distributions $\rho$ on $\mathcal{H}_{\mathcal{X}}\times \mathcal{Y}$ such that
\begin{enumerate}
\item A range space assumption is satisfied: $\exists G\in\mathcal{H}_{\Omega}$ s.t.
$
H_{\rho} =T^{\frac{c-1}{2}}G$ and $
\|G\|^2_{\mathcal{H}_{\Omega}}\leq \zeta
$
\item In the spectral decomposition $
T=\sum_{k=1}^{\infty} \lambda_k e_k \langle \cdot, e_k\rangle_{\mathcal{H}_{\Omega}}
$, where $\{e_k\}_{k=1}^{\infty}$ is a basis of $Ker(T)^{\perp}$, the eigenvalues satisfy $\alpha\leq k^b \lambda_k \leq \beta$ for some $\alpha,\beta>0$
\end{enumerate}
\end{hypothesis}
We define the power of operator $T$ with respect to its eigendecomposition; see Appendix~\ref{sec:cov_technical} for formal justification. The latter condition is interpretable as polynomial decay of eigenvalues: $\lambda_k=\Theta(k^{-b})$. Larger $b$ means faster decay of eigenvalues of the covariance operator $T$ and hence smaller effective input dimension. Larger $c$ corresponds to a smoother structural operator $H_{\rho}$ \cite{szabo2016learning}.
\textbf{Estimation and convergence:} The estimator has a closed form solution, as shown by \cite{szabo2015two,szabo2016learning} in the second stage of the distribution regression problem. We present the solution in notation similar to \cite{caponnetto2007optimal} to elucidate how the stage 1 and stage 2 estimators have the same structure.
\begin{theorem}\label{sol_2}
$\forall \xi>0$, the solution $H^{m}_{\xi}$ to $\mathcal{E}_{\xi}^m$
and the solution $\hat{H}^{m}_{\xi}$ to $\hat{\mathcal{E}}_{\xi}^m$ exist, are unique, and
\begin{align*}
H_{\xi}^{m}&=(\mathbf{T}+\xi)^{-1}\mathbf{g},\quad \mathbf{T}=\dfrac{1}{m}\sum_{i=1}^m T_{\mu(\tilde{z}_i)},\quad \mathbf{g}=\dfrac{1}{m}\sum_{i=1}^m \Omega_{\mu(\tilde{z}_i)}\tilde{y}_i \\
\hat{H}_{\xi}^{m}&=(\hat{\mathbf{T}}+\xi)^{-1}\hat{\mathbf{g}},\quad \hat{\mathbf{T}}=\dfrac{1}{m}\sum_{i=1}^m T_{\mu^n_{\lambda}(\tilde{z}_i)},\quad \hat{\mathbf{g}}=\dfrac{1}{m}\sum_{i=1}^m \Omega_{\mu^n_{\lambda}(\tilde{z}_i)}\tilde{y}_i
\end{align*}
\end{theorem}
We now present this paper's main theorem. In Appendix~\ref{sec:finalTheorems}, we provide a finite sample bound on the excess error of the estimator $\hat{H}^{m}_{\xi}$ with respect to its target $H_{\rho}$.
Adapting arguments by \cite{szabo2016learning}, we demonstrate that KIV is able to achieve the minimax optimal single-stage rate derived by \cite{caponnetto2007optimal}. In other words, our two-stage estimator is able to learn the causal relationship with confounded data equally well as single-stage RKHS regression is able to learn the causal relationship with unconfounded data.
\begin{theorem}\label{rate}
Assume Hypotheses 1-9. Choose $\lambda=n^{-\frac{1}{c_1+1}}$ and $n=m^{\frac{a(c_1+1)}{\iota(c_1-1)}}$ where $a>0$.
\begin{enumerate}
\item If $a\leq \frac{b(c+1)}{bc+1}$ then $\mathcal{E}(\hat{H}_{\xi}^{m})-\mathcal{E}(H_{\rho})=O_p(m^{-\frac{ac}{c+1}})$ with $\xi=m^{-\frac{a}{c+1}}$
\item If $a\geq \frac{b(c+1)}{bc+1}$ then $\mathcal{E}(\hat{H}_{\xi}^{m})-\mathcal{E}(H_{\rho})=O_p(m^{-\frac{bc}{bc+1}})$ with $\xi=m^{-\frac{b}{bc+1}}$
\end{enumerate}
\end{theorem}
At $a=\frac{b(c+1)}{bc+1}<2$, the convergence rate $m^{-\frac{bc}{bc+1}}$ is minimax optimal while requiring the fewest observations \cite{szabo2016learning}. This statistically efficient rate is calibrated by $b$, the effective input dimension, as well as $c$, the smoothness of structural operator $H_{\rho}$ \cite{caponnetto2007optimal}. The efficient ratio between stage 1 and stage 2 samples is $n=m^{\frac{b(c+1)}{bc+1}\cdot \frac{(c_1+1)}{\iota(c_1-1)}}$, implying $n>m$. As far as we know, asymmetric sample splitting is a novel prescription in the IV literature; previous analyses assume $n=m$ \cite{angrist1995split,hefny2015supervised}.
\section{Experiments}\label{sec:experiments}
We compare the empirical performance of KIV (\verb|KernelIV|) to four leading competitors: standard kernel ridge regression (\verb|KernelReg|) \cite{saunders1998ridge}, Nadaraya-Watson IV (\verb|SmoothIV|) \cite{carrasco2007linear,darolles2011nonparametric}, sieve IV (\verb|SieveIV|) \cite{newey2003instrumental,chen2018optimal}, and deep IV (\verb|DeepIV|) \cite{hartford2017deep}. To improve the performance of sieve IV, we impose Tikhonov regularization in both stages with KIV's tuning procedure. This adaptation exceeds the theoretical justification provided by \cite{chen2018optimal}. However, it is justified by our analysis insofar as sieve IV is a special case of KIV: set feature maps $\psi,\phi$ equal to the sieve bases.
\begin{wrapfigure}{R}{\textwidth/3}
\vspace{5pt}
\begin{center}
\includegraphics[width=\textwidth/3]{sigmoid.eps}
\vspace{-10pt}
\caption{Sigmoid design}
\label{sigmoid_main}
\end{center}
\vspace{-0pt}
\end{wrapfigure}
We implement each estimator on three designs. The \textit{linear} design \cite{chen2018optimal} involves learning counterfactual function $h(x)=4x-2$, given confounded observations of continuous variables $(X,Y)$ as well as continuous instrument $Z$. The \textit{sigmoid} design \cite{chen2018optimal} involves learning counterfactual function $h(x)=\ln (|16x-8|+1)\cdot sgn(x-0.5)$ under the same regime. The \textit{demand} design \cite{hartford2017deep} involves learning demand function $h(p,t,s)=100+(10+p)\cdot s\cdot \psi(t)-2p$ where $
\psi(t)$ is the complex nonlinear function in Figure~\ref{demand_dgp}. An observation consists of $(Y,P,T,S,C)$ where $Y$ is sales, $P$ is price, $T$ is time of year, $S$ is customer sentiment (a discrete variable), and $C$ is a supply cost shifter. The parameter $\rho\in\{0.9,0.75,0.5,0.25,0.1\}$ calibrates the extent to which price $P$ is confounded by supply-side market forces. In KIV notation, inputs are $X=(P,T,S)$ and instruments are $Z=(C,T,S)$.
For each algorithm, design, and sample size, we implement 40 simulations and calculate MSE with respect to the true structural function $h$. Figures~\ref{sigmoid_main},~\ref{demand}, and~\ref{uni} visualize results.
In the sigmoid design, \verb|KernelIV| performs best across sample sizes. In the demand design, \verb|SmoothIV| performs best for sample size $n+m=1000$. Like \cite{hartford2017deep}, we do not implement \verb|SmoothIV| for greater sample sizes due to running time. Among estimators that we are able to implement, \verb|KernelIV| performs best for sample sizes $n+m=5000$ and $n+m=10000$. \verb|KernelReg| ignores the instrument $Z$, and it is biased away from the structural function due to confounding noise $e$. This phenomenon can have counterintuitive consequences. Figure~\ref{demand} shows that in the highly nonlinear demand design, \verb|KernelReg| deviates further from the structural function as sample size increases because the algorithm is further misled by confounded data. Figure 2 of \cite{hartford2017deep} documents the same effect when a feedforward neural network is used. The remaining algorithms make use of the instrument $Z$ to overcome this issue.
\verb|KernelIV| improves on \verb|SieveIV| in the same way that kernel ridge regression improves on ridge regression: by using an infinite dictionary of implicit basis functions rather than a finite dictionary of explicit basis functions. \verb|KernelIV| improves on \verb|SmoothIV| by using kernel ridge regression in not only stage 2 but also stage 1, avoiding costly density estimation. Finally, it improves on \verb|DeepIV| by directly learning stage 1 mean embeddings, rather than performing costly density estimation and sampling from the estimated density. The experiments show that \verb|KernelIV| works particularly well when the true structural function $h$ is smooth, confirming the theoretical guarantee of Theorem~\ref{rate}. See Appendix~\ref{sec:sim} for representative plots, implementation details, and a robustness study.
\begin{figure}[H]
\centering
\vspace{-5pt}
\includegraphics[width=\textwidth]{demand.eps}
\vspace{-15pt}
\caption{Demand design}
\label{demand}
\vspace{-5pt}
\end{figure}
\section{Conclusion}
We introduce KIV, an algorithm for learning a nonlinear, causal relationship from confounded observational data. KIV is easily implemented and minimax optimal. As a contribution to the IV literature, we show how to estimate the stage 1 conditional expectation operator--an infinite by infinite dimensional object--by kernel ridge regression. As a contribution to the kernel methods literature, we show how the RKHS is well-suited to causal inference and ill-posed inverse problems.
In simulations, KIV outperforms state of the art algorithms for nonparametric IV regression. The success of KIV suggests RKHS methods may be an effective bridge between econometrics and machine learning.
\subsubsection*{Acknowledgments}
We are grateful to Alberto Abadie, Anish Agarwal, Michael Arbel, Victor Chernozhukov, Geoffrey Gordon, Jason Hartford, Motonobu Kanagawa, Anna Mikusheva, Whitney Newey, Nakul Singh, Bharath Sriperumbudur, and Suhas Vijaykumar. This project was made possible by the Marshall Aid Commemoration Commission.
\begin{thebibliography}{10}
\bibitem{anderson1981estimation}
Theodore~W Anderson and Cheng Hsiao.
\newblock Estimation of dynamic models with error components.
\newblock {\em Journal of the American Statistical Association},
76(375):598--606, 1981.
\bibitem{angrist1990lifetime}
Joshua~D Angrist.
\newblock Lifetime earnings and the {V}ietnam era draft lottery: Evidence from
{S}ocial {S}ecurity administrative records.
\newblock {\em The American Economic Review}, pages 313--336, 1990.
\bibitem{angrist1996identification}
Joshua~D Angrist, Guido~W Imbens, and Donald~B Rubin.
\newblock Identification of causal effects using instrumental variables.
\newblock {\em Journal of the American Statistical Association},
91(434):444--455, 1996.
\bibitem{angrist1995split}
Joshua~D Angrist and Alan~B Krueger.
\newblock Split-sample instrumental variables estimates of the return to
schooling.
\newblock {\em Journal of Business \& Economic Statistics}, 13(2):225--235,
1995.
\bibitem{arbel2018kernel}
Michael Arbel and Arthur Gretton.
\newblock Kernel conditional exponential family.
\newblock In {\em International Conference on Artificial Intelligence and
Statistics}, pages 1337--1346, 2018.
\bibitem{arellano1991some}
Manuel Arellano and Stephen Bond.
\newblock Some tests of specification for panel data: {M}onte {C}arlo evidence
and an application to employment equations.
\newblock {\em The Review of Economic Studies}, 58(2):277--297, 1991.
\bibitem{bell2016handout}
Jordan Bell.
\newblock Trace class operators and {H}ilbert-{S}chmidt operators.
\newblock Technical report, University of Toronto Department of Mathematics,
2016.
\bibitem{berlinet2011reproducing}
Alain Berlinet and Christine Thomas-Agnan.
\newblock {\em Reproducing Kernel Hilbert Spaces in Probability and
Statistics}.
\newblock Springer, 2011.
\bibitem{Bishop06}
Christopher Bishop.
\newblock {\em Pattern Recognition and Machine Learning}.
\newblock Springer, 2006.
\bibitem{blundell2007semi}
Richard Blundell, Xiaohong Chen, and Dennis Kristensen.
\newblock Semi-nonparametric {IV} estimation of shape-invariant {Engel} curves.
\newblock {\em Econometrica}, 75(6):1613--1669, 2007.
\bibitem{blundell2012measuring}
Richard Blundell, Joel~L Horowitz, and Matthias Parey.
\newblock Measuring the price responsiveness of gasoline demand: Economic shape
restrictions and nonparametric demand estimation.
\newblock {\em Quantitative Economics}, 3(1):29--51, 2012.
\bibitem{boots2013hilbert}
Byron Boots, Arthur Gretton, and Geoffrey~J Gordon.
\newblock Hilbert space embeddings of predictive state representations.
\newblock In {\em Uncertainty in Artificial Intelligence}, pages 92--101, 2013.
\bibitem{bound1995problems}
John Bound, David~A Jaeger, and Regina~M Baker.
\newblock Problems with instrumental variables estimation when the correlation
between the instruments and the endogenous explanatory variable is weak.
\newblock {\em Journal of the American Statistical Association},
90(430):443--450, 1995.
\bibitem{caponnetto2007optimal}
Andrea Caponnetto and Ernesto De~Vito.
\newblock Optimal rates for the regularized least-squares algorithm.
\newblock {\em Foundations of Computational Mathematics}, 7(3):331--368, 2007.
\bibitem{carmeli2006vector}
Claudio Carmeli, Ernesto De~Vito, and Alessandro Toigo.
\newblock Vector valued reproducing kernel {H}ilbert spaces of integrable
functions and {M}ercer theorem.
\newblock {\em Analysis and Applications}, 4(4):377--408, 2006.
\bibitem{carrasco2007linear}
Marine Carrasco, Jean-Pierre Florens, and Eric Renault.
\newblock Linear inverse problems in structural econometrics estimation based
on spectral decomposition and regularization.
\newblock {\em Handbook of Econometrics}, 6:5633--5751, 2007.
\bibitem{chen2018optimal}
Xiaohong Chen and Timothy~M Christensen.
\newblock Optimal sup-norm rates and uniform inference on nonlinear functionals
of nonparametric {IV} regression.
\newblock {\em Quantitative Economics}, 9(1):39--84, 2018.
\bibitem{chen1997shape}
Xiaohong Chen, Lars~P Hansen, and Jose Scheinkman.
\newblock Shape-preserving estimation of diffusions.
\newblock Technical report, University of Chicago Department of Economics,
1997.
\bibitem{chen2012estimation}
Xiaohong Chen and Demian Pouzo.
\newblock Estimation of nonparametric conditional moment models with possibly
nonsmooth generalized residuals.
\newblock {\em Econometrica}, 80(1):277--321, 2012.
\bibitem{ciliberto2016consistent}
Carlo Ciliberto, Lorenzo Rosasco, and Alessandro Rudi.
\newblock A consistent regularization approach for structured prediction.
\newblock In {\em Advances in Neural Information Processing Systems}, pages
4412--4420, 2016.
\bibitem{cortes2005general}
Corinna Cortes, Mehryar Mohri, and Jason Weston.
\newblock A general regression technique for learning transductions.
\newblock In {\em International Conference on Machine Learning}, pages
153--160, 2005.
\bibitem{cucker2002mathematical}
Felipe Cucker and Steve Smale.
\newblock On the mathematical foundations of learning.
\newblock {\em Bulletin of the American mathematical society}, 39(1):1--49,
2002.
\bibitem{darolles2011nonparametric}
Serge Darolles, Yanqin Fan, Jean-Pierre Florens, and Eric Renault.
\newblock Nonparametric instrumental regression.
\newblock {\em Econometrica}, 79(5):1541--1565, 2011.
\bibitem{darolles2004kernel}
Serge Darolles, Jean-Pierre Florens, and Christian Gourieroux.
\newblock Kernel-based nonlinear canonical analysis and time reversibility.
\newblock {\em Journal of Econometrics}, 119(2):323--353, 2004.
\bibitem{de2005risk}
Ernesto De~Vito and Andrea Caponnetto.
\newblock Risk bounds for regularized least-squares algorithm with
operator-value kernels.
\newblock Technical report, MIT CSAIL, 2005.
\bibitem{downey2017predictive}
Carlton Downey, Ahmed Hefny, Byron Boots, Geoffrey~J Gordon, and Boyue Li.
\newblock Predictive state recurrent neural networks.
\newblock In {\em Advances in Neural Information Processing Systems}, pages
6053--6064, 2017.
\bibitem{Dutordoir18Gaussian}
Vincent Dutordoir, Hugh Salimbeni, James Hensman, and Marc Deisenroth.
\newblock Gaussian process conditional density estimation.
\newblock In {\em Advances in Neural Information Processing Systems}, pages
2385--2395, 2018.
\bibitem{engl1996regularization}
Heinz~W Engl, Martin Hanke, and Andreas Neubauer.
\newblock {\em Regularization of Inverse Problems}, volume 375.
\newblock Springer Science \& Business Media, 1996.
\bibitem{florens2003inverse}
Jean-Pierre Florens.
\newblock Inverse problems and structural econometrics.
\newblock In {\em Advances in Economics and Econometrics: Theory and
Applications}, volume~2, pages 46--85, 2003.
\bibitem{fukumizu2004dimensionality}
Kenji Fukumizu, Francis~R Bach, and Michael~I Jordan.
\newblock Dimensionality reduction for supervised learning with reproducing
kernel {H}ilbert spaces.
\newblock {\em Journal of Machine Learning Research}, 5(Jan):73--99, 2004.
\bibitem{fukumizu2013kernel}
Kenji Fukumizu, Le~Song, and Arthur Gretton.
\newblock Kernel {B}ayes' rule: Bayesian inference with positive definite
kernels.
\newblock {\em Journal of Machine Learning Research}, 14(1):3753--3783, 2013.
\bibitem{gretton2018notes}
Arthur Gretton.
\newblock R{K}{H}{S} in machine learning: Testing statistical dependence.
\newblock Technical report, UCL Gatsby Unit, 2018.
\bibitem{grunewalder2013smooth}
Steffen Gr{\"u}new{\"a}lder, Arthur Gretton, and John Shawe-Taylor.
\newblock Smooth operators.
\newblock In {\em International Conference on Machine Learning}, pages
1184--1192, 2013.
\bibitem{grunewalder2012conditional}
Steffen Gr{\"u}new{\"a}lder, Guy Lever, Luca Baldassarre, Sam Patterson, Arthur
Gretton, and Massimiliano Pontil.
\newblock Conditional mean embeddings as regressors.
\newblock In {\em International Conference on Machine Learning}, volume~5,
2012.
\bibitem{grunewalder2012modelling}
Steffen Gr{\"u}new{\"a}lder, Guy Lever, Luca Baldassarre, Massimilano Pontil,
and Arthur Gretton.
\newblock Modelling transition dynamics in {MDP}s with {RKHS} embeddings.
\newblock In {\em International Conference on Machine Learning}, pages
1603--1610, 2012.
\bibitem{hartford2017deep}
Jason Hartford, Greg Lewis, Kevin Leyton-Brown, and Matt Taddy.
\newblock Deep {IV}: A flexible approach for counterfactual prediction.
\newblock In {\em International Conference on Machine Learning}, pages
1414--1423, 2017.
\bibitem{hefny2015supervised}
Ahmed Hefny, Carlton Downey, and Geoffrey~J Gordon.
\newblock Supervised learning for dynamical system learning.
\newblock In {\em Advances in Neural Information Processing Systems}, pages
1963--1971, 2015.
\bibitem{hernan2019causal}
Miguel~A Hernan and James~M Robins.
\newblock {\em Causal Inference}.
\newblock CRC Press, 2019.
\bibitem{Hsu2012tail}
Daniel Hsu, Sham Kakade, and Tong Zhang.
\newblock Tail inequalities for sums of random matrices that depend on the
intrinsic dimension.
\newblock {\em Electronic Communications in Probability}, 17(14):1--13, 2012.
\bibitem{imbens2009identification}
Guido~W Imbens and Whitney~K Newey.
\newblock Identification and estimation of triangular simultaneous equations
models without additivity.
\newblock {\em Econometrica}, 77(5):1481--1512, 2009.
\bibitem{kress1989linear}
Rainer Kress.
\newblock {\em Linear Integral Equations}, volume~3.
\newblock Springer, 1989.
\bibitem{micchelli2005learning}
Charles~A Micchelli and Massimiliano Pontil.
\newblock On learning vector-valued functions.
\newblock {\em Neural Computation}, 17(1):177--204, 2005.
\bibitem{muandet2017kernel}
Krikamol Muandet, Kenji Fukumizu, Bharath Sriperumbudur, Bernhard
Sch{\"o}lkopf, et~al.
\newblock Kernel mean embedding of distributions: A review and beyond.
\newblock {\em Foundations and Trends in Machine Learning}, 10(1-2):1--141,
2017.
\bibitem{murray2006avoiding}
Michael~P Murray.
\newblock Avoiding invalid instruments and coping with weak instruments.
\newblock {\em Journal of Economic Perspectives}, 20(4):111--132, 2006.
\bibitem{nashed1974convergence}
M~Zuhair Nashed and Grace Wahba.
\newblock Convergence rates of approximate least squares solutions of linear
integral and operator equations of the first kind.
\newblock {\em Mathematics of Computation}, 28(125):69--80, 1974.
\bibitem{nashed1974generalized}
M~Zuhair Nashed and Grace Wahba.
\newblock Generalized inverses in reproducing kernel spaces: An approach to
regularization of linear operator equations.
\newblock {\em SIAM Journal on Mathematical Analysis}, 5(6):974--987, 1974.
\bibitem{nashed1974regularization}
M~Zuhair Nashed and Grace Wahba.
\newblock Regularization and approximation of linear operator equations in
reproducing kernel spaces.
\newblock {\em Bulletin of the American Mathematical Society},
80(6):1213--1218, 1974.
\bibitem{newey2003instrumental}
Whitney~K Newey and James~L Powell.
\newblock Instrumental variable estimation of nonparametric models.
\newblock {\em Econometrica}, 71(5):1565--1578, 2003.
\bibitem{pearl2009causality}
Judea Pearl.
\newblock {\em Causality}.
\newblock Cambridge University Press, 2009.
\bibitem{saunders1998ridge}
Craig Saunders, Alexander Gammerman, and Volodya Vovk.
\newblock Ridge regression learning algorithm in dual variables.
\newblock In {\em International Conference on Machine Learning}, pages
515--521, 1998.
\bibitem{scholkopf2001generalized}
Bernhard Sch{\"o}lkopf, Ralf Herbrich, and Alex~J Smola.
\newblock A generalized representer theorem.
\newblock In {\em International Conference on Computational Learning Theory},
pages 416--426, 2001.
\bibitem{singh2019causal}
Rahul Singh.
\newblock Causal inference tutorial.
\newblock Technical report, MIT Economics, 2019.
\bibitem{smale2005shannon}
Steve Smale and Ding-Xuan Zhou.
\newblock Shannon sampling {II}: Connections to learning theory.
\newblock {\em Applied and Computational Harmonic Analysis}, 19(3):285--302,
2005.
\bibitem{smale2007learning}
Steve Smale and Ding-Xuan Zhou.
\newblock Learning theory estimates via integral operators and their
approximations.
\newblock {\em Constructive Approximation}, 26(2):153--172, 2007.
\bibitem{song2010nonparametric}
Le~Song, Arthur Gretton, and Carlos Guestrin.
\newblock Nonparametric tree graphical models.
\newblock In {\em International Conference on Artificial Intelligence and
Statistics}, pages 765--772, 2010.
\bibitem{song2009hilbert}
Le~Song, Jonathan Huang, Alex Smola, and Kenji Fukumizu.
\newblock Hilbert space embeddings of conditional distributions with
applications to dynamical systems.
\newblock In {\em International Conference on Machine Learning}, pages
961--968, 2009.
\bibitem{sriperumbudur2010relation}
Bharath Sriperumbudur, Kenji Fukumizu, and Gert Lanckriet.
\newblock On the relation between universality, characteristic kernels and
{RKHS} embedding of measures.
\newblock In {\em International Conference on Artificial Intelligence and
Statistics}, pages 773--780, 2010.
\bibitem{staiger1997instrumental}
Douglas Staiger and James~H Stock.
\newblock Instrumental variables regression with weak instruments.
\newblock {\em Econometrica}, 65(3):557--586, 1997.
\bibitem{steinwart2008support}
Ingo Steinwart and Andreas Christmann.
\newblock {\em Support Vector Machines}.
\newblock Springer, 2008.
\bibitem{stock2003retrospectives}
James~H Stock and Francesco Trebbi.
\newblock Retrospectives: Who invented instrumental variable regression?
\newblock {\em Journal of Economic Perspectives}, 17(3):177--194, 2003.
\bibitem{stock2002survey}
James~H Stock, Jonathan~H Wright, and Motohiro Yogo.
\newblock A survey of weak instruments and weak identification in generalized
method of moments.
\newblock {\em Journal of Business \& Economic Statistics}, 20(4):518--529,
2002.
\bibitem{Sugiyama2010leastsquares}
Masashi Sugiyama, Ichiro Takeuchi, Taiji Suzuki, Takafumi Kanamori, Hirotaka
Hachiya, and Daisuke Okanohara.
\newblock Least-squares conditional density estimation.
\newblock {\em {IEICE} Transactions}, 93-D(3):583--594, 2010.
\bibitem{sutherland2017fix}
Dougal Sutherland.
\newblock Fixing an error in {C}aponnetto and {D}e {V}ito.
\newblock Technical report, UCL Gatsby Unit, 2017.
\bibitem{szabo2015two}
Zolt{\'a}n Szab{\'o}, Arthur Gretton, Barnab{\'a}s P{\'o}czos, and Bharath
Sriperumbudur.
\newblock Two-stage sampled learning theory on distributions.
\newblock In {\em International Conference on Artificial Intelligence and
Statistics}, pages 948--957, 2015.
\bibitem{szabo2016learning}
Zolt{\'a}n Szab{\'o}, Bharath Sriperumbudur, Barnab{\'a}s P{\'o}czos, and
Arthur Gretton.
\newblock Learning theory for distribution regression.
\newblock {\em Journal of Machine Learning Research}, 17(152):1--40, 2016.
\bibitem{wahba1990spline}
Grace Wahba.
\newblock {\em Spline Models for Observational Data}, volume~59.
\newblock SIAM, 1990.
\bibitem{Wasserman06AllOF}
Larry Wasserman.
\newblock {\em All of Nonparametric Statistics}.
\newblock Springer, 2006.
\bibitem{wright1928tariff}
Philip~G Wright.
\newblock {\em Tariff on Animal and Vegetable Oils}.
\newblock Macmillan Company, 1928.
\end{thebibliography}
\newpage