# Nonlinear Causal Discovery via Kernel Anchor Regression

Wenqi Shi1 Wenkai Xu2

(1Department of Industrial Engineering, Tsinghua University
2Department of Statistics, University of Oxford
)

###### Abstract

Learning causal relationships is a fundamental problem in science. Anchor
regression has been developed to address this problem for a large class of
causal graphical models, though the relationships between the variables are
assumed to be linear. In this work, we tackle the nonlinear setting by
proposing kernel anchor regression (KAR). Beyond the natural formulation using
a classic two-stage least square estimator, we also study an improved variant
that involves nonparametric regression in three separate stages. We provide
convergence results for the proposed KAR estimators and the identifiability
conditions for KAR to learn the nonlinear structural equation models (SEM).
Experimental results demonstrate the superior performances of the proposed KAR
estimators over existing baselines.

##  1  Introduction

Causal relationships are concerned with consequences of actions or decisions;
thus, understanding these relationships can be the key ingredient in many
scientific studies. For instance, medical practitioners need to know whether a
treatment is effective to the target disease in clinical trials;
econometricians ask whether a particular purchasing behaviour drives a change
in Consumer Price Index (CPI); epidemiologists want to understand whether a
government intervention policy has a positive effect on the pandemic. While
the goal of revealing causal effects remains the same, the focus in causal
relationships can differ by applications. To describe different aspects of the
causal notion and design statistical procedures for inferring causal effects,
various frameworks have been developed including Rubin’s potential outcome
framework [Rubin,, 2004, 2005], counterfactual distributions [Chernozhukov et
al.,, 2013] and Pearl’s causal graphical models [Pearl et al.,, 2000, 2016]. A
succinct yet comprehensive introduction can be found in Peters et al., [2017].

Causality has also been an evolving field in machine learning community and
machine learning techniques have been considered to improve the statistical
procedures for causal discovery. In particular, nonparmetric independence
[Gretton et al.,, 2005] and conditional independence [Fukumizu et al.,, 2007]
measures have been exploited to infer causal graphical models [Colombo et
al.,, 2012; Mooij et al.,, 2009], especially with additive noise [Hoyer et
al.,, 2008; Peters et al.,, 2014]. Independent Component Analysis (ICA)
methods [Hyvärinen,, 2013; Hyvarinen and Morioka,, 2017] have been employed to
identify causal relationship in both linear [Hyvärinen et al.,, 2010; Shimizu
et al.,, 2006, 2011] and non-linear settings [Monti et al.,, 2020; Khemakhem
et al.,, 2021]. Score matching [Hyvärinen and Dayan,, 2005] has also been
considered [Rolland et al.,, 2022] for non-linear causal discovery. Moreover,
kernel methods, that utilize rich representation of reproducing kernel Hilbert
space (RKHS), have been applied to tackle nonparametric estimation [Muandet et
al.,, 2021; Singh et al.,, 2019] and regression [Singh et al.,, 2019; Zhu et
al.,, 2022] problems with causal implications. Deep neural networks have also
been attempted for learning treatment effect [Johansson et al.,, 2020;
Kallus,, 2020; Louizos et al.,, 2017] or useful causal representations
[Besserve et al.,, 2019; Schölkopf et al.,, 2021; Xu et al.,, 2020, 2021].

Recently, an elegant and statistically robust approach formulates causality as
an invariant risk minimization (IRM), see for example [Bühlmann,, 2018; Peters
et al.,, 2016]. The causal structure is thought to be invariant across the
environment and robust under intervention. The IRM learning procedure
[Arjovsky et al.,, 2019] on the observational data is then formulated as a
regularised empirical risk minimization (ERM) to achieve both in-distribution
performance and out-of-distribution generalization. In particular, anchor
regression [Rothenhäusler et al.,, 2018] has been developed under the IRM
framework to tackle a very general class of causal graphical models with the
confounders being partly (but not fully) observed. By choosing different
regularisation parameter, anchor regression is able to unify the ordinary
least square (OLS) regression, partialling out (PA) regression, and
instrumental variable (IV) regression. While existing works [Oberst et al.,,
2021; Rothenhäusler et al.,, 2018] mostly considered linear cases, we explore
the non-linear setting for anchor regression [Kook et al.,, 2022].
Specifically, we consider the nonparametric estimation to tackle non-linear
features via RKHS functions.

The paper is structured as follows. In Section 2, we review useful concepts
including instrumental variable (IV), anchor regression (AR), and reproducing
kernel Hilbert space (RKHS). Then we develop two versions of kernel anchor
regression (KAR) estimators in Section 3. Theoretical analysis on the
estimators and the causal interpretation with nonlinear SEM are provided in
Section 4. Experimental results for synthetic data and real-world applications
are shown in Section 5 followed by concluding discussion and future directions
in Section 6. The code for the experiments is available in at
https://github.com/Swq118/Kernel-Anchor-Regression.

##  2 Background

Directed Acyclic Graph (DAG) is a power class of graphical model for
characterising conditional dependency structures and has been widely used for
probabilistic modelling such as hidden Markov models [Rabiner and Juang,,
1986], latent variable models [Bishop,, 1998] and topic models [Blei,, 2012].
By enforcing certain Markov and faithfulness assumptions [Peters et al.,,
2011], as well as noise structures [Hoyer et al.,, 2008], DAG models the
causal relationships [Glymour et al.,, 2019; Spirtes et al.,, 2013] and the
learning procedures have been developed [Colombo et al.,, 2012; Spirtes et
al.,, 2000; Zhang et al.,, 2018].

### From Instrumental Variable to Anchor Regression

Instrumental variable (IV) has been developed to incorporate endogenous
explanatory variables in econometrics [Bowden and Turkington,, 1990] and then
applied for estimating causal effect [Angrist et al.,, 1996]. Consider the
linear regression problem $Y=X\beta+\epsilon$. OLS assumes independence
between noise $\epsilon$ and explanatory $X$ (the exogenous variable) and
$\beta$ is estimated via minimizing

$\beta^{OLS}=\operatorname*{arg\,min}_{\beta}\mathbb{E}_{train}[\|Y-X\beta\|^{2}].$
(1)

The IV setting assumes explicit dependency between $X$ and $\epsilon$ via
instrumental variable $Z$, i.e. $X=Z\theta+\varepsilon$ where
$Z\perp\varepsilon$. The two-stage least squares (2SLS) procedure, widely used
in economics, tackles the linear IV estimation by first regressing $Z$ over
$X$ to get conditional means $\bar{X}(z):=\mathbb{E}[X|Z=z]$ and secondly
regressing outputs $Y$ on these conditional means111For the second stage, by
writing
$Y=X\beta+\epsilon=\underset{\mathbb{E}[{X}|Z]}{\underbrace{Z\theta}}\beta+(\varepsilon\beta+\epsilon),$
then the regressor is independent of noise and the OLS estimator can then
apply.. This corresponds to minimizing the projected least square objective,

$\beta^{IV}=\operatorname*{arg\,min}_{\beta}\mathbb{E}_{train}[\|P_{Z}(Y-X\beta)\|^{2}].$
(2)

$P_{Z}$ denotes the projection to $Z$ where
$P_{Z=z}(X)=\mathbb{E}[X|Z=z]=\bar{X}(z)$. 2SLS works well when the underlying
assumptions hold. The corresponding DAG is shown in Figure 1 with only solid
lines. In practice, the relation between $Y$ and $X$ may not be linear, nor
may be the relation between $X$ and $Z$. Nonlinear IV has been explored
[Bennett et al.,, 2019; Centorrino et al.,, 2019; Hartford et al.,, 2017;
Singh et al.,, 2019; Xu et al.,, 2020; Zhu et al.,, 2022].

$Z$$X$$Y$$C$ Figure 1: IV regression (solid lines only) and anchor regression
(with dashed lines).

However, $Y$’s dependency on $Z$ may not be solely through $X$, i.e. as the
dashed lines from $Z$ to $Y$ in Figure 1 indicate, $Y$ may depend on $Z$
directly, even though the strength of such dependency may remain unknown. The
latent confounder $C$ may not be independent of $Z$, as indicated by dashed
line from $Z$ to $C$ in Figure 1. Incorporating such dependency structures
tackles a much more general class of DAG, e.g. IV is a special case. To
estimate $\beta$, anchor regression has been proposed [Rothenhäusler et al.,,
2018] that effectively combines Equation 1 and Equation 2. For chosen
regularization parameter $\gamma$ and identity operator $Id(Z):=Z$,

$\displaystyle\beta^{\gamma}=\operatorname*{arg\,min}_{\beta}\mathbb{E}_{train}[\|(Id-
P_{Z})(Y-X\beta)\|^{2}]$ (3)
$\displaystyle+\gamma\mathbb{E}_{train}[\|P_{Z}(Y-X\beta)\|^{2}].$ (4)

Here, $\gamma\geq 0$ can be thought of the level of dependencies of $Y$ from
$Z$ variable222The smaller $\gamma$ value dashed line, the stronger the
dependency, i.e. the more solid dashed line from $Z$ to $Y$.. By setting
different $\gamma$ values, anchor regression resembles classical settings,
i.e. $\gamma=1$ corresponds to OLS, $\beta^{1}=\beta^{OLS}$; $\gamma\to\infty$
corresponds to IV,
$\beta^{\to\infty}:=\lim_{\gamma\to\infty}\beta^{\gamma}=\beta^{IV}$;
$\gamma=0$ corresponds to the "partialling out" setting where only residuals
between regression of $Z$ to $X$ and $Y$ are of interest.

### Kernel-based Methods

Kernel methods employ functions in RKHS to tackle various statistical and
machine learning tasks with nonlinear features [Hofmann et al.,, 2008], e.g.
kernel ridge regression, support vector machine[Scholkopf and Smola,, 2018;
Steinwart and Christmann,, 2008], etc. Functions in RKHS have also been
developed to represent and characterize distributions, via kernel mean
embedding [Muandet et al.,, 2017]. For probability measure $p$, and kernel $k$
associated with RKHS ${\mathcal{H}}$, the mean embedding $\mu_{p}:=\int
k(x,\cdot)dp(x)\in{\mathcal{H}}$. This notion has been widely used to compare
distributions, e.g. via maximum-mean-discrepancy (MMD) [Gretton et al.,,
2012]. With $p$ being a conditional distribution, conditional mean embedding
[Song et al.,, 2009] has also been considered for learning and regression
problems. Various techniques have also been developed to formulate and learn
operators to manipulate conditional mean embeddings [Fukumizu et al.,, 2007;
Grünewälder et al.,, 2012]. With the rich representation of nonlinear
features, RKHS functions are also applicable of learning distribution directly
via distribution regression [Szabó et al.,, 2015, 2016].

##  3 Kernel Anchor Regression

To capture the non-linear features in the DAG, we kernelize the anchor
regression framework by utilizing the rich feature representation of RKHS
functions. The kernelizing procedure is inspired from kernel instrumental
variable (KIV) [Singh et al.,, 2019] where the operators are learned for
conditional mean embedding in two separate regression stages. The DAG
representation is illustrated in Figure 2.

Let
$k_{\mathcal{X}}:{\mathcal{X}}\times{\mathcal{X}}\rightarrow{{\mathbb{R}}}$,
$k_{\mathcal{Z}}:{\mathcal{Z}}\times{\mathcal{Z}}\rightarrow{{\mathbb{R}}}$ be
measurable positive definite kernels corresponding to RKHS
${\mathcal{H}}_{\mathcal{X}}$ and ${\mathcal{H}}_{\mathcal{Z}}$. Denote the
feature maps
$\psi:{\mathcal{X}}\rightarrow{\mathcal{H}}_{\mathcal{X}},x\rightarrow
k_{\mathcal{X}}(x,\cdot)$ and
$\phi:{\mathcal{Z}}\rightarrow{\mathcal{H}}_{\mathcal{Z}},z\rightarrow
k_{\mathcal{Z}}(z,\cdot).$ Let $P_{{\phi(Z)}}$ and $Id$ denote the
$L_{2}$-projection on the linear span from the components of $\phi(Z)$ and the
identity operator, respectively. Denote
$H:{\mathcal{H}}_{\mathcal{X}}\to{\mathcal{Y}}$ as the conditional operator we
aim to learn. Then for $\gamma\geq 0$, define the population-level kernel
anchor regression operator $H^{\gamma}$ as

$\displaystyle{H}^{\gamma}=\operatorname*{arg\,min}_{H}{\mathbb{E}}[\|(Id-
P_{{\phi(Z)}})(Y-H\psi(X))\|^{2}]$
$\displaystyle+\gamma{\mathbb{E}}[\|P_{{\phi(Z)}}(Y-H\psi(X))\|^{2}].$ (5)

To unravel $P_{\phi(Z)}$, both IV and AR estimators applied the two-stage
procedure, where the first stage is to estimate the projection operator
$P_{\phi(Z)}$ and the second stage is to perform the projection adjusted
regression.

$Z$$\mathcal{H}_{\mathcal{Z}}$$X$$\mathcal{H}_{\mathcal{X}}$$C$$Y$$\phi$$\psi$
Figure 2:  DAG representation for kernel anchor regression.

###  3.1 Projection Stage

The projection stage aims to tackle $P_{\phi(Z)}$ by transforming the problem
of learning $P_{\phi(Z)}\psi(X)$ and $P_{\phi(Z)}Y$ into two separate kernel
ridge regressions. Let operators
$E_{X}:{\mathcal{H}}_{\mathcal{Z}}\to{\mathcal{H}}_{\mathcal{X}}$ and
$E_{Y}:{\mathcal{H}}_{\mathcal{Z}}\to{\mathcal{Y}}$ be the projections to
learn333We note that due to the explicit dependency from $Z$ to $Y$,
$P_{\phi(Z)}Y$ needs to be treated separately from $P_{\phi(Z)}\psi(X)$. This
is different from the IV setting where $P_{\phi(Z)}Y=Y$ as the edge from $Z$
to $Y$ in the DAG is absent.; $\alpha_{1},\alpha_{2}>0$ be regularization
parameters. The objectives regularized by Hilbert-Schmidt (HS) norm are

