# Deep Proxy Causal Learning and its Application to Confounded Bandit Policy
Evaluation

Liyuan Xu
Gatsby Unit
<EMAIL_ADDRESS>
Heishiro Kanagawa
Gatsby Unit
<EMAIL_ADDRESS>
Arthur Gretton
Gatsby Unit
<EMAIL_ADDRESS>

###### Abstract

Proxy causal learning (PCL) is a method for estimating the causal effect of
treatments on outcomes in the presence of unobserved confounding, using
proxies (structured side information) for the confounder. This is achieved via
two-stage regression: in the first stage, we model relations among the
treatment and proxies; in the second stage, we use this model to learn the
effect of treatment on the outcome, given the context provided by the proxies.
PCL guarantees recovery of the true causal effect, subject to identifiability
conditions. We propose a novel method for PCL, the deep feature proxy variable
method (DFPV), to address the case where the proxies, treatments, and outcomes
are high-dimensional and have nonlinear complex relationships, as represented
by deep neural network features. We show that DFPV outperforms recent state-
of-the-art PCL methods on challenging synthetic benchmarks, including settings
involving high dimensional image data. Furthermore, we show that PCL can be
applied to off-policy evaluation for the confounded bandit problem, in which
DFPV also exhibits competitive performance.

##  1 Introduction

In causal learning, we aim to estimate the effect of our actions on the world.
For example, we may be interested in measuring the impact of flight ticket
prices on sales [34, 2], or the effect of grade retention on cognitive
development [4]. We refer to our action as a treatment, which results in a
particular outcome. It is often impossible to determine the effect of
treatment on outcome from observational data alone, since the observed joint
distribution of treatment and outcome can depend on a common confounder which
influences both, and which might not be observed. In our example on sales of
plane tickets given a price, the two might even be positively correlated in
some circumstances, such as the simultaneous increase in sales and prices
during the holiday season. This does not mean that raising the price causes
sales to increase. In this context, people’s desire to fly is a confounder,
since it affects both the number of ticket sales and the prices people are
willing to accept. Thus, we need to correct the bias caused by the confounder.

One common assumption to cope with confounding bias is to assume no unobserved
confounders exist [8], or more generally, the ignorable treatment assignment
assumption [28], which states that the treatment assignment is independent of
the potential outcomes caused by the treatment, given the background data
available. Although a number of methods are proposed based on this assumption
[7, 9, 36], it can be too restrictive, since it is often difficult to
determine how the confounder affects treatment assignments and outcomes.

A less restrictive assumption is that we have access to proxy variables, which
contain relevant side information on the confounder. In the flight tickets
example, we can use the number of views of the ticket reservation page as a
proxy variable, which reflects peoples’ desire for flights. Note that if we
can completely recover the confounder from proxy variables, the ignorable
treatment assignment assumption can be satisfied. Motivated by this, Lee et
al. [14] and Louizos et al. [17] aim to recover the distribution of
confounders from proxy variables using modern machine learning techniques such
as generative adversarial networks [5] or variational auto-encoders (VAE)
[11]. Although these methods exhibit powerful empirical performance, there is
little theory that guarantees the correct recovery of the causal effects.

Kuroki and Pearl [13] first considered the necessary conditions on proxy
variables to provably recover the underlying causal effect via direct recovery
of the hidden confounder. This work was in turn generalized by Miao et al.
[20]. In their work, it is shown that two types of proxy variables are
sufficient to recover the true causal effects without explicitly recovering
the confounder. One is an outcome-inducing proxy, which correlates with
confounders and causes the outcome, and the other is a treatment-inducing
proxy which is caused by confounders and correlates with the treatment. In the
flight ticket example, we can use the number of views of the ticket
reservation page as the outcome-inducing proxy, and the cost of fuel as the
treatment-inducing proxy. Given these proxy variables, Miao et al. [20] show
that the true causal effect can be recovered by solving a Fredholm integral
equation, which is referred to as the proxy causal learning (PCL) problem. The
PCL problem is known to have an interesting relation to the causal inference
with _multiple treatments_ , which uses a subset of treatments as proxy
variables in PCL [33, 22].

Although the PCL problem has a solid theoretical grounding, the question of
how to estimate the causal effect remains a practical challenge, in particular
when we consider nonlinear causal relationships or high-dimensional
treatments. In Miao et al. [20], the treatment and outcome are assumed to be
categorical variables. In a follow-up study, Miao et al. [21] show that we can
learn a linear causal effect by a method of moments. Deaner [3] models the
causal effect as a linear combination of nonlinear basis functions, which is
learned by solving two stage regression. These two methods are extended by
Mastouri and Zhu et al. [18], who learn the causal effect in a predefined
reproducing kernel Hilbert space (RKHS). We provide an overview of the PCL
problem and the two-stage regression in Section 2. Although these methods
enjoy desirable theoretical properties, the flexibility of the model is
limited, since all existing work uses pre-specified features.

In this paper, we propose a novel Deep Feature Proxy Variable (DFPV) method,
which is the first work to apply neural networks to the PCL problem. The
technique we employ builds on earlier work in instrumental variable (IV)
regression, which is a related causal inference setting to PCL. A range of
deep learning methods has recently been introduced for IV regression [6, 35,
1]. We propose to adopt the Deep Feature Instrumental Variable method [35],
which learns deep adaptive features within a two-stage regression framework.
Details of DFPV are given in Section 3. In Section 4, we empirically show that
DFPV outperforms other PCL methods in several examples. We further apply PCL
methods to the off-policy evaluation problem in a confounded bandit setting,
which aims to estimate the average reward of a new policy given data with
confounding bias. We discuss the setting in Section 3, and show the
superiority of DFPV in experiments in Section 4.

##  2 Preliminaries

In this section, we introduce the proxy causal learning (PCL) problem and
describe the existing two-stage regression methods to solve PCL.

#### Notation.

Throughout the paper, a capital letter (e.g. $A$) denotes a random variable,
and we denote the set where a random variable takes values by the
corresponding calligraphic letter (e.g. $\mathcal{A}$). The symbol
$\mathbb{P}\left(\cdot\right)$ denotes the probability distribution of a
random variable (e.g. $\mathbb{P}\left(A\right)$). We use a lowercase letter
to denote a realization of a random variable (e.g. $a$). We denote the
expectation over a random variable as $\mathbb{E}[\cdot]$ and
$\|f\|_{\mathbb{P}\left(\cdot\right)}$ as the $L^{2}$-norm of a function $f$
with respect to $\mathbb{P}\left(\cdot\right)$; i.e.
$\|f\|_{\mathbb{P}\left(A\right)}=\sqrt{\mathbb{E}_{A}\left[f^{2}(A)\right]}$.

###  2.1 Problem Setting for Proxy Causal Learning

$A$$Y$$Z$$W$$U$ Figure 1: Causal Graph.

We begin with a description of the PCL setting. We observe a treatment
$A\in\mathcal{A}$, where $\mathcal{A}\subset\mathbb{R}^{d_{A}}$, and the
corresponding outcome $Y\in\mathcal{Y}$, where $\mathcal{Y}\subset\mathbb{R}$.
We assume that there exists an unobserved confounder $U\in\mathcal{U}$ that
affects both $A$ and $Y$. The goal of PCL is to estimate the structural
function $f_{\mathrm{struct}}(X)$ defined as

$\displaystyle
f_{\mathrm{struct}}(a)=\mathbb{E}_{U}\left[\mathbb{E}_{Y}\left[Y|A=a,U\right]\right],$

which we assume to be continuous. This function is also known as the Average
Treatment Effect (ATE). The challenge of estimating $f_{\mathrm{struct}}$ is
that the confounder $U$ is not observable — we cannot estimate the structural
function from observations $A$ and $Y$ alone. To deal with this, we assume
access to a treatment-inducing proxy variable $Z$, and an outcome-inducing
proxy variable $W$, which satisfy the following structural assumption and
completeness assumption.

######  Assumption 1 (Structural Assumption [3, 18]).

We assume $Y\mathop{\perp\\!\\!\\!\perp}Z|A,U$, and
$W\mathop{\perp\\!\\!\\!\perp}(A,Z)|U$.

######  Assumption 2 (Completeness Assumption on Confounder [3, 18]).

Let $l:\mathcal{U}\to\mathbb{R}$ be any square integrable function
$\|l\|_{\mathbb{P}\left(U\right)}<\infty$. The following conditions hold:

$\displaystyle\mathbb{E}\left[l(U)\mid
A=a,W=w\right]=0\quad\forall(a,w)\in\mathcal{A}\times\mathcal{W}\quad\Leftrightarrow\quad
l(u)=0\quad\mathbb{P}\left(U\right)\text{-}\mathrm{a.e.}$
$\displaystyle\mathbb{E}\left[l(U)\mid
A=a,Z=z\right]=0\quad\forall(a,z)\in\mathcal{A}\times\mathcal{Z}\quad\Leftrightarrow\quad
l(u)=0\quad\mathbb{P}\left(U\right)\text{-}\mathrm{a.e.}$

Figure 1 shows the causal graph describing these relationships. In our
setting, we assume that there is no observable confounder, though this may be
easily included [18, 31] as presented in Appendix D. Here, the bidirectional
arrows mean that we allow both directions or even a common ancestor variable.
Given these assumptions, it is shown that the structural function can be
expressed using a solution of an integral equation.

######  Proposition 1 (Miao et al. [20]).

Let Assumptions 1, 2 and Assumptions 4, 5, 6 in Appendix B hold. Then there
exists at least one solution to the functional equation

$\displaystyle\mathbb{E}\left[Y|A=a,Z=z\right]=\int
h(a,W)\rho_{W}(w|A=a,Z=z)\mathrm{d}w,$ (1)

which holds for any $(a,z)\in\mathcal{A}\times\mathcal{Z}$. Here, we denote
$\rho_{W}(w|A=a,Z=z)$ as the density function of the conditional probability
$\mathbb{P}\left(W|A=a,Z=z\right)$. Let $h^{*}$ be a solution of (1). The
structural function $f_{\mathrm{struct}}$ is given as

$\displaystyle f_{\mathrm{struct}}(a)=\mathbb{E}_{W}\left[h^{*}(a,W)\right].$
(2)

For completeness, we present a proof in Appendix B (Lemma 2 and Corollary 1),
which is due to Miao et al. [21] and Deaner [3]. Following Miao et al. [21],
we call $h^{*}$ the bridge function. From Proposition 1, we can see that the
estimation of the structural function reduces to the estimation of the bridge
function, since once we obtain the bridge function, the structural function
directly follows from (2).

###  2.2 Two-stage Regression Scheme for Proximal Causal Learning with Fixed
Features

To obtain the bridge function, we need to solve the functional equation (1).
However, directly solving (1) in a rich function space can be ill-posed (see
discussion in Nashed and Wahba [25]). To address this, recent works [3, 18]
minimize the following regularized loss $\mathcal{L}_{\mathrm{PV}}$ to obtain
an estimate of the bridge function $\hat{h}$:

$\displaystyle\hat{h}=\operatorname*{arg\,min}_{h\in\mathcal{H}_{h}}\mathcal{L}_{\mathrm{PV}}(h),\quad\mathcal{L}_{\mathrm{PV}}(h)=\mathbb{E}_{Y,A,Z}\left[(Y-\mathbb{E}_{W|A,Z}\left[h(A,W)\right])^{2}\right]+\Omega(h),$
(3)

where $\mathcal{H}_{h}$ is an arbitrary space of continuous functions and
$\Omega(h)$ is a regularizer on $h$. Note that this loss can be interpreted as
the deviation of the r.h.s and the l.h.s in (1) measured in $L^{2}$-norm with
respect to the distribution $\mathbb{P}\left(Y,A,Z\right)$.

Deaner [3] and Mastouri and Zhu et al. [18]solve the minimization problem (3)
by successively solving two-stage regression problems. They model the bridge
function as

$h(a,w)=\bm{u}^{\top}(\bm{\psi}_{A}(a)\otimes\bm{\psi}_{W}(w))$ (4)

where $\bm{u}$ is a learnable weight vector,
$\bm{\psi}_{A}(a),\bm{\psi}_{W}(w)$ are vectors of fixed basis functions, and
$\otimes$ is a Kronecker product, defined as
$\bm{a}\otimes\bm{b}=\mathrm{vec}(\bm{a}\bm{b}^{\top})$ for any finite
dimensional vectors $\bm{a},\bm{b}$.111Throughout this paper, we assume the
number of basis functions to be finite. Mastouri and Zhu et al. [18] consider
an infinite number of basis function in a reproducing kernel Hilbert space,
and use the definitions of the inner and Kronecker products for that space. An
estimate $\hat{\bm{u}}$ is obtained by solving the successive regression
problems. In Stage 1, we estimate the conditional expectation
$\mathbb{E}_{W|A=a,Z=z}\left[\bm{\psi}_{W}(W)\right]$ as a function of $a,z$.
In Stage 2, we substitute the model (4) into the inner conditional expectation
in $\mathcal{L}_{\mathrm{PV}}$,

$\mathbb{E}_{W|A=a,Z=z}\left[h(a,W)\right]=\bm{u}^{\top}(\bm{\psi}_{A}(a)\otimes\mathbb{E}_{W|A=a,Z=z}\left[\bm{\psi}_{W}(W)\right]),$

and then minimize $\mathcal{L}_{\mathrm{PV}}$ with respect to $\bm{u}$ using
the estimate of $\mathbb{E}_{W|A=a,Z=z}\left[\bm{\psi}_{W}(W))\right]$ from
Stage 1.

The above idea can be implemented as follows. We model the conditional
expectation as

$\mathbb{E}_{W|A=a,Z=z}\left[\bm{\psi}(W)\right]=\bm{V}(\bm{\phi}_{A}(a)\otimes\bm{\phi}_{Z}(z)),$

where $\bm{\phi}_{A}(a),\bm{\phi}_{Z}(z)$ are another set of basis functions,
and $\bm{V}$ is a _matrix_ to be learned. Note that we can use different basis
functions for $\bm{\phi}_{A}(a)$ and $\bm{\psi}_{A}(a)$. In Stage 1, the
matrix $\bm{V}$ is learned by minimizing the following loss,

$\displaystyle\mathcal{L}_{1}(\bm{V})=\mathbb{E}_{W,A,Z}\left[\|\bm{\psi}_{W}(W)-\bm{V}(\bm{\phi}_{A}(A)\otimes\bm{\phi}_{Z}(Z))\|^{2}\right]+\lambda_{1}\|\bm{V}\|^{2},$
(5)

where $\lambda_{1}>0$ is a regularization parameter. This is a linear ridge
regression problem with multiple targets, which can be solved analytically. In
Stage 2, given the minimizer
$\hat{\bm{V}}=\operatorname*{arg\,min}_{\bm{V}}\mathcal{L}_{1}(\bm{V})$, we
can obtain $\hat{\bm{u}}$ by minimizing the loss

$\displaystyle\mathcal{L}_{2}(\bm{u})=\mathbb{E}_{Y,A,Z}\left[\|Y-\bm{u}^{\top}(\bm{\psi}_{A}(A)\otimes(\hat{\bm{V}}(\bm{\phi}_{A}(A)\otimes\bm{\phi}_{Z}(Z))))\|^{2}\right]+\lambda_{2}\|\bm{u}\|^{2},$
(6)

where $\lambda_{2}>0$ is another regularization parameter. Stage 2 corresponds
to another linear ridge regression from input
$\bm{\psi}_{A}(A)\otimes(\hat{\bm{V}}(\bm{\phi}_{A}(A)\otimes\bm{\phi}_{Z}(Z)))$
to target $Y$, and also enjoys a closed-form solution. Given the learned
weights
$\hat{\bm{u}}=\operatorname*{arg\,min}_{\bm{u}}\mathcal{L}_{2}(\bm{u})$, the
estimated structural function is
$\hat{f}_{\mathrm{struct}}(a)=\hat{\bm{u}}^{\top}(\bm{\psi}_{A}(a)\otimes\mathbb{E}_{W}\left[\bm{\psi}_{W}(W)\right])$.

When fixed feature dictionaries are used, this two-stage regression benefits
from strong theoretical guarantees [3, 18]. The use of pre-specified feature
maps, limits the scope and flexibility of the method, however, especially if
the treatment and proxies are high dimensional (e.g. images or text), and the
causal relations are nonlinear. To overcome this, we propose to use adaptive
features, expressed by neural nets, as described in the next section.

##  3 Deep Feature Proxy Causal Learning

In this section, we develop the DFPV algorithm, which learns adaptive features
modeled by neural nets using a technique similar to Xu et al. [35]. As in
Mastouri and Zhu et al. [18], we assume that we do not necessarily have access
to observations from the joint distribution of $(A,Y,Z,W)$. Instead, we are
given $m$ observations of $(A,Z,W)$ for Stage 1 and $n$ observations of
$(A,Z,Y)$ for Stage 2. We denote the stage 1 observations by
$(a_{i},z_{i},w_{i})$ and the stage 2 observations by
$(\tilde{a}_{i},\tilde{z}_{i},\tilde{y}_{i})$. If observations of $(A,Y,Z,W)$
are given for both stages, we can evaluate the out-of-sample loss of Stage 1
using Stage 2 data and vice versa, and these losses can be used for hyper-
parameter tuning of $\lambda_{1},\lambda_{2}$ (Appendix A). We first introduce
two-stage regression with adaptive feature maps and then describe the detailed
learning procedure of DFPV.

###  3.1 Two-stage regression with adaptive features

In DFPV, we consider the following models of the bridge function $h(a,w)$ and
conditional feature mean
$\mathbb{E}_{W|A,Z}\left[\bm{\psi}_{\theta_{W}}(W)\right]$:

$\displaystyle
h(a,w)=\bm{u}^{\top}(\bm{\psi}_{\theta_{A(2)}}(a)\otimes\bm{\psi}_{\theta_{W}}(w)),\quad\mathbb{E}_{W|A=a,Z=z}\left[\bm{\psi}_{\theta_{W}}(W)\right]=\bm{V}(\bm{\phi}_{\theta_{A(1)}}(a)\otimes\bm{\phi}_{\theta_{Z}}(z)),$

where $\bm{u}$ and $\bm{V}$ are parameters, and
$\bm{\phi}_{\theta_{A(1)}},\bm{\phi}_{\theta_{Z}},\bm{\psi}_{\theta_{A(2)}},\bm{\psi}_{\theta_{W}}$
are neural nets parametrized by
$\theta_{A(1)},\theta_{Z},\theta_{A(2)},\theta_{W}$, respectively. Again, we
may use different neural nets in the treatment features
$\bm{\phi}_{\theta_{A(1)}}$ and $\bm{\psi}_{\theta_{A(2)}}$.

As in the existing work [3, 18], we learn
$\mathbb{E}_{W|a,z}\left[\bm{\psi}_{\theta_{W}}(w)\right]$ in Stage 1 and
$h(a,w)$ in Stage 2, but in addition to the weights $\bm{u}$ and $\bm{V}$, we
also learn the parameters of the feature maps. Specifically, in Stage 1, we
learn $\bm{V}$ and parameters $\theta_{A(1)},\theta_{Z}$ by minimizing the
following empirical loss:

$\displaystyle\hat{\mathcal{L}}_{1}(\bm{V},\theta_{A(1)},\theta_{Z})=\frac{1}{m}\sum_{i=1}^{m}\left\|\bm{\psi}_{\theta_{W}}(w_{i})-\bm{V}\left(\bm{\phi}_{\theta_{A(1)}}(a_{i})\otimes\bm{\phi}_{\theta_{Z}}(z_{i})\right)\right\|^{2}+\lambda_{1}\|\bm{V}\|^{2}.$

Note that $\hat{\mathcal{L}}_{1}$ is an empirical estimate of
$\mathcal{L}_{1}$ in (5) with adaptive feature maps. Although
$\hat{\mathcal{L}}_{1}$ depends on $\theta_{W}$, at this stage, we do not
optimize $\theta_{W}$ with $\hat{\mathcal{L}}_{1}$ as
${\bm{\psi}_{\theta_{W}}(w)}$ is the “target variable” in Stage 1. Given the
minimizers
$(\hat{\bm{V}},\hat{\theta}_{A(1)},\hat{\theta}_{Z})=\operatorname*{arg\,min}\hat{\mathcal{L}}_{1}$,
we learn weights $\bm{u}$ and parameters $\theta_{W},\theta_{A(2)}$ by
minimizing the empirical stage 2 loss,

$\displaystyle\hat{\mathcal{L}}_{2}(\bm{u},\theta_{W},\theta_{A(2)})=\frac{1}{n}\sum_{i=1}^{n}\left(\tilde{y}_{i}-\bm{u}^{\top}\left(\bm{\psi}_{\theta_{A(2)}}(\tilde{a}_{i})\otimes\hat{\bm{V}}\left(\bm{\phi}_{\hat{\theta}_{A(1)}}(\tilde{a}_{i})\otimes\bm{\phi}_{\hat{\theta}_{Z}}(\tilde{z}_{i})\right)\right)\right)^{2}+\lambda_{2}\|\bm{u}\|^{2}.$

Again, the loss $\hat{\mathcal{L}}_{2}$ is an empirical estimate of
$\mathcal{L}_{2}$ in (6) with adaptive feature maps. Although the expression
of $\hat{\mathcal{L}}_{2}$ does not explicitly contain $\theta_{W}$, it
implicitly depends on $\theta_{W}$ as
$(\hat{\bm{V}},\hat{\theta}_{A(1)},\hat{\theta}_{Z})$ is the solution of a
minimization problem involving $\theta_{W}$. This implicit dependency makes it
challenging to update $\theta_{W}$, as we cannot directly obtain its gradient.
One possible solution is to use the implicit gradient method [16], but this
approach might be too computationally expensive. Instead, we use the method
proposed in Xu et al. [35], in which we ignore the dependency of $\theta_{W}$
on parameters $\hat{\theta}_{A(1)},\hat{\theta}_{Z}$, and compute the gradient
via the closed-form solution of $\hat{V}$.

###  3.2 Deep Feature Proxy Variable Method

We now describe the learning procedure for DFPV. First, we fix parameters in
the adaptive feature maps
$(\theta_{A(1)},\theta_{Z},\theta_{A(2)},\theta_{W})$. Given these parameters,
optimal weights $\hat{\bm{V}},\hat{\bm{u}}$ can be learned by minimizing the
empirical stage 1 loss $\hat{\mathcal{L}}_{1}$ and empirical stage 2 loss
$\hat{\mathcal{L}}_{2}$, respectively. These minimizations can be solved
analytically, where the solutions are

$\displaystyle\hat{\bm{V}}(\bm{\theta})=\Psi_{1}^{\top}\Phi_{1}(\Phi_{1}^{\top}\Phi_{1}+m\lambda_{1}I)^{-1},$
$\displaystyle\hat{\bm{u}}(\bm{\theta})=\left(\Phi_{2}^{\top}\Phi_{2}+n\lambda_{2}I\right)^{-1}\Phi_{2}^{\top}\bm{y}_{2},$
(7)

where we denote
$\bm{\theta}=(\theta_{A(1)},\theta_{Z},\theta_{A(2)},\theta_{W})$ and define
matrices as follows:

$\displaystyle\Psi_{1}=\left[\bm{\psi}_{\theta_{W}}(w_{1}),\dots,\bm{\psi}_{\theta_{W}}(w_{m})\right]^{\top},$
$\displaystyle\Phi_{1}=[\bm{v}_{1}(a_{1},z_{1}),\dots,\bm{v}_{1}(a_{m},z_{m})]^{\top},$
$\displaystyle\bm{y}_{2}=[\tilde{y}_{1},\dots,\tilde{y}_{n}]^{\top},$
$\displaystyle\Phi_{2}=[\bm{v}_{2}(\tilde{a}_{1},\tilde{z}_{1}),\dots,\bm{v}_{2}(\tilde{a}_{n},\tilde{z}_{n})]^{\top},$
$\displaystyle\bm{v}_{1}(a,z)=\bm{\phi}_{\theta_{A(1)}}(a)\otimes\bm{\phi}_{\theta_{Z}}(z),$
$\displaystyle\bm{v}_{2}(a,z)=\bm{\psi}_{\theta_{A(2)}}(a)\otimes\left(\hat{\bm{V}}(\bm{\theta})\left(\bm{\phi}_{\theta_{A(1)}}(a)\otimes\bm{\phi}_{\theta_{Z}}(z)\right)\right).$

Given these weights $\hat{\bm{u}}(\bm{\theta}),\hat{\bm{V}}(\bm{\theta})$, we
can update feature parameters by a gradient descent method with respect to the
residuals of the loss of each stage, while regrading $\hat{\bm{V}}$ and
$\hat{\bm{u}}$ as functions of parameters. Specifically, we take the gradient
of the losses

$\displaystyle\hat{\mathcal{L}}^{\mathrm{DFPV}}_{1}(\bm{\theta})=\frac{1}{m}\sum_{i=1}^{m}\left\|\bm{\psi}_{\theta_{W}}(w_{i})-\hat{\bm{V}}(\bm{\theta})\left(\bm{\phi}_{\theta_{A(1)}}(a_{i})\otimes\bm{\phi}_{\theta_{Z}}(z_{i})\right)\right\|^{2}+\lambda_{1}\|\hat{\bm{V}}(\bm{\theta})\|^{2},$
$\displaystyle\hat{\mathcal{L}}^{\mathrm{DFPV}}_{2}(\bm{\theta})=\frac{1}{n}\sum_{i=1}^{n}\left(\tilde{y}_{i}-\hat{\bm{u}}(\bm{\theta})^{\top}\left(\bm{\psi}_{\theta_{A(2)}}(\tilde{a}_{i})\otimes\hat{\bm{V}}(\bm{\theta})\left(\bm{\phi}_{\theta_{A(1)}}(\tilde{a}_{i})\otimes\bm{\phi}_{\theta_{Z}}(\tilde{z}_{i})\right)\right)\right)^{2}\\!+\\!\lambda_{2}\|\hat{\bm{u}}(\bm{\theta})\|^{2},$

where $\hat{\bm{V}}(\bm{\theta}),\hat{\bm{u}}(\bm{\theta})$ are given in (7).
Given these losses, $(\theta_{A(1)},\theta_{Z})$ are minimized with respect to
$\hat{\mathcal{L}}^{\mathrm{DFPV}}_{1}(\bm{\theta})$, and
$(\theta_{A(2)},\theta_{W})$ are minimized with respect to
$\hat{\mathcal{L}}^{\mathrm{DFPV}}_{2}(\bm{\theta})$. Finally, we take the
empirical mean of $\bm{\psi}_{\theta^{(t)}_{W}}$ based on additional output-
proxy data $S_{W}=\\{w^{\mathrm{extra}}_{i}\\}_{i=1}^{n_{W}}$, which is used
for estimating the structural function. Here, we assume access to $S_{W}$ for
proving consistency results, but empirically, we can use stage 1 data to
compute this mean.

The complete procedure is presented in Algorithm 1. Note that we may use any
sophisticated gradient-based learning method to optimize, such as Adam [10].
As reported in Xu et al. [35], we observe that the learning procedure is
stabilized by running several gradient descent steps on the stage 1 parameters
$(\theta_{A(1)},\theta_{Z})$ before updating the stage 2 features
$(\theta_{A(2)},\theta_{W})$. Furthermore, we may use mini-batch updates,
which sample subsets of the data at the beginning of each iteration and only
use these subsamples to update the parameters.

Algorithm 1 Deep Feature Proxy Causal Learning

0: Stage 1 data $S_{1}=\\{a_{i},z_{i},w_{i}\\}$, Stage 2 data
$S_{2}=\\{\tilde{a}_{i}\tilde{z}_{i},\tilde{y}_{i}\\}$, Additional outcome-
proxy data $S_{W}=\\{w^{\mathrm{extra}}_{i}\\}$, Regularization parameters
$(\lambda_{1},\lambda_{2})$. Initial values
$\bm{\theta}^{(0)}=(\theta^{(0)}_{A(1)},\theta^{(0)}_{Z},\theta^{(0)}_{A(2)},\theta^{(0)}_{W})$.
Learning rate $\alpha$.

0: Estimated structural function $\hat{f}_{\mathrm{struct}}(a)$

1: $t\leftarrow 0$

2: repeat

3: Compute $\hat{\bm{V}}(\bm{\theta}^{(t)}),\hat{\bm{u}}(\bm{\theta}^{(t)})$
in (7)

4: Update parameters in features
$\bm{\theta}^{(t+1)}\leftarrow(\theta^{(t+1)}_{A(1)},\theta^{(t+1)}_{Z},\theta^{(t+1)}_{A(2)},\theta^{(t+1)}_{W})$
as follows
$\displaystyle\theta^{(t+1)}_{A(1)}\leftarrow\theta^{(t)}_{A(1)}-\alpha\nabla_{\theta_{A(1)}}\hat{\mathcal{L}}^{\mathrm{DFPV}}_{1}(\bm{\theta})|_{\bm{\theta}=\bm{\theta}^{(t)}},$
$\displaystyle\theta^{(t+1)}_{Z}\leftarrow\theta^{(t)}_{Z}-\alpha\nabla_{\theta_{Z}}\hat{\mathcal{L}}^{\mathrm{DFPV}}_{1}(\bm{\theta})|_{\bm{\theta}=\bm{\theta}^{(t)}}$
$\displaystyle\theta^{(t+1)}_{A(2)}\leftarrow\theta^{(t)}_{A(2)}-\alpha\nabla_{\theta_{A(2)}}\hat{\mathcal{L}}^{\mathrm{DFPV}}_{2}(\bm{\theta})|_{\bm{\theta}=\bm{\theta}^{(t)}},$
$\displaystyle\theta^{(t+1)}_{W}\leftarrow\theta^{(t)}_{W}-\alpha\nabla_{\theta_{W}}\hat{\mathcal{L}}^{\mathrm{DFPV}}_{2}(\bm{\theta})|_{\bm{\theta}=\bm{\theta}^{(t)}}$

5: Increment counter $t\leftarrow t+1;$

6: until convergence

7: Compute $\hat{\bm{u}}(\bm{\theta}^{(t)})$ from (7)

8: Compute mean feature for $W$ using stage 1 dataset:
$\bm{\mu}_{\theta_{W}}\leftarrow\frac{1}{n}\sum\bm{\psi}_{\theta^{(t)}_{W}}(w^{\mathrm{extra}}_{i})$

9: return
$\hat{f}_{\mathrm{struct}}(a)=(\hat{\bm{u}}^{(t)})^{\top}\left(\bm{\psi}_{\hat{\theta}^{(t)}_{A(2)}}(a)\otimes\bm{\mu}_{\theta_{W}}\right)$

###  3.3 Consistency of DFPV

In this section, we show that the DFPV method yields a consistent estimate of
the bridge function. Our approach differs from the consistency result of
Deaner [3] and Mastouri and Zhu et al. [18], which depend on the assumption
that the true bridge function lies in a specific functional space (such as an
RKHS [18]) or satisfies certain smoothness properties [3]. Since our features
are adaptive, however, we instead build our result on a Rademacher complexity
argument [24], which holds for a wide variety of hypothesis spaces.
Furthermore, we do not assume that our hypothesis space contains the true
bridge function, which makes this theory of independent interest.

To derive the excess risk bound, we need another completeness assumption.

######  Assumption 3 (Completeness Assumption on Outcome-Inducing Proxy [3,
18]).

Let $l:\mathcal{W}\to\mathbb{R}$ be any square integrable function
$\|l\|_{\mathbb{P}\left(W\right)}<\infty$, We assume the following condition:

$\displaystyle\mathbb{E}\left[l(W)\mid
A=a,Z=z\right]=0\quad\forall(a,z)\in\mathcal{A}\times\mathcal{Z}\quad\Leftrightarrow\quad
l(w)=0\quad\mathbb{P}\left(W\right)\text{-}\mathrm{a.e.}$

This assumption is not necessary for identification (Proposition 1), but we
need it for connecting the two-stage loss and the deviation of the bridge
function. Given this assumption, we have the following consistency result.

######  Proposition 2.

Let Assumption 3 and Assumptions 7 and 8 in Appendix C hold. Given stage 1
data $S_{1}=\\{(w_{i},a_{i},z_{i})\\}_{i=1}^{m}$, stage 2 data
$S_{2}=\\{(\tilde{y}_{i},\tilde{a}_{i},\tilde{z}_{i})\\}_{i=1}^{n}$, and
additional output-proxy data
$S_{W}=\\{w_{i}^{\mathrm{extra}}\\}_{i=1}^{n_{W}}$, then for the minimizer of
$\hat{\mathcal{L}}_{2}$ denoted as
$(\hat{\bm{u}},\hat{\theta}_{A(2)},\hat{\theta}_{Z})$, we have

$\displaystyle\|f_{\mathrm{struct}}-\hat{f}_{\mathrm{struct}}\|_{\mathbb{P}\left(A\right)}\leq
O_{p}\left(\sqrt{\kappa_{1}+\hat{\mathfrak{R}}_{S_{1}}(\mathcal{H}_{1})+\frac{1}{\sqrt{m}}}+\kappa_{2}\\!+\\!\sqrt{\hat{\mathfrak{R}}_{S_{2}}(\mathcal{H}_{2})+\frac{1}{\sqrt{n}}}+\frac{1}{\sqrt{n_{W}}}\right),$

where $\mathcal{H}_{1},\mathcal{H}_{2}$ are the functional classes that are
considered in each stage, comprising both the predictions and their associated
losses, which are defined in (12) and (15), respectively; and
$\hat{\mathfrak{R}}_{S_{1}}(\mathcal{H}_{1}),\hat{\mathfrak{R}}_{S_{2}}(\mathcal{H}_{2})$
are their respective empirical Rademacher complexities; the quantities
$\kappa_{1},\kappa_{2}$ measure the misspecification in each regression, and
are defined in (13) and (16), respectively.

We formally restate Proposition 2 in Appendix C (Theorem 3), in which the
definitions of hypothesis spaces and constants can be found. From this, we can
see that if we correctly identify the hypothesis space
($\kappa_{1},\kappa_{2}=0$) and have vanishing Rademacher complexities
($\hat{\mathfrak{R}}_{S_{1}}(\mathcal{H}_{1}),\hat{\mathfrak{R}}_{S_{2}}(\mathcal{H}_{2})\to
0$ in probability), then the estimated bridge function converges in
probability to the true one. Note that the Rademacher complexity of certain
neural networks is known; e.g., ReLU two-layer networks are known to have the
complexity of order $O(\sqrt{1/n})$ [26]. However, it might be not
straightforward to derive the Rademacher complexities of the specific classes
of neural networks we use, because we are employing the outer product of two
neural networks in (5) and (6), whose Rademacher complexity remains a
challenging open problem. The analysis in Proposition 2 is similar to Xu et
al. [35], with two important distinctions: first, the target
$f_{\mathrm{struct}}$ is estimated using the bridge function $h$ marginalized
over one of its arguments $W$, but PCL learns the bridge function itself.
Second, our result is more general and includes the misspecified case, which
was not treated in the earlier work.

One limitation of our analysis is that we leave aside questions of
optimization. As we previously discussed, $\hat{\mathcal{L}}_{2}$ does not
explicitly include $\theta_{W}$, thus it is difficult to develop a consistency
result that includes the optimization. We emphasize that Algorithm 1 does not
guarantee to converge to the global minimizer, since we ignore the implicit
dependency of $\theta_{W}$ on $\theta_{A(1)},\theta_{Z}$ when calculating the
gradient. As we will see in the Section 4, however, the proposed method
outperforms competing methods with fixed feature dictionaries, which do not
have this issue.

###  3.4 DFPV for Policy Evaluation

PCL methods can estimate not only the structural function, but also other
causal quantities, such as Conditional Average Treatment Effect (CATE) [29]
and Average Treatment Effect on Treated (ATT) [3, 29]. Beyond the causal
context, [32] uses a PCL method to conduct the off-policy evaluation in a
partially observable Markov decision process with discrete action-state
spaces. Here, we present the application of PCL methods to bandit off-policy
evaluation with continuous actions, and propose a method to leverage adaptive
feature maps in this context.

In bandit off-policy evaluation, we regard the outcome $Y$ as the reward and
aim to evaluate the value function of a given policy. Denote policy
$\pi:\mathcal{C}\to\mathcal{A}$ and define the value function $v(\pi)$ as

$\displaystyle
v(\pi)=\mathbb{E}_{U,C}\left[\mathbb{E}\left[Y|A=\pi(C),U\right]\right],$

where $C\in\mathcal{C}$ is a variable on which the policy depends. In the
flight ticket sales prediction example, a company might be interested in
predicting the effect of discount offers. This requires us to solve the policy
evaluation with $C=A$ since a new price is determined based on the current
price. Alternatively, if the company is planning to introduce a new price
policy that depends on the fuel cost, we need to consider the policy
evaluation with $C=Z$, since the fuel cost can be regarded as the treatment-
inducing proxy variable. Here, we assume the policy $\pi$ to be deterministic,
but this can be easily generalized to a stochastic policy. Although the value
function $v$ contains the expectation over the latent confounder $U$, it can
be estimated through a bridge function.

######  Proposition 3.

Assume that the true bridge function
$h^{*}(a,w):\mathcal{A}\times\mathcal{W}\to\mathbb{R}$ is jointly measurable.
Suppose $C\mathop{\perp\\!\\!\\!\perp}W|U$. Under Assumptions 1, 2 and
Assumptions 4, 5, 6 in Appendix B, we have

$\displaystyle v(\pi)=\mathbb{E}_{C,W}\left[h^{*}(\pi(C),W)\right].$

The proof can be found in Appendix B (Corollary 2). From Proposition 3, we can
obtain the value function by taking the empirical average over the bridge
function $h(\pi(C),W)$. In DFPV, we estimate the bridge function as
$h(a,w)=(\hat{\bm{u}})^{\top}(\bm{\psi}_{\hat{\theta}_{A(2)}}(a)\otimes\bm{\psi}_{\hat{\theta}_{W}}(w))$
where $\hat{\bm{u}},\hat{\theta}_{A(2)},\hat{\theta}_{W}$ minimize
$\hat{\mathcal{L}}_{2}$. Hence, given $n^{\prime}$ observation of $(C,W)$
denoted as $\\{\check{c}_{i},\check{w}_{i}\\}_{i=1}^{n^{\prime}}$, we can
estimate the value function $v(\pi)$ as