$\displaystyle{\mathcal{E}}_{\alpha_{1}}(E_{X})={\mathbb{E}}\|\psi(X)-E_{X}\phi(Z)\|_{{\mathcal{H}}_{\mathcal{X}}}^{2}+\alpha_{1}\|E_{X}\|_{HS}^{2},$
(6)
$\displaystyle{\mathcal{E}}_{\alpha_{2}}(E_{Y})={\mathbb{E}}\|Y-E_{Y}\phi(Z)\|_{{{\mathcal{Y}}}}^{2}+\alpha_{2}\|E_{Y}\|_{HS}^{2}.$
(7)

Denote the optimal operators for the population risks as
$E^{p}_{\alpha_{1},X}=\operatorname*{arg\,min}_{E_{X}}{\mathcal{E}}_{\alpha_{1}}(E_{X})$,
and
$E^{p}_{\alpha_{2},Y}=\operatorname*{arg\,min}_{E_{Y}}{\mathcal{E}}_{\alpha_{2}}(E_{Y})$.
We then consider two variants of empirical risks and their corresponding
estimations.

####  3.1.1 Disjoint sample sets projection

Firstly, we treat two ridge regression in Equation 6 and Equation 7
independently, by using two disjoint sets of samples
$\mathbb{S}_{1}=\\{(x_{i},z_{i})\\}_{i\in[n_{1}]}$ and
$\mathbb{S}_{2}=\\{(y_{j},z_{j})\\}_{j\in[n_{2}]}$. The empirical forms for
Equation 6 and Equation 7 are

$\displaystyle\frac{1}{n_{1}}\sum_{i\in[n_{1}]}\|\psi(x_{i})-E_{X}\phi(z_{i})\|_{{\mathcal{H}}_{\mathcal{X}}}^{2}+\alpha_{1}\|E_{X}\|_{HS}^{2},$
(8)
$\displaystyle\frac{1}{n_{2}}\sum_{j\in[n_{2}]}\|y_{j}-E_{Y}\phi(z_{j})\|_{{\mathcal{Y}}}^{2}+\alpha_{2}\|E_{Y}\|_{HS}^{2}$
. (9)

Denote by $\Phi_{1,Z}=(\phi(z_{1}),\dots,\phi(z_{n_{1}}))$,
$\\{z_{i}\\}_{i\in[n_{1}]}\subset\mathbb{S}_{1}$;
$\Phi_{2,Z}=(\phi(z_{1}),\dots,\phi(z_{n_{2}}))$,
$\\{z_{j}\\}_{j\in[n_{2}]}\subset\mathbb{S}_{2}$; their corresponding gram
matrices $K_{1,ZZ}=\Phi_{1,Z}^{\top}\Phi_{1,Z}\in{{\mathbb{R}}}^{n_{1}\times
n_{1}}$ and
$K_{2,ZZ}=\Phi_{2,Z}^{\top}\Phi_{2,Z}\in{{\mathbb{R}}}^{n_{2}\times n_{2}}$.
Denote $\Psi_{1,X}=(\psi(x_{1}),\dots,\psi(x_{n_{1}}))$,
$\\{x_{i}\\}_{i\in[n_{1}]}\subset\mathbb{S}_{1}$ and
$Y_{2}=(y_{1},\dots,y_{n_{2}})$,
$\\{y_{j}\\}_{j\in[n_{2}]}\subset\mathbb{S}_{2}$. By the standard regression
formula, the optimal operators to minimize Equation 8 and Equation 9 are

$\displaystyle E_{\alpha_{1},X}^{n_{1}}$
$\displaystyle=\Psi_{1,X}(K_{1,ZZ}+n_{1}\alpha_{1}I)^{-1}\Phi_{1,Z}^{\top},$
(10) $\displaystyle E_{\alpha_{2},Y}^{n_{2}}$
$\displaystyle=Y_{2}(K_{2,ZZ}+n_{2}\alpha_{2}I)^{-1}\Phi_{2,Z}^{\top},$ (11)

where the superscripts $n_{1},n_{2}$ explicitly reveal sample sizes. We note
that the projections $P_{\phi(Z)}$ are estimated differently for
$P_{\phi(Z)}\psi(X)$ and $P_{\phi(Z)}Y$, through
$(K_{1,ZZ}+n_{1}\alpha_{1}I)^{-1}$ and $(K_{2,ZZ}+n_{2}\alpha_{2}I)^{-1}$,
respectively. $K_{1,ZZ}$ and $K_{2,ZZ}$ are independent due to the use of
disjoint i.i.d. sample sets of $Z$.

####  3.1.2 Joint sample set projection

On the other hand, we can also consider the projection analogous to
Rothenhäusler et al., [2018] where we jointly consider the samples used for
both projections, i.e. projecting onto the same $\phi(Z)$ subspace. Setting
$n=n_{1}+n_{2}$ and $\alpha_{1}=\alpha_{2}=\alpha$, we consider the joint
sample set $\mathbb{S}=\\{(x_{i},y_{i},z_{i})\\}_{i\in[n]}$ and the empirical
risks

$\displaystyle\frac{1}{n}\sum_{i\in[n]}\|\psi(x_{i})-E_{X}\phi(z_{i})\|_{{\mathcal{H}}_{\mathcal{X}}}^{2}+\alpha\|E_{X}\|_{HS}^{2},$
(12)
$\displaystyle\frac{1}{n}\sum_{i\in[n]}\frac{1}{n}\|y_{i}-E_{Y}\phi(z_{i})\|_{{\mathcal{Y}}}^{2}+\alpha\|E_{Y}\|_{HS}^{2}.$
(13)

Denote $K_{ZZ}\in{{\mathbb{R}}}^{n\times n}$ as the gram matrix from
$\\{z_{i}\\}_{i\in[n]}\subset\mathbb{S}$; $\Phi_{Z}=(\phi(z_{1}),\dots$,
$\phi(z_{n})),\\{z_{i}\\}_{i\in[n]}\subset\mathbb{S}$;
$\Psi_{X}=(\psi(x_{1}),\dots,\psi(x_{n})),\\{x_{i}\\}_{i\in[n]}\subset\mathbb{S}$
and $Y=(y_{1},\dots,y_{n}),y_{i}\in\mathbb{S}$. Then we have

$\displaystyle E_{\alpha,X}^{n}$ $\displaystyle=\Psi_{X}(K_{ZZ}+n\alpha
I)^{-1}\Phi_{Z}^{\top},$ (14) $\displaystyle E_{\alpha,Y}^{n}$
$\displaystyle=Y^{\top}(K_{ZZ}+n\alpha I)^{-1}\Phi_{Z}^{\top}.$ (15)

By setting the same level of regularisation, we can see that the $P_{\phi(Z)}$
projection, through $(K_{ZZ}+n\alpha I)^{-1}$, are the same for
$P_{\phi(Z)}\psi(X)$ and $P_{\phi(Z)}Y$.

###  3.2 Regression Stage

With the learned projections $P_{\phi(Z)}\psi(X)$ and $P_{\phi(Z)}Y$, we can
now tackle the overall objective in Equation 5.

Denote ${\mathcal{E}}(E_{X})$ and ${\mathcal{E}}(E_{Y})$ as the unregularized
version of Equation 6 and Equation 7; $E^{p}_{X}$ and $E^{p}_{Y}$ their
corresponding optimal operators, respectively. For given $\gamma$, define the
transformed input and output as

$\psi_{\gamma}(X)=\psi(X)-E^{p}_{X}\phi(Z)+\sqrt{\gamma}E^{p}_{X}\phi(Z)\in{\mathcal{H}}_{\mathcal{X}},$
(16)
$Y_{\gamma}=Y-E^{p}_{Y}\phi(Z)+\sqrt{\gamma}E^{p}_{Y}\phi(Z)\in{\mathcal{Y}}.$
(17)

######  Proposition 1 (Equivalence).

Let $H:{\mathcal{H}}_{\mathcal{X}}\to{\mathcal{Y}}$, and consider the
regression of transformed output in Equation 17 on transformed input in
Equation 16

${\mathcal{E}}^{\gamma}(H)={\mathbb{E}}_{(Z,X,Y)}\|Y_{\gamma}-H\psi_{\gamma}(X)\|_{\mathcal{Y}}^{2}.$
(18)

The solution to Equation 18 is equivalent to the KAR estimator in Equation 5,
i.e. $H^{\gamma}=\operatorname*{arg\,min}_{H}{\mathcal{E}}^{\gamma}(H).$

The proof is by expanding the projection operator $E^{p}_{X}$ and $E^{p}_{Y}$,
which is similar to the linear case in Rothenhäusler et al., [2018].

With regularization parameter $\xi\geq 0$, Equation 18 has the kernel ridge
regression form defined as

${\mathcal{E}}_{\xi}^{\gamma}(H)={\mathbb{E}}_{(Z,X,Y)}\|Y_{\gamma}-H\psi_{\gamma}(X)\|_{\mathcal{Y}}^{2}+\xi\|H\|_{HS}^{2}.$
(19)

The regression stage is formulated regardless how the projections are
estimated in Section 3.1. For the empirical version, we consider the estimated
operators $\widehat{E}_{X}\in\\{E_{\alpha_{1},X}^{n_{1}},E_{\alpha,X}^{n}\\}$
and $\widehat{E}_{Y}\in\\{E_{\alpha_{2},Y}^{n_{2}},E_{\alpha,Y}^{n}\\}$. With
samples $\mathbb{S}^{m}=\\{(x_{l},y_{l},z_{l})\\}_{l\in[m]}$, we compute the
transformed inputs and outputs as

$\widehat{\psi}_{\gamma,l}(x)=\psi(x_{l})+(\sqrt{\gamma}-1)\widehat{E}_{X}\phi(z_{l})\in{\mathcal{H}}_{\mathcal{X}},$
$\widehat{y}_{\gamma,l}=y_{l}+(\sqrt{\gamma}-1)\widehat{E}_{Y}\phi(z_{l})\in{\mathcal{Y}}.$

The empirical risk has the form

$\widehat{\mathcal{E}}_{\xi}^{\gamma,m}(H)=\frac{1}{m}\sum_{l\in[m]}\|\widehat{y}_{\gamma,l}-H\widehat{\psi}_{\gamma,l}(x)\|_{\mathcal{Y}}^{2}+\xi\|H\|_{HS}^{2},$
(20)
$\widehat{H}_{\xi}^{\gamma,m}=\operatorname*{arg\,min}\widehat{\mathcal{E}}_{\xi}^{\gamma,m}(H).$

Denote
$\widehat{\Psi}_{\gamma}=(\widehat{\psi}_{\gamma,1}(x),\dots,\widehat{\psi}_{\gamma,m}(x))$
and its gram matrix
$K_{\widehat{\Psi}_{\gamma}\widehat{\Psi}_{\gamma}}=\widehat{\Psi}_{\gamma}^{\top}\widehat{\Psi}_{\gamma}\in{{\mathbb{R}}}^{m\times
m}$;
$\widehat{Y}_{\gamma}=(\widehat{y}_{\gamma,1},\dots,\widehat{y}_{\gamma,m})$.
Again, by standard regression formula,

$\widehat{H}_{\xi}^{\gamma,m}=\widehat{Y}_{\gamma}(K_{\widehat{\Psi}_{\gamma}\widehat{\Psi}_{\gamma}}+m\xi
I)^{-1}\widehat{\Psi}_{\gamma}^{\top}.$ (21)

###  3.3 KAR Estimator

Given observational data of size $N$, $\\{(x_{i},y_{i},z_{i})\\}_{i\in[N]}$,
the KAR procedure can be performed in two ways based on the two variants in
the projection stage.

##### Three-stage KAR

To apply the disjoint sample sets projection in Section 3.1.1, we _randomly_
split the data set of size $N$ into three disjoint sets of sample size
$n_{1},n_{2},m$ where $N=n_{1}+n_{2}+m$ and re-index them from $\\{1:N\\}$.
The first two sets of data $\\{(x_{i},z_{i})\\}_{i\in\\{1:n_{1}\\}}$ and
$\\{(y_{j},z_{j})\\}_{j\in\\{n_{1}+1:n_{1}+n_{2}\\}}$ are used for learning
the projection operators in Equation 10 and Equation 11. We note that samples
$\\{y_{i}\\}_{i\in\\{1:n_{1}\\}}$ and
$\\{x_{j}\\}_{j\in\\{n_{1}+1:n_{1}+n_{2}\\}}$ are not used. The third set
$\\{(x_{l},y_{l},z_{l})\\}_{l\in\\{n_{1}+n_{2}+1:N\\}}$ is used for regression
stage to learn $\widehat{H}_{\xi}^{\gamma,m}$ in Equation 21. This procedure,
termed KAR, includes solving three different regression problems, which is
different from the two-stage settings used in linear anchor regression
[Rothenhäusler et al.,, 2018].

##### Two-stage KAR

For the joint sample set projection in Section 3.1.2, we only split the data
of size $N$ into two disjoint sets _randomly_ of size $n$ and $m$ where
$N=n+m$ and re-index them such that
$\\{(x_{i},y_{i},z_{i})\\}_{i\in\\{1:n\\}}$ and
$\\{(x_{l},y_{l},z_{l})\\}_{l\in\\{n+1:N\\}}$. The first set is then grouped
into $\\{(x_{i},z_{i})\\}_{i\in\\{1:n\\}}$ and
$\\{(y_{i},z_{i})\\}_{j\in\\{1:n\\}}$ to learn the projection operators in
Equation 10 and Equation 11. In this manner, $\\{z_{i}\\}_{i\in\\{1:n\\}}$ are
used twice. The second set $\\{(x_{l},y_{l},z_{l})\\}_{l\in\\{n+1:N\\}}$ is
used for regression stage to learn $\widehat{H}_{\xi}^{\gamma,m}$ in Equation
21, which is the same as the three-stage procedure above. This procedure,
termed KAR.2, replicates the 2SLS used in KIV [Singh et al.,, 2019] and linear
anchor regression [Rothenhäusler et al.,, 2018].

##  4  Analysis of KAR Estimators

###  4.1 Consistency

In this section, we first focus on the three-stage KAR procedure with disjoint
sample sets for projection in Section 3.1.1. The closed form solutions and
convergence rates of the estimators are extended from the analysis of 2SLS in
KIV [Singh et al.,, 2019]. We follow the integral operator notations in Singh
et al., [2019]. Define the projection stage operators as