$\displaystyle\hat{v}(\pi)=\frac{1}{n_{C}}\sum_{i=1}^{n^{\prime}}\hat{\bm{u}}^{\top}\left(\bm{\psi}_{\hat{\theta}_{A(2)}}(\pi(\check{c}_{i}))\otimes\bm{\psi}_{\hat{\theta}_{W}}(\check{w}_{i}))\right).$

We derive the following consistency result for policy value evaluation in a
bandit with confounding bias.

######  Proposition 4.

Let Assumption 3 and Assumptions 7, 8, 9 in Appendix C hold. Denote the
minimizers of $\hat{\mathcal{L}}_{2}$ as
$(\hat{u},\hat{\theta}_{A(2)},\hat{\theta}_{Z})$. Given stage 1 data
$S_{1}=\\{(w_{i},a_{i},z_{i})\\}_{i=1}^{m}$, stage 2 data
$S_{2}=\\{(\tilde{y}_{i},\tilde{a}_{i},\tilde{z}_{i})\\}_{i=1}^{n}$, and data
for policy evaluation
$S_{3}=\\{(\check{c}_{i},\check{w}_{i})\\}_{i=1}^{n^{\prime}}$, we have

$\displaystyle\\!|v(\pi)-\hat{v}(\pi)|\leq
O_{p}\left(\sqrt{\kappa_{1}+\hat{\mathfrak{R}}_{S_{1}}(\mathcal{H}_{1})+\frac{1}{\sqrt{m}}}+\kappa_{2}\\!+\\!\sqrt{\hat{\mathfrak{R}}_{S_{2}}(\mathcal{H}_{2})+\frac{1}{\sqrt{n}}}+\frac{1}{\sqrt{n^{\prime}}}\right)$

where $\mathcal{H}_{1},\mathcal{H}_{2}$ are the function classes that are
considered in each stage, comprising both the predictions and their associated
losses, which are defined in (12) and (15), respectively; and their empirical
Rademacher complexities are given by
$\hat{\mathfrak{R}}_{S_{1}}(\mathcal{H}_{1}),\hat{\mathfrak{R}}_{S_{2}}(\mathcal{H}_{2})$;
the quantities $\kappa_{1},\kappa_{2}$ measure the misspecification in each
regression, which are defined in (13) and (16), respectively.

The formal statement of Proposition 4 can be found in Appendix C (Theorem 3).
As in the structural function case, we can see that estimated value function
converges in probability to the true one, provided that we have vanishing
Rademacher complexities and there is no misspecification.

##  4 Experiments

In this section, we report the empirical performance of the DFPV method.
First, we present the results of estimating structural functions; we design
two experimental settings for low-dimensional treatments and high-dimensional
treatments, respectively. Then, we show the result of applying PCL methods to
the bandit off-policy evaluation problem with confounding. We include the
results for problems considered in prior work in Appendix E. The experiments
are implemented using PyTorch [27]. The code is included in the supplemental
material. All experiments can be run in a few minutes on Intel(R) Xeon(R) CPU
E5-2698 v4 2.20GHz.

#### Experiments for Structural Function

We present two structural function estimation experiments. One is a demand
design experiment based on a synthetic dataset introduced by Hartford et al.
[6], which is a standard benchmark for the instrumental variable regression.
Here, we modify the data generating process to provide a benchmark for PCL
methods. We consider the problem of predicting sales $Y$ from ticket price
$P$, where these are confounded by a potential demand $D\in[0,10]$. To correct
this confounding bias, we use the fuel price $(C_{1},C_{2})$ as the treatment-
inducing proxy, which has an impact on price $P$, and the number of views of
the ticket reservation page $V$ as the outcome-inducing proxy. Details of the
data generation process can be found in Appendix F.1.

Our second structural function estimation experiment considers high-
dimensional treatment variables. We test this using the dSprite dataset [19],
which is an image dataset described by five latent parameters (shape, scale,
rotation, posX and posY). The images are $64\times 64=4096$-dimensional. Based
on this, Xu et al. [35] introduced the causal experiment, where the treatment
is each figure, and the confounder is posY. Inspired by this, we consider the
PCL setting that learns the same structural functions with nonlinear
confounding, which is not possible to handle in the instrumental variable
setting. Specifically, the structural function $f_{\mathrm{struct}}$ and
outcome $Y$ are defined as

$\displaystyle
f_{\mathrm{struct}}(a)=\frac{\|\bm{B}a\|_{2}^{2}-5000}{1000},\quad
Y=\frac{1}{12}(\mathtt{posY}-0.5)f_{\mathrm{struct}}(A)+\varepsilon,\quad\varepsilon\sim\mathcal{N}(0,0.5),$

where each element of the matrix $\bm{B}\in\mathbb{R}^{10\times 4096}$ was
generated from $\mathrm{Unif}(0.0,1.0)$ and fixed throughout the experiment.
We fixed the shape parameter to heart and used other parameters as the
treatment-inducing proxy $Z$. We sampled another image that shared the same
posY as treatment $A$, which is used as output-inducing proxy $W$. Details of
data generation process can be found in Appendix F.2.

$A$$Y$$Q$$U$ Figure 2: Causal Graph in CEVAE

We compare the DFPV method to three competing methods, namely KPV [18], PMMR
[18], and an autoencoder approach derived from the CEVAE method [17]. In KPV,
the bridge function is estimated through the two-stage regression as described
in Section 2, where feature functions are fixed via their kernel functions.
PMMR also models the bridge function using kernel functions, but parameters
are learned by moment matching. CEVAE is not a PCL method, however, it
represents a state-of-the-art approach in correcting for hidden confounders
using observed proxies. The causal graph for CEVAE is shown in Figure 2, and
CEVAE uses a VAE [11] to recover the distribution of confounder $U$ from the
“proxy” $Q$. We make two modifications to CEVAE to apply it in our setting.
First, we include both the treatment-inducing proxy $Z$ and output-inducing
proxy $W$ as $Q$ in CEVAE (we emphasize that this does not follow the causal
graph in Figure 2, since there exist arrows from $Q$ to $A,Y$). Second, CEVAE
is originally used in the setting where $Q$ is conditioned on a particular
value, whereas we marginalize $Q$. See Appendix F.4 for the choice of the
network structure and hyper-parameters. We tuned the regularizers
$\lambda_{1},\lambda_{2}$ as discussed in Appendix A, with the data evenly
split for Stage 1 and Stage 2. We varied the dataset size and ran 20
simulations for each setting. Results are summarized in Figure 4.

Figure 3: Result of structural function experiment in demand design setting
(Left) and dSprite setting (Right).

Figure 4: Result of OPE experiment when the policy depends on the costs (Left)
and on the current price (Right).

In both experiments, DFPV consistently outperforms existing methods. This
suggests that DFPV is capable of learning complex structural functions by
taking the advantage of the flexibility of neural networks. KPV and PMMR
perform similarly in all settings, but KPV tends to perform slightly better
when the data size is large. This might be because KPV expresses its solution
using a larger number of parameters than PMMR. Although CEVAE also learns a
flexible model with a neural network, it’s highly unstable in the demand
design experiment and underperforms kernel methods in the dSprite experiment.
This is because CEVAE does not take advantage of the relations between proxy
variables and the structural function.

#### Experiments for Offline Policy Evaluation

We now describe the offline policy evaluation experiment based on the demand
design data. We set up synthetic experiments to evaluate two types of policy.
In the first case, a policy depends on costs $C=(C_{1},C_{2})$, to address the
question How much would we gain if we decide the price based on the fuel
costs. In the second case, the policy depends on current action $C=P$, which
addresses How much would we gain if we cut the price by 30 percent. The
detailed policies can be found in Appendix F.3.

We only compare PCL methods here, as it is not straightforward to apply CEVAE
to the off-policy evaluation problem. We evenly split the data for Stage 1,
Stage 2, and policy evaluation (i.e we set $n=m=n^{\prime}$).We ran 20
simulations for each setting. Results are summarized in Figure 4, in which
DFPV performs better than existing PCL methods. This is not surprising since,
as demonstrated in the structural function experiments, DFPV can estimate
complex bridge functions, which results in a more accurate estimation of the
value function.

##  5 Conclusion

We have proposed a novel approach for proxy causal learning, the Deep Feature
Proxy Variable (DFPV) method, which performs two-stage least squares
regression on flexible and expressive features. Motivated by the literature on
the instrumental variable problem, we showed how to learn these feature maps
adaptively with deep neural networks. We also showed that PCL learning can be
used for off-policy evaluation in the bandit setting with confounding, and
that DFIV performs competitively in this domain. This work thus brings
together research from the worlds of deep offline RL and causality from
observational data.

In future work, it would be interesting to explore different approaches to
learning deep models in the PCL problem. One direction would be to adapt the
method of moment matching, which has been studied extensively in the
instrumental variable literature [1]. Moreover, Liao et al. [15] recently
developed a novel adversarial method to solve a class of functional equations,
which would be also a promising approach in the PCL problem. In the RL
context, Tennenholtz et al. [32] shows an interesting connection between PCL
learning and off-policy evaluation in a partially observable Markov decision
process, and we believe that adapting DFPV to this setting will be of value.

## References

  * Bennett et al. [2019] A. Bennett, N. Kallus, and T. Schnabel.  Deep generalized method of moments for instrumental variable analysis.  In _Advances in Neural Information Processing Systems_ , volume 32, 2019.
  * Blundell et al. [2012] R. Blundell, J. Horowitz, and M. Parey.  Measuring the price responsiveness of gasoline demand: Economic shape restrictions and nonparametric demand estimation.  _Quantitative Economics_ , 3:29–51, 2012.
  * Deaner [2018] B. Deaner.  Proxy controls and panel data, 2018.
  * Fruehwirth et al. [2016] J. Fruehwirth, S. Navarro, and Y. Takahashi.  How the timing of grade retention affects outcomes: Identification and estimation of time-varying treatment effects.  _Journal of Labor Economics_ , 34(4):979 – 1021, 2016.
  * Goodfellow et al. [2014] I. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, and Y. Bengio.  Generative adversarial nets.  In _Advances in Neural Information Processing Systems_ , volume 27, 2014.
  * Hartford et al. [2017] J. Hartford, G. Lewis, K. Leyton-Brown, and M. Taddy.  Deep IV: A flexible approach for counterfactual prediction.  In _Proceedings of The 34th International Conference on Machine Learning_ , 2017.
  * Hill [2011] J. Hill.  Bayesian nonparametric modeling for causal inference.  _Journal of Computational and Graphical Statistics_ , 20:217–240, 2011.
  * Imbens [2003] G. W. Imbens.  Nonparametric estimation of average treatment effects under exogeneity: A review.  Working Paper 294, National Bureau of Economic Research, October 2003\.
  * Johansson et al. [2016] F. Johansson, U. Shalit, and D. Sontag.  Learning representations for counterfactual inference.  In _Proceedings of The 33rd International Conference on Machine Learning_ , 2016.
  * Kingma and Ba [2015] D. P. Kingma and J. Ba.  Adam: A method for stochastic optimization.  In _International Conference on Learning Representations_ , 2015.
  * Kingma and Welling [2014] D. P. Kingma and M. Welling.  Auto-encoding variational Bayes.  In _International Conference on Learning Representations_ , 2014.
  * Kress [1999] R. Kress.  _Linear Integral Equations_.  Applied Mathematical Sciences. Springer New York, 1999.
  * Kuroki and Pearl [2014] M. Kuroki and J. Pearl.  Measurement bias and effect restoration in causal inference.  _Biometrika_ , 101(2):423–437, 2014.
  * Lee et al. [2018] C. Lee, N. Mastronarde, and M. van der Schaar.  Estimation of individual treatment effect in latent confounder models via adversarial learning.  In _Advances in Neural Information Processing Systems_ , volume 32, 2018.
  * Liao et al. [2020] L. Liao, Y.-L. Chen, Z. Yang, B. Dai, M. Kolar, and Z. Wang.  Provably efficient neural estimation of structural equation models: An adversarial approach.  In _Advances in Neural Information Processing Systems_ , volume 33, 2020.
  * Lorraine et al. [2020] J. Lorraine, P. Vicol, and D. Duvenaud.  Optimizing millions of hyperparameters by implicit differentiation.  In _Proceedings of the 23rd International Conference on Artificial Intelligence and Statistics_ , 2020.
  * Louizos et al. [2017] C. Louizos, U. Shalit, J. M. Mooij, D. Sontag, R. S. Zemel, and M. Welling.  Causal effect inference with deep latent-variable models.  In _Advances in Neural Information Processing Systems_ , volume 31, 2017.
  * Mastouri et al. [2021] A. Mastouri, Y. Zhu, L. Gultchin, A. Korba, R. Silva, M. Kusner, A. Gretton, and K. Muandet.  Proximal causal learning with kernels: Two-stage estimation and moment restriction.  In _Proceedings of The 38th International Conference on Machine Learning_ , 2021.
  * Matthey et al. [2017] L. Matthey, I. Higgins, D. Hassabis, and A. Lerchner.  dSprites: Disentanglement testing sprites dataset, 2017.  URL https://github.com/deepmind/dsprites-dataset/.
  * Miao et al. [2018a] W. Miao, Z. Geng, and E. Tchetgen Tchetgen.  Identifying causal effects with proxy variables of an unmeasured confounder.  _Biometrika_ , 105(4):987—993, 2018a.
  * Miao et al. [2018b] W. Miao, X. Shi, and E. T. Tchetgen.  A confounding bridge approach for double negative control inference on causal effects, 2018b.
  * Miao et al. [2020] W. Miao, W. Hu, E. L. Ogburn, and X. Zhou.  Identifying effects of multiple treatments in the presence of unmeasured confounding, 2020.
  * Miyato et al. [2018] T. Miyato, T. Kataoka, M. Koyama, and Y. Yoshida.  Spectral normalization for generative adversarial networks.  In _International Conference on Learning Representations_ , 2018.
  * Mohri et al. [2012] M. Mohri, A. Rostamizadeh, and A. Talwalkar.  _Foundations of Machine Learning_.  MIT Press, 2012.
  * Nashed and Wahba [1974] M. Z. Nashed and G. Wahba.  Generalized inverses in reproducing kernel spaces: An approach to regularization of linear operator equations.  _SIAM Journal on Mathematical Analysis_ , 5(6):974–987, 1974.
  * Neyshabur et al. [2019] B. Neyshabur, Z. Li, S. Bhojanapalli, Y. LeCun, and N. Srebro.  The role of over-parametrization in generalization of neural networks.  In _International Conference on Learning Representations_ , 2019.  URL https://openreview.net/forum?id=BygfghAcYX.
  * Paszke et al. [2019] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Kopf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, and S. Chintala.  Pytorch: An imperative style, high-performance deep learning library.  In _Advances in Neural Information Processing Systems_ , volume 32, 2019.
  * Rosenbaum and Rubin [1983] P. R. Rosenbaum and D. B. Rubin.  The central role of the propensity score in observational studies for causal effects.  _Biometrika_ , 70(1):41–55, 1983.
  * Singh [2020] R. Singh.  Kernel methods for unobserved confounding: Negative controls, proxies, and instruments, 2020.
  * Singh et al. [2019] R. Singh, M. Sahani, and A. Gretton.  Kernel instrumental variable regression.  In _Advances in Neural Information Processing Systems_ , volume 32, 2019.
  * Tchetgen et al. [2020] E. J. T. Tchetgen, A. Ying, Y. Cui, X. Shi, and W. Miao.  An introduction to proximal causal learning, 2020.
  * Tennenholtz et al. [2020] G. Tennenholtz, U. Shalit, and S. Mannor.  Off-policy evaluation in partially observable environments.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 34, pages 10276–10283, 2020.
  * Wang and Blei [2021] Y. Wang and D. Blei.  A proxy variable view of shared confounding.  In M. Meila and T. Zhang, editors, _Proceedings of the 38th International Conference on Machine Learning_ , volume 139 of _Proceedings of Machine Learning Research_ , pages 10697–10707. PMLR, 18–24 Jul 2021.  URL https://proceedings.mlr.press/v139/wang21c.html.
  * Wright [1928] P. Wright.  _The Tariff on Animal and Vegetable Oils_.  Investigations in International Commercial Policies. Macmillan Company, 1928.
  * Xu et al. [2021] L. Xu, Y. Chen, S. Srinivasan, N. de Freitas, A. Doucet, and A. Gretton.  Learning deep features in instrumental variable regression.  In _International Conference on Learning Representations_ , 2021.
  * Yao et al. [2018] L. Yao, S. Li, Y. Li, M. Huai, J. Gao, and A. Zhang.  Representation learning for treatment effect estimation from observational data.  In _Advances in Neural Information Processing Systems_ , volume 31, 2018.

##  Appendix A Hyper-Parameter Tuning

If observations from the joint distribution of $(A,Y,Z,W)$ are available in
both stages, we can tune the regularization parameters
$\lambda_{1},\lambda_{2}$ using the approach proposed in Singh et al. [30], Xu
et al. [35]. Let the complete data of stage 1 and stage 2 be denoted as
$(a_{i},y_{i},z_{i},w_{i})$ and
$(\tilde{a}_{i},\tilde{y}_{i},\tilde{z}_{i},\tilde{w}_{i})$. Then, we can use
the data not used in each stage to evaluate the out-of-sample performance of
the other stage. Specifically, let Algorithm 1 converges at $t=T$, and the
regularizing parameters are given by

$\displaystyle\lambda_{1}^{*}=\operatorname*{arg\,min}\mathcal{L}_{1\text{-oos}},\quad\mathcal{L}_{1\text{-oos}}=\frac{1}{n}\sum_{j=1}^{n}\left\|\bm{\psi}_{\theta_{W}^{(T)}}(\tilde{w}_{j})-\hat{\bm{V}}^{(T)}\left(\bm{\phi}_{\theta^{(T)}_{A(1)}}(\tilde{a}_{i})\otimes\bm{\phi}_{\theta^{(T)}_{Z}}(\tilde{z}_{i})\right)\right\|^{2},$
$\displaystyle\lambda_{2}^{*}=\operatorname*{arg\,min}\mathcal{L}_{2\text{-oos}},$
$\displaystyle\quad\mathcal{L}_{2\text{-oos}}=\frac{1}{n}\sum_{i=1}^{n}\left(y_{i}-(\bm{u}^{(T)})^{\top}\left(\bm{\psi}_{\theta^{(T)}_{A(2)}}(a_{i})\otimes\left(\hat{\bm{V}}^{(T)}\left(\bm{\phi}_{\theta^{(T)}_{A(1)}}(a_{i})\otimes\bm{\phi}_{\theta^{(T)}_{Z}}(z_{i})\right)\right)\right)\right)^{2}$

where
$\theta_{Z}^{(T)},\theta_{A(1)}^{(T)},\theta_{W}^{(T)},\theta_{A(2)}^{(T)},\hat{\bm{V}}^{(T)},\bm{u}^{(T)}$
are the learned parameters by Algorithm 1.

##  Appendix B Identifiability

In this appendix, we prove propositions given in the main text. In the
following, we assume that the spaces $\mathcal{U}$, $\mathcal{A}$,
$\mathcal{Z},\mathcal{W}$ are separable and completely metrizable topological
spaces and equipped with Borel $\sigma$-algebras. In this section, we use the
notation $P_{A|Z=z}$ to express the distribution of a random variable $A$
given another variable $Z=z$.

###  B.1 Existence of bridge function

First, we discuss conditions to guarantee the existence of the bridge function
$h$. Let us consider the following operators:

$\displaystyle E_{a}:L^{2}(P_{W|A=a})\to L^{2}(P_{Z|A=a}),\ E_{a}f$
$\displaystyle\coloneqq\mathbb{E}\left[f(W)|A=a,Z=\cdot\right],$
$\displaystyle F_{a}:L^{2}(P_{Z|A=a})\to L^{2}(P_{W|A=a}),\ F_{a}g$
$\displaystyle\coloneqq\mathbb{E}\left[g(Z)|A=a,W=\cdot\right],$

where the conditional expectations are identified as equivalent classes with
natural inclusion into their corresponding $L^{2}$ spaces. Our goal is to show
that $\mathbb{E}\left[Y|A=a,Z=\cdot\right]$ is in the range of $E_{a},$ i.e.,
we seek a solution of the inverse problem defined by

$E_{a}h=\mathbb{E}\left[Y|A=a,Z=\cdot\right].$ (8)

This suffices to prove the existence of the function $h,$ for if there exists
a function $h^{*}_{a}$ for each $a\in{\cal A}$ such that

$\mathbb{E}\left[h^{*}_{a}|A=a,Z=\cdot\right]=\mathbb{E}\left[Y|A=a,Z=\cdot\right],$

we can define $h^{*}(a,w)\coloneqq h^{*}_{a}(w).$ The existence of the bridge
function corresponds to the results of Deaner [3, Lemma 1.1.b] and Miao et al.
[20, Propisition 1]. For simplicity, in constrast to Deaner [3], we assume
that there is no shared quantity between the proxy variables $W$ and $Z$.
Following the proofs of the previous works [3, Appendix C, p. 57] and [20,
Appendix 6], we aim to solve the integral equation in (8).

We will make use of the following theorem on the existence of a solution of a
linear integral equation [12, Theorem 15.18].

######  Proposition 5 ([12, Theorem 15.18]).

Let $\mathcal{X}$ and $\mathcal{Y}$ be Hilbert spaces. Let
$E:\mathcal{X}\to\mathcal{Y}$ be a compact linear operator with singular
system $\\{(\mu_{n},\varphi_{n},g_{n})\\}_{n=1}^{\infty}.$ The equation of the
first kind

$E\varphi=f$

is solvable if and only if $f\in N(E^{*})^{\perp}$ and

$\sum_{n=1}^{\infty}\frac{1}{\mu_{n}^{2}}\lvert\langle
f,g_{n}\rangle\rvert^{2}<\infty.$

Here, $N(E^{*})$ denotes the null space of the operator $E^{*}$. Then a
solution is given by

$\phi=\sum_{n=1}^{\infty}\frac{1}{\mu_{n}}\langle f,g_{n}\rangle\varphi_{n}.$

To apply Proposition 5, we make the following additional assumptions.

######  Assumption 4.

For each $a\in{\cal A},$ the operator $E_{a}$ is compact with singular system
$\\{(\mu_{a,n},\varphi_{a,n},g_{a,n})\\}_{n=1}^{\infty}.$

######  Assumption 5.

For each $a\in{\cal A},$ the conditional expectation
$f_{Y|a}\coloneqq\mathbb{E}\left[Y|A=a,Z=\cdot\right]$ satisfies

$\sum_{n=1}^{\infty}\frac{1}{\mu_{a,n}^{2}}\lvert\langle
f_{Y|a},g_{a,n}\rangle_{L^{2}(P_{Z|A=a})}\rvert^{2}<\infty,$

for a singular system $\\{(\mu_{a,n},\phi_{a,n},g_{a,n})\\}_{n=1}^{\infty}$
given in Assumption 4.

######  Remark 1.

Assumption 4 is a minimal requisite to apply Proposition 5. The existing works
[20, 3] assume stronger conditions; the operator $E_{a}$ is assumed to be
Hilbert-Schmidt, which implies the compactness. Deaner [3, Assumption A.1]
assumes that the joint distribution of $W$ and $Z$ conditioned on $a$ is
absolutely continuous with respect to the product measure of $P_{W|A=a}$ and
$P_{Z|A=a},$ and its density is square integrable. The density function serves
as the integral kernel of the operator $E_{a}$ whose $L^{2}(P_{W|A=a}\otimes
P_{Z|A=a})$-norm corresponds to the Hilbert-Schmidt norm of the operator,
where we use $\otimes$ to denote the product measure. On the other hand, in
the setting of Miao et al. [20], all probability distributions have densities
with respect to the Lebesgue measure. The operator $E_{a}$ (denoted by $K_{x}$
in [20]) is defined by the relevant densities accordingly (see the paragraph
after Lemma 2 of [20]). The compactness is, as in [3], established by the
square-integrability of the integral kernel, which implies that the operator
$E_{a}$ is Hilbert-Schmidt (see Condition A1 in [20]).

It is easy to see that Assumptions 4 and 5 are required for using Proposition
5. The remaining condition to show is that
$\mathbb{E}\left[Y|A=a,Z=\cdot\right]$ is in $N(E_{a}^{*})^{\perp}.$ We show
that the structural assumption (Assumption 1) and completeness assumption
(Assumption 2) imply the required condition. The result below closely follows
[3] in that we share the same completeness condition.

######  Lemma 1.

Under Assumptions 1 and 2, the conditional expectation
$\mathbb{E}\left[Y|A=a,Z=\cdot\right]$ is in the orthogonal complement of the
null space $N(E_{a}^{*}).$

###### Proof.

We first show that the adjoint of $E_{a}$ is given by $F_{a}.$ For the
operator $E_{a},$ any $f\in L_{2}(P_{W|A=a})$ and $g\in L_{2}(P_{Z|A=a}),$ we
have

$\displaystyle\left\langle E_{a}f,g\right\rangle_{L_{2}(P_{Z|A=a})}$
$\displaystyle=\mathbb{E}_{Z|A=a}\left[\mathbb{E}\left[f(W)|A=a,Z\right]g(Z)\right]$
$\displaystyle=\mathbb{E}_{Z|A=a}\left[\mathbb{E}_{U|A=a,Z}\left[\mathbb{E}\left[f(W)|A=a,Z,U\right]\right]g(Z)\right]$
$\displaystyle\overset{\mathrm{(a)}}{=}\mathbb{E}_{Z|A=a}\left[\mathbb{E}_{U|A=a,Z}\left[\mathbb{E}\left[f(W)|A=a,U\right]\right]g(Z)\right]$
$\displaystyle=\mathbb{E}_{U,Z|A=a}\left[\mathbb{E}\left[f(W)|A=a,U\right]g(Z)\right]$
$\displaystyle=\mathbb{E}_{U|A=a}\left[\mathbb{E}\left[f(W)|A=a,U\right]\mathbb{E}\left[g(Z)|A=a,U\right]\right]$

where (a) follows from $W\mathop{\perp\\!\\!\\!\perp}(A,Z)|U$, which is from
Assumption 1. Similarly,

$\displaystyle\left\langle
f,F_{a}g\right\rangle_{L^{2}(P_{W|A=a})}=\mathbb{E}_{W|A=a}\left[f(W)\mathbb{E}\left[g(Z)|A=a,W\right]\right]$
$\displaystyle=\mathbb{E}_{W|A=a}\left[f(W)\mathbb{E}_{U|A=a,W}\left[\mathbb{E}\left[g(Z)|A=a,W,U\right]\right]\right]$
$\displaystyle\overset{\mathrm{(b)}}{=}\mathbb{E}_{W|A=a}\left[f(W)\mathbb{E}_{U|A=a,W}\left[\mathbb{E}\left[g(Z)|A=a,U\right]\right]\right]$
$\displaystyle=\mathbb{E}_{W,U|A=a}\left[f(W)\mathbb{E}\left[g(Z)|A=a,U\right]\right]$
$\displaystyle=\mathbb{E}_{U|A=a}\left[\mathbb{E}\left[f(W)|A=a,U\right]\mathbb{E}\left[g(Z)|A=a,U\right]\right]$
$\displaystyle=\left\langle E_{a}f,g\right\rangle_{L_{2}(P_{Z|A=a})}.$

Again, (b) is given by $W\mathop{\perp\\!\\!\\!\perp}A,Z|U$ from Assumption 1.
For any $f^{*}\in N(E^{*}_{a})=N(F_{a}),$ by iterated expectations, we have

$\displaystyle 0$ $\displaystyle=\mathbb{E}\left[f^{*}(Z)|A=a,W=\cdot\right]$
$\displaystyle=\mathbb{E}_{U}\left[\mathbb{E}\left[f^{*}(Z)|A,U,W\right]|A=a,W=\cdot\right]$
$\displaystyle=\mathbb{E}_{U}\left[\mathbb{E}\left[f^{*}(Z)|A,U\right]|A=a,W=\cdot\right].$
(9)

From Assumption 2,

$\displaystyle\mathbb{E}\left[l(U)\mid
A=a,W=w\right]=0\quad\forall(a,w)\in\mathcal{A}\times\mathcal{W}\quad\Leftrightarrow\quad
l(u)=0\quad P_{U}\text{-}\mathrm{a.e.}$

for all functions $l\in L^{2}(P_{U|A=a})$. Note that
$\mathbb{E}\left[f^{*}(Z)|A,U\right]\in L^{2}(P_{U|A=a})$ since

$\displaystyle\mathbb{E}_{U|A}\left[\mathbb{E}[f^{*}(Z)|A,U]^{2}\right]$
$\displaystyle\leq\mathbb{E}_{U|A}\left[\mathbb{E}\left[f^{*}(Z)^{2}|A,U\right]\right]\quad(\because\text{Jensen's
Inequality})$ $\displaystyle=\mathbb{E}\left[f^{*}(Z)^{2}|A\right]<\infty$

Hence, (9) and Assumption 2 implies

$\mathbb{E}\left[f^{*}(Z)|A=a,U=\cdot\right]=0\quad
P_{U}\text{-}\mathrm{a.s.}$

Then, the inner product between $f^{*}$ and
$\mathbb{E}\left[Y|A=a,Z=\cdot\right]$ is given as follows:

$\displaystyle\left\langle
f^{*},\mathbb{E}\left[Y|A=a,Z=\cdot\right]\right\rangle_{L^{2}(P_{Z|A=a})}$
$\displaystyle=\mathbb{E}_{Z|A=a}\left[f^{*}(Z)\mathbb{E}\left[Y|A=a,Z\right]\right]$
$\displaystyle=\mathbb{E}_{Z|A=a}\left[f^{*}(Z)\mathbb{E}_{U|A=a,Z}\left[\mathbb{E}\left[Y|A=a,Z,U\right]\right]\right]$
$\displaystyle\overset{(c)}{=}\mathbb{E}_{Z|A=a}\left[f^{*}(Z)\mathbb{E}_{U|A=a,Z}\left[\mathbb{E}\left[Y|A=a,U\right]\right]\right]$
$\displaystyle=\mathbb{E}_{U,Z|A=a}\left[f^{*}(Z){\mathbb{E}\left[Y|A=a,U\right]}\right]$
$\displaystyle=\mathbb{E}_{U|A=a}\left[\mathbb{E}_{Z|U,A=a}\left[f^{*}(Z)\right]{\mathbb{E}\left[Y|A=a,U\right]}\right]$
$\displaystyle=0.$

Again (c) holds from $Y\mathop{\perp\\!\\!\\!\perp}Z|A,U$ in Assumption 1.
Hence, we have

$\mathbb{E}\left[Y|A=a,Z=\cdot\right]\in N(E_{a}^{*})^{\perp}.$

∎

######  Remark 2.

The difference between the first condition in Assumption 2 and Condition 3 in
[20] is in the approach to establishing that the conditional expectation
$\mathbb{E}[Y|A=a,Z=\cdot]$ belongs to $N(F_{a})^{\perp}.$ More specifically,
Condition 3 in [20] is equivalent to having $N(F_{a})=\\{0\\}$ and so any
nontrivial $L^{2}(P_{Z|A=a})$-function is in the orthogonal complement.
Assumption 2 does not require $N(F_{a})=\\{0\\}$ but only implies

$N(F_{a})\subset\\{f\in L^{2}(P_{Z|A=a}):\mathbb{E}[f(Z)|A=a,U=\cdot]=0\\}.$

For our identifiability proof, we use the weaker condition, Assumption 2.

Now, we are able to apply Proposition 5, leading to the following lemma.

######  Lemma 2.

Under Assumptions 1, 2, 4 and 5, for each $a\in{\cal A},$ there exists a
function $h^{*}_{a}\in L_{2}(P_{W|A=a})$ such that

$\mathbb{E}\left[Y|A=a,Z=\cdot\right]=\mathbb{E}\left[h^{*}_{a}(W)|A=a,Z=\cdot\right].$

###### Proof.

By Lemma 1, the regression function $\mathbb{E}\left[Y|A=a,Z=\cdot\right]$ is
in $N(E_{a}^{*})^{\perp}.$ Therefore, by Proposition 5, under the given
assumptions, there exists a solution to (8). Letting the solution be
$h^{*}_{a}$ completes the proof. ∎

###  B.2 Identifiability

Here, we show that the bridge function $h$ can be used to compute the various
causal quantities. In addition to the assumptions required for the existence
of bridge function listed in the previous section, we assume the conditional
expectation of outcome given treatment is square integrable.

######  Assumption 6.

For each $a\in{\cal A},$ under the observational distribution, we have

$\mathbb{E}\left[Y^{2}|A=a\right]<\infty.$

Given the assumption, we can prove the following theorem.

######  Theorem 1 (Identifiability).

Assume that there exists a function $h:{\cal A}\times{\cal W}\to\mathbb{R}$
such that for each $a\in{\cal A},$ the function $h(a,\cdot)$ is in
$L^{2}(P_{W|A=a})$ and satisfies

$\mathbb{E}\left[Y|A=a,Z=\cdot\right]=\mathbb{E}\left[h(a,W)|A=a,Z=\cdot\right]\quad
P_{Z|A=a}\text{-a.s.}$

Under Assumptions 1, 2 and 6, we have

$\mathbb{E}_{W|U}\left[h(a,W)\right]=\mathbb{E}\left[Y|A=a,U\right].$ (10)

###### Proof.

Note that we have

$\displaystyle\mathbb{E}\left[h(a,W)|A=a,Z\right]$
$\displaystyle=\mathbb{E}_{U|A=a,Z}\left[\mathbb{E}\left[h(a,W)|A=a,Z,U\right]\right]$
$\displaystyle=\mathbb{E}_{U|A=a,Z}\left[\mathbb{E}\left[h(a,W)|U\right]\right],$
$\displaystyle\mathbb{E}\left[Y|A=a,Z\right]$
$\displaystyle=\mathbb{E}_{U|A=a,Z}\left[\mathbb{E}\left[Y|A=a,Z,U\right]\right]$
$\displaystyle=\mathbb{E}_{U|A=a,Z}\left[\mathbb{E}\left[Y|A=a,U\right]\right],$

where the second line of each equation follows from Assumption 1. Moreover, by
assumption, we have

$\displaystyle\mathbb{E}_{U|A=a}\left[\mathbb{E}\left[h(a,W)|U\right]^{2}\right]$
$\displaystyle\leq\mathbb{E}_{W|A=a}\left[h(a,W)^{2}\right]<\infty,\
\text{and}$
$\displaystyle\mathbb{E}_{U|A=a}\left[\mathbb{E}\left[Y|A=a,U\right]^{2}\right]$
$\displaystyle\leq\mathbb{E}\left[Y^{2}|A=a\right]<\infty.$

Note from Assumption 2, we have

$\mathbb{E}\left[l(U)\mid
A=a,Z=z\right]=0\quad\forall(a,z)\in\mathcal{A}\times\mathcal{Z}\quad\Leftrightarrow\quad
l(u)=0\quad P_{U}\text{-a.e.}$

Therefore, by setting
$l(u)=\mathbb{E}\left[Y|A=a,U=u\right]-\mathbb{E}\left[h(a,W)|U=u\right]$, for
all $a\in\mathcal{A}$, we have

$\displaystyle\mathbb{E}\left[h(a,W)|A=a,Z=\cdot\right]=\mathbb{E}\left[Y|A=a,Z=\cdot\right]\quad
P_{Z|A=a}\text{-a.s.}$ $\displaystyle\Leftrightarrow$
$\displaystyle\mathbb{E}_{U|A=a,Z=\cdot}\left[\mathbb{E}\left[h(a,W)|U\right]-\mathbb{E}\left[Y|A=a,U\right]\right]=0\quad
P_{Z|A=a}\text{-a.s.}$ $\displaystyle\Leftrightarrow$
$\displaystyle\mathbb{E}_{W|U=\cdot}\left[h(a,W)\right]=\mathbb{E}\left[Y|A=a,U=\cdot\right]\quad
P_{U}\text{-a.s.}$

∎

Using Theorem 1, we can show following two corollaries, which are used in the
main body.

######  Corollary 1.

Let the assumptions in Theorem 1 hold. Given a bridge function $h^{*}$, we can
estimate structural function $f_{\mathrm{struct}}$ as

$\displaystyle f_{\mathrm{struct}}(a)$
$\displaystyle\coloneqq\mathbb{E}_{U}\left[\mathbb{E}\left[Y|A=a,U\right]\right]$
$\displaystyle=\mathbb{E}_{W}\left[h^{*}(a,W)\right]$

###### Proof.

From Theorem 1, we have

$\displaystyle\mathbb{E}_{U}\left[\mathbb{E}\left[Y|A=a,U\right]\right]=\mathbb{E}_{U}\left[\mathbb{E}_{W|U}\left[h^{*}(a,W)\right]\right]=\mathbb{E}_{W}\left[h^{*}(a,W)\right]$

∎

######  Remark 3.

The above corollary corresponds to the identifiablity results in obtained in
the previous works [20, 3]. We follow the proof of Theorem 1.1.a in [3,
Appendix B] (See also, Theorem 1 and Appendix 3 of [20]).

######  Corollary 2.

Assume we are given a bridge function $h^{*}(a,w)$ that is jointly measurable.
Suppose $C\mathop{\perp\\!\\!\\!\perp}W|U$. With the assumptions in Theorem 1,
we can write the value function of policy $\pi(C)$ as

$\displaystyle v(\pi)$
$\displaystyle\coloneqq\mathbb{E}_{C,U}\left[\mathbb{E}\left[Y|A=\pi(C),U\right]\right]$
$\displaystyle=\mathbb{E}_{C}\left[\mathbb{E}_{W|C}\left[h(\pi(C),W)\right]\right]$

###### Proof.

From Theorem 1, we have

$\displaystyle\mathbb{E}_{C,U}\left[\mathbb{E}\left[Y|A=\pi(C),U\right]\right]=\mathbb{E}_{C,U}\left[\mathbb{E}_{W|U}\left[h^{*}(\pi(C),W)\right]\right]=\mathbb{E}_{C}\left[\mathbb{E}_{W|C}\left[h(\pi(C),W)\right]\right],$