$\displaystyle S_{1}^{*}$ $\displaystyle:$
$\displaystyle{\mathcal{H}}_{\mathcal{Z}}\rightarrow
L^{2}({\mathcal{Z}},\rho_{\mathcal{Z}}),\quad l\rightarrow\langle
l,\phi(\cdot)\rangle_{{\mathcal{H}}_{\mathcal{Z}}},$ $\displaystyle S_{1}$
$\displaystyle:$ $\displaystyle
L^{2}({\mathcal{Z}},\rho_{\mathcal{Z}})\rightarrow{\mathcal{H}}_{\mathcal{Z}},\quad\tilde{l}\rightarrow\int\phi(z)\tilde{l}(z)d\rho_{\mathcal{Z}}(z),$

where $\rho$ denotes the joint distribution of $(Z,X,Y)$.
$L^{2}({\mathcal{Z}},\rho_{\mathcal{Z}})$ denotes the space of square
integrable functions from ${\mathcal{Z}}$ to ${\mathcal{Y}}$ with respect to
measure $\rho_{\mathcal{Z}}$, where $\rho_{\mathcal{Z}}$ is the restriction of
$\rho$ to ${\mathcal{Z}}$. $T_{1}=T_{2}=S_{1}\circ S_{1}^{*}$ are then
uncentered covariance operators. We define the power of operator $T_{1}$ with
respect to its eigendecomposition. Let
${\mathcal{H}}_{\Gamma}={\mathcal{H}}_{\mathcal{X}}\otimes{\mathcal{H}}_{\mathcal{Z}}$,
${\mathcal{H}}_{\Theta}={\mathcal{Y}}\otimes{\mathcal{H}}_{\mathcal{Z}}$ and
${\mathcal{H}}_{\Omega}={\mathcal{Y}}\otimes{\mathcal{H}}_{\mathcal{X}}$ be
the relevant tensor product spaces for the operators.

######  Condition 1.

(i) ${\mathcal{X}}$ and ${\mathcal{Z}}$ are Polish, i.e. separable and
completely metrizable topological spaces. (ii) $k_{\mathcal{X}}$ and
$k_{\mathcal{Z}}$ are continuous and bounded:
$\sup_{x\in{\mathcal{X}}}\|\psi(x)\|_{{\mathcal{H}}_{\mathcal{X}}}\leq Q_{1}$,
$\sup_{z\in{\mathcal{Z}}}\|\phi(z)\|_{{\mathcal{H}}_{\mathcal{Z}}}\leq\kappa$.
(iii) $\psi$ and $\phi$ are measurable. (iv) $k_{\mathcal{X}}$ is
characteristic. (v) $E_{X}^{p}\in{\mathcal{H}}_{\Gamma}$ s.t.
${\mathcal{E}}(E_{X}^{p})=\inf_{E_{X}\in{\mathcal{H}}_{\Gamma}}{\mathcal{E}}(E_{X})$.
(vi) Fix $\zeta_{1}<\infty$. For $c_{1}\in(1,2]$, define the prior
$\mathcal{P}(\zeta_{1},c_{1})$ as the set of probability distributions $\rho$
on ${\mathcal{X}}\times{\mathcal{Z}}$ s.t. $\exists
G_{1}\in{\mathcal{H}}_{\Gamma}$ s.t. $E_{X}^{p}=T_{1}^{(c_{1}-1)/2}\circ
G_{1}$ and $\|G_{1}\|_{{\mathcal{H}}_{\Gamma}}^{2}\leq\zeta_{1}$.

Condition 1 is adapted from Singh et al., [2019] to bound the approximation
error of the regularized estimator $E_{\alpha_{1},X}^{n_{1}}$. Parameter
$c_{1}$ suggests the smoothness of conditional operator
$E_{\alpha_{1},X}^{n_{1}}$. A larger $c_{1}$ corresponds to a smoother
operator.

######  Lemma 1.

$\forall\alpha_{1}>0$, the solution $E_{\alpha_{1},X}^{n_{1}}$ of the
regularized empirical objective in Equation 8 exists and is unique. With
$\mathbf{T}_{1}=\frac{1}{n_{1}}\sum_{i=1}^{n_{1}}\phi(z_{i})\otimes\phi(z_{i})$
and
$\mathbf{g}_{1}=\frac{1}{n_{1}}\sum_{i=1}^{n_{1}}\phi(z_{i})\otimes\psi(x_{i})$,
the estimator in Equation 10 has the form
$E_{\alpha_{1},X}^{n_{1}}=(\mathbf{T}_{1}+\alpha_{1})^{-1}\circ\mathbf{g}_{1}.$

Under Condition 1 and $\alpha_{1}=n_{1}^{-1/(c_{1}+1)}$, we have:

$\displaystyle\|E_{\alpha_{1},X}^{n_{1}}-E_{X}^{p}\|_{{\mathcal{H}}_{\Gamma}}=O_{p}(n_{1}^{-\frac{c_{1}-1}{2(c_{1}+1)}}).$

Lemma 1 follows from Singh et al., [2019], and shows that the efficient rate
of $\alpha_{1}$ is $n_{1}^{-1/(1+c_{1})}$. Note that the convergence rate of
$E_{\alpha_{1},X}^{n_{1}}$ is calibrated by $c_{1}$, which measures the
smoothness of the conditional expectation operator $E_{X}$.

For the disjoint set projection in Section 3.1.1, the closed form solution and
convergence rate for learning $P_{{\phi(Z)}}Y$ estimator is similar to that of
learning $P_{\phi(Z)}\psi(X)$ due to the independent estimation procedure and
further requires the following conditions.

######  Condition 2.

(i) ${\mathcal{Y}}$ is a Polish space. (ii) $Y$ is bounded:
$\sup_{y\in{\mathcal{Y}}}\|y\|_{{\mathcal{Y}}}\leq Q_{2}$. (iii)
$E_{Y}^{p}\in{\mathcal{H}}_{\Theta}$ s.t.
${\mathcal{E}}(E_{Y}^{p})=\inf_{E_{Y}\in{\mathcal{H}}_{\Theta}}{\mathcal{E}}(E_{Y})$.
(iv) Fix $\zeta_{2}<\infty$. For $c_{2}\in(1,2]$, define the prior
$\mathcal{P}(\zeta_{2},c_{2})$ as the set of probability distributions $\rho$
on ${\mathcal{Y}}\times{\mathcal{Z}}$ s.t. $\exists
G_{2}\in{\mathcal{H}}_{\Theta}$ s.t. $E_{Y}^{p}=T_{2}^{(c_{2}-1)/{2}}\circ
G_{2}$ and $\|G_{2}\|_{{\mathcal{H}}_{\Theta}}^{2}\leq\zeta_{2}$.

######  Lemma 2.

$\forall\alpha_{2}>0$, the solution $E_{\alpha_{2},Y}^{n_{2}}$ of the
regularized empirical objective in Equation 9 exists and is unique. With
$\mathbf{T}_{2}=\frac{1}{n_{2}}\sum_{j=1}^{n_{2}}\phi(z_{j})\otimes\phi(z_{j})$
and $\mathbf{g}_{2}=\frac{1}{n_{2}}\sum_{j=1}^{n_{2}}\phi(z_{j})y_{j}$, the
estimator in Equation 10 has the form
$E_{\alpha_{2},Y}^{n_{2}}=(\mathbf{T}_{2}+\alpha_{2})^{-1}\circ\mathbf{g}_{2}.$
Under Condition 1– 2 and $\alpha_{2}=n_{2}^{-1/(c_{2}+1)}$, we have:

$\displaystyle\|E_{\alpha_{2},Y}^{n_{2}}-E_{Y}^{p}\|_{{\mathcal{H}}_{\Theta}}=O_{p}(n_{2}^{-\frac{c_{2}-1}{2(c_{2}+1)}}).$

Similar to learning projection $P_{\phi(Z)}\psi(X)$, the efficient rate of
$\alpha_{2}$ is $n_{2}^{-1/(1+c_{2})}$, where $c_{2}$ measures the smoothness
of the conditional expectation operator $E_{Y}$.

Let $L^{2}({\mathcal{H}}_{\mathcal{X}},\rho_{{\mathcal{H}}_{\mathcal{X}}})$
denote the space of square integrable functions from
${\mathcal{H}}_{\mathcal{X}}$ to ${\mathcal{Y}}$ with respect to measure
$\rho_{{\mathcal{H}}_{\mathcal{X}}}$, where
$\rho_{{\mathcal{H}}_{\mathcal{X}}}$ is the extension of $\rho$ to
${\mathcal{H}}_{\mathcal{X}}$. Define the regression stage operator as

$\displaystyle S^{*}$ $\displaystyle:$
$\displaystyle{\mathcal{H}}_{\Omega}\rightarrow
L^{2}({\mathcal{H}}_{\mathcal{X}},\rho_{{\mathcal{H}}_{\mathcal{X}}}),\quad
H\rightarrow\Omega^{*}_{(\cdot)}H,$ $\displaystyle S$ $\displaystyle:$
$\displaystyle
L^{2}({\mathcal{H}}_{\mathcal{X}},\rho_{{\mathcal{H}}_{\mathcal{X}}})\rightarrow{\mathcal{H}}_{\Omega},$
$\displaystyle\tilde{H}\rightarrow\int\Omega_{\psi_{\gamma}}\circ\tilde{H}\psi_{\gamma}d\rho_{{\mathcal{H}}_{\mathcal{X}}}(\psi_{\gamma}),$

where $\Omega_{\psi_{\gamma}}:{\mathcal{Y}}\rightarrow{\mathcal{H}}_{\Omega}$
defined by $y\rightarrow\Omega(\cdot,\psi_{\gamma})y$ is the point evaluator
of Micchelli and Pontil, [2005]. Define
$T_{\psi_{\gamma}}=\Omega_{\psi_{\gamma}}\circ\Omega^{*}_{\psi_{\gamma}}$ and
covariance operator $T=S\circ S^{*}$. Define the power of operator $T$ with
respect to its eigendecomposition. Condition 3 below extends hypothesis 7–9 in
Singh et al., [2019], and is sufficient to bound the excess error of
$\widehat{H}_{\xi}^{\gamma,m}$ with the error propagated from the estimators
in the projection stage.

######  Condition 3.

  * (i)

The $\\{\Omega_{\psi_{\gamma}}\\}$ operator family is uniformly bounded in
Hilbert-Schmidt norm: $\exists B$ s.t. $\forall\psi_{\gamma}$,
$\|\Omega_{\psi_{\gamma}}\|_{L_{2}({\mathcal{Y}},{\mathcal{H}}_{\Omega})}^{2}=Tr(\Omega^{*}_{\psi_{\gamma}}\circ\Omega_{\psi_{\gamma}})\leq
B$.

  * (ii)

The $\\{\Omega_{\psi_{\gamma}}\\}$ operator family is is Hölder continuous in
operator norm: $\exists L>0,\iota\in(0,1]$ s.t.
$\forall\psi_{\gamma},\psi_{\gamma}^{\prime},\|\Omega_{\psi_{\gamma}}-\Omega_{\psi_{\gamma}^{\prime}}\|_{L({\mathcal{Y}},{\mathcal{H}}_{\Omega})}\leq
L\|\psi_{\gamma}-\psi_{\gamma}^{\prime}\|_{{\mathcal{H}}_{\mathcal{X}}}^{\iota}$.

  * (iii)

$H^{\gamma}\in{\mathcal{H}}_{\Omega}$, then
${\mathcal{E}}^{\gamma}(H^{\gamma})=\inf_{H\in{\mathcal{H}}_{\Omega}}{\mathcal{E}}^{\gamma}(H)$.

  * (iv)

$Y_{\gamma}$ is bounded, i.e. $\exists C<\infty$ s.t.
$\|Y_{\gamma}\|_{\mathcal{Y}}\leq C$.

  * (v)

Fix $\zeta<\infty$. For given $b_{\gamma}\in(1,\infty]$ and
$c_{\gamma}\in(1,2]$, define the prior
$\mathcal{P}(\zeta,b_{\gamma},c_{\gamma})$ as the set of probability
distributions $\rho$ on ${\mathcal{H}}_{\mathcal{X}}\times{\mathcal{Y}}$ s.t.

    * (a)

range space assumption is satisfied: $\exists G\in{\mathcal{H}}_{\Omega}$ s.t.
$H^{\gamma}=T^{\frac{(c_{\gamma}-1)}{2}}\circ G$ and
$\|G\|_{{\mathcal{H}}_{\Omega}}^{2}\leq\zeta$;

    * (b)

the eigenvalues from spectral decomposition
$T=\sum_{k=1}^{\infty}\lambda_{k}e_{k}\langle\cdot,e_{k}\rangle_{{\mathcal{H}}_{\Omega}}$,
where $\\{e_{k}\\}_{k=1}^{\infty}$ is a basis of $Ker(T)^{\bot}$ , satisfy
$\alpha\leq k^{b_{\gamma}}\lambda_{k}\leq\beta$ for some $\alpha,\beta>0$.

We note that all parameters mentioned in Condition 3 depend on $\gamma$,
though the function representations are not explicit. We set subscript
$\gamma$ especially for parameters $b_{\gamma}$ and $c_{\gamma}$ to emphasize
their dependency on $\gamma$. Parameter $b_{\gamma}$ measures the decay of
eigenvalues of the covariance operator $T$ and specifically, larger
$b_{\gamma}$ suggests smaller effective input dimension. A larger $c_{\gamma}$
corresponds to a smoother operator $H^{\gamma}$.

######  Lemma 3.

$\forall\xi>0$, the solution $\widehat{H}_{\xi}^{\gamma,m}$ to
$\widehat{\mathcal{E}}_{\xi}^{\gamma,m}$ exist, and is unique for each
$\gamma$. Let
$\widehat{\mathbf{T}}=\frac{1}{m}\sum_{l=1}^{m}T_{\widehat{\psi}_{\gamma,l}}$,
$\widehat{\mathbf{g}}=\frac{1}{m}\sum_{l=1}^{m}\Omega_{\widehat{\psi}_{\gamma,l}}\widehat{y}_{\gamma,l}$.
Equation 21 has the form

$\widehat{H}_{\xi}^{\gamma,m}=(\widehat{\mathbf{T}}+\xi)^{-1}\circ\widehat{\mathbf{g}}.$

######  Condition 4.

For $c_{1},c_{2}$ set in Conditions 1– 2 and $\iota$ satisfying Condition 3,
assume $n_{2}\geq n_{1}^{\frac{\iota(c_{1}-1)(c_{2}+1)}{(c_{1}+1)(c_{2}-1)}}$.

######  Remark 1.

Condition 4 is sufficient but not necessary to ensure that the error
propagates to regression stage from estimating $E_{Y}^{p}$ is smaller than
that from estimating $E_{X}^{p}$ in disjoint sample sets projection.

The main challenge of extending the convergence rate of KIV estimator [Singh
et al.,, 2019] to KAR estimator is that in our case, the excess error depends
not only on the accuracy of $E_{X}^{p}$ estimator but also on the accuracy of
$E_{Y}^{p}$ estimator. However, by assuming Condition 4, we ensure the error
from estimating $E_{Y}^{p}$ is dominated by that of $E_{X}^{p}$, and manage to
illustrate the optimal convergence rate for KAR as shown below in Thereom 1.

######  Theorem 1.

Under Condition 1– 4, let $d_{1},d_{2}>0$ and choose
$\alpha_{1}=n_{1}^{-\frac{1}{c_{1}+1}}$,
$\alpha_{2}=n_{2}^{-\frac{1}{c_{2}+1}}$,
$n_{1}=m^{\frac{d_{1}(c_{1}+1)}{\iota(c_{1}-1)}}$,
$n_{2}=m^{\frac{d_{2}(c_{2}+1)}{\iota(c_{2}-1)}}$, We have:

  * (i)

If $d_{1}\leq\frac{b_{\gamma}(c_{\gamma}+1)}{b_{\gamma}c_{\gamma}+1}$, then
${\mathcal{E}}^{\gamma}(\widehat{H}_{\xi}^{\gamma,m})-{\mathcal{E}}^{\gamma}(H^{\gamma})=O_{p}(m^{-\frac{d_{1}c_{\gamma}}{c_{\gamma}+1}})$
with $\xi=m^{-\frac{d_{1}}{c_{\gamma}+1}}$.

  * (ii)

If $d_{1}>\frac{b_{\gamma}(c_{\gamma}+1)}{b_{\gamma}c_{\gamma}+1}$, then
${\mathcal{E}}^{\gamma}(\widehat{H}_{\xi}^{\gamma,m})-{\mathcal{E}}^{\gamma}(H^{\gamma})=O_{p}(m^{-\frac{b_{\gamma}c_{\gamma}}{b_{\gamma}c_{\gamma}+1}})$
with $\xi=m^{-\frac{b_{\gamma}}{b_{\gamma}c_{\gamma}+1}}$.

At $d_{1}=b_{\gamma}(c_{\gamma}+1)/(b_{\gamma}c_{\gamma}+1)<2$, the
convergence rate of KAR estimator
$m^{-b_{\gamma}c_{\gamma}/(b_{\gamma}c_{\gamma}+1)}$ is optimal. This
statistically efficient rate is calibrated by $b_{\gamma}$, the effective
input dimension, together with $c_{\gamma}$, the smoothness of the operator
$H^{\gamma}$. The condition
$d_{1}=b_{\gamma}(c_{\gamma}+1)/(b_{\gamma}c_{\gamma}+1)<2$ also suggests that
$n_{1}>m$.

Additional results and discussions including the two-stage approach are
included in the Appendix.

###  4.2 Causal effect and target KAR estimate

In this section, we discuss the scenarios assuming that the data are generated
from a structural causal model with nonlinear features as shown below,

$\displaystyle\begin{pmatrix}C\\\ \psi(X)\\\
Y\end{pmatrix}=B\begin{pmatrix}\phi(Z)\\\ C\\\ \psi(X)\\\
Y\end{pmatrix}+\begin{pmatrix}\epsilon_{C}\\\ \epsilon_{X}\\\
\epsilon_{Y}\end{pmatrix},$ (22)

where we write operator $B$ in the following matrix form

$B=\begin{pmatrix}B_{CZ}&0&0&0\\\ B_{XZ}&B_{XC}&0&0\\\
B_{YZ}&B_{YC}&B_{YX}&0\end{pmatrix}.$

We note that each operator $B_{\triangle\square}$ represents an operator that
takes an element from $\square$-related space to $\triangle$-related space,
e.g. $B_{XZ}:{\mathcal{H}}_{\mathcal{Z}}\to{\mathcal{H}}_{\mathcal{X}}$ and
$B_{YZ}:{\mathcal{H}}_{\mathcal{Z}}\to{\mathcal{Y}}$. The noise variables
$\epsilon_{Z}$, $\epsilon_{C}$, $\epsilon_{X}$ and $\epsilon_{Y}$ are
independent of each other. Let $\Sigma_{Z}$, $\Sigma_{C}$, $\Sigma_{X}$ and
$\Sigma_{Y}$ denote the covariance of $\epsilon_{Z}$, $\epsilon_{C}$,
$\epsilon_{X}$ and $\epsilon_{Y}$, respectively. Here each operator in $B$
represents a line in the model shown in Figure 2. For instance, $B_{CZ}$
stands for the line from ${\mathcal{H}}_{\mathcal{Z}}$ to $C$; $B_{YX}$
corresponds to the line from ${\mathcal{H}}_{\mathcal{X}}$ to $Y$. $B_{YX}$ in
Equation 22 reflects the causal effect we are interested in. We study the
identifiability scenarios where operator $B_{YX}$ can be learned via KAR
estimator $H^{\gamma}$.

######  Theorem 2.

An operator $B_{XZ}$ is a zero operator written by $B_{XZ}=0$, if
$\langle\psi(x),B_{XZ}\phi(z)\rangle_{{\mathcal{H}}_{\mathcal{X}}}=0$,
$\forall\psi(X)\in{\mathcal{H}}_{\mathcal{X}},\phi(z)\in{\mathcal{H}}_{\mathcal{Z}}$.
Operator $B_{CZ}=0$ if $c^{\top}B_{CZ}\phi(z)=0$, $\forall
c\in\mathcal{C},\phi(z)\in{\mathcal{H}}_{\mathcal{Z}}$. A matrix-valued
operator, e.g. $B_{YC}=0$ if all entries are $0$. For data generation process
following Equation 22, we have $H^{\gamma}=B_{YX}$ in following cases.

  * (i)

$B_{YC}=0$ and $\gamma=0$, i.e. no latent confounder.

  * (ii)

$B_{YZ}+B_{YC}B_{CZ}=0$ and $\gamma=\infty$, where kernel IV is a special
case, i.e. both $B_{YZ}=0$ and $B_{CZ}=0$.

  * (iii)

$B_{YC}=0$, $B_{YZ}+B_{YC}B_{CZ}=0$, and $\gamma\geq 0$.

  * (iv)

$\Sigma_{YX}^{||}=-a\Sigma_{YX}^{\bot}$ for some $a>0$, and $\gamma=1/a$.

$\Sigma_{YX}^{||}=(B_{YZ}+B_{YC}B_{CZ})\Sigma_{Z}(B_{ZX}+B_{ZC}B_{CX})$
denotes the covariance between $\psi(X)$ and $Y$ projected on the linear span
from the components of $\phi(Z)$; and
$\Sigma_{YX}^{\bot}=B_{YC}\Sigma_{C}B_{CX}$ denote the covariance between the
residuals of $\psi(X)$ and $Y$.

Thereom 2 (i) suggests that KPA is optimal when there is no unobserved
confounder; (ii) is a generalized condition including KIV; (iii) shows the KAR
estimator identifies the causal relation from $X$ to $Y$ regardless of
$\gamma$ with generalized KIV condition in (ii) and no latent confounder in
(i); (iv) shows the KAR identifiability condition with appropriate choice of
$\gamma$ when $\Sigma_{XY}^{||}$ and $\Sigma_{XY}^{\bot}$ are in the flipped
direction. In the next section, we show the empirical results for KAR and
relevant baseline methods.

##  5  Empirical Results

###  5.1 Synthetic experiments

We consider the data generating process of the following nonlinear structural
equation,

$Y=0.75C-0.25Z+\ln(|16X-8|+1)sgn(X-0.5),$

where $sgn(x)\in\\{-1,0,+1\\}$ denotes the sign of $x$. The explanatory
variables $X,Z,C$ are generated from

$\displaystyle\begin{pmatrix}C\\\ V\\\ W\end{pmatrix}\sim
N\left(\begin{pmatrix}0\\\ 0\\\ 0\end{pmatrix},\begin{pmatrix}1,0.3,0.2\\\
0.3,1,0\\\ 0.2,0,1\end{pmatrix}\right),$ $\displaystyle
X=F\left(\frac{W+V}{\sqrt{2}}\right),\quad Z=F(W)-0.5,$

where $F$ denote the c.d.f of standard normal distribution.

We generate $\\{(x_{i},y_{i},z_{i})\\}_{i\in[N]}$ samples with $N=700$. To
perform the data-splitting procedures described in Section 3, we set
$n_{1}=n_{2}=250$ and $n=500$ for a fair comparison in the projection stage
(Section 3.1); and $m=200$ in the regression stage (Section 3.2). We set
regularizers as $\alpha_{1}=1.5n_{1}^{-0.5}$, $\alpha_{2}=1.5n_{2}^{-0.5}$,
$\alpha=1.5n^{-0.5}$ and $\xi=1.5m^{-0.5}$.

Figure 3:  Restuls for the synthetic example: fitted (top left) nonlinear
models; and (top right) linear models; (bottom): log MSE.

##### Fitting methods

We consider estimations via the three-stage kernel anchor regression with
disjoint data set projection (KAR) and two-stage kernel anchor regression with
joint data set projection (KAR.2). The baseline approaches include the kernel-
based nonlinear methods: kernel instrument variable regression (KIV), kernel
partialling out regression (KPA), kernel ridge regression (KReg); and the
linear models: linear anchor regression (AR), linear instrument variable
regression (IV), linear partialling out regression (PA) and ordinary least
square (OLS). We use Gaussian kernel for all kernel methods, where the median
heuristic is used for choosing the bandwidth [Gretton et al.,, 2012]. For the
synthetic example, we set $\gamma=2$ for all anchor regressions (KAR, KAR.2
and AR).