where the last equality holds by the conditional independence
$C\mathop{\perp\\!\\!\\!\perp}W|U$. ∎

Note that in the existence claim, the bridge function $h^{*}(a,w)$ is
constructed by aggregating over $\\{h^{*}_{a}\\}_{a\in\mathcal{A}}$. This
construction does not guarantee that the function is measurable with respect
to $a$; the lack of measurability renders its expectation undefined. Thus, we
have additionally assumed the joint measurability of the bridge function in
Corollary 2. Validating this assumption is crucial theoretically, and we leave
it for future work.

##  Appendix C Consistency of DFPV algorithm

In this appendix, we prove consistency of the DFPV approach. Following Xu et
al. [35], we establish consistency of the end-to-end procedure incorporating
Stages 1 and 2. We establish the result by first showing a Stage 1 consistency
result (Lemma 3), and then establishing the consistency of Stage 2 with the
empirical Stage 1 solution used as input (Lemma 4). The desired result then
follows in Theorem 2. Here, we assume the bridge function $h^{*}(a,w)$ is
jointly measurable so that the expectation of $h^{*}$ is defined.

Consistency results will be expressed in terms of the complexity of the
function classes used in Stages 1 and 2, as encoded in the Rademacher
complexity of the function space of these functions (see Proposition 6 below).
Consistency for particular function classes can then be shown by establishing
that the respective Rademacher complexities vanish. We leave for future work
the task of demonstrating this property for individual function classes of
interest.

###  C.1 Operator view of DFPV

The goal of DFPV is to learn a bridge function $h^{*}$, which satisfies

$\displaystyle\mathbb{E}_{Y|A=a,Z=z}\left[Y\right]=\mathbb{E}_{W|A=a,Z=z}\left[h^{*}(a,W)\right].$
(11)

We use the model

$h(a,w)=\bm{u}^{\top}(\bm{\psi}_{\theta_{A(2)}}(a)\otimes\bm{\psi}_{\theta_{W}}(w))$

and denote the hypothesis spaces for $\bm{\psi}_{\theta_{W}}$ and $h$ as
$\mathcal{H}_{\bm{\psi}_{\theta_{W}}}:\mathcal{W}\to\mathbb{R}^{d_{W}}$ and
$\mathcal{H}_{h}:\mathcal{A}\times\mathcal{W}\to\mathbb{R}$, respectively. To
learn the parameters, we minimize the following stage 2 loss:

$\displaystyle\hat{\bm{u}},\hat{\theta}_{W},\hat{\theta}_{A(2)}=\operatorname*{arg\,min}_{\bm{u},\theta_{W},\theta_{A(2)}}\hat{\mathcal{L}}_{2}(\bm{u},\theta_{X},\theta_{A(2)}),$
$\displaystyle\hat{\mathcal{L}}_{2}=\frac{1}{n}\sum_{i=1}^{n}(\tilde{y}_{i}-\bm{u}^{\top}(\bm{\psi}_{\theta_{A(2)}}(\tilde{a}_{i})\otimes\hat{\mathbb{E}}_{W|A,Z}\left[\bm{\psi}_{\theta_{W}}(W)\right](\tilde{a}_{i},\tilde{z}_{i})))^{2}.$

We denote the resulting estimated bridge function as

$\hat{h}(a,w)=(\hat{\bm{u}})^{\top}(\bm{\psi}_{\hat{\theta}_{A(2)}}(a)\otimes\bm{\psi}_{\hat{\theta}_{W}}(w)).$

For simplicity, we set all regularization terms to zero. Here,
$\hat{\mathbb{E}}_{W|A,Z}\left[\bm{\psi}_{\theta_{W}}(W)\right]$ is the
empirical conditional expectation operator, which maps an element of
$\mathcal{H}_{\bm{\psi}_{\theta_{W}}}$ to some function $\bm{g}(a,z)$ which is
defined as

$\displaystyle\hat{\mathbb{E}}_{W|A,Z}\left[\bm{\psi}_{\theta_{W}}(W)\right]=\operatorname*{arg\,min}_{\bm{g}\in\mathcal{G}}\hat{\mathcal{L}}_{1}(\bm{g};\bm{\psi}_{{\theta}_{W}}),$
$\displaystyle\hat{\mathcal{L}}_{1}=\frac{1}{n}\sum_{i=1}^{m}\|\bm{\psi}_{{\theta}_{W}}(w_{i})-\bm{g}(a_{i},z_{i})\|^{2},$

where $\|\cdot\|$ is the $\ell_{2}$-norm, and $\mathcal{G}$ is an arbitrary
function space. In DFPV, we specify $\mathcal{G}$ as the set consisting of
functions $\bm{g}$ of the form

$\bm{g}={\bm{V}}(\bm{\phi}_{{\theta}_{A(1)}}(a_{i})\otimes\bm{\phi}_{{\theta}_{Z}}(z_{i})).$

Note that this formulation is equivalent to the one introduced in Section 3.
With a slight abuse of notation, for
$h(a,w)=\bm{u}^{\top}\bm{\psi}_{\theta_{A(2)}}(a)\otimes\bm{\psi}_{\theta_{W}}(w)\in\mathcal{H}_{h}$,
we define $\hat{\mathbb{E}}_{W|A,Z}\left[h\right]$ to be

$\displaystyle\hat{\mathbb{E}}_{W|A,Z}\left[h(A,W)\right](a,z)=\bm{u}^{\top}\left(\bm{\psi}_{\theta_{A(2)}}(a)\otimes\hat{\mathbb{E}}_{W|A,Z}\left[\bm{\psi}_{\theta_{W}}\right](a,z)\right)$

since this is the empirical estimate of
$\mathbb{E}_{W|A,Z}\left[h(A,W)\right]$.

###  C.2 Generalization errors for regression

Here, we bound the generalization errors of both stages using Rademacher
complexity bounds [24].

######  Proposition 6.

[Theorem 3.3 24, with slight modification] Let $\mathcal{S}$ be a measurable
space and $\mathcal{H}$ be a family of functions mapping from $\mathcal{S}$ to
$[0,M]$. Given fixed dataset $S=(s_{1},s_{2},\dots,s_{n})\in\mathcal{S}^{n}$,
the empirical Rademacher complexity is given by

$\displaystyle\hat{\mathfrak{R}}_{S}(\mathcal{H})=\mathbb{E}_{\bm{\sigma}}\left[\frac{1}{n}\sup_{h\in\mathcal{H}}\sum_{i=1}^{n}\sigma_{i}h(s_{i})\right],$

where $\bm{\sigma}=(\sigma_{1},\dots,\sigma_{n})$, with $\sigma_{i}$
independent random variables taking values in $\\{-1,+1\\}$ with equal
probability. Then, for any $\delta>0$, with probability at least $1-\delta$
over the draw of an i.i.d sample $S$ of size $n$, each of following holds for
all $h\in\mathcal{H}$:

$\displaystyle\mathbb{E}\left[h(s)\right]\leq\frac{1}{n}\sum_{i=1}^{n}h(s_{i})+2\hat{\mathfrak{R}}_{S}+3M\sqrt{\frac{\log
2/\delta}{2n}},$
$\displaystyle\frac{1}{n}\sum_{i=1}^{n}h(s_{i})\leq\mathbb{E}\left[h(s)\right]+2\hat{\mathfrak{R}}_{S}+3M\sqrt{\frac{\log
2/\delta}{2n}}.$

We list the assumptions below.

######  Assumption 7.

The following hold:

  1. 1.

Bounded outcome variable $|Y|\leq M$.

  2. 2.

Bounded stage 1 hypothesis space: $\forall
a\in\mathcal{A},z\in\mathcal{Z},\|\bm{g}(a,z)\|\leq 1$.

  3. 3.

Bounded stage 2 feature map $\forall
w\in\mathcal{W},\|\bm{\psi}_{\theta_{W}}(w)\|\leq 1$.

  4. 4.

Bounded stage 2 weight: $\forall\|\bm{v}\|\leq 1,\forall
a\in\mathcal{A},|\bm{u}^{\top}(\bm{\psi}_{\theta_{A(2)}}(a)\otimes\bm{v})|\leq
M$.

Note that unlike Xu et al. [35], we consider the case where the bridge
function $h^{*}$ and the conditional expectation
$\mathbb{E}_{W|A,Z}\left[\bm{\psi}_{\theta_{W}}\right]$ might not be included
in the hypothesis spaces $\mathcal{H}_{h}$, $\mathcal{G}$, respectively.

As in Xu et al. [35], we leave aside questions of optimization. Thus, we
assume that the optimization procedure over
$(\theta_{A(1)},\theta_{Z},\bm{V})$ is sufficient to recover
$\hat{\mathbb{E}}_{W|A,Z}\left[\bm{\psi}_{\theta_{W}}(W)\right],$ and that the
optimization procedure over $(\theta_{A(2)},\theta_{W},\bm{u})$ is sufficient
to recover $\hat{h}$ (which requires, in turn, the correct
$\hat{\mathbb{E}}_{W|A,Z},$ for this $\bm{\psi}_{\theta_{W}}$). We emphasize
that Algorithm 1 does not guarantee these properties. Based on these
assumptions, we derive the generalization error in terms of $L_{2}$-norm with
respect to $\mathbb{P}\left(A,W\right)$.

Following the discussion in Xu et al. [35, Lemma 1], we can show the following
lemma that provides a generalization error bound for Stage 1 regression.

######  Lemma 3.

Under Assumption 7, and given stage 1 data
$S_{1}=\\{(a_{i},z_{i},w_{i})\\}_{i=1}^{m}$, for any $\delta>0$, with at least
probability $1-2\delta$, we have

$\displaystyle\left\|\hat{\mathbb{E}}_{W|A,Z}\left[h(A,W)\right]-\mathbb{E}_{W|A,Z}\left[h(A,W)\right]\right\|_{P(A,Z)}\leq
M\sqrt{\kappa_{1}+4\hat{\mathfrak{R}}_{S_{1}}(\mathcal{H}_{1})+24\sqrt{\frac{\log
2/\delta}{2m}}}$

for any $\bm{\psi}_{\theta_{W}}\in\mathcal{H}_{\bm{\psi}_{\theta_{W}}}$, where
hypothesis space $\mathcal{H}_{1}$ is defined as

$\displaystyle\mathcal{H}_{1}=\\{(w,a,z)\in\mathcal{W}\times\mathcal{A}\times\mathcal{Z}\mapsto\|\bm{\psi}_{\theta_{W}}(w)-\bm{g}(a,z)\|^{2}\in\mathbb{R}~{}|~{}\bm{g}\in\mathcal{G},\bm{\psi}_{\theta_{W}}\in\mathcal{H}_{\bm{\psi}_{\theta_{W}}}\\},$
(12)

and $\kappa_{1}$ is the misspecificaiton error in the stage 1 regression
defined as

$\displaystyle\kappa_{1}=\max_{\bm{\psi}_{\theta_{W}}\in\mathcal{H}_{\bm{\psi}_{\theta_{W}}}}\min_{\bm{g}\in\mathcal{G}}\mathbb{E}_{W,A,Z}\left[\|\bm{\psi}_{\theta_{W}}(W)-\bm{g}(A,Z)\|^{2}\right].$
(13)

###### Proof.

From Assumption 7, we have

$\displaystyle\left\|\hat{\mathbb{E}}_{W|A,Z}\left[h(A,W)\right]-\mathbb{E}_{W|A,Z}\left[h(A,W)\right]\right\|_{P(A,Z)}$
$\displaystyle=\left\|\bm{u}^{\top}\left(\bm{\psi}_{\theta_{A(2)}}(a)\otimes\left(\mathbb{E}_{W|A,Z}\left[\bm{\psi}_{\theta_{W}}\right](a,z)-\hat{\mathbb{E}}_{W|A,Z}\left[\bm{\psi}_{\theta_{W}}\right](a,z)\right)\right)\right\|_{P(A,Z)}$
$\displaystyle=\sqrt{\mathbb{E}_{A,Z}\left[\left|\bm{u}^{\top}\left(\bm{\psi}_{\theta_{A(2)}}(a)\otimes\left(\mathbb{E}_{W|A,Z}\left[\bm{\psi}_{\theta_{W}}\right](a,z)-\hat{\mathbb{E}}_{W|A,Z}\left[\bm{\psi}_{\theta_{W}}\right](a,z)\right)\right)\right|^{2}\right]}$
$\displaystyle\leq\sqrt{\mathbb{E}_{A,Z}\left[M^{2}\left\|\mathbb{E}_{W|A,Z}\left[\bm{\psi}_{\theta_{W}}\right](a,z)-\hat{\mathbb{E}}_{W|A,Z}\left[\bm{\psi}_{\theta_{W}}\right](a,z)\right\|^{2}\right]}$
$\displaystyle=M\left\|\mathbb{E}_{W|A,Z}\left[\bm{\psi}_{\theta_{W}}\right](a,z)-\hat{\mathbb{E}}_{W|A,Z}\left[\bm{\psi}_{\theta_{W}}\right](a,z)\right\|_{P(A,Z)}$

From Proposition 6, we have

$\displaystyle\mathbb{E}_{W,A,Z}\left[\|\bm{\psi}_{\theta_{W}}(W)-\hat{\mathbb{E}}_{W|A,Z}\left[\bm{\psi}_{\theta_{W}}(W)\right](A,Z)\|^{2}\right]$
$\displaystyle\quad\leq\hat{\mathcal{L}}_{1}(\hat{\mathbb{E}}_{W|A,Z}\left[\bm{\psi}_{\theta_{W}}(W)\right])+2\hat{\mathfrak{R}}_{S_{1}}(\mathcal{H}_{1})+12\sqrt{\frac{\log
2/\delta}{2m}}$

with probability $1-\delta$, since all functions $f\in\mathcal{H}_{1}$ satisfy
$\|f\|\leq 4$ from $\|\bm{\psi}_{\theta_{W}}\|\leq 1,\|\bm{g}\|\leq 1.$ Now,
let $\bm{g}^{*}_{\theta_{W}}$ be

$\displaystyle\bm{g}^{*}_{\theta_{W}}=\operatorname*{arg\,min}_{\bm{g}\in\mathcal{G}}\mathbb{E}_{W,A,Z}\left[\|\bm{\psi}_{\theta_{W}}(W)-\bm{g}(A,Z)\|^{2}\right].$

Then, again from Proposition 6, we have

$\displaystyle\hat{\mathcal{L}}_{1}(\bm{g}^{*}_{\theta_{W}})\leq\mathbb{E}_{W,A,Z}\left[\|\bm{\psi}_{\theta_{W}}(W)-\bm{g}^{*}_{\theta_{W}}(A,Z)\|^{2}\right]+2\hat{\mathfrak{R}}_{S_{1}}(\mathcal{H}_{1})+12\sqrt{\frac{\log
2/\delta}{2m}}$

From the optimality of
$\hat{\mathbb{E}}_{W|A,Z}\left[\bm{\psi}_{\theta_{W}}(W)\right]$, we have
$\hat{\mathcal{L}}_{1}(\bm{g}^{*}_{\theta_{W}})\geq\hat{\mathcal{L}}_{1}(\hat{\mathbb{E}}_{W|A,Z}\left[\bm{\psi}_{\theta_{W}}(W)\right])$,
thus

$\displaystyle\mathbb{E}_{W,A,Z}\left[\|\bm{\psi}_{\theta_{W}}(W)-\hat{\mathbb{E}}_{W|A,Z}\left[\bm{\psi}_{\theta_{W}}(W)\right](A,Z)\|^{2}\right]$
$\displaystyle\quad\leq\mathbb{E}_{W,A,Z}\left[\|\bm{\psi}_{\theta_{W}}(W)-\bm{g}^{*}_{\theta_{W}}(A,Z)\|^{2}\right]+4\hat{\mathfrak{R}}_{S_{1}}(\mathcal{H}_{1})+24\sqrt{\frac{\log
2/\delta}{2m}}.$ (14)

holds with probability $1-2\delta$. Since we have

$\displaystyle\mathbb{E}_{W,A,Z}\left[\|\bm{\psi}_{\theta_{W}}(W)-g(A,Z)\|^{2}\right]$
$\displaystyle=\mathbb{E}_{W,A,Z}\left[\|\bm{\psi}_{\theta_{W}}(W)-\mathbb{E}_{W|A,Z}\left[\bm{\psi}_{\theta_{W}}(W)\right]\|^{2}\right]+\mathbb{E}_{A,Z}\left[\|g(A,Z)-\mathbb{E}_{W|A,Z}\left[\bm{\psi}_{\theta_{W}}(W)\right]\|^{2}\right]$

for all $\bm{g}\in\mathcal{G}$, by subtracting
$\mathbb{E}_{W,A,Z}\left[\|\bm{\psi}_{\theta_{W}}(W)-\mathbb{E}_{W|A,Z}\left[\bm{\psi}_{\theta_{W}}(W)\right]\|^{2}\right]$
from both sides of (14), we have

$\displaystyle\mathbb{E}_{A,Z}\left[\|\mathbb{E}_{W|A,Z}\left[\bm{\psi}_{\theta_{W}}(W)\right]-\hat{\mathbb{E}}_{W|A,Z}\left[\bm{\psi}_{\theta_{W}}(W)\right](A,Z)\|^{2}\right]$
$\displaystyle\quad\leq\mathbb{E}_{A,Z}\left[\|\mathbb{E}_{W|A,Z}\left[\bm{\psi}_{\theta_{W}}(W)\right]-\bm{g}^{*}_{\theta_{W}}(A,Z)\|^{2}\right]+4\hat{\mathfrak{R}}_{S_{1}}(\mathcal{H}_{1})+24\sqrt{\frac{\log
2/\delta}{2m}}$
$\displaystyle\quad\leq\kappa_{1}+4\hat{\mathfrak{R}}_{S_{1}}(\mathcal{H}_{1})+24\sqrt{\frac{\log
2/\delta}{2m}}$

By taking the square root of both sides, we have

$\displaystyle\left\|\hat{\mathbb{E}}_{W|A,Z}\left[\bm{\psi}_{\theta_{W}}(W)\right]-\mathbb{E}_{W|A,Z}\left[\bm{\psi}_{\theta_{W}}(W)\right]\right\|_{P(A,Z)}\leq
M\sqrt{\kappa_{1}+4\hat{\mathfrak{R}}_{S_{1}}(\mathcal{H}_{1})+24\sqrt{\frac{\log
2/\delta}{2m}}}$

∎

The generalization error for Stage 2 can be shown by similar reasoning as in
Xu et al. [35, Lemma 2].

######  Lemma 4.

Under Assumption 7, given stage 1 data
$S_{1}=\\{(a_{i},z_{i},w_{i})\\}_{i=1}^{m},$ stage 2 data
$S_{2}=\\{(\tilde{a}_{i},\tilde{z}_{i},\tilde{y}_{i})\\}_{i=1}^{n},$ and the
estimated structural function
$\hat{h}(a,z)=(\hat{\bm{u}})^{\top}(\bm{\psi}_{\hat{\theta}_{A(2)}}(a)\otimes\bm{\psi}_{\hat{\theta}_{W}}(w))$,
then for any $\delta>0$, with at least probability $1-4\delta$, we have

$\displaystyle\left\|\mathbb{E}_{Y|A,Z}\left[Y\right]-\hat{\mathbb{E}}_{W|A,Z}\left[\hat{h}(A,W)\right]\right\|_{P(A,Z)}$
$\displaystyle\quad\leq\kappa_{2}+M\sqrt{\kappa_{1}+4\hat{\mathfrak{R}}_{S_{1}}(\mathcal{H}_{2})+24\sqrt{\frac{\log
2/\delta}{2m}}}+\sqrt{4\hat{\mathfrak{R}}_{S_{2}}(\mathcal{H}_{2})+24M^{2}\sqrt{\frac{\log
2/\delta}{2n}}},$

where $\mathcal{H}_{1}$ is defined in (12), and $\mathcal{H}_{2}$ is defined
as

$\displaystyle\mathcal{H}_{2}=\\{(y,a,z)$
$\displaystyle\in\mathcal{Y}\times\mathcal{A}\times\mathcal{Z}$
$\displaystyle\mapsto(y-\bm{u}^{\top}(\bm{\psi}_{\theta_{A(2)}}(a)\otimes\bm{g}(a,z)))^{2}\in\mathbb{R}~{}|~{}\bm{g}\in\mathcal{G},\bm{u},\theta_{A(2)}\\},$
(15)

and $\kappa_{2}$ is the misspecification error in Stage 2 defined as

$\displaystyle\kappa_{2}=\min_{h\in\mathcal{H}_{h}}\left\|\mathbb{E}_{W|A,Z}\left[h(A,W)\right](A,Z)-\mathbb{E}_{Y|A,Z}\left[Y\right](A,Z)\right\|_{P(A,Z)}.$
(16)

###### Proof.

From Proposition 6, we have

$\displaystyle\mathbb{E}_{Y,A,Z}\left[\left|Y-\hat{\mathbb{E}}_{W|A,Z}\left[\hat{h}(A,W)\right](A,Z)\right|^{2}\right]\leq\hat{\mathcal{L}}_{2}(\hat{h})+2\hat{\mathfrak{R}}_{S_{2}}(\mathcal{H}_{2})+12M^{2}\sqrt{\frac{\log
2/\delta}{2n}}$

with probabiltiy $1-\delta$, since all functions $f\in\mathcal{H}_{2}$ satisfy
$\|f\|\leq 4M^{2}$ from $\|Y\|\leq
M,|\bm{u}^{\top}(\bm{\psi}_{\theta_{A(2)}}(a)\otimes\bm{g})|\leq M$. Let
$\tilde{h}$ be

$\displaystyle\tilde{h}=\operatorname*{arg\,min}_{h\in\mathcal{H}_{h}}\left\|\mathbb{E}_{W|A,Z}\left[h(A,W)\right](A,Z)-\mathbb{E}_{Y|A,Z}\left[Y\right](A,Z)\right\|_{P(A,Z)}.$

Again from Proposition 6, we have

$\displaystyle\hat{\mathcal{L}}_{2}(\tilde{h})\leq\mathbb{E}_{Y,A,Z}\left[\left|Y-\hat{\mathbb{E}}_{W|A,Z}\left[\tilde{h}(A,W)\right](A,Z)\right|^{2}\right]+2\hat{\mathfrak{R}}_{S_{2}}(\mathcal{H}_{2})+12M^{2}\sqrt{\frac{\log
2/\delta}{2n}}.$

with probabiltiy $1-\delta$. From the optimality of $\hat{h}$, we have
$\hat{\mathcal{L}}_{2}(\tilde{h})\geq\hat{\mathcal{L}}_{2}(\hat{h})$, hence

$\displaystyle\mathbb{E}_{Y,A,Z}\left[\left|Y-\hat{\mathbb{E}}_{W|A,Z}\left[\hat{h}(A,W)\right](A,Z)\right|^{2}\right]$
$\displaystyle\leq\mathbb{E}_{Y,A,Z}\left[\left|Y-\hat{\mathbb{E}}_{W|A,Z}\left[\tilde{h}(A,W)\right](A,Z)\right|^{2}\right]+4\hat{\mathfrak{R}}_{S_{2}}(\mathcal{H}_{2})+24M^{2}\sqrt{\frac{\log
2/\delta}{2n}}$

By subtracting
$\mathbb{E}_{Y,A,Z}\left[\|Y-\mathbb{E}_{Y|A,Z}\left[Y\right]\|^{2}\right]$
from both sides, we have

$\displaystyle\mathbb{E}_{A,Z}\left[\left|\mathbb{E}_{Y|A,Z}\left[Y\right]-\hat{\mathbb{E}}_{W|A,Z}\left[\hat{h}(A,W)\right](A,Z)\right|^{2}\right]$
$\displaystyle\leq\mathbb{E}_{A,Z}\left[\left|\mathbb{E}_{Y|A,Z}\left[Y\right]-\hat{\mathbb{E}}_{W|A,Z}\left[\tilde{h}(A,W)\right](A,Z)\right|^{2}\right]+4\hat{\mathfrak{R}}_{S_{2}}(\mathcal{H}_{2})+24M^{2}\sqrt{\frac{\log
2/\delta}{2n}}$

with probability $1-2\delta$. By taking the square root of both sides, with
probability $1-2\delta$, we have

$\displaystyle\|\mathbb{E}_{Y|A,Z}\left[Y\right]-\hat{\mathbb{E}}_{W|A,Z}\left[\tilde{h}(A,W)\right](A,Z)\|_{P(A,Z)}$
$\displaystyle\leq\sqrt{\mathbb{E}_{A,Z}\left[\left|\mathbb{E}_{Y|A,Z}\left[Y\right]-\hat{\mathbb{E}}_{W|A,Z}\left[\tilde{h}(A,W)\right](A,Z)\right|^{2}\right]+4\hat{\mathfrak{R}}_{S_{2}}(\mathcal{H}_{2})+24M^{2}\sqrt{\frac{\log
2/\delta}{2n}}}$
$\displaystyle\leq\|\mathbb{E}_{Y|A,Z}\left[Y\right]-\hat{\mathbb{E}}_{W|A,Z}\left[\tilde{h}(A,W)\right](A,Z)\|_{P(A,Z)}+\sqrt{4\hat{\mathfrak{R}}_{S_{2}}(\mathcal{H}_{2})+24M^{2}\sqrt{\frac{\log
2/\delta}{2n}}}$
$\displaystyle\overset{(a)}{\leq}\left\|\mathbb{E}_{W|A,Z}\left[\tilde{h}(A,W)\right](A,Z)-\hat{\mathbb{E}}_{W|A,Z}\left[\tilde{h}(A,W)\right](A,Z)\right\|_{P(A,Z)}$
$\displaystyle\quad\quad+\left\|\mathbb{E}_{W|A,Z}\left[\tilde{h}(A,W)\right](A,Z)-\mathbb{E}_{Y|A,Z}\left[Y\right](A,Z)\right\|_{P(A,Z)}+\sqrt{4\hat{\mathfrak{R}}_{S_{2}}(\mathcal{H}_{2})+24M^{2}\sqrt{\frac{\log
2/\delta}{2n}}}$
$\displaystyle\overset{(b)}{\leq}\kappa_{2}+M\sqrt{\kappa_{1}+4\hat{\mathfrak{R}}_{S_{1}}(\mathcal{H}_{2})+24\sqrt{\frac{\log
2/\delta}{2m}}}+\sqrt{4\hat{\mathfrak{R}}_{S_{2}}(\mathcal{H}_{2})+24M^{2}\sqrt{\frac{\log
2/\delta}{2n}}}$

where (a) holds from the triangular inequality and (b) holds from Lemma 3. ∎

Given the generalization errors in both stages, we can bound the error in (11)
for the estimated bridge function $\hat{h}$. We need Assumption 3 to connect
the error in (11) and the error
$\|\hat{h}-h^{*}\|_{\mathbb{P}\left(A,W\right)}$. Let us restate Assumption 3.

See 3 Given this assumption, we can consider the following constant
$\tau_{a}$.

$\displaystyle\tau_{a}=\max_{h\in L^{2}(\mathbb{P}\left(W|A=a\right)),h\neq
h^{*}}\frac{\|h^{*}(a,W)-h(a,W)\|_{\mathbb{P}\left(W|A=a\right)}}{\|\mathbb{E}_{W|A=a,Z}\left[h^{*}(a,W)\right]-\mathbb{E}_{W|A=a,Z}\left[h(a,W)\right]\|_{\mathbb{P}\left(Z|A=a\right)}}$

Note that Assumption 3 ensures $\tau_{a}<\infty$. We can bound the error using
the supremum of this constant, $\tau=\sup_{a}\tau_{a}$.

######  Theorem 2.

Let Assumptions 3 and 7 hold. Given stage 1 data
$S_{1}=\\{(w_{i},a_{i},z_{i})\\}_{i=1}^{m}$ and stage 2 data
$S_{2}=\\{(\tilde{y}_{i},\tilde{a}_{i},\tilde{z}_{i})\\}_{i=1}^{n}$, for any
$\delta>0$, with at least probability of $1-6\delta$, we have

$\displaystyle\|h^{*}(A,W)-\hat{h}(A,W)\|_{P(A,W)}$
$\displaystyle\quad\leq\tau\left(\kappa_{2}+2M\sqrt{\kappa_{1}+4\hat{\mathfrak{R}}_{S_{1}}(\mathcal{H}_{1})+24\sqrt{\frac{\log
2/\delta}{2m}}}+\sqrt{4\hat{\mathfrak{R}}_{S_{2}}(\mathcal{H}_{2})+24M^{2}\sqrt{\frac{\log
2/\delta}{2n}}}\right),$

where $\tau=\sup_{a}\tau_{a}$ and
$\mathcal{H}_{1},\kappa_{1},\mathcal{H}_{2},\kappa_{2}$ are defined in (12),
(13), (15), (16), respectively.

###### Proof.

$\displaystyle\|h^{*}(A,W)-\hat{h}(A,W)\|_{P(A,W)}$
$\displaystyle\leq\sqrt{\mathbb{E}_{A}\left[\|h^{*}(a,W)-\hat{h}(a,W)\|^{2}_{\mathbb{P}\left(W|A=a\right)}\right]}$
$\displaystyle\leq\sqrt{\mathbb{E}_{A}\left[\tau_{a}^{2}\left\|\mathbb{E}_{W|A=a,Z}\left[h^{*}(a,W)-\hat{h}(a,W)\right]\right\|^{2}_{\mathbb{P}\left(Z|A=a\right)}\right]}$
$\displaystyle\leq\tau\sqrt{\mathbb{E}_{A}\left[\left\|\mathbb{E}_{W|A=a,Z}\left[h^{*}(a,W)-\hat{h}(a,W)\right]\right\|^{2}_{\mathbb{P}\left(Z|A=a\right)}\right]}$
$\displaystyle\leq\tau\left\|\mathbb{E}_{Y|A,Z}\left[Y\right]-\mathbb{E}_{W|A,Z}\left[\hat{h}(A,W)\right]\right\|_{\mathbb{P}\left(A,Z\right)}$
$\displaystyle\leq\tau\left\|\mathbb{E}_{Y|A,Z}\left[Y\right]-\hat{\mathbb{E}}_{W|A,Z}\left[\hat{h}(A,W)\right]\right\|_{\mathbb{P}\left(A,Z\right)}+\tau\left\|\mathbb{E}_{W|A,Z}\left[\hat{h}(A,W)\right]-\hat{\mathbb{E}}_{W|A,Z}\left[\hat{h}(A,W)\right]\right\|_{\mathbb{P}\left(A,Z\right)}$

Using Lemmas 3 and 4, the result thus follows. ∎

From this result, we obtain the following corollary.

######  Corollary 3.

Let Assumption 7 hold and $\kappa_{1},\kappa_{2}=0$. If
$\hat{\mathfrak{R}}_{S_{1}}(\mathcal{H}_{1})\to 0$ and
$\hat{\mathfrak{R}}_{S_{2}}(\mathcal{H}_{2})\to 0$ in probability as the
dataset size increases, $\hat{h}$ converges to ${h^{*}}$ in probability with
respect to $\|\cdot\|_{\mathbb{P}\left(A,W\right)}$.

###  C.3 Consistency Result for Causal Parameters

In this section, we develop a consistency result for causal parameters
discussed in the main body. First, we consider the structural function. Given
estimated bridge function
$\hat{h}(a,w)=\hat{\bm{u}}^{\top}(\bm{\psi}_{\hat{\theta}_{A(2)}}(a)\otimes\bm{\psi}_{\hat{\theta}_{W}}(w))$,
we estimate the structural function by taking the empirical mean over $W$. To
make the discussion simple, we assume the access to an additional data sample
$S_{W}=\\{w^{\mathrm{extra}}_{i}\\}_{i=1}^{n_{W}}$, such that the estimated
structural function is given as

$\displaystyle\hat{f}_{\mathrm{struct}}(a)=\hat{\bm{u}}^{\top}(\bm{\psi}_{\hat{\theta}_{A(2)}}(a)\otimes\bm{\mu}_{\hat{\theta}_{W}}),$
(17)

where

$\displaystyle\bm{\mu}_{\hat{\theta}_{W}}=\sum_{i=1}^{n_{W}}\bm{\psi}_{\hat{\theta}_{W}}(w^{\mathrm{extra}}_{i}).$

Note that empirically, we can use outcome-proxy data in $S_{1}$ instead of
$S_{W}$.

To bound the deviation from the true structural function, we need the
following assumption.

######  Assumption 8.

Let $\rho_{W}(w),\rho_{W|A}(w|a)$ be the respective density functions of
probability distributions
$\mathbb{P}\left(W\right),\mathbb{P}\left(W|A\right)$. The densities satisfy

$\displaystyle\eta_{a}\coloneqq\left\|\frac{\rho_{W}(W)}{\rho_{W|A}(W|a)}\right\|_{\mathbb{P}\left(W|A=a\right)}<\infty.$

Given Assumption 8, we can bound the error in structural function estimation.
Before stating the theorem, let us introduce a useful concentration inequality
for multi-dimensional random variables. 222Lemma 5 is discussed in
MathOverflow (https://mathoverflow.net/questions/186097/hoeffdings-inequality-
for-vector-valued-random-variables, accessed on June 2nd 2021)

######  Lemma 5.

Let $\bm{x}_{1},\dots,\bm{x}_{n}\in[-1,1]^{d}$ be independent random
variables. Then, with probability at least $1-\delta$, we have

$\displaystyle\left\|\frac{1}{n}\sum_{i=1}^{n}\bm{x}_{i}-\bm{\mu}\right\|\leq\sqrt{\frac{2d\log
2d/\delta}{n}}$

where $\bm{\mu}=\mathbb{E}\left[\frac{1}{n}\sum_{i=1}^{n}\bm{x}_{i}\right]$.

###### Proof.

Let $j$-th coordinate of $\bm{x}_{i}$ be denoted as $(\bm{x}_{i})_{j}$. Then,

$\displaystyle\mathbb{P}\left(\left\|\frac{1}{n}\sum_{i=1}^{n}\bm{x}_{i}-\bm{\mu}\right\|\leq\varepsilon\right)$
$\displaystyle=\mathbb{P}\left(\sum_{j=1}^{d}\left(\frac{1}{n}\sum_{i=1}^{n}(\bm{x}_{i})_{j}-(\bm{\mu})_{j}\right)^{2}\leq\varepsilon^{2}\right)$
$\displaystyle\leq\mathbb{P}\left(\bigcap_{j=1}^{d}\left|\frac{1}{n}\sum_{i=1}^{n}(\bm{x}_{i})_{j}-(\bm{\mu})_{j}\right|\leq\frac{\varepsilon}{\sqrt{d}}\right)$
$\displaystyle\leq 1-2d\exp\left(-\frac{n\varepsilon^{2}}{2d}\right),$

Where the last inequality holds from Hoeffding’s inequality. The claim follows
by solving $\delta=d\exp\left(-\frac{n\varepsilon^{2}}{2d}\right)$. ∎

Given this lemma, we can prove the following result.

######  Theorem 3.

Let Assumptions 3, 7 and 8 hold. Let $\eta=\sup_{a\in\mathcal{A}}\eta_{a}$.
Given stage 1 data $S_{1}=\\{(w_{i},a_{i},z_{i})\\}_{i=1}^{m}$, stage 2 data
$S_{2}=\\{(\tilde{y}_{i},\tilde{a}_{i},\tilde{z}_{i})\\}_{i=1}^{n}$,
additional outcome-proxy variable data
$S_{W}=\\{w^{\mathrm{extra}}_{i}\\}_{i=1}^{n_{W}}$ , then with probability at
least $1-7\delta$, we have

$\displaystyle\|f^{*}_{\mathrm{struct}}-\hat{f}_{\mathrm{struct}}\|_{\mathbb{P}\left(A\right)}$
$\displaystyle\leq\sqrt{\frac{2d_{W}\log
2d_{W}/\delta}{n_{W}}}+\eta\tau\left(\kappa_{2}+2M\sqrt{\kappa_{1}+4\hat{\mathfrak{R}}_{S_{1}}(\mathcal{H}_{1})+24\sqrt{\frac{\log
2/\delta}{2m}}}+\sqrt{4\hat{\mathfrak{R}}_{S_{2}}(\mathcal{H}_{2})+24M^{2}\sqrt{\frac{\log
2/\delta}{2n}}}\right),$

where $\hat{f}_{\mathrm{struct}}$ is given in (17) and $d_{W}$ is the
dimension of $\bm{\psi}_{\theta_{W}}$.

###### Proof.

From the relationship between structural function and bridge function, we have

$\displaystyle\|f^{*}_{\mathrm{struct}}-\hat{f}_{\mathrm{struct}}\|_{\mathbb{P}\left(A\right)}$
$\displaystyle=\|\mathbb{E}_{W}\left[h^{*}(A,W)\right]-\hat{f}_{\mathrm{struct}}\|_{\mathbb{P}\left(A\right)}$
$\displaystyle\leq\left\|\mathbb{E}_{W}\left[h^{*}(A,W)\right]-\mathbb{E}_{W}\left[\hat{h}(A,W)\right]\right\|_{\mathbb{P}\left(A\right)}+\left\|\mathbb{E}_{W}\left[\hat{h}(A,W)\right]-\hat{f}_{\mathrm{struct}}\right\|_{\mathbb{P}\left(A\right)}.$

We can bound each term as follows. For the first term, we have