For each algorithm, we implement 50 trials and calculate the mean squared
error (MSE) with respect to the true causal model
${\mathbb{E}}(Y|do(x))$444Setting a particular value $X=x$ while ignoring
other variables that may potentially changing the distribution of $y$,
$p(y|X=x)$ is noted as $p(y|do(x)$ [Pearl,, 2009; Peters et al.,, 2016].
${\mathbb{E}}[Y|do(x)]$ is set us the mean over $p(y|do(x))$ averaging out
different $Z$ values in this case., which can be computed from the structural
model. A trial is shown in Figure 3 as a visual example. We can see that the
KAR produces a closest estimation to the true model among all other methods
and outperforms KAR.2. The comparison with linear models are also shown. IV
model fits better than other linear models. We report $\log_{10}(\text{MSE})$
in the bottom of Figure 3, which shows that both KAR methods have smaller
errors than others. KAR performs slightly better than KAR.2 in this case. To
check the robustness of KAR estimators, we study a less smooth variant of the
generating process and show the results in Section B.2.

Figure 4:  Effects of different $\gamma$ choices on MSE.

##### The effect of $\gamma$ choices

To investigate how the change of $\gamma$ affects the estimator, we consider
KIV as our baseline as the IV setting corresponds to $\gamma\to\infty$. We
consider the data generating process used in the KIV paper [Singh et al.,,
2019]. The $\log_{10}$(MSE) results of KAR and KAR.2, in comparison with KIV,
are shown in Figure 4. For the simulation, we set $N=1000$, $n_{1}=n_{2}=200$,
$n=n_{1}+n_{2}=400$ and $m=600$. From the result, we can see that both KAR and
KAR.2 achieves smaller error when choosing $\gamma=2$. Data generation and
model implementation details are included in Section B.1.

Figure 5:  Prediction error with distributional intervention.

##### Intervention and Generalization

To evaluate the robustness and generalization performance of both KAR
estimators under distribution shift, as discussed in Rothenhäusler et al.,
[2018], we intervene the anchor variable $Z$. We train the model on a
subpopulation of samples with $Z<0$ and test on the samples with $Z\geq 0$.
The performance is measured by prediction error (PE) of fitted model with
respect to ${\mathbb{E}}(Y|X=x,Z\geq 0)$, where the true conditional model is
not known in closed form but estimated from samples.

We also exchange the training set and the testing set. As shown in Figure 5,
our KAR estimator has the lowest PE among others, showing better out-of-
distribution generalization performance. More importantly, by checking the two
(flipped) scenarios, i.e. train on $Z<0$ v.s. train on $Z\geq 0$, we also see
that KAR is the most invariant in terms of PE. On the contrary, linear version
of AR and IV achieves very different PE in both cases. Variances of PE for KPA
are also very different in the two cases. Despite KReg achieves a relatively
low PE in both cases, the distributions of PE can be found very different.

Figure 6:  Fitted models (top) and prediction errors (bottom) when training on
male subjects and testing on female subjects.

###  5.2 Real-world application

We consider the smoking dataset extracted from National Medical Expenditure
Survey (NMES) [Johnson et al.,, 2003] to study the effect of smoking amount on
medical expenditure [Imai and Van Dyk,, 2004]555The dataset is accessible
through using the R package for “estimating causal dose response function”
causaldrf [Galagate,, 2016]
https://cran.r-project.org/web/packages/causaldrf/index.html..

The treatment variable $X$ is the $\log$ of smoking amount, the outcome $Y$ is
the $\log$ of medical expenditure, and the anchor $Z$ is set to be the last
age for smoking. We use $1000$ samples, randomly selected from $9708$
available samples, to fit the model. We set $n_{1}=n_{2}=300$, $n=600$ and
$m=400$. We also set $\gamma=2.9$ and apply Gaussian kernel with median
heuristic bandwidth [Gretton et al.,, 2012] for all kernel methods. As shown
in the upper part of Figure 6, KAR estimators show that the effect of $X$ on
$Y$ is more significant when $X\in[-2,1]$ compared to $X\in[1,4]$. Our method
can also be used in complement with the approaches finding causal directions,
e.g. [Peters et al.,, 2016]666Implementation with R package CAM can be found
at https://rdrr.io/cran/CAM/man/CAM.html. We run the CAM to ensure that there
is a causal effect in the direction from $X$ to $Y$ and KAR procedure further
learns the specific function representing such effect. However, existing work
such as propensity score approaches [Imai and Van Dyk,, 2004] did not manage
to extract such causal relationship between smoking and medical expenditure.

To strengthen our finding, we quantify the performance of the estimators.
Since we do not know the real generating process of the data, we cannot
compare the MSE as Figure 3 and 4. Instead, it’s feasible to evaluate the
performance of estimators under distribution perturbation via PE, similar to
Figure 5. We train models on male subjects and compute the prediction accuracy
of fitted model on female subjects. The results are shown in the bottom of
Figure 6. From the result, we see that both KAR approaches outperform other
kernel-based approaches as well as the linear version of AR, suggesting a
better learned effect from the smoking amount to medical expenditure.

##  6  Conclusion

In this work, we consider learning a more general class of causal DAG in a
nonlinear setting using kernelized anchor regression. By considering different
data splitting strategies to estimate the projection operators, we show that
the three-stage approach not only performs better empirically than baseline
approaches as well as the 2SLS approach, but also achieves optimal rate under
given conditions. Identifiability results are provided and are shown to
generalize KIV and “no latent confounder” scenarios. For the future, data
adaptive choice of $\gamma$ can be an interesting direction to explore.

### Acknowledgements

The authors thank Ana Korba for helpful discussions. W.X. acknowledges the
support from EPSRC grant EP/T018445/1.

## References

  * Angrist et al., [1996] Angrist, J. D., Imbens, G. W., and Rubin, D. B. (1996).  Identification of causal effects using instrumental variables.  Journal of the American statistical Association, 91(434):444–455.
  * Arjovsky et al., [2019] Arjovsky, M., Bottou, L., Gulrajani, I., and Lopez-Paz, D. (2019).  Invariant risk minimization.  arXiv preprint arXiv:1907.02893.
  * Bennett et al., [2019] Bennett, A., Kallus, N., and Schnabel, T. (2019).  Deep generalized method of moments for instrumental variable analysis.  Advances in neural information processing systems, 32.
  * Besserve et al., [2019] Besserve, M., Mehrjou, A., Sun, R., and Schölkopf, B. (2019).  Counterfactuals uncover the modular structure of deep generative models.  In International Conference on Learning Representations.
  * Bishop, [1998] Bishop, C. M. (1998).  Latent variable models.  In Learning in graphical models, pages 371–403. Springer.
  * Blei, [2012] Blei, D. M. (2012).  Probabilistic topic models.  Communications of the ACM, 55(4):77–84.
  * Bowden and Turkington, [1990] Bowden, R. J. and Turkington, D. A. (1990).  Instrumental variables.  Number 8. Cambridge university press.
  * Bühlmann, [2018] Bühlmann, P. (2018).  Invariance, causality and robustness.  arXiv preprint arXiv:1812.08233.
  * Centorrino et al., [2019] Centorrino, S., Fève, F., and Florens, J.-P. (2019).  Nonparametric instrumental regressions with (potentially discrete) instruments independent of the error term.  arXiv preprint arXiv:1905.07812.
  * Chen and Christensen, [2018] Chen, X. and Christensen, T. M. (2018).  Optimal sup-norm rates and uniform inference on nonlinear functionals of nonparametric iv regression.  Quantitative Economics, 9(1):39–84.
  * Chernozhukov et al., [2013] Chernozhukov, V., Fernández-Val, I., and Melly, B. (2013).  Inference on counterfactual distributions.  Econometrica, 81(6):2205–2268.
  * Colombo et al., [2012] Colombo, D., Maathuis, M. H., Kalisch, M., and Richardson, T. S. (2012).  Learning high-dimensional directed acyclic graphs with latent and selection variables.  The Annals of Statistics, pages 294–321.
  * Fukumizu et al., [2007] Fukumizu, K., Gretton, A., Sun, X., and Schölkopf, B. (2007).  Kernel measures of conditional dependence.  Advances in neural information processing systems, 20.
  * Galagate, [2016] Galagate, D. (2016).  Causal inference with a continuous treatment and outcome: Alternative estimators for parametric dose-response functions with applications.  PhD thesis, University of Maryland, College Park.
  * Glymour et al., [2019] Glymour, C., Zhang, K., and Spirtes, P. (2019).  Review of causal discovery methods based on graphical models.  Frontiers in genetics, 10:524.
  * Gretton et al., [2012] Gretton, A., Borgwardt, K. M., Rasch, M. J., Schölkopf, B., and Smola, A. (2012).  A kernel two-sample test.  The Journal of Machine Learning Research, 13(1):723–773.
  * Gretton et al., [2005] Gretton, A., Bousquet, O., Smola, A., and Schölkopf, B. (2005).  Measuring statistical dependence with hilbert-schmidt norms.  In International conference on algorithmic learning theory, pages 63–77. Springer.
  * Grünewälder et al., [2012] Grünewälder, S., Lever, G., Baldassarre, L., Patterson, S., Gretton, A., and Pontil, M. (2012).  Conditional mean embeddings as regressors.  In Proceedings of the 29th International Coference on International Conference on Machine Learning, pages 1803–1810.
  * Hartford et al., [2017] Hartford, J., Lewis, G., Leyton-Brown, K., and Taddy, M. (2017).  Deep IV: A flexible approach for counterfactual prediction.  In International Conference on Machine Learning, pages 1414–1423. PMLR.
  * Hofmann et al., [2008] Hofmann, T., Schölkopf, B., and Smola, A. J. (2008).  Kernel methods in machine learning.  The annals of statistics, 36(3):1171–1220.
  * Hoyer et al., [2008] Hoyer, P., Janzing, D., Mooij, J. M., Peters, J., and Schölkopf, B. (2008).  Nonlinear causal discovery with additive noise models.  Advances in neural information processing systems, 21.
  * Hyvärinen, [2013] Hyvärinen, A. (2013).  Independent component analysis: recent advances.  Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 371(1984):20110534.
  * Hyvärinen and Dayan, [2005] Hyvärinen, A. and Dayan, P. (2005).  Estimation of non-normalized statistical models by score matching.  Journal of Machine Learning Research, 6(4).
  * Hyvarinen and Morioka, [2017] Hyvarinen, A. and Morioka, H. (2017).  Nonlinear ica of temporally dependent stationary sources.  In Artificial Intelligence and Statistics, pages 460–469. PMLR.
  * Hyvärinen et al., [2010] Hyvärinen, A., Zhang, K., Shimizu, S., and Hoyer, P. O. (2010).  Estimation of a structural vector autoregression model using non-gaussianity.  Journal of Machine Learning Research, 11(5).
  * Imai and Van Dyk, [2004] Imai, K. and Van Dyk, D. A. (2004).  Causal inference with general treatment regimes: Generalizing the propensity score.  Journal of the American Statistical Association, 99(467):854–866.
  * Johansson et al., [2020] Johansson, F. D., Shalit, U., Kallus, N., and Sontag, D. (2020).  Generalization bounds and representation learning for estimation of potential outcomes and causal effects.  arXiv preprint arXiv:2001.07426.
  * Johnson et al., [2003] Johnson, E., Dominici, F., Griswold, M., and Zeger, S. L. (2003).  Disease cases and their medical costs attributable to smoking: an analysis of the national medical expenditure survey.  Journal of Econometrics, 112(1):135–151.
  * Kallus, [2020] Kallus, N. (2020).  Deepmatch: Balancing deep covariate representations for causal inference using adversarial training.  In International Conference on Machine Learning, pages 5067–5077. PMLR.
  * Khemakhem et al., [2021] Khemakhem, I., Monti, R., Leech, R., and Hyvarinen, A. (2021).  Causal autoregressive flows.  In International conference on artificial intelligence and statistics, pages 3520–3528. PMLR.
  * Kook et al., [2022] Kook, L., Sick, B., and Bühlmann, P. (2022).  Distributional anchor regression.  Statistics and Computing, 32(3):1–19.
  * Louizos et al., [2017] Louizos, C., Shalit, U., Mooij, J. M., Sontag, D., Zemel, R., and Welling, M. (2017).  Causal effect inference with deep latent-variable models.  Advances in neural information processing systems, 30.
  * Micchelli and Pontil, [2005] Micchelli, C. A. and Pontil, M. (2005).  On learning vector-valued functions.  Neural computation, 17(1):177–204.
  * Monti et al., [2020] Monti, R. P., Zhang, K., and Hyvärinen, A. (2020).  Causal discovery with general non-linear relationships using non-linear ica.  In Uncertainty in artificial intelligence, pages 186–195. PMLR.
  * Mooij et al., [2009] Mooij, J., Janzing, D., Peters, J., and Schölkopf, B. (2009).  Regression by dependence minimization and its application to causal inference in additive noise models.  In Proceedings of the 26th annual international conference on machine learning, pages 745–752.
  * Muandet et al., [2017] Muandet, K., Fukumizu, K., Sriperumbudur, B., Schölkopf, B., et al. (2017).  Kernel mean embedding of distributions: A review and beyond.  Foundations and Trends® in Machine Learning, 10(1-2):1–141.
  * Muandet et al., [2021] Muandet, K., Kanagawa, M., Saengkyongam, S., and Marukatat, S. (2021).  Counterfactual mean embeddings.  Journal of Machine Learning Research, 22(162):1–71.
  * Oberst et al., [2021] Oberst, M., Thams, N., Peters, J., and Sontag, D. (2021).  Regularizing towards causal invariance: Linear models with proxies.  In International Conference on Machine Learning, pages 8260–8270. PMLR.
  * Pearl, [2009] Pearl, J. (2009).  Causality.  Cambridge university press.
  * Pearl et al., [2000] Pearl, J. et al. (2000).  Models, reasoning and inference.  Cambridge, UK: CambridgeUniversityPress, 19(2).
  * Pearl et al., [2016] Pearl, J., Glymour, M., and Jewell, N. P. (2016).  Causal Inference in Statistics: A Primer.  John Wiley & Sons.
  * Peters et al., [2016] Peters, J., Bühlmann, P., and Meinshausen, N. (2016).  Causal inference by using invariant prediction: identification and confidence intervals.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 78(5):947–1012.
  * Peters et al., [2017] Peters, J., Janzing, D., and Schölkopf, B. (2017).  Elements of causal inference: foundations and learning algorithms.  MIT press.
  * Peters et al., [2011] Peters, J., Mooij, J. M., Janzing, D., and Schölkopf, B. (2011).  Identifiability of causal graphs using functional models.  In Proceedings of the Twenty-Seventh Conference on Uncertainty in Artificial Intelligence, pages 589–598.
  * Peters et al., [2014] Peters, J., Mooij, J. M., Janzing, D., and Schölkopf, B. (2014).  Causal discovery with continuous additive noise models.  Journal of Machine Learning Research, 15:2009–2053.
  * Rabiner and Juang, [1986] Rabiner, L. and Juang, B. (1986).  An introduction to hidden markov models.  ieee assp magazine, 3(1):4–16.
  * Rolland et al., [2022] Rolland, P., Cevher, V., Kleindessner, M., Russell, C., Janzing, D., Schölkopf, B., and Locatello, F. (2022).  Score matching enables causal discovery of nonlinear additive noise models.  In International Conference on Machine Learning, pages 18741–18753. PMLR.
  * Rothenhäusler et al., [2018] Rothenhäusler, D., Meinshausen, N., Bühlmann, P., and Peters, J. (2018).  Anchor regression: heterogeneous data meets causality.  arXiv preprint arXiv:1801.06229.
  * Rubin, [2004] Rubin, D. B. (2004).  Direct and indirect causal effects via potential outcomes.  Scandinavian Journal of Statistics, 31(2):161–170.
  * Rubin, [2005] Rubin, D. B. (2005).  Causal inference using potential outcomes: Design, modeling, decisions.  Journal of the American Statistical Association, 100(469):322–331.
  * Schölkopf et al., [2021] Schölkopf, B., Locatello, F., Bauer, S., Ke, N. R., Kalchbrenner, N., Goyal, A., and Bengio, Y. (2021).  Toward causal representation learning.  Proceedings of the IEEE, 109(5):612–634.
  * Scholkopf and Smola, [2018] Scholkopf, B. and Smola, A. J. (2018).  Learning with kernels: support vector machines, regularization, optimization, and beyond.  MIT press.
  * Shimizu et al., [2006] Shimizu, S., Hoyer, P. O., Hyvärinen, A., Kerminen, A., and Jordan, M. (2006).  A linear non-gaussian acyclic model for causal discovery.  Journal of Machine Learning Research, 7(10).
  * Shimizu et al., [2011] Shimizu, S., Inazumi, T., Sogawa, Y., Hyvärinen, A., Kawahara, Y., Washio, T., Hoyer, P. O., and Bollen, K. (2011).  Directlingam: A direct method for learning a linear non-gaussian structural equation model.  The Journal of Machine Learning Research, 12:1225–1248.
  * Singh et al., [2019] Singh, R., Sahani, M., and Gretton, A. (2019).  Kernel instrumental variable regression.  arXiv preprint arXiv:1906.00232.
  * Song et al., [2009] Song, L., Huang, J., Smola, A., and Fukumizu, K. (2009).  Hilbert space embeddings of conditional distributions with applications to dynamical systems.  In Proceedings of the 26th Annual International Conference on Machine Learning, pages 961–968.
  * Spirtes et al., [2000] Spirtes, P., Glymour, C. N., Scheines, R., and Heckerman, D. (2000).  Causation, prediction, and search.  MIT press.
  * Spirtes et al., [2013] Spirtes, P. L., Meek, C., and Richardson, T. S. (2013).  Causal inference in the presence of latent variables and selection bias.  arXiv preprint arXiv:1302.4983.
  * Steinwart and Christmann, [2008] Steinwart, I. and Christmann, A. (2008).  Support vector machines.  Springer Science & Business Media.
  * Szabó et al., [2015] Szabó, Z., Gretton, A., Póczos, B., and Sriperumbudur, B. (2015).  Two-stage sampled learning theory on distributions.  In Artificial Intelligence and Statistics, pages 948–957. PMLR.
  * Szabó et al., [2016] Szabó, Z., Sriperumbudur, B. K., Póczos, B., and Gretton, A. (2016).  Learning theory for distribution regression.  The Journal of Machine Learning Research, 17(1):5272–5311.
  * Xu et al., [2020] Xu, L., Chen, Y., Srinivasan, S., de Freitas, N., Doucet, A., and Gretton, A. (2020).  Learning deep features in instrumental variable regression.  arXiv preprint arXiv:2010.07154.
  * Xu et al., [2021] Xu, L., Kanagawa, H., and Gretton, A. (2021).  Deep proxy causal learning and its application to confounded bandit policy evaluation.  Advances in Neural Information Processing Systems, 34:26264–26275.
  * Zhang et al., [2018] Zhang, K., Schölkopf, B., Spirtes, P., and Glymour, C. (2018).  Learning causality and causality-related learning: some recent progress.  National science review, 5(1):26–29.
  * Zhu et al., [2022] Zhu, Y., Gultchin, L., Gretton, A., Kusner, M. J., and Silva, R. (2022).  Causal inference with treatment measurement error: a nonparametric instrumental variable approach.  In Uncertainty in Artificial Intelligence, pages 2414–2424. PMLR.

Supplementary Material for Nonlinear Causal Discovery
via Kernel Anchor Regression

##  Appendix A Proofs and derivations

###  A.1 Proof of Theorem 1

Before proving Theorem 1, we introduce the exact bounds of the approximation
errors for estimating $E_{X}^{p}$ and $E_{Y}^{p}$ in the disjoint sample sets
projection stage. Lemma A1 and A2 below are adapted from Theorem 2 in Singh et
al., [2019].

######  Lemma A1.

Under Condition 1, $\forall\delta\in(0,1)$, the following holds w.p.
$1-\delta$:

$\displaystyle\|E_{\alpha_{1},X}^{n_{1}}-E_{X}^{p}\|_{{\mathcal{H}}_{\Gamma}}\leq
r_{E_{1}}(\delta,n_{1},c_{1}):=\frac{\sqrt{\zeta_{1}}(c_{1}+1)}{4^{\frac{1}{c_{1}+1}}}\left(\frac{4\kappa(Q_{1}+\kappa\|E_{X}^{p}\|_{{\mathcal{H}}_{\Gamma}}\ln(2/\delta)}{\sqrt{n_{1}\zeta_{1}}(c_{1}-1)}\right)^{\frac{c_{1}-1}{c_{1}+1}},$
$\displaystyle\alpha_{1}=\left(\frac{8\kappa(Q_{1}+\kappa\|E_{X}^{p}\|_{{\mathcal{H}}_{\Gamma}}\ln(2/\delta)}{\sqrt{n_{1}\zeta_{1}}(c_{1}-1)}\right)^{\frac{2}{c_{1}+1}}.$

######  Lemma A2.

Under Condition 1 and Condition 2, $\forall\epsilon\in(0,1)$, the following
holds w.p. $1-\epsilon$:

$\displaystyle\|E_{\alpha_{2},Y}^{n_{2}}-E_{Y}^{p}\|_{{\mathcal{H}}_{\Theta}}\leq
r_{E_{2}}(\epsilon,n_{2},c_{2}):=\frac{\sqrt{\zeta_{2}}(c_{2}+1)}{4^{\frac{1}{c_{2}+1}}}\left(\frac{4\kappa(Q_{2}+\kappa\|E_{Y}^{p}\|_{{\mathcal{H}}_{\Theta}}\ln(2/\epsilon)}{\sqrt{n_{2}\zeta_{2}}(c_{2}-1)}\right)^{\frac{c_{2}-1}{c_{2}+1}},$
$\displaystyle\alpha_{2}=\left(\frac{8\kappa(Q_{2}+\kappa\|E_{Y}^{p}\|_{{\mathcal{H}}_{\Theta}}\ln(2/\epsilon)}{\sqrt{n_{2}\zeta_{2}}(c_{2}-1)}\right)^{\frac{2}{c_{2}+1}}.$

Recall that we define the population-level risk for the regression stage
${\mathcal{E}}^{\gamma}(H)$, population-level risk with regularization
${\mathcal{E}}^{\gamma}_{\xi}(H)$, and the empirical risk
$\widehat{\mathcal{E}}^{\gamma,m}_{\xi}(H)$ with $E_{X}^{p}$ and $E_{Y}^{p}$
being replaced by $E_{\alpha_{1},X}^{n_{1}}$ and $E_{\alpha_{2},Y}^{n_{2}}$,
respectively. Denote the optimal operator to ${\mathcal{E}}^{\gamma}_{\xi}(H)$
as
$H_{\xi}^{\gamma}=\operatorname*{arg\,min}_{H}{\mathcal{E}}_{\xi}^{\gamma}(H)$.
We now define the empirical risk ${\mathcal{E}}_{\gamma,m}^{\xi}(H)$ with true
$E_{X}^{p}$ and $E_{Y}^{p}$, and the corresponding optimal operator.

$\displaystyle{\mathcal{E}}_{\xi}^{\gamma,m}(H)=\frac{1}{m}\sum_{l=1}^{m}\|y_{\gamma,l}-H\psi_{\gamma,l}\|_{{\mathcal{Y}}}^{2}+\xi\|H\|_{{\mathcal{H}}_{\Omega}}^{2},\quad
H_{\xi}^{\gamma,m}=\operatorname*{arg\,min}_{H}{\mathcal{E}}_{\xi}^{\gamma,m}(H),$

where the true transformed inputs and outputs are given by

$\psi_{\gamma,l}=\psi(x_{l})+(\sqrt{\gamma}-1)E_{X}^{p}\phi(z_{l})\in{\mathcal{H}}_{\mathcal{X}},\quad
y_{\gamma,l}=y_{l}+(\sqrt{\gamma}-1)E_{Y}^{p}\phi(z_{l})\in{\mathcal{Y}}.$

The closed form solution of $H_{\xi}^{\gamma,m}$ is given by Lemma A3 below,
and it’s adapted from Theorem 3 in Singh et al., [2019]

######  Lemma A3.

$\forall\xi>0$, the solution $H_{\xi}^{\gamma,m}$ to
${\mathcal{E}}_{\xi}^{\gamma,m}$ exists, is unique, and

$\displaystyle\mathbf{T}=\frac{1}{m}\sum_{l=1}^{m}T_{\psi_{\gamma,l}},\quad\mathbf{g}=\frac{1}{m}\sum_{l=1}^{m}\Omega_{\psi_{\gamma,l}}y_{\gamma,l},$
$\displaystyle H_{\xi}^{\gamma,m}=(\mathbf{T}+\xi)^{-1}\circ\mathbf{g}.$

We then define the following terms.

######  Definition 1.

Fix $\eta\in(0,1)$ and define the following constants

$\displaystyle C_{\eta}=96\ln^{2}(6/\eta),\quad
M=2(C+\|H^{\gamma}\|_{{\mathcal{H}}_{\Omega}}\sqrt{B}),\quad\Sigma=\frac{M}{2}.$

For the excess error of KAR estimator $\widehat{H}_{\xi}^{\gamma,m}$, we can
bound it by five terms according to Proposition 32 in Singh et al., [2019].

######  Lemma A4.

The excess error can be bounded as follows

$\displaystyle{\mathcal{E}}^{\gamma}(\widehat{H}^{\gamma,m}_{\xi})-{\mathcal{E}}^{\gamma}(H^{\gamma})\leq
5[S_{-1}+S_{0}+{\mathcal{A}}(\xi)+S_{1}+S_{2}],$

where

$\displaystyle S_{-1}$ $\displaystyle=$
$\displaystyle\|\sqrt{T}\circ(\widehat{\mathbf{T}}+\xi)^{-1}(\widehat{\mathbf{g}}-\mathbf{g})\|_{{\mathcal{H}}_{\Omega}}^{2},$
$\displaystyle S_{0}$ $\displaystyle=$
$\displaystyle\|\sqrt{T}\circ(\widehat{\mathbf{T}}+\xi)^{-1}(\mathbf{T}-\widehat{\mathbf{T}})H^{\gamma,m}_{\xi}\|_{{\mathcal{H}}_{\Omega}}^{2},$
$\displaystyle S_{1}$ $\displaystyle=$
$\displaystyle\|\sqrt{T}\circ(\widehat{\mathbf{T}}+\xi)^{-1}(\mathbf{g}-\mathbf{T}H^{\gamma})\|_{{\mathcal{H}}_{\Omega}}^{2},$
$\displaystyle S_{2}$ $\displaystyle=$
$\displaystyle\|\sqrt{T}\circ(\widehat{\mathbf{T}}+\xi)^{-1}(T-\mathbf{T})(H^{\gamma}_{\xi}-H^{\gamma})\|_{{\mathcal{H}}_{\Omega}}^{2},$
$\displaystyle{\mathcal{A}}(\xi)$ $\displaystyle=$
$\displaystyle\|\sqrt{T}(H_{\xi}^{\gamma}-H^{\gamma})\|_{{\mathcal{H}}_{\Omega}}^{2}.$

For all five terms above, only $\widehat{\mathbf{g}}-\mathbf{g}$ in $S_{-1}$
depends on the approximation error of $E_{Y}^{p}$. The bounds for other four
terms are same to the KIV case. Below we introduce without proof the bond of
$S_{0}$, $S_{1}$, $S_{2}$ and ${\mathcal{A}}(\xi)$ according to Theorem 7 in
Singh et al., [2019].

######  Lemma A5.

Under Condition 1–3, if $m$ is large enough and
$\xi\leq\|T\|_{L({\mathcal{H}}_{\Omega})}$ then $\forall\delta,\eta\in(0,1)$,
the following holds up w.p. $1-\eta-\delta$:

$\displaystyle S_{0}$ $\displaystyle\leq$
$\displaystyle\frac{4}{\xi}4BL^{2}r_{x}^{2\iota}\|H^{\gamma,m}_{\xi}\|^{2}_{{\mathcal{H}}_{\Omega}},$
$\displaystyle S_{1}$ $\displaystyle\leq$ $\displaystyle
32\ln^{2}(6\eta)\left[\frac{BM^{2}}{m^{2}\xi}+\frac{\Sigma^{2}}{m}\beta^{1/b_{\gamma}}\frac{\pi/b_{\gamma}}{\sin(\pi\b{)}\xi^{-1/b_{\gamma}}}\right],$
$\displaystyle S_{2}$ $\displaystyle\leq$ $\displaystyle
8\ln^{2}(6/\eta)\left[\frac{4B^{2}\zeta\xi^{c_{\gamma}-1}}{m^{2}\xi}+\frac{B\zeta\xi^{c_{\gamma}}}{m\xi}\right],$
$\displaystyle{\mathcal{A}}(\xi)$ $\displaystyle\leq$
$\displaystyle\zeta\xi^{c_{\gamma}}.$

To extend the convergence rate of KIV estimator to KAR estimator. We then
illustrate the bound for $S_{-1}$. To begin with, the bound of term
$\sqrt{T}\circ(\widehat{\mathbf{T}}+\xi)^{-1}$ in $S_{-1}$ is given by
Proposition 39 in Singh et al., [2019].

######  Lemma A6.

If
$\|\widehat{\psi}_{\gamma}-\psi_{\gamma}\|_{{\mathcal{H}}_{\mathcal{X}}}\leq
r_{x}$ w.p. $1-\delta$, $\xi\leq\|T\|_{L({\mathcal{H}}_{\Omega})}$, $m$ is
sufficiently large and Condition 3 holds, then w.p. $1-\eta/3-\delta$

$\|\sqrt{T}\circ(\widehat{\mathbf{T}}+\xi)^{-1}\|_{L({\mathcal{H}}_{\Omega})}\leq\frac{2}{\sqrt{\xi}}.$

With the the error propagated from the estimators in the projection stage, we
can bound $\widehat{\psi}_{\gamma}-\psi_{\gamma}$ and
$\widehat{y}_{\gamma}-y_{\gamma}$ as shown in Lemma A7–A8.

######  Lemma A7.

Under Condition 1, $\forall\delta\in(0,1)$, the following statement holds w.p.
$1-\delta$: $\forall z\in{\mathcal{Z}},x\in{\mathcal{X}}$,

$\|\widehat{\psi}_{\gamma}-\psi_{\gamma}\|_{{\mathcal{H}}_{\mathcal{X}}}\leq
r_{x}(\gamma,\delta,n_{1},c_{1}):=|\sqrt{\gamma}-1|\kappa
r_{E_{1}}(\delta,n_{1},c_{1}).$

###### Proof.

By definition, we have

$\displaystyle\|\widehat{\psi}_{\gamma}-\psi_{\gamma}\|_{{\mathcal{H}}_{\mathcal{X}}}$
$\displaystyle=$
$\displaystyle\|\left(\sqrt{\gamma}-1\right)\left(E^{n_{1}}_{\alpha_{1},X}-E_{X}^{p}\right)\phi(z)\|_{{\mathcal{H}}_{\mathcal{X}}}$
$\displaystyle\leq$
$\displaystyle|\sqrt{\gamma}-1|\|E^{n_{1}}_{\alpha_{1},X}-E_{X}^{p}\|_{{\mathcal{H}}_{\Gamma}}\|\phi(z)\|_{{\mathcal{H}}_{\mathcal{Z}}}.$

This, together with Lemma A1 and Condition 1, ensures that w.p. $1-\delta$

$\|\widehat{\psi}_{\gamma}-\psi_{\gamma}\|_{{\mathcal{H}}_{\mathcal{X}}}\leq
r_{x}(\gamma,\delta,n_{1},c_{1}):=|\sqrt{\gamma}-1|\kappa
r_{E_{1}}(\delta,n_{1},c_{1}).$

∎

######  Remark A1.

Corollary 1 in Singh et al., [2019] is a special case of Lemma A7 with
$\gamma=0$.

######  Lemma A8.

Under Condition 1– 2, $\forall\epsilon\in(0,1)$, the following statement holds
w.p. $1-\epsilon$: $\forall z\in{\mathcal{Z}},y\in{\mathcal{Y}}$,

$\|\widehat{y}_{\gamma}-y_{\gamma}\|_{{\mathcal{H}}_{\mathcal{Y}}}\leq
r_{y}(\gamma,\epsilon,n_{2},c_{2}):=|\sqrt{\gamma}-1|\kappa
r_{E_{2}}(\epsilon,n_{2},c_{2}).$

###### Proof.

Lemma A8 is analogous to Lemma A7 by replacing $\psi_{\gamma}$ with
$y_{\gamma}$. The proof is thus omitted. ∎

Combining Lemma A6\- A8, we can derive the bound of
$\widehat{\mathbf{g}}-\mathbf{g}$ and then the bound of $S_{-1}$.

######  Lemma A9.

If
$\|\widehat{\psi}_{\gamma}-\psi_{\gamma}\|_{{\mathcal{H}}_{\mathcal{X}}}\leq
r_{x}$ w.p. $1-\delta$ and
$\|\widehat{y}_{\gamma}-y_{\gamma}\|_{{\mathcal{Y}}}\leq r_{y}$ w.p.
$1-\epsilon$, then w.p. $1-\delta-\epsilon$

$\|\widehat{\mathbf{g}}-\mathbf{g}\|_{{\mathcal{H}}_{\Omega}}^{2}\leq
3(L^{2}r_{x}^{2\iota}r_{y}^{2}+B^{2}r_{y}^{2}+L^{2}r_{x}^{2\iota}C^{2}).$

###### Proof.

By definition, we have

$\displaystyle\widehat{\mathbf{g}}-\mathbf{g}$ $\displaystyle=$
$\displaystyle\frac{1}{m}\sum_{l=1}^{m}\left(\Omega_{\widehat{\psi}_{\gamma,l}}\widehat{y}_{\gamma,l}-\Omega_{\psi_{\gamma,l}(x)}y_{\gamma,l}\right)$
$\displaystyle=$
$\displaystyle\frac{1}{m}\sum_{l=1}^{m}\left\\{\Omega_{\widehat{\psi}_{\gamma,l}}-\Omega_{\psi_{\gamma,l}}\right\\}\left\\{\widehat{y}_{\gamma,l}-y_{\gamma,l}\right\\}+\Omega_{\widehat{\psi}_{\gamma,l}}\left\\{\widehat{y}_{\gamma,l}-y_{\gamma,l}\right\\}+\left\\{\Omega_{\widehat{\psi}_{\gamma,l}}-\Omega_{\psi_{\gamma,l}}\right\\}y_{\gamma,l}.$

We then have

$\displaystyle\|\widehat{\mathbf{g}}-\mathbf{g}\|_{{\mathcal{H}}_{\Omega}}^{2}$
$\displaystyle\leq$
$\displaystyle\frac{3m}{m^{2}}\sum_{l=1}^{m}\|\left\\{\Omega_{\widehat{\psi}_{\gamma,l}}-\Omega_{\psi_{\gamma,l}}\right\\}\left\\{\widehat{y}_{\gamma,l}-y_{\gamma,l}\right\\}\|_{{\mathcal{H}}_{\Omega}}^{2}+\|\Omega_{\widehat{\psi}_{\gamma,l}}\left\\{\widehat{y}_{\gamma,l}-y_{\gamma,l}\right\\}\|_{{\mathcal{H}}_{\Omega}}^{2}$
$\displaystyle+\|\left\\{\Omega_{\widehat{\psi}_{\gamma,l}}-\Omega_{\psi_{\gamma,l}}\right\\}y_{\gamma,l}\|_{{\mathcal{H}}_{\Omega}}^{2}$
$\displaystyle\leq$
$\displaystyle\frac{3}{m}\sum_{l=1}^{m}\|\Omega_{\widehat{\psi}_{\gamma,l}}-\Omega_{\psi_{\gamma,l}}\|_{L({\mathcal{Y}},{\mathcal{H}}_{\Omega})}^{2}\|\widehat{y}_{\gamma,l}-y_{\gamma,l}\|_{{\mathcal{Y}}}^{2}+\|\Omega_{\psi_{\gamma,l}}\|_{L({\mathcal{Y}},{\mathcal{H}}_{\Omega})}^{2}\|\widehat{y}_{\gamma,l}-y_{\gamma,l}\|_{{\mathcal{Y}}}^{2}$
$\displaystyle+\|\Omega_{\widehat{\psi}_{\gamma,l}}-\Omega_{\psi_{\gamma,l}}\|_{L({\mathcal{Y}},{\mathcal{H}}_{\Omega})}^{2}\|y_{\gamma,l}\|_{{\mathcal{Y}}}^{2}.$

By the boundedness and the Hölder property in Condition 3, we obtain that w.p.
$1-\delta-\epsilon$,

$\displaystyle\|\widehat{\mathbf{g}}-\mathbf{g}\|_{{\mathcal{H}}_{\Omega}}^{2}$
$\displaystyle\leq$
$\displaystyle\frac{3}{m}\sum_{l=1}^{m}L^{2}\|\widehat{\psi}_{\gamma,l}-\psi_{\gamma,l}\|_{{\mathcal{H}}_{\mathcal{X}}}^{2\iota}\|\widehat{y}_{\gamma,l}-y_{\gamma,l}\|_{{\mathcal{Y}}}^{2}+\|\Omega_{\psi_{\gamma,l}}\|_{L({\mathcal{Y}},{\mathcal{H}}_{\Omega})}^{2}\|\widehat{y}_{\gamma,l}-y_{\gamma,l}\|_{{\mathcal{Y}}}^{2}$
$\displaystyle+L^{2}\|\widehat{\psi}_{\gamma,l}-\psi_{\gamma,l}\|_{{\mathcal{H}}_{\mathcal{X}}}^{2\iota}\|y_{\gamma,l}\|_{{\mathcal{Y}}}^{2}$
$\displaystyle\leq$ $\displaystyle
3(L^{2}r_{x}^{2\iota}r_{y}^{2}+B^{2}r_{y}^{2}+L^{2}r_{x}^{2\iota}C^{2}).$

∎

######  Lemma A10.

Under Condition 1–3, then w.p. $1-\delta-\epsilon$

$S_{-1}\leq\frac{4}{\xi}3(L^{2}r_{x}^{2\iota}r_{y}^{2}+B^{2}r_{y}^{2}+L^{2}r_{x}^{2\iota}C^{2}).$

###### Proof.

We can derive from the definition of $S_{-1}$ that

$S_{-1}\leq\|\sqrt{T}\circ(\widehat{\mathbf{T}}+\xi)^{-1}\|_{L({\mathcal{H}}_{\Omega})}^{2}\|\widehat{\mathbf{g}}-\mathbf{g}\|_{{\mathcal{H}}_{\Omega}}^{2}.$

This, together with Lemma A6 and Lemma A9, ensures

$S_{-1}\leq\frac{4}{\xi}3(L^{2}r_{x}^{2\iota}r_{y}^{2}+B^{2}r_{y}^{2}+L^{2}r_{x}^{2\iota}C^{2}).$

∎

We then show the order of the sum $S_{0}+S_{1}+S_{2}+{\mathcal{A}}(\xi)$,
which is adapted from Theorem 4 in Singh et al., [2019].

######  Lemma A11.

Under Condition 1– 3, choose $\alpha_{1}=n_{1}^{-\frac{1}{c_{1}+1}}$,
$n_{1}=m^{\frac{d_{1}(c_{1}+1)}{\iota(c_{1}-1)}}$, where $d_{1}>0$. Let

$f(m)=\frac{1}{m^{2+{d_{1}}}\xi^{3}}+\frac{1}{m^{1+{d_{1}}}\xi^{2+1/b_{\gamma}}}+\frac{1}{m^{d_{1}}\xi}+\xi^{c_{\gamma}}+\frac{1}{m^{2}\xi}+\frac{1}{m\xi^{1/b_{\gamma}}},$

we then have

$O_{p}(S_{0}+{\mathcal{A}}(\xi)+S_{1}+S_{2})=O(f(m)).$

  * (i)

If ${d_{1}}\leq\frac{b_{\gamma}({c_{\gamma}}+1)}{{b_{\gamma}c_{\gamma}}+1}$
then $O(f(m))=O(m^{-\frac{{d_{1}}{c_{\gamma}}}{{c_{\gamma}}+1}})$ with
$\xi=m^{-\frac{{d_{1}}}{{c_{\gamma}}+1}}$;

  * (ii)

If ${d_{1}}>\frac{b_{\gamma}({c_{\gamma}}+1)}{{b_{\gamma}c_{\gamma}}+1}$ then
$O(f(m))=O(m^{-\frac{b_{\gamma}{c_{\gamma}}}{{b_{\gamma}c_{\gamma}}+1}})$ with
$\xi=m^{-\frac{b_{\gamma}}{{b_{\gamma}c_{\gamma}}+1}}$.

###### Proof of Theorem 1.

The choices of $\alpha_{1},\alpha_{2}$ and $n_{1},n_{2}$ in the statement of
Theorem 1 ensure that

$\displaystyle
r_{x}^{2}=O([(n_{1}^{-\frac{1}{2}})^{\frac{2}{c_{1}+1}}]^{2\iota})=O(m^{-{d_{1}}}),\quad
r_{y}^{2}=O([(n_{2}^{-\frac{1}{2}})^{\frac{2}{c_{2}+1}}]^{2})=O(m^{-{d_{2}}}).$

Thus, by Lemma A10, we have

$O_{p}(S_{-1})=O_{p}(1/{\xi}(r_{x}^{2\iota}r_{y}^{2}+r_{y}^{2}+r_{x}^{2\iota}))=O_{p}(1/{\xi}\left\\{m^{-{d_{1}}}+m^{-{d_{2}}}+m^{-{d_{1}}-{d_{2}}}\right\\}).$

Since ${d_{1}},{d_{2}}>0$, and ${d_{1}}\leq{d_{2}}$ by Condition 4,
$m^{-{d_{1}}}/{\xi}$ then dominates two other terms in $S_{-1}$.

Note that $f(m)$ in Lemma A11 also includes $m^{-{d_{1}}}/{\xi}$. Therefore,
given Condition 4, the sum of four terms
$S_{0}+{\mathcal{A}}(\xi)+S_{1}+S_{2}$ dominates $S_{-1}$, which suggests that
the approximation error of $E_{Y}^{p}$ is dominated by that of $E_{X}^{p}$. We
can then derive the result from Lemma A11. ∎

###  A.2 Proof of Theorem 2

###### Proof of Theorem 2.

Under the kernel structural equation model, simple calculation gives

$\displaystyle C=$ $\displaystyle B_{CZ}\Phi(Z)+\epsilon_{C},$ (23)
$\displaystyle\Psi(X)=$
$\displaystyle(B_{XZ}+B_{XC}B_{CZ})\Phi(Z)+B_{XC}\epsilon_{C}+\epsilon_{X},$
(24) $\displaystyle Y=$
$\displaystyle[B_{YZ}+B_{YC}B_{CZ}+B_{YX}(B_{XZ}+B_{XC}B_{CZ})]\Phi(Z)$
$\displaystyle+(B_{YC}+B_{YX}B_{XC})\epsilon_{C}+B_{YX}\epsilon_{X}+\epsilon_{Y}.$
(25)

We denote $B_{\square\triangle}$ as the adjoint operator of
$B_{\triangle\square}$, $B_{\square\triangle}=B_{\triangle\square}^{*}$. When
no ambiguity arise, we use the transpose matrix notation
$B_{\square\triangle}=B_{\triangle\square}^{\top}$. For instance,
$B_{XZ}=B_{ZX}^{\top}$, $B_{YC}=B_{CY}^{\top}$. Recall that the transformed
input and output in Equation 16 and Equation 17 has the form

$\psi_{\gamma}(X)=\psi(X)-E^{p}_{X}\phi(Z)+\sqrt{\gamma}E^{p}_{X}\phi(Z),$

and

$Y_{\gamma}=Y-E^{p}_{Y}\phi(Z)+\sqrt{\gamma}E^{p}_{Y}\phi(Z).$

In the SEM case, the projections $E_{X}^{p}$ and $E_{Y}^{p}$ into $\phi(Z)$
are noted by the (composition of) operators in Equation 24 and Equation 25,
where

$E_{X}^{p}=(B_{XZ}+B_{XC}B_{CZ}),$

and

$E_{Y}^{p}=[B_{YZ}+B_{YC}B_{CZ}+B_{YX}(B_{XZ}+B_{XC}B_{CZ})].$

As such, the transformed input and output has the form

$\psi_{\gamma}(x)=B_{XC}\epsilon_{C}+\epsilon_{X}+\sqrt{\gamma}(B_{XZ}+B_{XC}B_{CZ})\phi(Z),$
(26)

and

$y_{\gamma}=(B_{YC}+B_{YX}B_{XC})\epsilon_{C}+B_{YX}\epsilon_{X}+\epsilon_{Y}+\gamma[B_{YZ}+B_{YC}B_{CZ}+B_{YX}(B_{XZ}+B_{XC}B_{CZ})]\phi(Z).$
(27)

Define relevant covariance matrix/operators as
$\Sigma_{C}={\mathbb{E}}[\epsilon_{C}\epsilon_{C}^{\top}]$,
$\Sigma_{X}={\mathbb{E}}[\epsilon_{X}\otimes\epsilon_{X}]$ and
$\Sigma_{Z}={\mathbb{E}}[\phi(Z)\otimes\phi(Z)]$, where $\otimes$ denotes the
tensor outer product. Then the solution for the least square objective on the
transformed input output can be written as

$H^{\gamma}={\mathbb{E}}[Y_{\gamma}\psi_{\gamma}(X)]({\mathbb{E}}[\psi_{\gamma}(X)\otimes\psi_{\gamma}(X)])^{-1}.$

Plug in the transformed terms in the form of Equation 26 and Equation 27, we
have

$\displaystyle{\mathbb{E}}[\psi_{\gamma}(X)\otimes\psi_{\gamma}(X)]$
$\displaystyle={\mathbb{E}}[(B_{XC}\epsilon_{C}+\epsilon_{X}+\sqrt{\gamma}(B_{XZ}+B_{XC}B_{CZ})\phi(Z))(B_{XC}\epsilon_{C}+\epsilon_{X}$
$\displaystyle\qquad+\sqrt{\gamma}(B_{XZ}+B_{XC}B_{CZ})\phi(Z))^{\top}]$
$\displaystyle=B_{XC}{\mathbb{E}}[\epsilon_{C}\epsilon_{C}^{\top}]B_{CX}+{\mathbb{E}}[\epsilon_{X}\otimes\epsilon_{X}]$
$\displaystyle\quad+\gamma(B_{XZ}+B_{XC}B_{CZ}){\mathbb{E}}[\phi(Z)\otimes\phi(Z)](B_{ZX}+B_{ZC}B_{CX})$
$\displaystyle=B_{XC}\Sigma_{C}B_{CX}+\Sigma_{X}+\gamma(B_{XZ}+B_{XC}B_{CZ})\Sigma_{Z}(B_{ZX}+B_{ZC}B_{CX}).$

Moreover, ${\mathbb{E}}[Y_{\gamma}\psi_{\gamma}(X)]$ has the form

$\displaystyle(B_{YC}+B_{YX}B_{XC}){\mathbb{E}}[\epsilon_{C}\epsilon_{C}^{\top}]B_{CX}+B_{YX}{\mathbb{E}}[\epsilon_{X}\otimes\epsilon_{X}]+$
$\displaystyle\gamma[B_{YZ}+B_{YC}B_{CZ}+B_{YX}(B_{XZ}+B_{XC}B_{CZ})]{\mathbb{E}}[\phi(Z)\otimes\phi(Z)](B_{ZX}+B_{ZC}B_{CX})$
$\displaystyle=$
$\displaystyle(B_{YC}+B_{YX}B_{XC})\Sigma_{C}B_{CX}+B_{YX}\Sigma_{X}+$
$\displaystyle\gamma[B_{YZ}+B_{YC}B_{CZ}+B_{YX}(B_{XZ}+B_{XC}B_{CZ})]\Sigma_{Z}(B_{ZX}+B_{ZC}B_{CX})$

as $\epsilon_{C}$, $\epsilon_{X}$ and $\epsilon_{Y}$ are independent
variables, which are also independent of $Z$. In overall, we have

$\displaystyle H^{\gamma}=$
$\displaystyle[(B_{YC}+B_{YX}B_{XC})\Sigma_{C}B_{CX}+B_{YX}\Sigma_{X}$
$\displaystyle+\gamma[B_{YZ}+B_{YC}B_{CZ}+B_{YX}(B_{XZ}+B_{XC}B_{CZ})]\Sigma_{Z}(B_{ZX}+B_{ZC}B_{CX})]$
$\displaystyle\left[B_{XC}\Sigma_{C}B_{CX}+\Sigma_{X}+\gamma(B_{XZ}+B_{XC}B_{CZ})\Sigma_{Z}(B_{ZX}+B_{ZC}B_{CX})\right]^{-1}$

The bias of the target KAR estimator is then given by

$\displaystyle H^{\gamma}-B_{YX}=$
$\displaystyle\Big{[}(B_{YC}+B_{YX}B_{XC})\Sigma_{C}B_{CX}+B_{YX}\Sigma_{X}$
$\displaystyle+\gamma[B_{YZ}+B_{YC}B_{CZ}+B_{YX}(B_{XZ}+B_{XC}B_{CZ})]\Sigma_{Z}(B_{ZX}+B_{ZC}B_{CX})\Big{]}$
$\displaystyle\Big{[}B_{XC}\Sigma_{C}B_{CX}+\Sigma_{X}+\gamma(B_{XZ}+B_{XC}B_{CZ})\Sigma_{Z}(B_{ZX}+B_{ZC}B_{CX})\Big{]}^{-1}-B_{YX}=$
$\displaystyle\Big{[}(B_{YC}+B_{YX}B_{XC})\Sigma_{C}B_{CX}+B_{YX}\Sigma_{X}$
$\displaystyle+\gamma[B_{YZ}+B_{YC}B_{CZ}+B_{YX}(B_{XZ}+B_{XC}B_{CZ})]\Sigma_{Z}(B_{ZX}+B_{ZC}B_{CX})$
$\displaystyle\qquad-
B_{YX}(B_{XC}\Sigma_{C}B_{CX}+\Sigma_{X}+\gamma(B_{XZ}+B_{XC}B_{CZ})\Sigma_{Z}(B_{ZX}+B_{ZC}B_{CX}))\Big{]}$
$\displaystyle\Big{[}B_{XC}\Sigma_{C}B_{CX}+\Sigma_{X}+\gamma(B_{XZ}+B_{XC}B_{CZ})\Sigma_{Z}(B_{ZX}+B_{ZC}B_{CX})\Big{]}^{-1}$

Collecting all the common terms we get

$\displaystyle H^{\gamma}-B_{YX}=$
$\displaystyle\Big{[}\underset{\Sigma_{YX}^{{\bot}}}{\underbrace{B_{YC}\Sigma_{C}B_{CX}}}+\gamma\underset{\Sigma_{YX}^{{||}}}{\underbrace{(B_{YZ}+B_{YC}B_{CZ})\Sigma_{Z}(B_{ZX}+B_{ZC}B_{CX})}}\Big{]}$
$\displaystyle\Big{[}B_{XC}\Sigma_{C}B_{CX}+\Sigma_{X}+\gamma(B_{XZ}+B_{XC}B_{CZ})\Sigma_{Z}(B_{ZX}+B_{ZC}B_{CX})\Big{]}^{-1}$

Thus, $\forall x\in{\mathcal{X}},y\in{\mathcal{Y}}$, consider the inner
product $y^{\top}(H^{\gamma}-B_{YX})\psi(x)=0$ when the following holds: (i)
$B_{YC}=0$ and $\gamma=0$, or (ii) $B_{YZ}+B_{YC}B_{CZ}=0$ and
$\gamma=\infty$, or (iii) $B_{YC}=0$, $B_{YZ}+B_{YC}B_{CZ}=0$ and $\gamma\geq
0$, or (iv) $\Sigma_{YX}^{||}=a\Sigma_{YX}^{\bot}$ for some $a>0$, and
$\gamma=\infty$, or (v) $\Sigma_{XY}^{||}=-a\Sigma_{XY}^{\bot}$ for some
$a>0$, and $\gamma=1/c$. As such, we conclude $H^{\gamma}=B_{XY}$.

∎

###  A.3 Convergence rate for KAR.2 estimator

In this section, we will further discuss the convergence rate of KAR.2
estimator, and show that the rate does not improve upon the convergence rate
of KAR estimator.

In the three-stage KAR procedure, we approximate $E_{X}^{p}$ and $E_{Y}^{p}$
by $E_{\alpha_{1},X}^{n_{1}}$ and $E_{\alpha_{2},Y}^{n_{2}}$, respectively. In
the two-stage KAR procedure, instead, we approximate the two operators by
$E_{\alpha,X}^{n}$ and $E_{\alpha,Y}^{n}$, respectively. Note that the
estimated operators $E_{\alpha,X}^{n}$ and $E_{\alpha,Y}^{n}$ use the same
$\alpha$. The shared $\alpha$ may fail to ensure the optimal approximation
error for $E_{\alpha,X}^{n}$ and $E_{\alpha,Y}^{n}$ at the same time.

######  Lemma A12.

Under Condition 1, $\forall\delta\in(0,1)$, the following holds w.p.
$1-\delta$:

$\displaystyle\|E_{\alpha,X}^{n}-E_{X}^{p}\|_{{\mathcal{H}}_{\Gamma}}\leq
r_{1}(\alpha):=\frac{4\kappa(Q_{1}+\kappa\|E_{X}^{p}\|_{{\mathcal{H}}_{\Gamma}})\ln(2/\delta)}{\sqrt{n}\alpha}+\alpha^{\frac{c_{1}-1}{2}}\sqrt{\zeta_{1}}.$

Under Condition 1 and Condition 2, $\forall\epsilon\in(0,1)$, the following
holds w.p. $1-\epsilon$:

$\displaystyle\|E_{\alpha,Y}^{n}-E_{Y}^{p}\|_{{\mathcal{H}}_{\Theta}}\leq
r_{2}(\alpha):=\frac{4\kappa(Q_{2}+\kappa\|E_{Y}^{p}\|_{{\mathcal{H}}_{\Theta}})\ln(2/\epsilon)}{\sqrt{n}\alpha}+\alpha^{\frac{c_{2}-1}{2}}\sqrt{\zeta_{2}}.$

Approximation error bound $r_{1}(\alpha)$ for $E_{\alpha,X}^{n}$ achieves its
minimum at rate $O(n^{-\frac{c_{1}-1}{2(c_{1}+1)}})$ when

$\alpha=\left(\frac{8\kappa(Q_{1}+\kappa\|E_{X}^{p}\|_{{\mathcal{H}}_{\Gamma}})\ln(2/\delta)}{\sqrt{n\zeta_{1}}(c_{1}-1)}\right)^{\frac{2}{c_{1}+1}}=O(n^{\frac{-1}{c_{1}+1}});$

and approximation error bound $r_{2}(\alpha)$ for $E_{\alpha,Y}^{n}$ achieves
its minimum at rate $O(n^{-\frac{c_{2}-1}{2(c_{2}+1)}})$ when

$\alpha=\left(\frac{8\kappa(Q_{2}+\kappa\|E_{Y}^{p}\|_{{\mathcal{H}}_{\Theta}})\ln(2/\epsilon)}{\sqrt{n\zeta_{2}}(c_{2}-1)}\right)^{\frac{2}{c_{2}+1}}=O(n^{\frac{-1}{c_{2}+1}}).$

Lemma A12 above provides the upper bounds of the approximation errors for
$E_{\alpha,X}^{n}$ and $E_{\alpha,Y}^{n}$, and it’s adapted from Theorem 2 in
Singh et al., [2019]. We can see that if $c_{1}\neq c_{2}$, we cannot claim
the optimal convergence rate for $E_{\alpha,X}^{n}$ and $E_{\alpha,Y}^{n}$ at
the same time, which disjoint sample sets projection estimators can guarantee
by setting different $\alpha_{1}$ and $\alpha_{2}$ as shown in Lemma 1 and 2.
In other words, in KAR.2 procedure, the error propagated to the final stage,
which are caused by using $E_{\alpha,X}^{n}$ and $E_{\alpha,Y}^{n}$, can have
larger order than using $E_{\alpha_{1},X}^{n_{1}}$ and
$E_{\alpha_{2},Y}^{n_{2}}$ separately in the KAR procedure. Therefore, we
cannot ensure a same or improved convergence rate for KAR.2 estimator compared
to KAR estimator.

##  Appendix B Additional simulation details and results

###  B.1 Synthetic example in KIV setting

In this section, we show the data generating process and implementation
details for the example used in the KIV [Singh et al.,, 2019] that follows the
simulation case of learning counterfactual functions [Chen and Christensen,,
2018] studied in Singh et al., [2019]. The structural model is set as follows,

$Y=C+\ln(|16X-8|+1)sgn(X-0.5).$

The explanatory variables are generated from

$\displaystyle\left(\begin{array}[]{c}C\\\ V\\\ W\end{array}\right)$
$\displaystyle\sim$ $\displaystyle N\left(\left(\begin{array}[]{c}0\\\ 0\\\
0\end{array}\right),\left(\begin{array}[]{ccc}1,0.5,0\\\ 0.5,1,0\\\
0,0,1\end{array}\right)\right),$ $\displaystyle X$ $\displaystyle=$
$\displaystyle F\left(\frac{W+V}{\sqrt{2}}\right),$ $\displaystyle Z$
$\displaystyle=$ $\displaystyle F(W),$

where $F$ denote the c.d.f of standard normal distribution. This structural
model ensures that anchor $Z$ is a valid instrumental variable, so that KIV is
supposed to perform well in this case. We conduct kernel anchor regression
with three-stage algorithm (KAR), kernel anchor regression with two-stage
algorithm (KAR.2) and multiple $\gamma$s and kernel instrument variable
regression (KIV). Set $n_{1}=200$, $n_{2}=200$, $m=600$, $n=n_{1}+n_{2}=400$.
For KAR and KAR.2, we set $\gamma$ to be 0, 0.5, 1, 2, 5, 10, and 100. We set
$\alpha_{1}=c_{\alpha}n_{1}^{-0.5}$, $\alpha_{2}=c_{\alpha}n_{2}^{-0.5}$,
$\alpha=c_{\alpha}n^{-0.5}$, and $\xi=1m^{-0.5}$, where $c_{\alpha}>0$ is a
constant chosen from $\\{0.01,0.05,0.1,0.5,0.8,1,2,3\\}$ for each estimator
separately to minimise the corresponding MSE. We use Gaussian kernel for all
kernel methods, where the lengthscales are set according to median heuristic
[Gretton et al.,, 2012].

For each algorithm, we then implement 50 simulations and calculate MSE with
respect to the true causal model ${\mathbb{E}}(Y|do(x))$, which can be
computed from the structural model. As shown in Figure 2(a), though KIV
performs better than most KAR and KAR.2 estimators, KAR and KAR.2 with
$\gamma=2$ defeat KIV in the KIV setting. The parameters $c_{\alpha}$s are
chosen to be 1, 0.1, 3, 0.8, 3, 3, 3, 1, 0.1, 3, 1, 3, 3, 3 and 2 for KAR with
$\gamma$ being 0, 0.5, 1, 2, 5, 10, 100, KAR.2 with same $\gamma$ series and
KIV, respectively.

Figure B1: Variant synthetic example: fitted nonlinear (left) and linear
(right) methods.

###  B.2 A variant of the synthetic data example

We also consider a variant case where the structural equation is same to the
case in Section 5.1 in the main text

$Y=0.75C-0.25Z+\ln(|16X-8|+1)sgn(X-0.5),$

and the explanatory variables are generated as

$\displaystyle\begin{pmatrix}C\\\ V\\\ W\end{pmatrix}\sim
N\left(\begin{pmatrix}0\\\ 0\\\ 0\end{pmatrix},\begin{pmatrix}1,0.3,0.2\\\
0.3,1,0\\\ 0.2,0,1\end{pmatrix}\right).$

Instead, $X$ and $Z$ are set via the following transformation.

$\displaystyle X=F\left(\frac{\left|W\right|+V}{\sqrt{2}}\right),\quad
Z=F(\left|W\right|)-0.5.$

(a) MSE results in the KIV setting

(b) MSE results of all estimators in the variant case.

Figure B2: Experimental results for additional experiments.

The fitted result of nonlinear and linear methods is shown in Figure B1. The
MSE averaged over 50 simulations is shown in Figure 2(b). From the result, we
can also see that the proposed kernel anchor regression estimators still
performs the best among others under the variant case.