$\displaystyle\left\|\mathbb{E}_{W}\left[h^{*}(\cdot,W)\right]-\mathbb{E}_{W}\left[\hat{h}(\cdot,W)\right]\right\|_{\mathbb{P}\left(A\right)}$
$\displaystyle=\sqrt{\int\left(\mathbb{E}_{W}\left[h^{*}(a,W)-\hat{h}(a,W)\right]\right)^{2}\rho_{A}(a)\mathrm{d}a}$
$\displaystyle=\sqrt{\int\left(\mathbb{E}_{W|A=a}\left[(h^{*}(a,W)-\hat{h}(a,W))\frac{\rho_{W}(W)}{\rho_{W|A}(W|a)}\right]\right)^{2}\rho_{A}(a)\mathrm{d}a}$
$\displaystyle\leq\sqrt{\int\left\|h^{*}(a,W)-\hat{h}(a,W)\right\|^{2}_{\mathbb{P}\left(W|A=a\right)}\left\|\frac{\rho_{W}(W)}{\rho_{W|A}(W|a)}\right\|^{2}_{\mathbb{P}\left(W|A=a\right)}\rho_{A}(a)\mathrm{d}a}\quad\because\text{~{}Cauchy\textendash
Schwarz inequality}$
$\displaystyle\leq\sqrt{\int\eta^{2}_{a}\mathbb{E}_{W|A=a}\left[(h^{*}(a,W)-\hat{h}(a,W))^{2}\right]\rho_{A}(a)\mathrm{d}a}\quad\because\text{~{}Assumption~{}\ref{assume-
good-support}}$
$\displaystyle\leq\eta\sqrt{\mathbb{E}_{W,A}\left[(h^{*}(A,W)-\hat{h}(A,W))^{2}\right]}$
$\displaystyle=\eta\|h^{*}-\hat{h}\|_{\mathbb{P}\left(A,W\right)}.$

From Theorem 2, with probability at least $1-6\delta$, we have

$\displaystyle\left\|\mathbb{E}_{W}\left[h^{*}(\cdot,W)\right]-\mathbb{E}_{W}\left[\hat{h}(\cdot,W)\right]\right\|_{\mathbb{P}\left(A\right)}$
$\displaystyle\leq\eta\tau\left(\kappa_{2}+2M\sqrt{\kappa_{1}+4\hat{\mathfrak{R}}_{S_{1}}(\mathcal{H}_{1})+24\sqrt{\frac{\log
2/\delta}{2m}}}+\sqrt{4\hat{\mathfrak{R}}_{S_{2}}(\mathcal{H}_{2})+24M^{2}\sqrt{\frac{\log
2/\delta}{2n}}}\right).$

For the second term, from Assumption 7, we have with probability at least
$1-\delta$,

$\displaystyle\left\|\mathbb{E}_{W}\left[\hat{h}(A,W)\right]-\hat{f}_{\mathrm{struct}}\right\|_{\mathbb{P}\left(A\right)}$
$\displaystyle\leq
M\left\|\bm{\mu}_{\hat{\theta}_{W}}-\mathbb{E}_{W}\left[\bm{\psi}_{\hat{\theta}_{W}}\right]\right\|$
$\displaystyle\leq\sqrt{\frac{2d_{W}\log 2d_{W}/\delta}{n_{W}}}$

The last inequality holds from Lemma 5. Using the uniform inequality, we have
shown the claim. ∎

We can evaluate the error in estimating a value function as well. Given
$S_{3}=\\{\check{w}_{i},\check{c}_{i}\\}_{i=1}^{n^{\prime}}$, we estimate the
value function as

$\displaystyle\hat{v}(\pi)=\frac{1}{n^{\prime}}\sum_{i=1}^{n}\hat{\bm{u}}^{\top}(\bm{\psi}_{\hat{\theta}_{A(2)}}(\pi(\check{c}_{i}))\otimes\bm{\psi}_{\hat{\theta}_{W}}(\check{w}_{i})),$

Furthermore, we assume the following relationship between distributions of $A$
and $C$

######  Assumption 9.

There exists a constant $\sigma$ such that

$\displaystyle\|l(\pi(C),W)\|_{\mathbb{P}\left(C,W\right)}\leq\sigma\|l(A,W)\|_{\mathbb{P}\left(A,W\right)}$

for all square integrable functions
$l:\mathcal{A}\times\mathcal{W}\to\mathbb{R},\|l\|<\infty$.

Given these assumptions, we have the following theorem.

######  Theorem 4.

Let Assumptions 3, 7, 8, 9 hold. Given stage 1 data
$S_{1}=\\{(w_{i},a_{i},z_{i})\\}_{i=1}^{m}$, stage 2 data
$S_{2}=\\{(\tilde{y}_{i},\tilde{a}_{i},\tilde{z}_{i})\\}_{i=1}^{n}$, stage 3
data $S_{3}=\\{\check{w}_{i},\check{c}_{i}\\}_{i=1}^{n^{\prime}}$, with at
least probability $1-7\delta$, we have

$\displaystyle|v(\pi)-\hat{v}(\pi)|$
$\displaystyle\leq\sigma\tau\left(\kappa_{2}+2M\sqrt{\kappa_{1}+4\hat{\mathfrak{R}}_{S_{1}}(\mathcal{H}_{1})+24\sqrt{\frac{\log
2/\delta}{2m}}}+\sqrt{4\hat{\mathfrak{R}}_{S_{2}}(\mathcal{H}_{2})+24M^{2}\sqrt{\frac{\log
2/\delta}{2n}}}\right)$ $\displaystyle\quad\quad+\sqrt{\frac{M^{2}\log
2/\delta}{n^{\prime}}}.$

###### Proof.

We have

$\displaystyle|v(\pi)-\hat{v}(\pi)|$
$\displaystyle=\left|\mathbb{E}_{W|C}\left[h^{*}(\pi(C),W)\right]-\hat{v}(\pi)\right|$
$\displaystyle\leq\left|\mathbb{E}_{W,C}\left[h(\pi(C),W)\right]-\mathbb{E}_{W,C}\left[\hat{h}(\pi(C),W)\right]\right|$
$\displaystyle\quad\quad+\left|\mathbb{E}_{W,C}\left[\hat{h}(\pi(C),W)\right]-\hat{v}(\pi)\right|.$

We bound each term as follows. For the first term, we have

$\displaystyle\left|\mathbb{E}_{W,C}\left[h^{*}(\pi(C),W)\right]-\mathbb{E}_{W,C}\left[\hat{h}(\pi(C),W)\right]\right|$
$\displaystyle\leq\|h^{*}(\pi(C),W)-\hat{h}(\pi(C),W)\|_{\mathbb{P}\left(C,W\right)}\quad\because\text{Jensen's
inequality}$
$\displaystyle\leq\sigma\|h^{*}(A,W)-\hat{h}(A,W)\|_{\mathbb{P}\left(A,W\right)}\quad\because\text{Assumption~{}\ref{assume-
policy-l2}}$

Hence, from Theorem 2, we have

$\displaystyle\left|\mathbb{E}_{W,C}\left[h^{*}(\pi(C),W)\right]-\mathbb{E}_{W,C}\left[\hat{h}(\pi(C),W)\right]\right|$
$\displaystyle\leq\sigma\tau\left(\kappa_{2}+2M\sqrt{\kappa_{1}+4\hat{\mathfrak{R}}_{S_{1}}(\mathcal{H}_{1})+24\sqrt{\frac{\log
2/\delta}{2m}}}+\sqrt{4\hat{\mathfrak{R}}_{S_{2}}(\mathcal{H}_{2})+24M^{2}\sqrt{\frac{\log
2/\delta}{2n}}}\right)$

with probability at least $1-6\delta$.

For the second term, we have

$\displaystyle\left|\mathbb{E}_{W,C}\left[\hat{\bm{u}}^{\top}\left(\bm{\psi}_{\hat{\theta}_{A(2)}}(\pi(C))\otimes\bm{\psi}_{\hat{\theta}_{W}}(W)\right)\right]-\hat{v}(\pi)\right|\leq\sqrt{\frac{M^{2}\log
2/\delta}{n^{\prime}}}$

with probability at least $1-\delta$, from Hoeffding’s inequality. The result
is obtained by taking the union bound. ∎

##  Appendix D DFPV algorithm with observable confounders

In this appendix, we formulate the DFPCL method when observable confounders
are present, building on Tchetgen et al. [31] and Mastouri and Zhu et al.
[18]. Here, we consider the causal graph given in Figure 5. In addition to
variables $(A,Y,Z,W)$, we have an observable confounder $X\in\mathcal{X}$. The
structural function $f_{\mathrm{struct}}$ we aim to learn is

$\displaystyle
f_{\mathrm{struct}}(a)=\mathbb{E}_{X,U}\left[\mathbb{E}\left[Y\mid
A=a,X,U\right]\right].$

$A$$Y$$Z$$W$$U$$X$ Figure 5: Causal graph with observable confounder

The structural assumption and completeness assumption including observable
confounders are given as follows.

######  Assumption 10 (Structural Assumption [18]).

We assume $Y\mathop{\perp\\!\\!\\!\perp}Z|A,U,X$, and
$W\mathop{\perp\\!\\!\\!\perp}(A,Z)|U,X$.

######  Assumption 11 (Completeness Assumption [18]).

Let $l:\mathcal{U}\to\mathbb{R}$ be any square integrable function
$\|l\|_{\mathbb{P}\left(U\right)}$. We assume the following:

$\displaystyle\mathbb{E}\left[l(U)\mid
A=a,Z=z,X=x\right]=0\quad\forall(a,z,x)\in\mathcal{A}\times\mathcal{Z}\times\mathcal{X}\quad\Leftrightarrow\quad
l(u)=0~{}\mathrm{a.s.}$ $\displaystyle\mathbb{E}\left[l(U)\mid
A=a,Z=z,X=x\right]=0\quad\forall(a,z,x)\in\mathcal{A}\times\mathcal{Z}\times\mathcal{X}\quad\Leftrightarrow\quad
l(u)=0~{}\mathrm{a.s.}$

Following similar reasoning as in Section 2, we can estimate the bridge
function $\hat{h}:\mathcal{A}\times\mathcal{X}\times\mathcal{W}\to\mathbb{R}$
by minimizing the following loss:

$\displaystyle\hat{h}=\operatorname*{arg\,min}_{h\in\mathcal{H}_{h}}\tilde{\mathcal{L}}(h),\quad\tilde{\mathcal{L}}(h)=\mathbb{E}_{Y,A,Z,X}\left[(Y-\mathbb{E}_{W|Z,A,X}\left[h(A,X,W)\right])^{2}\right]+\Omega(h).$

Given bridge function, we can estimate the structural function by

$\displaystyle
f_{\mathrm{struct}}(a)=\mathbb{E}_{X,W}\left[h^{*}(a,X,W)\right].$

Similar to Mastouri and Zhu et al. [18], we model

$\displaystyle\mathbb{E}_{W|a,z}\left[\bm{\psi}_{\theta_{W}}(w)\right]=\bm{V}(\bm{\phi}_{\theta_{A(1)}}(A)\otimes\bm{\phi}_{\theta_{Z}}(Z))\otimes\bm{\phi}_{\theta_{X(1)}}(X))$
$\displaystyle
h(a,x,w)=\bm{u}^{\top}(\bm{\psi}_{\theta_{A(2)}}(a)\otimes\bm{\psi}_{\theta_{X(2)}}(x)\otimes\bm{\psi}_{\theta_{W}}(w)),$

where $\bm{\phi}_{\theta_{X(1)}}(X),\bm{\psi}_{\theta_{X(2)}}(X)$ are the
feature maps of $X$ parameterized by $\theta_{X(1)},\theta_{X(2)}$,
respectively. Then, in stage 1, we learn
$(\bm{V},\theta_{A(1)},\theta_{Z},\theta_{X(1)})$ by minimizing

$\displaystyle\hat{\mathcal{L}}_{1}(\bm{V},\theta_{A(1)},\theta_{Z},\theta_{X(1)})=\frac{1}{m}\sum_{i=1}^{m}\left\|\bm{\psi}_{\theta_{W}}(w_{i})-\bm{V}\left(\bm{\phi}_{\theta_{A(1)}}(a_{i})\otimes\bm{\phi}_{\theta_{Z}}(z_{i})\otimes\bm{\phi}_{\theta_{X}}(x_{i})\right)\right\|^{2}+\lambda_{1}\|\bm{V}\|^{2},$
(18)

which estimates the conditional expectation
$\mathbb{E}_{W}\left[\bm{\psi}_{W}(W)|A=a,Z=z,X=x\right]$. Let
$(\hat{\bm{V}},\hat{\theta}_{A(1)},\hat{\theta}_{Z},\hat{\theta}_{X(1)})$ be
the minimizer of $\hat{\mathcal{L}}_{1}$. Then, in stage 2, we can learn
$(\bm{u},\theta_{A(2)},\theta_{Z},\theta_{X(2)})$ using

$\displaystyle\hat{\mathcal{L}}_{2}(\bm{u},\theta_{A(2)},\theta_{W},\theta_{X(2)})=\frac{1}{n}\sum_{i=1}^{n}\left(\tilde{y}_{i}-\bm{u}^{\top}\left(\bm{\psi}_{\theta_{A(2)}}(\tilde{a}_{i})\otimes\bm{\psi}_{\theta_{X(2)}}(\tilde{x}_{i})\otimes\hat{\bm{V}}\bm{v}_{1}(\tilde{a}_{i},\tilde{x}_{i},\tilde{z}_{i})\right)\right)^{2}+\lambda_{2}\|\bm{u}\|^{2},$
(19)

where we denote
$\bm{v}_{1}(a,x,z)=\left(\bm{\phi}_{\hat{\theta}_{A(1)}}(\tilde{a}_{i})\otimes\bm{\phi}_{\hat{\theta}_{X}}(\tilde{x}_{i})\otimes\bm{\phi}_{\hat{\theta}_{Z}}(\tilde{z}_{i})\right)$.

In DFPV, we first fix parameters
$\bm{\theta}=(\theta_{A(1)},\theta_{Z},\theta_{X(1)},\theta_{A(2)},\theta_{Z},\theta_{X(2)})$
and obtain weights. This is given as

$\displaystyle\hat{\bm{V}}(\bm{\theta})=\Psi_{1}^{\top}\Phi_{1}(\Phi_{1}^{\top}\Phi_{1}+m\lambda_{1}I)^{-1},$
$\displaystyle\hat{\bm{u}}(\bm{\theta})=\left(\Phi_{2}^{\top}\Phi_{2}+n\lambda_{2}I\right)^{-1}\Phi_{2}^{\top}\bm{y}_{2},$
(20)

where we denote
$\bm{\theta}=(\theta_{A(1)},\theta_{Z},\theta_{A(2)},\theta_{W})$ and define
matrices as follows:

$\displaystyle\Psi_{1}=\left[\bm{\psi}_{\theta_{W}}(w_{1}),\dots,\bm{\psi}_{\theta_{W}}(w_{m})\right]^{\top},$
$\displaystyle\Phi_{1}=[\bm{v}_{1}(a_{1},z_{1}),\dots,\bm{v}_{1}(a_{m},z_{m})]^{\top},$
$\displaystyle\bm{y}_{2}=[\tilde{y}_{1},\dots,\tilde{y}_{n}]^{\top},$
$\displaystyle\Phi_{2}=[\bm{v}_{2}(\tilde{a}_{1},\tilde{z}_{1}),\dots,\bm{v}_{2}(\tilde{a}_{n},\tilde{z}_{n})]^{\top},$
$\displaystyle\bm{v}_{1}(a,x,z)=\bm{\phi}_{\theta_{A(1)}}(a)\otimes\bm{\phi}_{\theta_{X(1)}}(x)\otimes\bm{\phi}_{\theta_{Z}}(z),$
$\displaystyle\bm{v}_{2}(a,z)=\bm{\psi}_{\theta_{A(2)}}(a)\otimes\bm{\psi}_{\theta_{X(2)}}(x)\otimes\left(\hat{\bm{V}}(\bm{\theta})\bm{v}_{1}(a,x,z)\right).$

We learn the parameters by minimizing the following:

$\displaystyle\hat{\mathcal{L}}^{\mathrm{DFPV}}_{1}(\bm{\theta})=\frac{1}{m}\sum_{i=1}^{m}\left\|\bm{\psi}_{\theta_{W}}(w_{i})-\hat{\bm{V}}(\bm{\theta})\bm{v}_{1}(a_{i},x_{i},z_{i})\right\|^{2}+\lambda_{1}\|\hat{\bm{V}}(\bm{\theta})\|^{2},$
$\displaystyle\hat{\mathcal{L}}^{\mathrm{DFPV}}_{2}(\bm{\theta})=\frac{1}{n}\sum_{i=1}^{n}\left(\tilde{y}_{i}-\hat{\bm{u}}(\bm{\theta})^{\top}\bm{v}_{2}(\tilde{a}_{i},\tilde{x}_{i},\tilde{z}_{i})\right)^{2}\\!+\\!\lambda_{2}\|\hat{\bm{u}}(\bm{\theta})\|^{2}.$

The algorithm is given in Algorithn 2.

Algorithm 2 Deep Feature Instrumental Variable with Observable Confounder

0: Stage 1 data $(a_{i},z_{i},w_{i},x_{i})$, Stage 2 data
$(\tilde{a}_{i}\tilde{z}_{i},\tilde{y}_{i},\tilde{x}_{i})$, Regularization
parameters $(\lambda_{1},\lambda_{2})$. Initial values
$\bm{\theta}^{(0)}=(\theta^{(0)}_{A(1)},\theta^{(0)}_{Z},\theta^{(0)}_{A(2)},\theta^{(0)}_{W},\theta^{(0)}_{X(1)},\theta^{(0)}_{X(2)})$.
Learning rate $\alpha$, additional data
$(x_{i}^{\mathrm{extra}},w_{i}^{\mathrm{extra}})$

0: Estimated structural function $\hat{f}_{\mathrm{struct}}(a)$

1: $t\leftarrow 0$

2: repeat

3: Compute $\hat{\bm{V}}(\bm{\theta}),\hat{\bm{u}}(\bm{\theta})$ in (20)

4: Update parameters in features
$\bm{\theta}^{(t+1)}\leftarrow(\theta^{(t+1)}_{A(1)},\theta^{(t+1)}_{Z},\theta^{(t+1)}_{A(2)},\theta^{(t+1)}_{W},\theta^{(t+1)}_{X(1)},\theta^{(t+1)}_{X(2)})$
by
$\displaystyle\theta^{(t+1)}_{A(1)}\leftarrow\theta^{(t)}_{A(1)}-\alpha\nabla_{\theta_{A(1)}}\hat{\mathcal{L}}^{\mathrm{DFPV}}_{1}(\bm{\theta})|_{\bm{\theta}=\bm{\theta}^{(t)}},$
$\displaystyle\theta^{(t+1)}_{Z}\leftarrow\theta^{(t)}_{Z}-\alpha\nabla_{\theta_{Z}}\hat{\mathcal{L}}^{\mathrm{DFPV}}_{1}(\bm{\theta})|_{\bm{\theta}=\bm{\theta}^{(t)}}$
$\displaystyle\theta^{(t+1)}_{X(1)}\leftarrow\theta^{(t)}_{X(1)}-\alpha\nabla_{\theta_{X(1)}}\hat{\mathcal{L}}^{\mathrm{DFPV}}_{1}(\bm{\theta})|_{\bm{\theta}=\bm{\theta}^{(t)}},$
$\displaystyle\theta^{(t+1)}_{X(2)}\leftarrow\theta^{(t)}_{X(2)}-\alpha\nabla_{\theta_{X(2)}}\hat{\mathcal{L}}^{\mathrm{DFPV}}_{2}(\bm{\theta})|_{\bm{\theta}=\bm{\theta}^{(t)}}$
$\displaystyle\theta^{(t+1)}_{A(2)}\leftarrow\theta^{(t)}_{A(2)}-\alpha\nabla_{\theta_{A(2)}}\hat{\mathcal{L}}^{\mathrm{DFPV}}_{2}(\bm{\theta})|_{\bm{\theta}=\bm{\theta}^{(t)}},$
$\displaystyle\theta^{(t+1)}_{W}\leftarrow\theta^{(t)}_{W}-\alpha\nabla_{\theta_{W}}\hat{\mathcal{L}}^{\mathrm{DFPV}}_{2}(\bm{\theta})|_{\bm{\theta}=\bm{\theta}^{(t)}}$

5: Increment counter $t\leftarrow t+1;$

6: until convergence

7: Compute $\hat{\bm{u}}(\bm{\theta}^{(t)})$ from (20)

8: Compute mean feature for $W$ using stage 1 dataset
$\displaystyle\bm{\mu}_{\theta_{X(2)}\otimes\theta_{W}}\leftarrow\frac{1}{n}\sum\bm{\psi}_{\theta^{(t)}_{X(2)}}(x^{\mathrm{extra}}_{i})\otimes\bm{\psi}_{\theta^{(t)}_{W}}(w^{\mathrm{extra}}_{i})$

9: return
$\hat{f}_{\mathrm{struct}}(a)=(\hat{\bm{u}}^{(t)})^{\top}\left(\bm{\psi}_{\hat{\theta}^{(t)}_{A(2)}}(a)\otimes\bm{\mu}_{\theta_{X(2)}\otimes\theta_{W}}\right)$

##  Appendix E Additional Experiments

In this appendix, we report the results of two additional experiments. One is
a synthetic setting introduced in Mastouri et al. [18], which has a simpler
data generating process. The other is based on the real-world setting
introduced by Deaner [3]. In both setting, DFPV performs similarly to or
better than existing methods.

###  E.1 Experiments with Simpler Data Generating Process

Here, we show the result for the synthetic setting proposed in Mastouri et al.
[18]. The data generating process for each variable is given as follows:

$\displaystyle U:=[U_{1},U_{2}],\quad U_{2}\sim\mathrm{Unif}[-1,2]\quad
U_{1},\sim\mathrm{Unif}[0,1]-\mathbbm{1}[0\leq U_{2}\leq 1]$ $\displaystyle
Z:=[U_{1}+\mathrm{Unif}[-1,1],~{}U_{2}+\mathcal{N}(0,3)]$ $\displaystyle
W:=[U_{1}+\mathcal{N}(0,3),~{}U_{2}+\mathrm{Unif}[-1,1]]$ $\displaystyle
A:=U_{2}+\mathcal{N}(0,0.05)$ $\displaystyle Y:=U_{2}\cos(2(A+0.3U_{1}+0.2))$

From observations of $(Y,W,Z,A)$, we estimate $\hat{f}_{\mathrm{struct}}$ by
PCL. For each estimated $\hat{f}_{\mathrm{struct}}$, we measure out-of-sample
error as the mean square error of $\hat{f}$ versus true $f_{\mathrm{struct}}$
obtained from Monte-Carlo simulation. Specifically, we consider 20 evenly
spaced values of $A\in[0.0,1.0]$ as the test data. The results with data size
$n=m=\\{500,1000\\}$ are shown in Figure 6.

Figure 6: Result of structural function experiment in the setting in Mastouri
et al. [18]

From Figure 6, we can see that DFPV and CEVAE methods perform worse and have
larger variances than KPV and PMMR methods. This is not surprising, since DFPV
tends to require more data than KPV and PMMR, as needed to learn the neural
net feature maps (rather than using fixed pre-defined kernel features). Hence,
we can say that we should favor KPV and PMMR over DFPV when the data is low-
dimensional and the relations between the variables are smooth. We would like
to note, however, DFPV outperforms CEVAE, which shows that the proxy setting
is still required.

###  E.2 Experiments using Grade Retention dataset

To test the performance of DFPV in a more realistic setting, we conducted the
experiment on the Grade Retention dataset introduced by Deaner [3]. This aims
to estimate the effect of grade retention based on the score of math and
reading on the long-term cognitive outcomes, in which we use scores in
elementary school as a treatment-inducing proxy (Z) and cognitive test scores
from Kindergarten as the an outcome-inducing proxy (W). Following Mastouri et
al. [18], we generate a synthetic "ground truth" by fitting a generalized
additive model to learn a structured causal model (SCM), and a Gaussian
mixture model to learn unmeasured confounder based on the learned SCM. Note,
this is needed since for real-world data there is no measured ground truth.

| DFPV | CEVAE | KPV | PMMR
---|---|---|---|---
Math | 0.023(0.001) | 0.054(0.007) | 0.043(0.000) | 0.032(0.001)
Reading | 0.027(0.002) | 0.082(0.007) | 0.028(0.000) | 0.022(0.000)
Table 1: Results of grade retension dataset

Table 1 shows the result for this dataset. In this setting, the performance of
DFPV matches KPV and PMMR. As in the experiment described in the revious
section, the setting is low-dimensional (one-dim treatment variable, three-dim
treatment-inducing proxy, four-dim outcome-inducing proxy) and the generative
model is smooth (the "ground truth" being a generalized additive model and a
Gaussian mixture model). For these reasons, we might again expect this data to
favor kernel methods, such as KPV and PMMR; nonetheless, our method matches
them. DFPV again outperforms CEVAE in this setting.

##  Appendix F Experiment Details

In this section, we present the data generation process of experiments and the
detailed settings of hyper-parameters.

###  F.1 Demand Design Experiment

Here, we introduce the details of demand design experiments. The observations
are generated from the following causal model,

$\displaystyle Y=P\left(\exp\left(\frac{V-P}{10}\right)\wedge
5\right)-5g(D)+\varepsilon,\quad\varepsilon\sim\mathcal{N}(0,1),$

where $Y$ represents sales, $P$ is the treatment variable (price), and these
are confounded by potential demand $D$. Here we denote $a\wedge b=\min(a,b)$,
and the function $g$ as

$\displaystyle
g(d)=2\left(\frac{(d-5)^{4}}{600}+\exp(-4(d-5)^{2})+\frac{d}{10}-2\right).$

To correct this confounding bias, we introduce cost-shifter $C_{1},C_{2}$ as a
treatment-inducing proxy, and views $V$ of the reservation page as the
outcome-inducing proxy. Data is sampled as

$\displaystyle D\sim\mathrm{Unif}[0,10]$ $\displaystyle C_{1}\sim
2\sin(2D\pi/10)+\varepsilon_{1}$ $\displaystyle C_{2}\sim
2\cos(2D\pi/10)+\varepsilon_{2}$ $\displaystyle V\sim
7g(D)+45+\varepsilon_{3}$ $\displaystyle
P=35+(C_{1}+3)g(D)+C_{2}+\varepsilon_{4}$

where
$\varepsilon_{1},\varepsilon_{2},\varepsilon_{3},\varepsilon_{4}\sim\mathcal{N}(0,1)$.
From observations of $(Y,P,C_{1},C_{2},V)$, we estimate
$\hat{f}_{\mathrm{struct}}$ by PCL. For each estimated
$\hat{f}_{\mathrm{struct}}$, we measure out-of-sample error as the mean square
error of $\hat{f}$ versus true $f_{\mathrm{struct}}$ obtained from Monte-Carlo
simulation. Specifically, we consider 10 evenly spaced values of $p\in[10,30]$
as the test data.

###  F.2 dSprite Experiment

Here, we describe the data generation process for the dSprites dataset
experiment. This is an image dataset parametrized via five latent variables
(shape, scale, rotation, posX and posY). The images are $64\times
64=4096$-dimensional. In this experiment, we fixed the shape parameter to
heart, i.e. we only used the heart-shaped images. The other latent parameters
take values of
$\mathtt{scale}\in[0.5,1],~{}\mathtt{rotation}\in[0,2\pi],~{}\mathtt{posX}\in[0,1],~{}\mathtt{posY}\in[0,1]$.

From this dataset, we generate the treatment variable $A$ and outcome $Y$ as
follows:

  1. 1.

Uniformly samples latent parameters
$(\mathtt{scale},\mathtt{rotation},\mathtt{posX},\mathtt{posY})$.

  2. 2.

Generate treatment variable $A$ as

$\displaystyle
A=\mathtt{Fig}(\mathtt{scale},\mathtt{rotation},\mathtt{posX},\mathtt{posY})+\bm{\eta}_{A}.$

  3. 3.

Generate outcome variable $Y$ as

$\displaystyle
Y=\frac{1}{12}(\mathtt{posY}-0.5)\frac{\|BA\|_{2}^{2}-5000}{1000}+\varepsilon,\quad\varepsilon\sim\mathcal{N}(0,0.5).$

Here, function $\mathtt{Fig}$ returns the corresponding image for the latent
parameters, and $\bm{\eta},\varepsilon$ are noise variables generated from
$\bm{\eta}_{A}\sim\mathcal{N}(0.0,0.1I)$ and
$\varepsilon\sim\mathcal{N}(0.0,0.5)$. Each element of the matrix
$B\in\mathbb{R}^{10\times 4096}$ is generated from $\mathrm{Unif}(0.0,1.0)$
and fixed throughout the experiment. From the data generation process, we can
see that $A$ and $Y$ are confounded by posY. Treatment variable $A$ is given
as a figure corrupted with Gaussian random noise. The variable $\mathtt{posY}$
is not revealed to the model, and there is no observable confounder. The
structural function for this setting is

$\displaystyle f_{\mathrm{struct}}(A)=\frac{\|BA\|_{2}^{2}-5000}{1000}.$

To correct this confounding bias, we set up the following PCL setting. We
define the treatment-inducing variable
$Z=(\mathtt{scale},\mathtt{rotation},\mathtt{posX})\in\mathbb{R}^{3}$, and the
outcome-inducing variable by another figure that shares the same
$\mathtt{posY}$, with the remaining latent parameters fixed as follows:

$\displaystyle W=\mathtt{Fig}(0.8,0,0.5,\mathtt{posY})+\bm{\eta}_{W},$

where $\bm{\eta}_{W}\sim\mathcal{N}(0.0,0.1I)$.

We use 588 test points for measuring out-of-sample error, which are generated
from the grid points of latent variables. The grids consist of 7 evenly spaced
values for $\mathtt{posX},\mathtt{posY}$, 3 evenly spaced values for
$\mathtt{scale}$, and 4 evenly spaced values for $\mathtt{orientation}$.

###  F.3 Policy Evaluation Experiments

We use the same data $(Y,P,C_{1},C_{2},V)$ in demand design for policy
evaluation experiments. We consider two policies. One is a policy depends on
costs $C_{1},C_{2}$ which is

$\displaystyle\pi_{C_{1},C_{2}}(C_{1},C_{2})=23+C_{1}C_{2}.$

To conduct offline-policy evaluation, we use data $(C_{1},C_{2},V)$ to compute
the empirical average of $h(\pi_{C_{1},C_{2}}(C_{1},C_{2}),V)$. In our second
experiment, the policy depends on current price $P$, which is given as

$\displaystyle\pi_{P}(P)=\max(0.7P,10).$

Again, we use data $(P,V)$ to compute the empirical average of
$h(\pi_{P}(P),V)$.

###  F.4 Hyper-parameters and network architectures

Here, we describe the network architecture and hyper-parameters of all
experiments.

For KPV and PMMR method, we used the Gaussian kernel where the bandwidth is
determined by the median trick. We follow the procedure for hyper-parameter
tuning proposed in Mastouri et al. [18] in selecting the regularizers
$\lambda_{1},\lambda_{2}$.

For DFPV, we optimize the model using Adam [10] with learning rate = 0.001,
$\beta_{1}=0.9$, $\beta_{2}=0.999$ and $\varepsilon=10^{-8}$. Regularizers
$\lambda_{1},\lambda_{2}$ are both set to 0.1 as a result of the tuning
procedure described in Appendix A. Network structure is given in Tables 2, 3.

In CEVAE, we attempt to reconstruct the latent variable $L$ from $(A,Z,W,Y)$
using a VAE. Here, we use a 20-dim latent variable $L$, whose the conditional
distribution is specified as follows:

$\displaystyle
q(L|A,Z,W,Y)=\mathcal{N}(\bm{V}_{1}\bm{\psi}_{q}(A,Z,W,Y),\mathrm{diag}(\bm{V}_{2}\bm{\psi}_{q}(A,Z,W,Y)\vee
0.1))$

where $\bm{\psi}_{q}$ is a neural net and $\bm{V}_{1},\bm{V}_{2}$ are matrices
to be learned, and we denote $\bm{a}\vee 0.1=(\max(a_{i},0.1))_{i}$.
Furthermore, we specify the likelihood distribution as follows:

$\displaystyle
p(W,Z|L)=\mathcal{N}(\bm{V}_{3}\bm{\psi}_{p(W,Z|L)}(L),\mathrm{diag}(\bm{V}_{4}\bm{\psi}_{p(W,Z|L)}(L)\vee
0.1))$ $\displaystyle
p(A|L)=\mathcal{N}(\bm{V}_{5}\bm{\psi}_{p(A|L)}(L),\mathrm{diag}(\bm{V}_{6}\bm{\psi}_{p(A|L))}(L)\vee
0.1))$ $\displaystyle p(Y|A,L)=\mathcal{N}(\bm{\mu}_{p(Y|A,L)}(A,L),0.5)$

Here, $\bm{\psi}_{p(W,Z|L)},\bm{\psi}_{p(A|L)},\bm{\mu}_{p(Y|A,L)}$ are neural
networks. We provide the structure of neural nets in Table 4 and 5. Following
the orignal work [17], we train all neural nets by Adamax [10] with a learning
rate of 0.01, which was annealed with an exponential decay schedule. We also
performed early stopping according to the lower bound on a validation set. To
predict structural function, we obtain $q(L)$ by marginalizing $q(L|A,Z,W,Y)$
by observed data $A,Z,W,Y$. We then output
$\hat{f}_{\mathrm{struct}}(a)=\mathbb{E}_{L\sim q(L)}\left[\mathbb{E}_{Y\sim
p(Y|A=a,L)}\left[Y\right]\right]$.

Table 2: Network structures of DFPV for demand design experiments. For the
input layer, we provide the input variable. For the fully-connected layers
(FC), we provide the input and output dimensions.

Stage 1 Treatment Feature $\bm{\phi}_{\theta_{A(1)}}$
---
Layer | Configuration
1 | Input($P$)
2 | FC(1, 32), ReLU
3 | FC(32, 16), ReLU
4 | FC(16, 8), ReLU

Treatment-inducing Proxy Feature $\bm{\phi}_{\theta_{Z}}$
---
Layer | Configuration
1 | Input($C_{1}$, $C_{2}$)
2 | FC(2, 32), ReLU
3 | FC(32, 16), ReLU
4 | FC(16, 8), ReLU

Stage 2 Treatment Feature $\bm{\psi}_{\theta_{A(2)}}$
---
Layer | Configuration
1 | Input($P$)
2 | FC(1, 32), ReLU
3 | FC(32, 16), ReLU
4 | FC(16, 8), ReLU

Outcome-inducing Proxy Feature $\bm{\psi}_{\theta_{W}}$
---
Layer | Configuration
1 | Input($V$)
2 | FC(1, 32), ReLU
3 | FC(32, 16), ReLU
4 | FC(16, 8), ReLU

Table 3: Network structures of DFPV for dSprite experiments. For the input
layer, we provide the input variable. For the fully-connected layers (FC), we
provide the input and output dimensions. SN denotes Spectral Normalization
[23]. BN denotes Batch Normalization.

Stage 1 Treatment Feature $\bm{\phi}_{\theta_{A(1)}}$
---
Layer | Configuration
1 | Input($A$)
2 | FC(4096, 1024), SN, ReLU
3 | FC(1024, 512), SN, ReLU, BN
4 | FC(512, 128), SN, ReLU
5 | FC(128, 32), SN, ReLU

Treatment-inducing Proxy Feature $\bm{\phi}_{\theta_{Z}}$
---
Layer | Configuration
1 | Input($Z$)
2 | FC(3, 128), ReLU
3 | FC(128, 64), ReLU
4 | FC(64, 32), ReLU

Stage 2 Treatment Feature $\bm{\psi}_{\theta_{A(2)}}$
---
Layer | Configuration
1 | Input($A$)
2 | FC(4096, 1024), SN, ReLU
3 | FC(1024, 512), SN, ReLU, BN
4 | FC(512, 128), SN, ReLU
5 | FC(128, 32), SN, ReLU

Outcome-inducing Proxy Feature $\bm{\psi}_{\theta_{W}}$
---
Layer | Configuration
1 | Input($W$)
2 | FC(4096, 1024), SN, ReLU
3 | FC(1024, 512), SN, ReLU, BN
4 | FC(512, 128), SN, ReLU
5 | FC(128, 32), SN, ReLU

Table 4: Network structures of CEVAE for demand design experiment. For the
input layer, we provide the input variable. For the fully-connected layers
(FC), we provide the input and output dimensions.

Structure of $\bm{\psi}_{q}$
---
Layer | Configuration
1 | Input($P,Y,C_{1},C_{2},V$)
2 | FC(5, 128), ReLU
3 | FC(128, 64), ReLU
4 | FC(64, 32), ReLU

Structure of $\bm{\psi}_{p(W,Z|L)}$
---
Layer | Configuration
1 | Input($L$)
2 | FC(20, 64), ReLU
3 | FC(64, 32), ReLU
4 | FC(32, 16), ReLU

Structure of $\bm{\psi}_{p(A|L)}$
---
Layer | Configuration
1 | Input($L$)
2 | FC(20, 64), ReLU
3 | FC(64, 32), ReLU
4 | FC(32, 16), ReLU

Structure of $\bm{\mu}_{p(Y|A,L)}$
---
Layer | Configuration
1 | Input($L,P$)
2 | FC(21, 64), ReLU
3 | FC(64, 32), ReLU
4 | FC(32, 16), ReLU
5 | FC(16, 1)

Table 5: Network structures of CEVAE for dSprite experiment. For the input
layer, we provide the input variable. For the fully-connected layers (FC), we
provide the input and output dimensions. SN denotes Spectral Normalization
[23]. BN denotes Batch Normalization.

Structure of $\bm{\psi}_{q}$
---
Layer | Configuration
1 | Input($W,A,Z,Y$)
2 | FC(8196, 1024), SN, ReLU
3 | FC(1024, 512), SN, ReLU, BN
4 | FC(512, 128), SN, ReLU
5 | FC(128, 32), SN, ReLU

Structure of $\bm{\psi}_{p(W,Z|L)}$
---
Layer | Configuration
1 | Input($L$)
2 | FC(20, 64), ReLU
3 | FC(64, 128), ReLU
4 | FC(128, 256), ReLU

Structure of $\bm{\psi}_{p(A|L)}$
---
Layer | Configuration
1 | Input($L$)
2 | FC(20, 64), ReLU
3 | FC(64, 128), ReLU
4 | FC(128, 256), ReLU

Structure of $\bm{\mu}_{p(Y|A,L)}$
---
Layer | Configuration
1 | Input($L,A$)
2 | FC(4116, 1024), SN, ReLU, BN
3 | FC(1024, 512), SN, ReLU, BN
4 | FC(512, 128), SN, ReLU
5 | FC(128, 32), SN, ReLU
6 | FC(32, 1)