# Fast Instrument Learning with Faster Rates

Ziyu Wang, Yuhao Zhou, Jun Zhu
Dept. of Comp. Sci. and Tech., BNRist Center, State Key Lab for Intell. Tech.
& Sys.,
Institute for AI, Tsinghua-Bosch Joint Center for ML, Tsinghua University
<EMAIL_ADDRESS><EMAIL_ADDRESS>JZ is the corresponding
author. Code available at https://github.com/meta-inf/fil.

###### Abstract

We investigate nonlinear instrumental variable (IV) regression given high-
dimensional instruments. We propose a simple algorithm which combines
kernelized IV methods and an arbitrary, adaptive regression algorithm,
accessed as a black box. Our algorithm enjoys faster-rate convergence and
adapts to the dimensionality of informative latent features, while avoiding an
expensive minimax optimization procedure, which has been necessary to
establish similar guarantees. It further brings the benefit of flexible
machine learning models to quasi-Bayesian uncertainty quantification,
likelihood-based model selection, and model averaging. Simulation studies
demonstrate the competitive performance of our method.

###  1 Introduction

Instrumental variable (IV) analysis is widely used for causal inference [1, 2,
3]. Given confounded observational data, IV analysis identifies the causal
effect through the use of instruments. Nonlinear IV regression is typically
defined by the following conditional moment restrictions (CMRs):

$\mathbb{E}(\mathbf{y}-f_{0}(\mathbf{x})\mid\mathbf{z})=0\quad a.s.~{}[P(dz)]$
(1)

where $f_{0}$ is the causal effect function of interest, and
$\mathbf{x},\mathbf{y},\mathbf{z}$ denote the observed treatment, response and
instrument, respectively. Similar CMR problems also appear in other
applications of causal statistics and machine learning [see, e.g., 4, 5, for
examples].

Starting from [6], recent works have demonstrated great promise in applying
flexible machine learning (ML) methods to IV regression. Modern ML methods are
appealing due to their adaptivity to the informative latent structure in data
[7]: they may adapt to the low dimensionality of the informative latent
features even if the observed input is high-dimensional and its signal-to-
noise ratio is low. Sample complexity gaps have been established in such
settings, between deep models based on neural networks (NNs) or trees, and
linear models such as fixed-form kernels [8, 9, 10]. In IV regression, such
adaptivity will be highly desirable when the observed instruments are high-
dimensional, which is prevalent in applications such as genomics [11], and may
also arise from the general desire to use structured data as instrument.
Following previous work in the parametric setting [12, 13], we refer to this
problem of learning informative latents in instruments as instrument learning.
It generalizes the classical problem of instrument selection [14, 15, 16].

Comparing with standard supervised learning, IV regression is more
challenging, due to the need to estimate a conditional expectation operator
which defines (1). Consequently, establishing adaptivity guarantees becomes
more difficult. While many recent works have demonstrated promising empirical
results using deep models, they are often used as heuristics [e.g., 6], or
justified with crude slow-rate analyses, which establish convergence rates
that saturate at $\Omega(n^{-1/4})$ [e.g., 17]. This is in contrast to faster
rates which approach $n^{-1/2}$ as the regularity of model improves. The only
exception is a minimax formulation of IV estimation [5, 18, 19, 20]. [19]
establish faster rate convergence for this formulation, for models with local
Rademacher complexity bounds. Though local Rademacher analysis covers many
adaptive ML procedures [e.g., 10, 21], it still does not fully explain the
success of modern ML approaches, with prominent alternatives including
implicit regularization [22] and PAC-Bayesian analyses [23]. From a practical
perspective, minimax optimization is computationally expensive, yet it cannot
be avoided in the framework of [19], unless we instantiate their method with
the less flexible kernel models. It requires additional hyperparameter tuning,
which can be challenging in causal problems where validation is indirect and
more difficult. It also prevents the use of such flexible models for
uncertainty quantification, or inference, for which reliable methods have only
been developed for linear nonparametric models [24, 25, 26].

This work bridges the gap between sharp theoretical guarantees and robust,
practical implementation. We assume our prior knowledge about the causal
effect function $f_{0}$ is characterized by a reproducing kernel Hilbert space
(RKHS) $\mathcal{H}$, and focus on the flexibility of conditional expectation
estimation. This is often possible, because the treatment variable is
determined by the problem at hand, and thus has a fixed dimensionality.111Our
method can also be applied when $\mathbf{x}$ contains high-dimensional
exogenous covariates; see Appendix G. We then present a surprisingly simple
algorithm, with faster rate guarantees which in many cases match the best
known in literature. The algorithm defines the conditional expectation
estimates using a learned kernel, the basis of which is defined by applying
adaptive regression algorithms to random draws from a Gaussian process (GP)
prior. Given this learned first-stage kernel and $\mathcal{H}$, we can
estimate $f_{0}$ using kernelized IV methods [19, 20, 27], which have closed-
form solutions and can be efficiently approximated (e.g., with Nyström [19]).
Our method allows easy hyperparameter tuning, and exhibits competitive
performance in simulations. It accesses the regression algorithm as a black
box, thus allowing for the use of any ML methods and benefits from their
established theoretical guarantees. It also enables fast quasi-Bayesian
uncertainty quantification [24, 28, 29] with improved flexibility.

Our algorithm connects to many ideas in literature. Most notably, it can be
viewed as an infinite-dimensional generalization of [12, 13], which consider a
linear outcome model with fixed dimensionality, and use ML methods to learn
the optimal instruments [30, 31]. Our setting requires different analyses, for
defining an infinite-dimensional estimation target and quantifying errors with
its intrinsic complexity. Additionally, analysis of the resultant IV estimator
is complicated by the ill-posedness of infinite dimensional IV models [32]. A
more subtle distinction is in the choice of basis: while for finite-
dimensional function spaces we can pick any set of basis (i.e., features) and
apply the black-box regressor separately, in our case seemingly obvious
choices of basis lead to inferior results (Appendix C.5.3). Our analysis also
connects to the multi-task learning literature, to which we make technical
contributions. Section 6 discusses related work in detail.

The remaining of the paper is organized as follows: Section 2 reviews
background knowledge. Section 3 introduces the instrument learning problem,
and reduces it to a general kernel learning problem. We solve the latter
problem in Section 4, and return to IV in Section 5 with our main results. We
review related work in Section 6, and present numerical experiments in Section
7.

###  2 Notations and Setup

###### Notations

Denote the joint data distribution as $P(dz\times dx\times dy)$, its marginal
distributions as $P(dx),P(dz)$, etc., and their support as
$\mathcal{X},\mathcal{Z}$. For functions of observed variables (e.g.,
$\mathbf{x}$ or $\mathbf{z}$), $\|\cdot\|_{2}$ denotes the $L_{2}$ norm w.r.t.
the respective marginal data distribution. $\|\cdot\|_{\infty}$ denotes the
$L_{\infty}$ norm. We use the notation $[m]:=\\{1,\ldots,m\\}$. Boldface
($\mathbf{x},\mathbf{y},\mathbf{z}$) emphasizes the denotation of random
variables. For any kernel $k$, $\mathcal{GP}(0,k)$ refers to the “standard
Gaussian process” [33] with zero mean, and covariance defined by $k$.
$\lesssim,\gtrsim,\asymp$ represent (in)equalities up to constants; the hidden
constants will not depend on any sample size. $\tilde{\mathcal{O}}(\cdot)$
denotes inequality up to logarithm factors.

###### Problem Setup

Nonparametric IV regression (NPIV) is formulated as (1) [34, 35]. Introduce
the conditional expectation operator $E:L_{2}(P(dx))\to
L_{2}(P(dz)),f\mapsto\mathbb{E}(f(\mathbf{x})\mid\mathbf{z}=\cdot)$, and
define $g_{0}(z):=\mathbb{E}(\mathbf{y}\mid\mathbf{z}=z)$. We can then express
(1) as a linear inverse problem:

$Ef_{0}=g_{0},$ (2)

where we observe $g_{0}$ up to regression error. NPIV deviates from standard
inverse problems in its need to estimate both $f_{0}$ and $E$. Following
conventions in the two stage least square method [1], we refer to the modeling
of $f_{0}$ as the _second stage_ , and that of $E$ – or equivalently, that of
$Ef$ for all $f$ in a hypothesis space – as the _first stage_. The following
assumption describes the setup in full detail:

######  Assumption 2.1 (NPIV).

(i) The data variables $\mathbf{x},\mathbf{y},\mathbf{z}$ satisfy (1), and
$\mathbf{y}$ is bounded by $B$. (ii) We observe two sets of i.i.d. samples,
with matching (marginal) distributions:
$\mathcal{D}_{s2}^{(n_{2})}:=\\{(z_{i},x_{i},y_{i}):i\in[n_{2}]\\}$,
$\mathcal{D}_{s1}^{(n_{1})}:=\\{(\tilde{z}_{i},\tilde{x}_{i}):i\in[n_{1}]\\}$.

We impose (ii) since such additional samples are sometimes available, as
discussed in e.g. [27]. If only $n$ samples from the joint distribution are
available, we can set $n_{1}=n_{2}=n/2$.

In the main text, we assume that our prior knowledge about $f_{0}$ is fully
characterized by an RKHS $\mathcal{H}$, in the following sense.

######  Assumption 2.2 (second stage RKHS).

(i) $\mathcal{X}$ is a bounded subset of $\mathbb{R}^{d_{x}}$; the reproducing
kernel $k_{x}$ of $\mathcal{H}$ is bounded and continuous. (ii) The integral
operator $T_{x}:f\mapsto\int k_{x}(x,\cdot)f(x)P(dx)$ has eigenvalues
$\lambda_{i}(T_{x})\lesssim i^{-(b+1)}$, for some $b>0$. (iii) One of the
following holds true:

  1. (iii).a

(“kernel scheme”): $f_{0}\in\mathcal{H}$.

  2. (iii).b

(“GP scheme”): $b>1$; and for all $n$, $\exists
f^{\dagger}_{n}\in\mathcal{H}$ s.t.
$\|f^{\dagger}_{n}-f_{0}\|_{\mathcal{H}}\lesssim
n^{\frac{1/2}{b+1}},~{}\|f^{\dagger}_{n}-f_{0}\|_{2}\lesssim
n^{-\frac{b/2}{b+1}}.$

In the above, (i) and (ii) are common technical assumptions: (i) ensures the
existence of Mercer’s representation, and (ii) is a complexity measure, with a
larger value of $b$ indicating a smaller hypothesis space. (iii) requires
$\mathcal{H}$ is correctly specified for regression; its two cases cover the
different assumptions in standard RKHS-based estimation and GP modeling.
(iii).a is intuitive. (iii).b is standard in the posterior contraction
literature [36]; it roughly requires $f_{0}$ to be (at least) as regular as
“typical” samples from $\mathcal{GP}(0,k_{x})$, in the sense of [36, Theorem
2.1]. This is different from (iii).a because when $\mathcal{H}$ is infinite
dimensional, almost all GP samples fall out of $\mathcal{H}$ [37, 38]. Our
algorithm applies to both settings, but analysis of the GP scheme requires
additional effort. It is useful as it allows for quasi-Bayesian uncertainty
quantification using a $\mathcal{GP}(0,k_{x})$ prior [26].

The RKHS assumption has been employed in a thread of recent work [27, 20, 39,
26], and generalizes the sieve method in literature [34, 40, 4]. It is most
reasonable when $\mathbf{x}$ has moderate dimensions; for example, when
$f_{0}$ satisfies certain $L_{2}$-Sobolev regularity conditions, we can set
$k_{x}$ to be a suitable Matérn kernel (Example A.1-A.2). Appendix G studies a
more general setting, where $\mathbf{x}$ and $\mathbf{z}$ include additional,
high-dimensional exogenous covariates. Nonetheless, the assumption will be
less reasonable when the treatment variable is high-dimensional and variable
selection is needed for it.

NPIV is typically an ill-posed inverse problem [35]. We now quantify the
degree of ill-posedness:

######  Assumption 2.3.

The operator $E$ is compact, with singular values $s_{i}(E)\asymp i^{-p}$,
where $p>0$.

Such _mildly ill-posed_ settings [41] match our polynomial eigendecay
assumption for the kernel. In the _severely ill-posed_ setting where the decay
of $s_{i}(E)$ is exponential, kernels with a similar eigendecay should be
used. While the analyses of the two settings share many ideas, the Bayesian
inverse problem literature typically restricts to the former for technical
reasons [42, 43].

###  3 From Instrument to Kernel Learning

As we assume $\mathcal{H}$ is a correctly specified second-stage model, it
remains to determine the first stage. In this section, we show that an ideal
first stage model can be defined using another RKHS $\mathcal{I}$, determined
by $\mathcal{H}$ and $E$. Although its kernel $k_{z}$ has an unknown form, we
demonstrate that we can access noisy samples from $\mathcal{GP}(0,k_{z})$,
which, as Section 4 below shows, enable efficient learning of $\mathcal{I}$.
This can be viewed as instrument learning, as $\mathcal{I}$ will only depend
on the informative features in $\mathbf{z}$ (Example 3.1).

Let us first consider the GP scheme (Assumption 2.2 (iii).b) which roughly
requires $f_{0}$ to be similar to typical samples from
$\mathcal{GP}(0,k_{x})$. From a Bayesian perspective, an ideal prior for
$Ef_{0}$ should match the distribution of $Ef$, for
$f\sim\mathcal{GP}(0,k_{x})$. This distribution is “almost equivalent” to
another GP:

######  Lemma 3.1 (proof in Appendix B).

Denote by $[g]_{\sim}$ the $L_{2}$ equivalence class of $g$.222Recall the
$L_{2}$ space is not a function space, and consists of equivalence classes of
functions. Note that for readability, we may occasionally ignore this
distinction in the main text, and use (a version of) $E$ to also denote the
corresponding map between function spaces. All such denotations can be made
unambiguous (Remark B.1), and all null set ambiguities in this section can be
removed under mild additional assumptions (Lemma B.5). Under Assumptions 2.1,
2.2, there exists a kernel $k_{z}$, with integral operator
$T_{z}=ET_{x}E^{\top}$, s.t. for
$f\sim\mathcal{GP}(0,k_{x}),g\sim\mathcal{GP}(0,k_{z})$, $[g]_{\sim}$ has the
same distribution as $E[f]_{\sim}$.

Informally, the lemma shows that $\mathcal{GP}(0,k_{z})$ matches the
distribution of $Ef$. It is thus intuitive that $k_{z}$ could be a good choice
for the first stage. The following lemma further motivates its use in the
kernel scheme: its (i) shows that $\mathcal{I}$ fulfills the conditions in
previous work [19, 27]: the restriction of $E$ on $\mathcal{H}$ has image
contained in $\mathcal{I}$, and is a bounded linear map to $\mathcal{I}$. (ii)
shows that $\mathcal{I}$, as a set of functions, cannot be made smaller while
maintaining (i) .

######  Lemma 3.2 (proof in Appendix B).

Let $\mathcal{I}$ be the RKHS defined by $k_{z}$. Under Assumptions 2.1, 2.2,
(i) for any $f\in\mathcal{H}$, there exists $g\in\mathcal{I}$ s.t.
$[g]_{\sim}=E[f]_{\sim}$; (ii) for any $g\in\mathcal{I}$, there exists
$f\in\mathcal{H}$ satisfying the above. In both cases, we have
$\|f\|_{\mathcal{H}}=\|g\|_{\mathcal{I}}$.

We now demonstrate that $\mathcal{I}$ only depends on the informative latent
features.

######  Example 3.1 (informative latent structure).

Let $\Phi:\mathcal{Z}\to\bar{\mathcal{Z}}$ be a feature extractor that maps
the observed instruments $\mathbf{z}$ to latent features
$\bar{\mathbf{z}}:=\Phi(\mathbf{z})$, s.t.
$\mathbb{E}(f(\mathbf{x})\mid\mathbf{z})=\mathbb{E}(f(\mathbf{x})\mid\Phi(\mathbf{z}))$
for all $L_{2}$-integrable $f$. Then we can apply Lemma 3.1, with $E$ replaced
by $\bar{E}:f\mapsto\mathbb{E}(f(\mathbf{x})\mid\bar{\mathbf{z}})\in
L_{2}(P(d\bar{z}))$, leading to a latent-space RKHS $\bar{\mathcal{I}}$ with
kernel $\bar{k}_{z}$. $\bar{\mathcal{I}}$ induces the input-space RKHS

$\mathcal{I}:=\\{g=\bar{g}\circ\Phi:\bar{g}\in\bar{\mathcal{I}}\\},~{}~{}\|\bar{g}\circ\Phi\|_{\mathcal{I}}:=\|\bar{g}\|_{\bar{\mathcal{I}}};~{}~{}k_{z}(z,z^{\prime})=\bar{k}_{z}(\Phi(z),\Phi(z^{\prime})).$

The above $k_{z}$ satisfies Lemma 3.1-3.2.333There may be multiple kernels
satisfying Lemma 3.1, but they are equivalent up to null sets (Claim B.2); the
ambiguity can be removed under mild assumptions (Lemma B.5). Observe that
$\mathcal{I}$ perfectly approximates $\\{Ef:f\in\mathcal{H}\\}$, but its
complexity only depends on $\bar{\mathcal{I}}$. In particular, $k_{z}$ has the
same Mercer eigenvalues as $\bar{k}_{z}$ (Claim B.2), the decay of which is a
standard complexity measure [e.g., 44, Ch. 7].

While $k_{z}$ has ideal properties, it cannot be used directly as it involves
the unknown operator $E$. Instead, we need to construct an approximation from
data. Our main insight is that we can effectively draw noisy samples from
$\mathcal{GP}(0,k_{z})$; as we develop in Section 4, such samples enable the
approximation of $k_{z}$. To see how the noisy samples are obtained, consider
$f\sim\mathcal{GP}(0,k_{x})$. By Lemma 3.1, $g=Ef$ is $L_{2}$-equivalent to
clean samples from $\mathcal{GP}(0,k_{z})$; and we have
$f(\mathbf{x})=g(\mathbf{z})+(f(\mathbf{x})-(Ef)(\mathbf{z}))$, where the
latter term is unpredictable given $\mathbf{z}$, and from this perspective can
be viewed as noise. Thus, if we apply any regression algorithm to
$f\sim\mathcal{GP}(0,k_{x})$, with $\mathbf{z}$ as input, we will recover a
“denoised” sample from $\mathcal{GP}(0,k_{z})$, up to regression errors.

In the informative latent feature setting, optimal regression error can only
be achieved by methods that adapt to such structures [8, 9, 10]. Approximating
$\mathcal{I}$ with such “denoised” samples can then be viewed as a knowledge
distillation procedure, which results in a compact representation of the
adaptive regression algorithm. This is particularly beneficial in the NPIV
setting: as discussed in the introduction, using a learned kernel eliminates
the need of minimax optimization in estimation, and allows the (indirect) use
of adaptive methods for uncertainty quantification.

###  4 Black-Box Kernel Learning

In this section, we address the problem of kernel learning given noisy GP
samples. As our results apply to more general settings, we first state the
assumptions with full generality.

######  Assumption 4.1 (RKHS).

There exist a continuous function $\Phi:\mathcal{Z}\to\bar{\mathcal{Z}}$, and
a reproducing kernel $\bar{k}_{z}$ over $\bar{\mathcal{Z}}$, s.t. (i) the
random variable $\bar{\mathbf{z}}=\Phi(\mathbf{z})$ is supported on a bounded
subset of $\mathbb{R}^{d_{l}}$; $\bar{k}_{z}$ is bounded. (ii) The eigenvalues
of the integral operator
$T_{\bar{z}}:\bar{g}\mapsto\int\bar{k}_{z}(\bar{z},\cdot)g(\bar{z})P(d\bar{z})$
satisfy $\lambda_{i}(T_{\bar{z}})\lesssim i^{-({\bar{b}}+1)}$, for some
${\bar{b}}>0$. (iii) $\bar{g}\sim\mathcal{GP}(0,\bar{k}_{z})$ have finite sup
norm with probability 1.

The above assumption applies to a latent-space kernel $\bar{k}_{z}$. As shown
in Example 3.1, $\Phi$ and $\bar{k}_{z}$ induce an input-space kernel $k_{z}$,
and RKHS $\mathcal{I}$, which inherit the assumed regularity conditions. Our
goal is to estimate $k_{z}$. This is harder than the estimation of
$\bar{k}_{z}$, as it also involves $\Phi$.

All conditions for $\bar{k}_{z}$ are satisfied by Matérn kernels with a
sutiable order; see Appendix A. Appendix B.1 discusses its applicability in
the IV setting, where $\mathcal{I}$ is defined as in Section 3. Briefly, (ii)
always holds for ${\bar{b}}\geq\max\\{b,2p-1\\}$, and if $\mathcal{H}$ is
further correctly specified in the sense of Assumption E.1, ${\bar{b}}=b+2p$.
(i) and (iii) hold under mild technical assumptions.

We will “denoise” noisy $\mathcal{GP}(0,k_{z})$ samples using a regression
oracle, which is specified below:

######  Assumption 4.2 (regression oracle).

Let $\mathcal{D}^{(n_{1})}:=\\{(\tilde{z}_{i},g(\tilde{z}_{i})+e_{i})\\}$ be
$n_{1}$ iid replications of the rvs $(\mathbf{z},g(\mathbf{z})+\mathbf{e})$,
s.t. $\mathbb{E}(\mathbf{e}\mid\mathbf{z})=0$ and $g(\mathbf{z})+\mathbf{e}$
has a $1$-subgaussian distribution. Then the oracle returns estimator
$\hat{g}_{u,n_{1}}$ s.t.
$\mathbb{E}_{g\sim\mathcal{GP}(0,k_{z})}\mathbb{E}_{\mathcal{D}^{(n_{1})}}\|\hat{g}_{u,n_{1}}-g\|_{2}^{2}\leq\xi_{n_{1}}^{2},$
for some $\xi_{n_{1}}\to 0.$

In the IV setting, we have $g=Ef\sim\mathcal{GP}(0,k_{z})$, and
$\mathbb{E}(\mathbf{e}\mid\mathbf{z})=\mathbb{E}(f(\mathbf{x})-(Ef)(\mathbf{z})\mid\mathbf{z})=0$;
the subgaussian condition is verified by Lemma A.4.

To provide some intuition on adaptivity, we instantiate the assumption with
the DNN model in [10], and compare the resulted $\xi_{n}$ with fixed-form
kernels:

######  Example 4.1 (adaptivity of DNN oracles).

Let $\mathcal{Z}\subset\mathbb{R}^{d_{z}}$,
$\Phi:\mathcal{Z}\to\bar{\mathcal{Z}}$ be $\beta_{1}$-Hölder regular,
$\bar{\mathcal{I}}$ be a Matérn-$\beta_{2}$ RKHS, and $\beta_{1},\beta_{2}\geq
1$. Let the regression oracle return a $\epsilon_{opt}^{2}$-approximate
empirical risk minimizer for the model in [10]. Then for any $\epsilon>0$, it
holds that (see Appendix C.5.1 for derivations)

$\xi_{n}=\tilde{\mathcal{O}}\Big{(}n^{-\frac{\beta_{1}}{2\beta_{1}+d_{z}}}+n^{-\frac{\beta_{2}-\epsilon}{2\beta_{2}+d_{l}}}+\epsilon_{opt}\Big{)}=:\tilde{\mathcal{O}}\Big{(}\epsilon_{fea,n}+n^{-\frac{\beta_{2}-\epsilon}{2\beta_{2}+d_{l}}}+\epsilon_{opt}\Big{)}$
(3)

In the above, $\epsilon_{fea,n}$ characterizes the hardness of feature
learning, i.e., learning $\Phi$. The second term characterizes that of
kernelized regression given the optimal features: it matches the optimal
regression rate if we had full knowledge about $\Phi$, or equivalently,
$\mathcal{I}$, and would be attainable by kernel ridge regression (KRR) using
$\mathcal{I}$.

As long as $\epsilon_{opt}$ is small, (3) will match the minimax rate up to
logarithms. When
$\nicefrac{{\beta_{1}}}{{d_{z}}}<\nicefrac{{\beta_{2}}}{{d_{l}}}$, the minimax
rate is $\epsilon_{fea,n}\gg n^{-\beta_{2}/(2\beta_{2}+d_{l})},$ meaning that
the hardness of feature learning cannot be overlooked. Otherwise, the rate
$\xi_{n}$ nearly matches the rate given full knowledge of the unknown
$\mathcal{I}$, up to the infinitesimal $\epsilon>0$; this is realistic when,
e.g., $d_{z}\gg d_{l}$ and $\Phi$ is linear ($\beta_{1}=\infty$).

We are interested in the high-dimensional regime where $d_{z}\gg d_{l}$. In
this case, fixed-form Matérn or RBF kernels could only attain the rate of
$\mathcal{O}\big{(}n^{-\frac{\min\\{\beta_{1},\beta_{2}\\}}{2\min\\{\beta_{1},\beta_{2}\\}+d_{z}}}\big{)}$,
which can always be much worse than (3), regardless of the hardness of feature
learning. This comparison suggests that fixed-form kernels cannot adapt to the
latent feature structure to avoid the curse of dimensionality.444[10]
establishes formal lower bounds. Also, for small $\beta_{2}$, we can replace
$d_{z}$ with a manifold dimensionality of $\mathcal{Z}$, but it can still be
much larger than $d_{l}$.

We now define the approximate RKHS. Let $\\{g^{(j)}:j\in[m]\\}$ be i.i.d.
samples from the GP prior, and $\hat{g}^{(j)}_{u,n_{1}}$ be the respective
estimate returned by the regression oracle, constructed from the shared
dataset
$\mathcal{D}^{(n_{1})}=\\{(\tilde{z}_{i},g^{(j)}(\tilde{z}_{i})+\epsilon^{(j)}_{i}):i\in[n_{1}],j\in[m]\\}$
where $\epsilon^{(j)}_{i}$ are subgaussian, mean-zero noise. Let
$\hat{g}^{(j)}_{n_{1}}:=\min\\{C_{k}\log n,\hat{g}^{(j)}_{u,n_{1}}(\cdot)\\}$,
where $C_{k}$ is a constant determined by $\mathcal{I}$. Define
$\hat{G}_{n}(z):=(\hat{g}^{(1)}_{n_{1}}(z),\ldots,\hat{g}^{(m)}_{n_{1}}(z))$.
Our approximate RKHS is defined as

$\tilde{\mathcal{I}}:=\big{\\{}g(z)=\theta^{\top}\hat{G}_{n_{1}}(z)\text{ for
some }\theta\in\mathbb{R}^{m}\big{\\}},~{}~{}\text{with norm
}\|g\|_{\tilde{\mathcal{I}}}:=\sqrt{m}\|\theta\|_{2}.$ (4)

As $\tilde{\mathcal{I}}$ is a finite-dimensional linear space, it is an RKHS.
We can check that $\|g\|_{\infty}\leq C_{k}\log
n_{1}\|g\|_{\tilde{\mathcal{I}}}$.

###### Theoretical Results

Under a given model, regression error is decomposed into approximation and
estimation (i.e., generalization) errors. We first present the approximation
error bound:

######  Theorem 4.1 (proof in Appendix C.2).

Under Assumptions 4.1, 4.2, there exists a universal constant $c_{r}>0$, and
an event $E_{n_{1}}$ determined by $g^{(1\ldots m)}$ and
$\mathcal{D}^{(n_{1})}$ with $\mathbb{P}_{\mathcal{D}^{(n_{1})}}E_{n_{1}}\to
1$, on which for any $g^{*}\in L_{2}(P(dz))$, there exists
$\tilde{g}^{*}\in\tilde{\mathcal{I}}$ s.t.

$\displaystyle\|\tilde{g}^{*}\|_{\tilde{\mathcal{I}}}$ $\displaystyle\leq
c_{r}\|\mathrm{Proj}_{m^{\prime}}{g^{*}}\|_{\mathcal{I}},$ (5)
$\displaystyle\|\tilde{g}^{*}-g^{*}\|_{2}$ $\displaystyle\leq
c_{r}\|\mathrm{Proj}_{m^{\prime}}{g^{*}}\|_{\mathcal{I}}(\xi_{n_{1}}+m^{-({\bar{b}}+1)/2})\sqrt{\log
n_{1}}+\|g^{*}-\mathrm{Proj}_{m^{\prime}}{g^{*}}\|_{2},$ (6)

where $m^{\prime}=[m/2]$, and $\mathrm{Proj}_{m^{\prime}}$ denotes the
projection onto the top $m^{\prime}$ Mercer eigenfunctions of $k_{z}$.

We will use $\tilde{\mathcal{I}}$ to estimate functions on a separate dataset
with $n_{2}$ samples. For a single regression task, the estimation error can
be simply bounded as $\tilde{\mathcal{O}}(\sqrt{m/n_{2}})$ [45]. However, our
analysis of IV estimation will require quantifying the intrinsic complexity of
$\tilde{\mathcal{I}}$, which will also allow the use of a larger $m$ in
practice. The following proposition provides one such result; it will be used
in Section 5, to analyze IV estimation in the kernel scheme (Assumption 2.2
(iii).a).

######  Proposition 4.2 (proof in Appendix C.3).

Let $\tilde{\mathcal{I}},\mathcal{D}^{(n_{1})}$ be defined as above, and
$\delta_{n_{2}}$ be the critical radius of the local Rademacher complexity of
the norm ball $\tilde{\mathcal{I}}_{1}$ [46, Ch. 14]. On the event defined in
Theorem 4.1, we have
$\delta_{n_{2}}=\tilde{O}(n_{2}^{-({\bar{b}}+1)/2({\bar{b}}+2)}+m^{-({\bar{b}}+1)/2}+\xi_{n_{1}}).$

IV estimation in the GP scheme is more delicate, and requires additional
analysis of $\tilde{\mathcal{I}}$, which is deferred to App. D. Before we
proceed, however, we illustrate the results on a simple regression task:

######  Example 4.2 (Example 4.1, cont’d).

Let $\Phi,\mathcal{I}$ be defined as before, and $\hat{g}^{(j)}_{n_{1}}$ be
estimated by the DNN oracle. Suppose $\epsilon_{opt}$ is not greater than the
other terms. Then

  1. i.

Let $m=\lceil n_{1}^{{\bar{b}}/({\bar{b}}+1)^{2}}\rceil$. On the event in
Theorem 4.1, for any $g^{*}\in\mathcal{I}$, there exists
$\tilde{g}^{*}\in\tilde{\mathcal{I}}$ s.t.
$\|\tilde{g}^{*}\|_{\tilde{\mathcal{I}}}\leq
c_{r}\|g^{*}\|_{\mathcal{I}},\|\tilde{g}^{*}-g^{*}\|_{2}=\tilde{O}(\|g^{*}\|_{\mathcal{I}}\xi_{n_{1}})$.

  2. ii.

Let $g^{*}\sim\mathcal{GP}(0,k_{z})$. A refined analysis, based on Corollary
C.6, shows that when $m=\lceil n_{1}^{1/({\bar{b}}+1)}\rceil$, there exists
$\tilde{g}$ s.t. $\|\tilde{g}\|_{\tilde{\mathcal{I}}}\lesssim
n_{1}^{1/2({\bar{b}}+1)},\mathbb{E}_{g^{*}\sim\mathcal{GP}(0,k_{z})}\|\tilde{g}-g^{*}\|_{2}=\tilde{\mathcal{O}}(\xi_{n_{1}}).$

(See Appendix C.5.2 for derivations, and another high-probability bound in the
GP scheme.)

Let $\hat{g}^{*}_{n}$ be the truncated OLS estimate using
$\tilde{\mathcal{I}}$, on a dataset
$\\{(z_{i},g^{*}(z_{i})+e_{i}):i\in[n_{2}]\\}$ where $\mathbb{E}(e_{i}\mid
z_{i})=0,\mathrm{Var}(e_{i})\leq 1,\|g^{*}\|_{\infty}\leq B$. Then
$\mathbb{E}\|\hat{g}^{*}_{n}-g^{*}\|_{2}=\tilde{\mathcal{O}}(\|\tilde{g}^{*}-g^{*}\|_{2}+B\sqrt{m/n_{2}})$
[45, Thm. 11.3]. When $n_{1}=n_{2}$, the latter term is $\ll\xi_{n_{2}}$, and
case (ii) above always matches the DNN rate. Case (i) matches the DNN rate
when feature learning becomes harder ($\xi_{n_{1}}^{2}\gtrsim
n_{1}^{-{\bar{b}}/({\bar{b}}+1)}$); otherwise the rate may be slightly
inferior, but still approaches $n_{1}^{-1/2}$ as the regularity ${\bar{b}}$
improves.

As discussed in Example 4.1, when $d_{z}>d_{l}$, the DNN rate can outperform
fixed-form kernels by a large margin. The above example demonstrates a similar
superiority of the learned kernel.

###  5 Results for IV Regression

We shall use the approximate first stage $\tilde{\mathcal{I}}$ for IV
regression, by plugging $\tilde{\mathcal{I}}$ and $\mathcal{H}$ to the
kernelized estimators in [19, 26]; see Algorithm 1. We analyze the resulted
estimators in this section, while deferring implementation details, including
hyperparameter selection, to Appendix F.

Algorithm 1 Kernelized IV with learned instruments.

0: $\mathcal{D}_{s1}^{(n_{1})},\mathcal{D}_{s2}^{(n_{2})}$; regression
algorithm $\mathsf{Regress}$; second-stage kernel $k_{x}$; $m\in\mathbb{N}$

1: for $j\leftarrow 1$ to $m$ do

2: Sample $f^{(j)}\sim\mathcal{GP}(0,k_{x})$

3:
$\hat{g}^{(j)}_{u,n_{1}}\leftarrow\mathsf{Regress}(\\{(\tilde{z}_{i},f^{(j)}(\tilde{x}_{i})):i\in[n_{1}]\\})$

4: end for {the $m$ invocations of $\mathsf{Regress}$ may be replaced with a
single vector-valued regression}

5: Define
$\tilde{k}_{z}(z,z^{\prime}):=\frac{1}{m}\sum_{j=1}^{m}\hat{g}^{(j)}_{n_{1}}(z)\hat{g}^{(j)}_{n_{1}}(z^{\prime}),$where
$\hat{g}^{(j)}_{n_{1}}:=\min\\{\hat{g}^{(j)}_{u,n_{1}}(\cdot),C\log m\\}$.

6: return
$\mathsf{KernelizedIV}(\mathcal{D}_{s2}^{(n_{2})},\tilde{k}_{z},k_{x})$ {See
(7) below, or Appendix F for the closed-form solution}

Both [19] and the posterior mean estimator of [26] have the form

$\arg\min_{f\in\mathcal{H}}\ell_{n_{2}}(f)+\mu\|f\|_{\mathcal{H}}^{2}:=\arg\min_{f\in\mathcal{H}}\max_{g\in\tilde{\mathcal{I}}}\frac{1}{n_{2}}\sum_{i=1}^{n_{2}}(y_{i}-f(x_{i})-\kappa
g(z_{i}))g(z_{i})-\lambda\|g\|_{\tilde{\mathcal{I}}}^{2}+\mu\|f\|_{\mathcal{H}}^{2}.$
(7)

Their difference lies in the regularization scaling, which arises from the
different assumptions about $f_{0}$ and $\mathcal{H}$ (Assumption 2.2). Thus,
we analyze the resulted two estimators separately, in Section 5.1 and Section
5.2 below. In the setting of [26] we are also able to justify the use of
likelihood-based model comparison and (quasi-)Bayesian model averaging (BMA).

####  5.1 Estimation in the Kernel Scheme

[19] establish faster rate convergence of the point estimator under simple
assumptions. We now provide corresponding results using our learned
$\tilde{\mathcal{I}}$, by plugging in the results in Section 4.

######  Proposition 5.1 (proof in Appendix E.1).

Assume Assumptions 2.1, 2.2 (kernel scheme), 4.1 and 4.2. Let
$\tilde{\mathcal{I}}$ be defined by $\tilde{k}_{z}$ in Algorithm 1, and
$\hat{f}_{n_{2}}$ be defined by (7), with $\kappa,\lambda,\nu$ set as in
Appendix E.1. On the event defined in Theorem 4.1, we have

$\\!\|E(\hat{f}_{n_{2}}-f_{0})\|_{2}=\tilde{\mathcal{O}}_{p}\big{(}\big{(}\xi_{n_{1}}+n_{2}^{-\frac{b+1}{2(b+2)}}\big{)}(1+\|f_{0}\|_{\mathcal{H}}^{2})\big{)}.\\!\\!$
(8)

Let us compare the result with [19] in the setting of Example 4.1. Suppose
$n_{1}=n_{2}$. [19] establishes the rate of
$\mathcal{O}_{p}((\xi^{\prime}_{n_{2}}+n_{2}^{-(b+1)/(b+2)})(1+\|f_{0}\|_{\mathcal{H}}^{2}))$,
where $\xi^{\prime}_{n_{2}}$ is comparable with a first-stage regression rate
established from local Rademacher analysis. For the DNN model in Example 4.1,
we have
$\xi^{\prime}_{n}=\tilde{\Theta}(\xi_{n})=\tilde{\mathcal{O}}(\epsilon_{fea,n}+n^{-{\bar{b}}/2({\bar{b}}+1)})$
in the general case,555With some abuse of notation, we also use
$\tilde{\mathcal{O}}$ to hide the infinitesimal deterioration of the
polynomial order.or
$\xi_{n}^{\prime}=\tilde{\mathcal{O}}(\epsilon_{fea,n}+n^{-({\bar{b}}+1)/2({\bar{b}}+2)})$
assuming additional regularity (Remark C.3). Thus, the two rates are
equivalent if $\epsilon_{fea}$ is sufficiently large, meaning that the
difficulty of feature learning cannot be ignored. Otherwise, [19] may be
better if ${\bar{b}}<b+1$; this is a somewhat narrow range, as
${\bar{b}}\geq\max\\{2p-1,b\\}$. Recall that our method is more appealing
computationally: directly instantiating [19] with DNNs requires solving a
minimax problem similar to (7), while for our learned kernel (7) can be
evaluated in closed form.

We can also compare (8) with kernelized IV using a fixed-form first stage. In
the above setting, its best rate is also provided by [19], and is dominated by
the kernel regression error in Ex. 4.1 which, as we discussed, can be much
worse than $\xi_{n}$. Our improved rate has been made possible by the fact
that we are approximating a first-stage model with optimal adaptivity (Ex.
3.1), at a rate that is also adaptive to the informative latent structure (Ex.
4.1).

In summary, our algorithm combines the best of both worlds: it maintains the
sharp guarantees of adaptive models, and the simplicity of kernel methods.

####  5.2 Quasi-Bayesian Estimation and Uncertainty Quantification

Quasi-Bayesian analysis enables efficient uncertainty quantification for NPIV,
without introducing extra risks of model misspecification [24, 28]. [26]
studies a quasi-Bayesian posterior constructed from (7) and a
$\mathcal{GP}(0,k_{x})$ prior. It is defined through the Radon-Nikodym
derivative w.r.t. the prior,
$\big{(}\nicefrac{{d\Pi(\cdot\mid\mathcal{D}_{s2}^{(n_{2})})}}{{d\Pi}}\big{)}(f)\propto
e^{-n_{2}\ell_{n_{2}}(f)}$. For kernel first-stage models, the quasi-posterior
can be evaluated in closed form (App. F). For general models, however, it is
entirely unclear if approximate inference can be possible, since for any
parameter $f$, evaluation of $\ell_{n}(f)$ involves solving a separate
optimization problem. Our kernel learning algorithm enables the (indirect) use
of such models.

Analysis of (quasi-)Bayesian procedures is more challenging, partly because of
the weaker regularization. Thus, [26] introduced additional assumptions. Our
analysis is further complicated by a different assumption on $\mathcal{I}$,
and approximation errors in $\tilde{\mathcal{I}}$, which necessitate further
assumptions. App. E.2 discusses these assumptions in detail. For simplicity,
we state the result in a “rate-optimal” case:666 Classical NPIV lower bounds
continue to hold given full knowledge of $E$ [47], so the rate
$n_{2}^{-\nicefrac{{b}}{{2(b+2p+1)}}}$ is minimax optimal irrespective of
$n_{1}$. In our setting, it is certainly desirable to improve the dependency
on $n_{1}$, and our restriction is only employed to simplify proof. In
simulations we find the choice of $n_{1}=n_{2}$ works well.

######  Theorem 5.2 (posterior contraction; proof in App. E.3).

Assume Asms. 2.1, 2.2 (GP scheme), 2.3, 4.1, 4.2, D.1, D.2, E.1, E.2. Let
$n_{1}$ be s.t. $\xi_{n_{1}}^{2}\log n_{2}+n_{1}^{-(b+2p)/(b+2p+1)}\lesssim
n_{2}^{-1}$, and $m\asymp n_{1}^{1/(b+2p+1)}$. Let
$\Pi_{n_{1}}(\cdot\mid\mathcal{D}_{s2}^{(n_{2})})$ be defined in (52) in
appendix. Then, with $\mathcal{D}_{s1}^{(n_{1})}$-probability $\to 1$, we have

$\displaystyle\mathbb{E}_{\mathcal{D}_{s2}^{(n_{2})}}\Pi_{n_{1}}(\\{f:\|f-f_{0}\|_{2}\geq
M\bar{\epsilon}_{n_{2}}\\}\mid\mathcal{D}_{s2}^{(n_{2})})$ $\displaystyle\to
0,$
$\displaystyle\mathbb{E}_{\mathcal{D}_{s2}^{(n_{2})}}\Pi_{n_{1}}(\\{f:\|E(f-f_{0})\|_{2}\geq
M\bar{\delta}_{n_{2}}\\}\mid\mathcal{D}_{s2}^{(n_{2})})$ $\displaystyle\to 0,$

where $\bar{\delta}_{n_{2}}=\tilde{\mathcal{O}}(n_{2}^{-(b+2p)/2(b+2p+1)})$,
$\bar{\epsilon}_{n_{2}}=\tilde{\mathcal{O}}(n_{2}^{-b/2(b+2p+1)})$.

Theorem 5.2 immediately implies Theorems 5, 6 in [26] for our
$\tilde{\mathcal{I}}$, with the extra logarithms, as their proofs do not
involve the first stage. Those results establish Sobolev norm rates, and
justify uncertainty quantification by lower bounding the magnitude of
posterior spread.

In the nonparametric Bayes literature, contraction results like Theorem 5.2
often lead to the justification of marginal likelihood-based model selection
and averaging. This is also the case here. The key ingredient is the following
marginal quasi-likelihood bound:

######  Corollary 5.3 (proof in Appendix E.4).

In the setting of Theorem 5.2, for some $C>0$ we have

$\displaystyle\mathbb{P}_{\mathcal{D}_{s2}^{(n_{2})}}\big{(}C^{-1}n_{2}^{\frac{1}{b+2p+1}}\log^{-\frac{6}{b}}n_{2}$
$\displaystyle\leq-\log\Pi_{n_{1}}(\mathcal{D}_{s2}^{(n_{2})})\leq
Cn_{2}^{\frac{1}{b+2p+1}}\log^{2}n_{2}\big{)}\to 1.$

This result allows the comparison of a finite number of second-stage RKHSes.
Of particular interest is the comparison between power RKHSes (Defn. A.1),
which often have intuitive interpretations: e.g., for a Matérn RKHS
$\mathcal{H}$ and $\gamma\in(2/(b+1),1)$, the power RKHS
$\mathcal{H}^{\gamma}$ is equivalent to lower-order Matérn RKHSes ([48]; Ex.
A.3). We can verify that such $\mathcal{H}^{\gamma}$ fulfills the assumptions
about $\mathcal{H}$. Thus, provided the other assumptions continue to hold,
Corollary 5.3 will hold for all such $\mathcal{H}^{\gamma}$, with $b+1$
replaced by $\gamma(b+1)$, showing the marginal likelihood has a different
asymptotics. Consequently, it establishes asymptotically valid comparison
between such models, and justifies the use of BMA.

Analysis of more general settings requires additional effort: NPIV is an
inverse problem, and we anticipate the subtleties of model selection in
nonparametric inverse problems. For example, analyses are usually restricted
to the selection of $\gamma$ [43, 49, 50], and the $\gamma>1$ case requires
additional assumptions [49].777We do not cover it here for brevity, noting
that it is well-understood in inverse problem settings [43, 49]. In the IV
setting, it should also be noted that valid model comparison requires a good
approximation to $E|_{\mathcal{H}}$, since otherwise the quasi-likelihood
becomes less meaningful at any finite sample size. The same intuition applies
to other model selection procedures [18, 20, 27] based on the estimated
violation of (1). When the approximation cannot be guaranteed, it could be
preferable to stick to the prior knowledge and fix a conservative choice for
$\mathcal{H}$.

###  6 Related Work

###### Multi-Task Learning

Our Example 4.2 can also be viewed as quantifying sample efficiency
improvements in multi-task learning, if we view the GP prior draws as the
labeling functions for a handful of diverse training tasks, which share the
representation $\Phi$. This general idea is not new: starting from [51, 52], a
line of recent work establishes similar results. Most related is [52, Sec. 5],
which assumes a fixed-dimensional linear model for $\bar{g}$, and an adaptive
$\Phi$ with metric entropy bounds. We assume more general models for both
components, and do not require different training tasks to have separate
inputs. On the flip side, [52] allows for non-iid training tasks. [52, Sec. 6]
investigated infinite-dimensional $\bar{g}$, but established a slow rate. We
are unaware of any work that established fast-rate convergence for infinite-
dimensional top-level models, or used ML models as a black box. Both aspects
may be interesting for multi-task learning, and are necessary for instrument
learning.

###### Causal Statistics

The double machine learning framework [53] also uses black-box ML models to
estimate certain nuisance parameters in the model. While the operator $E$ can
be viewed as a nuisance parameter, the structure of the NPIV problem is quite
different: [54, p. 8] noted that it is very unclear if such a view can be
helpful for NPIV estimation; consistent with their remarks, we have also been
unable to cast our problem into the double ML framework. Note that double ML
has been applied to semiparametric estimation and inference for IV [53, 55,
56, 57], which are orthogonal to our goal.

It has long been known [30, 31] that under a linear outcome model
$f_{0}(\mathbf{x})=\theta^{\top}\mathbf{x}$, using
$\mathbb{E}(\mathbf{x}\mid\mathbf{z})$ as instrument leads to
$\sqrt{n}$-consistent estimates. Our Section 3 can be viewed as an infinite-
dimensional generalization of this observation.888As noted in [27], when
$f_{0}\in\mathcal{H}$ for some RKHS $\mathcal{H}$, the first stage should
model $\mathbb{E}(f(\mathbf{x})\mid\mathbf{z})$ for all $f\in\mathcal{H}$, as
opposed to merely modeling $\mathbb{E}(\mathbf{x}\mid\mathbf{z})$. Note that
[27] did not study the optimal choice of the first stage.  Given high-
dimensional instruments and a parametric outcome model, there is a large body
of literature on efficient inference; see [12] for a review. As we move to
nonparametric models, we focus on estimation which becomes much more
challenging, in the spirit of [54]. Still, we have provided qualitative
characterization for uncertainty estimates in Section 5.2.

For the use of ML for nonlinear IV, [6] studied a heuristic application of
NNs. We discussed the minimax formulation in introduction. [5, 26] justified
the use of NNs with the respective neural tangent kernels (NTKs) which, like
other fixed-form kernels, cannot adapt to the informative latent structure [8,
9]. [39, 58] investigated the combination of an NN-based second stage and a
linear first stage, which could be useful in complementary scenarios. [17]
considered feature learning in both stages, but only established a slow rate;
as the authors noted, it is also unclear if their algorithm reliably minimizes
the empirical risk.

For model selection in the setting of Section 5.2, [59] prove the validity of
bootstrap-based selection for the sieve estimator [34]. [39, 60] investigate
the use of marginal likelihood for two different kernel-based IV estimators:
[60] establish a crude $-\nicefrac{{1}}{{4}}\log n$ upper bound for the log
marginal likelihood, and [39] connect it to the empirical leave-one-out
validation error. Neither result fully justifies model selection as our
Corollary 5.3. For kernelized IV models, [20, 27, 61] proposed validation
statistics for comparing a finite number of first stage models.

Figure 1: Predictive performance: test MSE vs sample size $n_{1}=n_{2}$ for
all method, and $D\in\\{2,100\\}$. Full results are in App. H.2.

###  7 Simulation Study

Our main simulation setup is adapted from [18, 19]; Appendix H.4 presents
additional experiment on the demand dataset [6, 17]. In [18, 19], the observed
$\mathbf{z},\mathbf{x},\mathbf{y}$ are generated by

$\displaystyle\bar{\mathbf{z}}$
$\displaystyle\sim\mathrm{Unif}[-3,3]^{\lfloor\frac{D}{2}\rfloor},~{}\mathbf{z}=h(\bar{\mathbf{z}}),~{}\mathbf{u}\sim\mathcal{N}(0,1),~{}\mathbf{x}:=\bar{\mathbf{z}}_{1}+\mathbf{u}+\mathbf{e}_{x},~{}\mathbf{y}:=({f_{0}(\mathbf{x})+\mathbf{u}+\mathbf{e}_{y}-\mu})/\sigma,$

where $\mathbf{u}$ is the confounder,
$\mathbf{e}_{x},\mathbf{e}_{y}\sim\mathcal{N}(0,0.1^{2})$ are independent
noise, and the constants $\mu,\sigma$ standardize $\mathbf{y}$. We consider
three choices for $h$: (i) $D=2$, $h$ is the identity function; this recovers
the setup in previous work, and quantifies the hardness of the NPIV problem
given true instruments. (ii) $\mathrm{dim}\>\mathbf{z}=D\in\\{40,100\\}$, $h$
is a three-layer DNN; this simulates a feature learning scenario, and ensures
the observation has a low signal-to-noise ratio ($O(1/D)$). (iii) $h$ maps
$\bar{\mathbf{z}}_{1}$ to a MNIST [62] or CIFAR-10 [63] image with matching
label; the MNIST setting also appeared in previous work.

Table 1: Runtime results for all methods in the predictive experiment, for $N=2500,D=100$. Method | AGMM-Tree | AGMM-NN | AGMM-RBF | Proposed
---|---|---|---|---
Runtime / s | $1374\pm 418$ | $303\pm 16$ | $6.7\pm 0.1$ | $25.9\pm 5.6$
Table 2: Test MSE, radius and estimated coverage rate of the $90\%$ $L_{2}$ credible ball (CB), and the average coverage of pointwise $90\%$ credible interval (CI), for $f_{0}\sim\mathcal{GP},D=100$. For the CB coverage rate estimate, we report its $95\%$ Wilson score interval [64]. Full results are in App. H.3.  Method | $n_{1}=n_{2}$ | Test MSE | 90% CB. Rad. | 90% CB. Cvg. | 90% CI. Cvg.
---|---|---|---|---|---
Proposed | $500$ |  $.097$ $\pm.065$ |  $.201$ $\pm.025$ |  $.923$ $[.888,.948]$ |  $.915$ $\pm.123$
$2500$ |  $.035$ $\pm.024$ |  $.074$ $\pm.008$ |  $.917$ $[.880,.943]$ |  $.908$ $\pm.127$
$5000$ |  $.024$ $\pm.016$ |  $.049$ $\pm.004$ |  $.920$ $[.884,.946]$ |  $.905$ $\pm.134$
RBF | $500$ |  $.431$ $\pm.192$ |  $.240$ $\pm.036$ |  $.187$ $[.147,.235]$ |  $.640$ $\pm.191$
$2500$ |  $.176$ $\pm.089$ |  $.175$ $\pm.023$ |  $.517$ $[.460,.573]$ |  $.822$ $\pm.136$
$5000$ |  $.126$ $\pm.072$ |  $.156$ $\pm.019$ |  $.660$ $[.605,.711]$ |  $.855$ $\pm.143$

We consider two choices for $f_{0}$: (i) a widely used collection of functions
(e.g., $\sin,\mathrm{abs}$) in [18]. (ii) $f_{0}\sim\mathcal{GP}(0,k_{x})$.
(ii) ensures the correct specification of $\mathcal{H}$ and allows us to focus
on the first stage.

We use a DNN as the black-box learner, and a RBF kernel for $\mathcal{H}$,
with bandwidth determined by marginal likelihood (F). We set
$N_{1}=N_{2}\in\\{500,2500,5000\\}$. We defer setup details and full results
to Appendix H, and summarize the findings below:

###### Hyperparameter Selection (App. H.1)

We first study hyperparameter selection in instrument learning. We set
$f_{0}\sim\mathcal{GP}$, $D\in\\{2,40,100\\}$. We find our validation
statistics (71) always correlates with the counterfactual MSE
$\|\hat{f}_{n}-f_{0}\|_{2}^{2}$, and that across a large hyperparameter space,
trained DNNs always outperform first-stage models based on RBF kernels, or
randomly initialized DNNs.

###### Predictive Performance (App. H.2)

For $h$ defined as in (i-ii), we compare our algorithm with [19, AGMM],
instantiated with kernel, tree and NN models. As shown in [19], the baselines
have competitive performance on this setup; the latter two models also enjoy
adaptivity guarantees. A representative subset of results are plotted in Fig.
1: our method has stable performance as we move to high dimensions,
demonstrating excellent adaptivity. In contrast, fixed-form kernels fail to
identify the informative features. AGMM-tree and AGMM-NN also have
deteriorated performance as $D$ increases, despite their theoretical
guarantees, presumably due to the challenges in optimization. Table 1 reports
the run time of all methods in this experiment. As we can see, our method is
more efficient than both adaptive baselines.

For image-based $h$, we compare with AGMM-NN and [39], which report the best
results in the MNIST setting. Our method outperforms both baselines.

###### Uncertainty Quantification (App. H.3)

Table 2 presents a subset of results for $f_{0}\sim\mathcal{GP}$. Comparing
with a fixed-form RBF first stage, our method produces sharper credible
intervals, which also have better coverage. For $f_{0}$ specified as in [18],
we experiment with BMA over a grid of RBF kernels, and present visualizations
in Appendix H.3. We find that when the model is more correctly specified, BMA
produces conservative uncertainty estimates which are nonetheless informative.
However, when all models are severely misspecified (e.g., when $f_{0}$ is a
step function), we cannot expect model-based uncertainty estimates to have
ideal coverage.

###### Exogenous Covariates (App. H.4)

We evaluate the extended algorithm in Appendix G on the demand dataset [6],
which is a widely used simulation design with high-dimensional exogenous
covariates. As shown in the appendix, our extended algorithm has competitive
performance.

### Acknowledgements

This work was supported by NSFC Projects (Nos. 62061136001, 62076145,
62076147, U19B2034, U1811461, U19A2081, 61972224), Beijing NSF Project (No.
JQ19016), BNRist (BNR2022RC01006), Tsinghua Institute for Guo Qiang, and the
High Performance Computing Center, Tsinghua University. J.Z is also supported
by the XPlorer Prize.

### References

  * Angrist and Pischke [2008] J. D. Angrist and J.-S. Pischke, _Mostly harmless econometrics_. Princeton university press, 2008.
  * Brookhart et al. [2010] M. A. Brookhart, J. A. Rassen, and S. Schneeweiss, “Instrumental variable methods in comparative safety and effectiveness research,” _Pharmacoepidemiology and drug safety_ , vol. 19, no. 6, pp. 537–554, 2010\.
  * Burgess et al. [2017] S. Burgess, D. S. Small, and S. G. Thompson, “A review of instrumental variable estimators for mendelian randomization,” _Statistical methods in medical research_ , vol. 26, no. 5, pp. 2333–2355, 2017.
  * Chen and Pouzo [2012] X. Chen and D. Pouzo, “Estimation of Nonparametric Conditional Moment Models With Possibly Nonsmooth Generalized Residuals,” _Econometrica_ , vol. 80, no. 1, pp. 277–321, 2012.
  * Liao et al. [2020] L. Liao, Y.-L. Chen, Z. Yang, B. Dai, Z. Wang, and M. Kolar, “Provably Efficient Neural Estimation of Structural Equation Model: An Adversarial Approach,” _arXiv:2007.01290 [cs, stat]_ , 2020, arXiv: 2007.01290.
  * Hartford et al. [2017] J. Hartford, G. Lewis, K. Leyton-Brown, and M. Taddy, “Deep IV: A flexible approach for counterfactual prediction,” in _International Conference on Machine Learning_ , 2017, pp. 1414–1423.
  * Bach [2017] F. Bach, “Breaking the curse of dimensionality with convex neural networks,” _The Journal of Machine Learning Research_ , vol. 18, no. 1, pp. 629–681, 2017.
  * Wei et al. [2019] C. Wei, J. D. Lee, Q. Liu, and T. Ma, “Regularization matters: Generalization and optimization of neural nets v.s. their induced kernel,” in _Advances in Neural Information Processing Systems_ , H. Wallach, H. Larochelle, A. Beygelzimer, F. d’Alché Buc, E. Fox, and R. Garnett, Eds., vol. 32, 2019.
  * Ghorbani et al. [2019] B. Ghorbani, S. Mei, T. Misiakiewicz, and A. Montanari, “Limitations of lazy training of two-layers neural networks,” _arXiv preprint arXiv:1906.08899_ , 2019.
  * Schmidt-Hieber [2020] J. Schmidt-Hieber, “Nonparametric regression using deep neural networks with ReLU activation function,” _The Annals of Statistics_ , vol. 48, no. 4, 2020, arXiv: 1708.06633.
  * Lin et al. [2015] W. Lin, R. Feng, and H. Li, “Regularization Methods for High-Dimensional Instrumental Variables Regression With an Application to Genetical Genomics,” _Journal of the American Statistical Association_ , vol. 110, no. 509, pp. 270–288, 2015.
  * Singh et al. [2020a] A. Singh, K. Hosanagar, and A. Gandhi, “Machine learning instrument variables for causal inference,” in _Proceedings of the 21st ACM Conference on Economics and Computation_ , 2020, pp. 835–836.
  * Chen et al. [2021a] J. Chen, D. L. Chen, and G. Lewis, “Mostly Harmless Machine Learning: Learning Optimal Instruments in Linear IV Models,” _arXiv:2011.06158 [cs, econ, stat]_ , 2021, arXiv: 2011.06158.
  * Okui [2011] R. Okui, “Instrumental variable estimation in the presence of many moment conditions,” _Journal of Econometrics_ , vol. 165, no. 1, pp. 70–86, 2011\.
  * Carrasco [2012] M. Carrasco, “A regularization approach to the many instruments problem,” _Journal of Econometrics_ , vol. 170, no. 2, pp. 383–398, 2012.
  * Belloni et al. [2012] A. Belloni, D. Chen, V. Chernozhukov, and C. Hansen, “Sparse models and methods for optimal instruments with an application to eminent domain,” _Econometrica_ , vol. 80, no. 6, pp. 2369–2429, 2012.
  * Xu et al. [2020] L. Xu, Y. Chen, S. Srinivasan, N. de Freitas, A. Doucet, and A. Gretton, “Learning deep features in instrumental variable regression,” _arXiv preprint arXiv:2010.07154_ , 2020.
  * 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_ , 2019.
  * Dikkala et al. [2020] N. Dikkala, G. Lewis, L. Mackey, and V. Syrgkanis, “Minimax estimation of conditional moment models,” in _Advances in Neural Information Processing Systems_ , vol. 33, 2020, pp. 12 248–12 262.
  * Muandet et al. [2020] K. Muandet, A. Mehrjou, S. K. Lee, and A. Raj, “Dual instrumental variable regression,” _arXiv preprint arXiv:1910.12358_ , 2020.
  * Syrgkanis and Zampetakis [2020] V. Syrgkanis and M. Zampetakis, “Estimation and Inference with Trees and Forests in High Dimensions,” _arXiv:2007.03210 [cs, math, stat]_ , 2020, arXiv: 2007.03210.
  * Neyshabur et al. [2014] B. Neyshabur, R. Tomioka, and N. Srebro, “In search of the real inductive bias: On the role of implicit regularization in deep learning,” _arXiv preprint arXiv:1412.6614_ , 2014.
  * Alquier [2021] P. Alquier, “User-friendly introduction to PAC-Bayes bounds,” _arXiv preprint arXiv:2110.11216_ , 2021.
  * Kato [2013] K. Kato, “Quasi-Bayesian analysis of nonparametric instrumental variables models,” _The Annals of Statistics_ , vol. 41, no. 5, pp. 2359–2390, 2013.
  * Chen and Christensen [2015] X. Chen and T. M. Christensen, “Optimal sup-norm rates and uniform inference on nonlinear functionals of nonparametric IV regression,” _arXiv preprint arXiv:1508.03365_ , 2015.
  * Wang et al. [2021] Z. Wang, Y. Zhou, T. Ren, and J. Zhu, “Quasi-Bayesian dual instrumental variable regression,” _arXiv preprint arXiv:2106.08750_ , 2021.
  * Singh et al. [2019] R. Singh, M. Sahani, and A. Gretton, “Kernel Instrumental variable regression,” _arXiv preprint arXiv:1906.00232_ , 2019.
  * Chernozhukov and Hong [2003] V. Chernozhukov and H. Hong, “An MCMC approach to classical estimation,” _Journal of Econometrics_ , vol. 115, no. 2, pp. 293–346, 2003.
  * Florens and Simoni [2012] J.-P. Florens and A. Simoni, “Nonparametric estimation of an instrumental regression: A quasi-Bayesian approach based on regularized posterior,” _Journal of Econometrics_ , vol. 170, no. 2, pp. 458–475, 2012.
  * Kelejian [1971] H. H. Kelejian, “Two-stage least squares and econometric systems linear in parameters but nonlinear in the endogenous variables,” _Journal of the American Statistical Association_ , vol. 66, no. 334, pp. 373–374, 1971.
  * Chamberlain [1987] G. Chamberlain, “Asymptotic efficiency in estimation with conditional moment restrictions,” _Journal of Econometrics_ , vol. 34, no. 3, pp. 305–334, 1987.
  * Horowitz [2014] J. L. Horowitz, “Ill-Posed Inverse Problems in Economics,” _Annual Review of Economics_ , vol. 6, no. 1, pp. 21–51, 2014.
  * van Zanten and van der Vaart [2008] J. van Zanten and A. van der Vaart, “Reproducing kernel hilbert spaces of gaussian priors,” in _Pushing the limits of contemporary statistics: contributions in honor of Jayanta K. Ghosh_. Institute of Mathematical Statistics, 2008, pp. 200–222.
  * Newey and Powell [2003] W. K. Newey and J. L. Powell, “Instrumental variable estimation of nonparametric models,” _Econometrica_ , vol. 71, no. 5, pp. 1565–1578, 2003\.
  * Horowitz [2011] J. L. Horowitz, “Applied nonparametric instrumental variables estimation,” _Econometrica_ , vol. 79, no. 2, pp. 347–394, 2011.
  * van der Vaart and van Zanten [2008] A. W. van der Vaart and J. H. van Zanten, “Rates of contraction of posterior distributions based on Gaussian process priors,” _The Annals of Statistics_ , vol. 36, no. 3, 2008\.
  * Driscoll [1973] M. F. Driscoll, “The reproducing kernel hilbert space structure of the sample paths of a gaussian process,” _Zeitschrift für Wahrscheinlichkeitstheorie und verwandte Gebiete_ , vol. 26, no. 4, pp. 309–316, 1973.
  * van der Vaart and van Zanten [2011] A. W. van der Vaart and J. H. van Zanten, “Information rates of nonparametric Gaussian process regression,” _Journal of Machine Learning Research_ , vol. 12, pp. 2095–2119, 2011.
  * Zhang et al. [2020] R. Zhang, M. Imaizumi, B. Schölkopf, and K. Muandet, “Maximum moment restriction for instrumental variable regression,” _arXiv preprint arXiv:2010.07684_ , 2020.
  * Blundell et al. [2007] R. Blundell, X. Chen, and D. Kristensen, “Semi-Nonparametric IV Estimation of Shape-Invariant Engel Curves,” _Econometrica_ , vol. 75, no. 6, pp. 1613–1669, 2007, publisher: [Wiley, Econometric Society].
  * Cavalier [2008] L. Cavalier, “Nonparametric statistical inverse problems,” _Inverse Problems_ , vol. 24, no. 3, p. 034004, 2008.
  * Knapik et al. [2011] B. T. Knapik, A. W. van der Vaart, and J. H. van Zanten, “Bayesian inverse problems with Gaussian priors,” _The Annals of Statistics_ , vol. 39, no. 5, pp. 2626–2657, 2011.
  * Knapik et al. [2016] B. T. Knapik, B. T. Szabó, A. W. Van Der Vaart, and J. H. van Zanten, “Bayes procedures for adaptive inference in inverse problems for the white noise model,” _Probability Theory and Related Fields_ , vol. 164, no. 3-4, pp. 771–813, 2016.
  * Steinwart and Christmann [2008] I. Steinwart and A. Christmann, _Support Vector Machines_. Springer Science & Business Media, 2008.
  * Györfi et al. [2002] L. Györfi, M. Kohler, A. Krzyzak, H. Walk _et al._ , _A distribution-free theory of nonparametric regression_. Springer, 2002, vol. 1.
  * Wainwright [2019] M. J. Wainwright, _High-Dimensional Statistics: A Non-Asymptotic Viewpoint_. Cambridge University Press, 2019, vol. 48.
  * Chen and Reiss [2007] X. Chen and M. Reiss, “On Rate Optimality For Ill-posed Inverse Problems In Econometrics,” _Econometric Theory_ , no. 3, pp. 497–521, 2007\.
  * Fischer and Steinwart [2020] S. Fischer and I. Steinwart, “Sobolev norm learning rates for regularized least-squares algorithms,” _Journal of Machine Learning Research_ , vol. 21, pp. 205–1, 2020.
  * Szabó et al. [2015] B. Szabó, A. W. Van Der Vaart, and J. van Zanten, “Frequentist coverage of adaptive nonparametric Bayesian credible sets,” _The Annals of Statistics_ , vol. 43, no. 4, pp. 1391–1428, 2015.
  * Jia et al. [2021] J. Jia, J. Peng, and J. Gao, “Posterior contraction for empirical bayesian approach to inverse problems under non-diagonal assumption,” _Inverse Problems & Imaging_, vol. 15, no. 2, pp. 201–228, 2021.
  * Tripuraneni et al. [2020] N. Tripuraneni, M. I. Jordan, and C. Jin, “On the Theory of Transfer Learning: The Importance of Task Diversity,” _arXiv:2006.11650 [cs, stat]_ , 2020, arXiv: 2006.11650.
  * Du et al. [2021] S. S. Du, W. Hu, S. M. Kakade, J. D. Lee, and Q. Lei, “Few-Shot Learning via Learning the Representation, Provably,” _arXiv:2002.09434 [cs, math, stat]_ , 2021, arXiv: 2002.09434.
  * Chernozhukov et al. [2018] V. Chernozhukov, D. Chetverikov, M. Demirer, E. Duflo, C. Hansen, W. Newey, and J. Robins, “Double/debiased machine learning for treatment and structural parameters,” _The Econometrics Journal_ , vol. 21, no. 1, pp. C1–C68, 2018.
  * Foster and Syrgkanis [2019] D. J. Foster and V. Syrgkanis, “Orthogonal statistical learning,” _arXiv preprint arXiv:1901.09036_ , 2019.
  * Syrgkanis et al. [2019] V. Syrgkanis, V. Lei, M. Oprescu, M. Hei, K. Battocchi, and G. Lewis, “Machine learning estimation of heterogeneous treatment effects with instruments,” _arXiv preprint arXiv:1905.10176_ , 2019.
  * Singh et al. [2020b] R. Singh, L. Xu, and A. Gretton, “Generalized kernel ridge regression for nonparametric structural functions and semiparametric treatment effects,” _arXiv preprint arXiv:2010.04855_ , 2020.
  * Jung et al. [2021] Y. Jung, J. Tian, and E. Bareinboim, “Double machine learning density estimation for local treatment effects with instruments,” in _Advances in Neural Information Processing Systems_ , A. Beygelzimer, Y. Dauphin, P. Liang, and J. W. Vaughan, Eds., 2021.
  * Chen et al. [2021b] J. Chen, X. Chen, and E. Tamer, “Efficient Estimation in NPIV Models: A Comparison of Various Neural Networks-Based Estimators,” _arXiv:2110.06763 [cs, econ, math, stat]_ , 2021, arXiv: 2110.06763.
  * Chen et al. [2021c] X. Chen, T. M. Christensen, and S. Kankanala, “Adaptive estimation and uniform confidence bands for nonparametric iv,” Cowles Foundation for Research in Economics, Yale University, Tech. Rep., 2021.
  * Hsu and Ramos [2019] K. Hsu and F. Ramos, “Bayesian deconditional kernel mean embeddings,” in _International Conference on Machine Learning_. PMLR, 2019, pp. 2830–2838.
  * Zhang et al. [2021] R. Zhang, K. Muandet, B. Schölkopf, and M. Imaizumi, “Instrument Space Selection for Kernel Maximum Moment Restriction,” _arXiv:2106.03340 [cs]_ , 2021, arXiv: 2106.03340.
  * LeCun [1998] Y. LeCun, “The mnist database of handwritten digits,” _http://yann. lecun. com/exdb/mnist/_ , 1998.
  * Krizhevsky et al. [2009] A. Krizhevsky, G. Hinton _et al._ , “Learning multiple layers of features from tiny images,” 2009.
  * Wilson [1927] E. B. Wilson, “Probable inference, the law of succession, and statistical inference,” _Journal of the American Statistical Association_ , vol. 22, no. 158, pp. 209–212, 1927.
  * Steinwart and Scovel [2012] I. Steinwart and C. Scovel, “Mercer’s theorem on general domains: On the interaction between measures, kernels, and RKHSs,” _Constructive Approximation_ , vol. 35, no. 3, pp. 363–417, 2012.
  * Steinwart [2019] I. Steinwart, “Convergence types and rates in generic Karhunen-Loève expansions with applications to sample path properties,” _Potential Analysis_ , vol. 51, no. 3, pp. 361–395, 2019.
  * Kanagawa et al. [2018] M. Kanagawa, P. Hennig, D. Sejdinovic, and B. K. Sriperumbudur, “Gaussian processes and kernel methods: A review on connections and equivalences,” _arXiv preprint arXiv:1807.02582_ , 2018.
  * Steinwart et al. [2009] I. Steinwart, D. R. Hush, C. Scovel _et al._ , “Optimal rates for regularized least squares regression,” in _COLT_ , 2009, pp. 79–93.
  * Adams and Fournier [2003] R. A. Adams and J. J. Fournier, _Sobolev Spaces_. Elsevier, 2003.
  * Eldredge [2016] N. Eldredge, “Analysis and probability on infinite-dimensional spaces,” _arXiv preprint arXiv:1607.03591_ , 2016.
  * Ghosal and Van der Vaart [2017] S. Ghosal and A. Van der Vaart, _Fundamentals of Nonparametric Bayesian Inference_. Cambridge University Press, 2017, vol. 44.
  * Bhatia [2013] R. Bhatia, _Matrix analysis_. Springer Science & Business Media, 2013, vol. 169.
  * Kallenberg and Kallenberg [1997] O. Kallenberg and O. Kallenberg, _Foundations of modern probability_. Springer, 1997, vol. 2.
  * Stuart [2010] A. M. Stuart, “Inverse problems: a Bayesian perspective,” _Acta Numerica_ , vol. 19, pp. 451–559, 2010.
  * Simon [2015] B. Simon, _Operator theory_. American Mathematical Soc., 2015, vol. 4.
  * Edelman [1988] A. Edelman, “Eigenvalues and condition numbers of random matrices,” _SIAM journal on matrix analysis and applications_ , vol. 9, no. 4, pp. 543–560, 1988.
  * Vershynin [2018] R. Vershynin, _High-Dimensional Probability: An Introduction with Applications in Data Science_. Cambridge university press, 2018, vol. 47.
  * Wielandt [1955] H. Wielandt, “An extremum property of sums of eigenvalues,” _Proceedings of the American Mathematical Society_ , vol. 6, no. 1, pp. 106–110, 1955.
  * van der Vaart and van Zanten [2009] A. W. van der Vaart and J. H. van Zanten, “Adaptive Bayesian estimation using a Gaussian random field with inverse Gamma bandwidth,” _The Annals of Statistics_ , vol. 37, no. 5B, 2009.
  * Rudi et al. [2015] A. Rudi, R. Camoriano, and L. Rosasco, “Less is more: Nyström computational regularization.” in _NIPS_ , 2015, pp. 1657–1665.
  * Rudi and Rosasco [2016] A. Rudi and L. Rosasco, “Generalization properties of learning with random features,” _arXiv preprint arXiv:1602.04474_ , 2016.
  * Li and Linde [1999] W. V. Li and W. Linde, “Approximation, metric entropy and small ball estimates for Gaussian measures,” _The Annals of Probability_ , vol. 27, no. 3, 1999\.
  * Pisier [1979] G. Pisier, “Sur les espaces de banach $k$-convexes,” _Séminaire Analyse fonctionnelle (dit" Maurey-Schwartz")_ , pp. 1–15, 1979.
  * Pajor and Tomczak-Jaegermann [1989] A. Pajor and N. Tomczak-Jaegermann, “Volume ratio and other s-numbers of operators related to local properties of Banach spaces,” _Journal of Functional Analysis_ , vol. 87, no. 2, pp. 273–293, 1989.
  * Shawe-Taylor et al. [2002] J. Shawe-Taylor, C. Williams, N. Cristianini, and J. Kandola, “On the eigenspectrum of the gram matrix and its relationship to the operator eigenspectrum,” in _Alt 2002, Lnai 2533_. Springer-Verlag, 2002, pp. 23–40.
  * Bhattacharya et al. [2019] A. Bhattacharya, D. Pati, and Y. Yang, “Bayesian fractional posteriors,” _The Annals of Statistics_ , vol. 47, no. 1, 2019.
  * Hall and Horowitz [2005] P. Hall and J. L. Horowitz, “Nonparametric methods for inference in the presence of instrumental variables,” _The Annals of Statistics_ , vol. 33, no. 6, 2005\.
  * Rousseau and Szabo [2017] J. Rousseau and B. Szabo, “Asymptotic behaviour of the empirical Bayes posteriors associated to maximum marginal likelihood estimator,” _The Annals of Statistics_ , vol. 45, no. 2, 2017.
  * Särkkä [2013] S. Särkkä, “State space representation of gaussian processes,” _Aalto University, Espoo, Finland. Pg_ , pp. 22–23, 2013.
  * Horowitz [2012] J. L. Horowitz, “Specification testing in nonparametric instrumental variable estimation,” _Journal of Econometrics_ , vol. 167, no. 2, pp. 383–396, Apr. 2012.
  * Liang et al. [2018] F. Liang, Q. Li, and L. Zhou, “Bayesian neural networks for selection of drug sensitive genes,” _Journal of the American Statistical Association_ , vol. 113, no. 523, pp. 955–972, 2018, pMID: 31354179.
  * Bauer and Kohler [2019] B. Bauer and M. Kohler, “On deep learning as a remedy for the curse of dimensionality in nonparametric regression,” _The Annals of Statistics_ , vol. 47, no. 4, 2019\.
  * Davies et al. [2015] N. M. Davies, S. von Hinke Kessler Scholder, H. Farbmacher, S. Burgess, F. Windmeijer, and G. D. Smith, “The many weak instruments problem and Mendelian randomization,” _Statistics in Medicine_ , vol. 34, no. 3, pp. 454–468, 2015.

## Appendix

###  Appendix A Background and Technical Lemmas

###### Kernels

The two following lemma applies to our $\mathcal{H}$ satisfying Assumption
2.2, but we will also apply them to the RKHS $\bar{\mathcal{I}}$ defined in
Section 4.999$\mathcal{I}$ may not necessarily satisfy the requirement in
Lemma A.1 (iv), but a weaker version always holds; see Appendix B. In the
latter case, $x,\mathcal{X}$ should be replaced by
$\bar{z},\bar{\mathcal{Z}}$; and as discussed in the main text, the results
will transfer to the RKHS
$\mathcal{I}=\\{\bar{g}\circ\Phi:\bar{g}\in\mathcal{I}\\}$.

######  Lemma A.1 (Mercer’s representation).

Let $\mathcal{H}$ be any RKHS with kernel $k_{x}$ s.t. $\int
P(dx)k_{x}(x,x)<\infty.$ Then

  1. i.

$\mathcal{H}$ can be embedded into $L_{2}(P(dx))$, and the natural inclusion
operator $\iota_{x}:\mathcal{H}\to L_{2}(P(dx))$ and $\iota_{x}^{\top}$ are
Hilbert-Schmidt; the map $T_{x}:f\mapsto\int P(dx)k_{x}(x,\cdot)f(x)$ defines
a positive, self-adjoint and trace-class operator;
$T_{x}=\iota_{x}\iota_{x}^{\top}$.

  2. ii.

$T_{x}$ has the decomposition

$T_{x}f=\sum_{i\in I}\mu_{i}\langle\bar{e}_{i},f\rangle_{2}\bar{e}_{i},$

where the index set $I\subset\mathbb{N}$ is at most countable, and
$\\{\bar{e}_{i}\\}$ is an orthonormal system in $L_{2}(P(dx))$.

  3. iii.

There exists an orthogonal system $\\{e_{i}:i\in I\\}$ of $\mathcal{H}$ s.t.
$[e_{i}]_{\sim}=\sqrt{\lambda_{i}}\bar{e}_{i}$.

  4. iv.

If $k_{x}$ is additionally bounded and continuous, $\\{e_{i}:i\in I\\}$ will
define a Mercer’s representation whose convergence is absolute and uniform.

###### Proof.

[65, Lemma 2.3, 2.2 (for i), 2.12 (for ii-iii), Corollary 3.5 (for iv)]. ∎

The following material on power spaces are adapted from [26], which collected
them from [65, 48].

######  Definition A.1 (power space, embedding property).

Let $\mathcal{H}$ be an RKHS with Mercer’s representation
$\\{(\lambda_{i},\varphi_{i}):i\in\mathbb{N}\\}$. For $\gamma\geq 1$, the
_power space_ $[\mathcal{H}]^{\gamma}\subset L_{2}(P(dx))$ is defined as

$[\mathcal{H}]^{\gamma}=\left\\{[f]_{\sim}:=\sum_{i=1}^{\infty}a_{i}[\varphi_{i}]_{\sim}:\|[f]_{\sim}\|_{[\mathcal{H}]^{\gamma}}^{2}:=\sum_{i=1}^{\infty}\lambda_{i}^{-\gamma}a_{i}^{2}<\infty\right\\}.$

We say $\mathcal{H}$ satisfies an _embedding property_ with order $\gamma$ if
$[\mathcal{H}]^{\gamma}$ is continuously embedded into $L_{\infty}(P(dx))$,
denoted as

$[\mathcal{H}]^{\gamma}\hookrightarrow L_{\infty}(P(dx)).$ (EMB)

Clearly, $\mathcal{I}$ and $\bar{\mathcal{I}}$ will satisfy (EMB) with the
same order.

######  Lemma A.2.

Under (EMB), (i) the function space

$\mathcal{H}^{\gamma}=\left\\{f:=\sum_{i=1}^{\infty}a_{i}\varphi_{i}:(\lambda_{i}^{-\gamma/2}a_{i})_{i\in\mathbb{N}}\in\ell_{2}(\mathbb{N})\right\\},$

with a similarly defined norm, will be an RKHS (a “power RKHS”) with a bounded
kernel; (ii) the kernel has a pointwise convergent Mercer representation
$\\{(\lambda_{i}^{\gamma},\varphi_{i}):i\in\mathbb{N}\\}$. (iii) We have the
interpolation inequality

$\|f\|_{\infty}\lesssim\|f\|_{\mathcal{H}}^{\gamma}\|f\|_{2}^{1-\gamma},\quad\forall
f\in\mathcal{H}.$ (9)

###### Proof.

[66, Theorem 5.5 (for i)], [65, Theorem 3.1 (for ii), Thm. 5.3 (for iii)]. ∎

The embedding property is stronger when $\gamma$ can be chosen to be smaller.
The following example shows that for Matérn kernels, we can choose the best
posible $\gamma$:

######  Example A.1 (Matérn kernels and Sobolev regularity).

Let $\mathcal{X}$ be a bounded open set in $\mathbb{R}^{d}$ with a smooth
boundary, $P(dx)$ have its Lebesgue density bounded from both sides, and
$\mathcal{H}$ be the Matérn-$\alpha$ RKHS. Then

  1. i.

$\mathcal{H}$ is norm-equivalent to the $L_{2}$-Sobolev space
$W^{\alpha+\frac{d}{2},2}$ [67, Example 2.6].

  2. ii.

Its Mercer eigenvalues decay at $\lambda_{i}\asymp
i^{-(1+\frac{2\alpha}{d})}$, and it satisfies (EMB) for all
$\gamma>(1+\frac{2\alpha}{d})^{-1}$ [48, Section 4].

Among kernels with the same eigendecay, this is the best possible $\gamma$
[68].

We now provide some intuition on the “GP scheme” approximation condition,
Assumption 2.2 (iii).b:

######  Example A.2.

For any $\mathcal{H}$ satisfying the eigendecay assumption, simple calculation
shows that any $f_{0}\in[\mathcal{H}]^{b/(b+1)}$ satisfies Assumption 2.2
(iii).b [26, Lemma 23]. If we are further in the setting of Example A.1,
$f_{0}$ will satisfy Assumption 2.2 (iii).b if $f_{0}\in W^{\alpha,2}$ [69,
Chapter 7].

######  Lemma A.3.

Let $\mathcal{I}$ satisfy Assumption 4.1. Then for all
$g\in\mathcal{I},m\in\mathbb{N}$,

$\|\mathrm{Proj}_{m}{g}\|_{\mathcal{I}}\leq\|g\|_{\mathcal{I}},~{}~{}\|g-\mathrm{Proj}_{m}{g}\|_{2}^{2}\lesssim\|g\|_{\mathcal{I}}^{2}m^{-({\bar{b}}+1)},$
(10)

where the constant hidden in $\lesssim$ only depends on $\mathcal{I}$.

###### Proof.

101010Similar result has been stated in [26]. We restate the proof to drop
some unnecessary assumptions.

Let $\\{(\lambda_{i},\varphi_{i})\\}$ be the Mercer eigendecomposition, so
that $\\{\sqrt{\lambda_{i}}\varphi_{i}\\}$ constintute a countable ONB for
$\mathcal{I}$ (Lemma A.1). Thus, the RKHS norm bound holds, and

$\displaystyle\|g-\mathrm{Proj}_{m}{g}\|_{2}^{2}=\sum_{j=m+1}^{\infty}\langle
f,\sqrt{\lambda_{i}}\varphi_{i}\rangle_{\mathcal{I}}^{2}\|\sqrt{\lambda_{i}}\varphi_{i}\|_{2}^{2}\leq\|g-\mathrm{Proj}_{m}{g}\|_{\mathcal{I}}^{2}\cdot\lambda_{m}\|\varphi_{m}\|_{2}^{2}\lesssim\|g\|_{\mathcal{I}}^{2}m^{-({\bar{b}}+1)}.$

∎

###### Gaussian Measure and Gaussian Process

######  Definition A.2 (Gaussian measure, [70]).

Let $(\mathbb{B},\|\cdot\|)$ be a Banach space, $W\sim\mu$ be a Borel
measurable map. $\mu$ is a _Gaussian measure_ if for any
$b^{*}\in\mathbb{B}^{*}$, the pushforward measure $b^{*}_{\\#\mu}$ is normally
distributed.

A Gaussian measure defines a bilinear form on $\mathbb{B}^{*}$:
$q(f,g)=\mathbb{E}f(W)g(W)$. When $\mathbb{B}$ is additionally a Hilbert
space, $q$ will correspond to a bilinear form on $\mathbb{B}$, denoted as
$\Lambda$. We then introduce the notation $W\sim N(0,\Lambda)$, meaning that
for all $l\in\mathbb{B}$, the random variable $\langle
l,W\rangle_{H}\sim\mathcal{N}(0,\langle l,\Lambda l\rangle_{\mathbb{B}})$.

######  Lemma A.4 (Borell-TIS, [71], Proposition I.8).

Let $W$ be any mean-zero Gaussian process defined on a Banach space
$\mathbb{B}$, and $\|\cdot\|$ denote any Banach norm. If $\|W\|<\infty$ a.s.,
it will hold that

$\mathbb{P}(|\|W\|-E\|W\||>x)\leq 2e^{-x^{2}/(2\sigma^{2}(W))}\quad\forall
x>0,$

where
$\sigma(W):=\sup_{b^{*}\in\mathbb{B}^{*}:\|b^{*}\|=1}\sqrt{\mathbb{E}_{W}[b^{*}(W)^{2}]}$
is less than the median of $\|W\|$.

In the following, $C^{\beta}$ denotes the Hölder space of order $\beta$ on
$\mathcal{X}$.

######  Example A.3 (Matérn processes).

Let $k_{x}$ be a Matérn-$\alpha$ kernel, $\mathcal{X}$ be as in Example A.1,
$\underline{\alpha}<\alpha$ be any positive number. Then there exists a
modification of $\mathcal{GP}(0,k_{x})$ which always has finite
$C^{\underline{\alpha}}$ norm [38, p. 2104].

###### Miscellaneous Results

######  Definition A.3 (entropy number).

Let $H,J$ be Banach spaces, $A\subset J$ be a bounded set. For all
$i\in\mathbb{N}$, the $i$-th entropy number is defined as

$e_{i}(A,J)=\inf\\{\epsilon>0:N(A,\|\cdot\|_{J},\epsilon)\leq 2^{i}\\},$

where $N$ denotes the covering number. Further, let $T:H\to J$ be any bounded
linear operator. Then the $i$-th entropy number of the operator $T$ is defined
using the image of the unit-norm ball $H_{1}$ under $T$:

$e_{i}(T)=e_{i}(T(H_{1}),J).$

The following singular value inequality will be frequently used, both to
$s_{j}(AB)$ and the $j$-th largest eigenvalue
$\lambda_{j}(ABB^{\top}A^{\top})=s_{j}(AB)^{2}$:

######  Lemma A.5 (72, Problem III.6.2).

Let $A,B$ be any two operators, $\|\cdot\|$ denote the operator norm, and
$s_{j}$ denote the $j$-th largest singular value. Then

$s_{j}(AB)\leq\min\\{\|B\|s_{j}(A),\|A\|s_{j}(B)\\}.$ (11)

###  Appendix B Deferred Proofs: Function Spaces

######  Remark B.1 (versions of $E$).

Conditional expectations are only defined up to $P(dz)$-null sets. As
$\mathbf{x}$ is supported on a bounded open subset of $\mathbb{R}^{d_{x}}$,
there exists a regular conditional probability $\mu$, which defines a version
of conditional expectation [73]

$\text{for all square integrable
}f,~{}~{}\mathbb{E}(f(\mathbf{x})\mid\mathbf{z})=\int\mu(dx,\mathbf{z})f(x)~{}a.s.~{}[P(dz)].$

Throughout the work, we work with the above version of conditional
expectation.111111The choice of $\mu$ is only unique up to a $P(dz)$-null set;
we fix an arbitrary version to define $E_{r}$. What matters to us is the fact
that $E_{r}$ is defined with a regular conditional probability, so that (12)
always holds.  It represents a linear operator between spaces of functions,
denoted as

$(E_{r}f)(z):=\int\mu(dx,z)f(x).$

As $\mu(\cdot,z)$ is a probability measure for all $z\in\mathcal{Z}$, we now
have

$\|E_{r}f\|_{\infty}\leq\|f\|_{\infty}.$ (12)

Our focus in this work is in estimation; thus, in the main text and other
sections of the appendix, we will abuse notation, and use $E$ to also refer to
$E_{r}$ for readability. In this section, however, we make the distinction
clear for full clarity.

The following claim is well-known. Note that by requiring $[f]_{\sim}$ to be a
Gaussian measure in $L_{2}(P(dx))$, we are requiring our Gaussian process to
possess a possibly richer $\sigma$-algebra than e.g., the version returned by
the Kolmogorov extension theorem. However, they will induce the same marginal
distributions for $f(X)$, and the resulted estimators. See Definition A.2, and
e.g. van Zanten and van der Vaart [33], Stuart [74] for an accessible review
of related issues.

######  Claim B.1.

Let $\mathcal{X}$ be a bounded subset of $\mathbb{R}^{d}$, $k_{x}$ be a
reproducing kernel on $\mathcal{X}$, s.t.
$\mathbb{E}_{P(dx)}k_{x}(x,x)<\infty$. Let $[f]_{\sim}\sim N(0,C)$ be a
Gaussian random element on $L_{2}(P(dx))$ (Definition A.2), s.t. the marginal
distributions of $f$ match $\mathcal{GP}(0,k_{x})$ s.t. $[f]_{\sim}$. Then $C$
equals the integral operator of $k_{x}$.

###### Proof.

Stuart [74, p. 538]. ∎

Note that for our $k_{x}$, the integral operator
$T_{x}=\iota_{x}\iota_{x}^{\top}$ (Lemma A.1).

###### Proof for Lemma 3.1.

By definition of Gaussian measure (A.2) and Claim B.1, we have
$E[f]_{\sim}\sim N(0,ET_{x}E^{\top})$, so it suffices to construct a $k_{z}$
with integral operator $ET_{x}E^{\top}$.

By Lemma A.1 (i) and boundedness of $E:L_{2}(P(dx))\to L_{2}(P(dz))$, the
operator $E\iota_{x}:\mathcal{H}\to L_{2}(P(dz))$ is Hilbert-Schmidt. Thus,
the operator $E\iota_{x}\iota_{x}^{\top}E^{\top}$ is trace-class, and we can
invoke Steinwart and Scovel [65, Theorem 3.10], which shows the existence of
an RKHS $\mathcal{I}$ with a measurable reproducing kernel $k_{0}$, such that

  1. i.

The integral operator of $k_{0}$ equals $E\iota_{x}\iota_{x}^{\top}E^{\top}$.

  2. ii.

For appropriate choices of $e^{z}_{i}$ s.t. $[e^{z}_{i}]_{\sim}$ diagonalizes
$E\iota_{x}\iota_{x}^{\top}E^{\top}$,
$\\{\sqrt{\lambda_{i}^{z}}e_{i}^{z}:i\in\mathbb{N},\lambda_{i}^{z}>0\\}$ form
an ONB of $\mathcal{I}$.

  3. iii.

$k_{0}(z,z^{\prime})=\sum_{i\in
I}\lambda^{z}_{i}e^{z}_{i}(z)e^{z}_{i}(z^{\prime}),~{}~{}\mathbb{E}_{P(dz)}\>k_{0}(z,z)<\infty$.

Combining (i, iii) and Claim B.1 above completes the proof. ∎

Observe the statement (ii) above shows that, the RKHS $\mathcal{I}$ satisfies

$\mathcal{I}=\Big{\\{}\sum_{i\in
I}b_{i}\sqrt{\lambda_{i}^{z}}e_{i}^{z}:~{}(b_{i})\in\ell_{2}(I)\Big{\\}},\quad\Big{\|}\sum_{i\in
I}b_{i}\sqrt{\lambda_{i}^{z}}e_{i}^{z}\Big{\|}_{\mathcal{I}}=\|(b_{i})\|_{\ell_{2}(I)},$
(13)

where $I\subset\mathbb{N}$ denotes an index set which is at most countable.

###### Proof for Lemma 3.2.

The operator $\iota_{x}^{\top}E^{\top}E\iota_{x}$ is also trace-class. Let
$\\{(\lambda^{z}_{i},e_{i}^{z}):i\in I\\}$ be defined as in the proof of Lemma
3.1, so that
$\\{e_{i}^{x}:=(\lambda^{z}_{i})^{-1/2}\iota_{x}^{\top}E^{\top}[e_{i}^{z}]_{\sim}:i\in
I\\}\subset\mathcal{H}$ diagonalizes $\iota_{x}^{\top}E^{\top}E\iota_{x}$.
Then for any $f\in\mathcal{H}$, it holds that

$\displaystyle\infty$ $\displaystyle>\sum_{i\in I}\langle
f,e_{i}^{x}\rangle_{\mathcal{H}}^{2}=\sum_{i\in
I}(\lambda_{i}^{z})^{-2}\langle
f,\iota_{x}^{\top}E^{\top}E\iota_{x}e_{i}^{x}\rangle_{\mathcal{H}}^{2}=\sum_{i\in
I}(\lambda_{i}^{z})^{-1}\langle
E\iota_{x}f,[e_{i}^{z}]_{\sim}\rangle_{2}^{2}.$

Comparing with (13), we can see that for any function $g$ s.t.
$[g]_{\sim}=E\iota_{x}f$, the RHS shows that $g\in\mathcal{I}$, and equals
$\|g\|_{\mathcal{I}}^{2}$.

Conversely, for any $g\in\mathcal{I}$, the sequence
$\\{(\lambda_{i}^{z})^{-1/2}\langle[g]_{\sim},[e_{i}^{z}]_{\sim}\rangle_{2}:i\in
I\\}$ must be in $\ell_{2}$. Additionally, $\\{e_{i}^{x}\\}_{i\in I}$ is an
ONS of $\mathcal{H}$, so the limit

$\sum_{i\in
I}(\lambda_{i}^{z})^{-1/2}\langle[g]_{\sim},[e_{i}^{z}]_{\sim}\rangle_{2}e_{i}^{x}=:f$

must exists in $\mathcal{H}$, and $\|f\|_{\mathcal{H}}=\|g\|_{\mathcal{I}}$
holds by the above display. Similarly, it holds that

$E[f]_{\sim}=\sum_{i\in
I}(\lambda_{i}^{z})^{-1/2}\langle[g]_{\sim},[e_{i}^{z}]_{\sim}\rangle_{2}E\iota_{x}e_{i}^{x}=\sum_{i\in
I}\langle[g]_{\sim},[e_{i}^{z}]_{\sim}\rangle_{2}[e_{i}^{z}]_{\sim}=[g]_{\sim}.$

This completes the proof. ∎

We prove the following claim, made in Example 3.1.

######  Claim B.2.

(i) Let $k_{z},k_{z}^{\prime}$ satisfy Lemma 3.1. Then
$\mathbb{E}_{\mathbf{z},\mathbf{z}^{\prime}\sim
P}(k_{z}(\mathbf{z},\mathbf{z}^{\prime})-k_{z}^{\prime}(\mathbf{z},\mathbf{z}^{\prime}))^{2}=0$.
(ii) Let $\bar{k}_{z},k_{z}$ be defined as in the example. Then the non-zero
Mercer eigenvalues of $\bar{k}_{z},k_{z}$ coincide.

###### Proof.

(i) Both kernels are $L_{2}(P(dz)\otimes P(dz))$ bounded; the claim thus
follows from the isometry between $L_{2}$-bounded kernels and their (Hilbert-
Schmidt) integral operators, and the fact that both kernels have the same
integral operator. (ii) Let $\\{(\lambda_{i},[e_{i}^{\bar{z}}]_{\sim})\\}$
denote the eigendecomposition of the integral operator $T_{\bar{k}}$.
Following the definitions we find that
$\\{[e_{i}^{\bar{z}}\circ\Phi]_{\sim}\\}$ are eigenfunctions of $T_{z}$, with
the same eigenvalues; and it is not possible for $T_{z}$ to have additional
non-zero eigenpairs. ∎

####  B.1 Further Regularity Properties of $\mathcal{I}$

###### Eigendecay

Recall the proof of Lemma 3.1 invokes [65], and leads to the following
results:

######  Claim B.3.

The kernel $k_{z}$ is measurable, satisfies
$\mathbb{E}_{P(dz)}k_{z}(z,z)<\infty$, and has integral operator equal to
$E\iota_{x}\iota_{x}^{\top}E^{\top}$.

The last statement bounds the decay of the Mercer eigenvalues: using
Assumption 2.2, 2.3, and (11), we immediately find
$\lambda^{z}_{i}=\lambda_{i}(E\iota_{x}\iota_{x}^{\top}E^{\top})\lesssim
i^{-\max\\{b+1,2p\\}}.$ We further have the following:

######  Claim B.4.

Under Assumption E.1, it holds that $\lambda^{z}_{i}\lesssim i^{-(b+2p+1)}$.

###### Proof.

Let $\\{[\psi_{i}]_{\sim}:i\in\mathbb{N}\\}$ be the Mercer eigenfunctions of
$\mathcal{H}$, s.t.
$\\{\sqrt{\lambda_{i}(\iota_{x}\iota_{x}^{\top})}\psi_{i}:i\in\mathbb{N}\\}$
form an ONB of $\mathcal{H}$ (Lemma A.1). By the min-max theorem for
eigenvalues [e.g., 75, Theorem 3.2.4],

$\displaystyle\lambda_{i}(E\iota_{x}\iota_{x}^{\top}E^{\top})=\lambda_{i}(\iota_{x}^{\top}E^{\top}E\iota_{x})$
$\displaystyle=\inf_{V\subset\mathcal{H},\mathrm{dim}V=i-1}\sup_{e\perp
V,\|e\|_{\mathcal{H}}=1}e^{\top}\iota_{x}^{\top}E^{\top}E\iota_{x}e$
$\displaystyle\leq\|E\iota_{x}(\sqrt{\lambda_{i}(\iota_{x}\iota_{x}^{\top})}\psi_{i})\|_{2}^{2}\lesssim
i^{-(b+1)}\|E[\psi_{i}]_{\sim}\|_{2}^{2}\lesssim i^{-(b+2p+1)}.$

The last inequality follows by the link condition. ∎

###### Bounded Kernel and GP Prior Draws

To establish boundedness of the kernel $k_{z}$ and (a version of)
$g\sim\mathcal{GP}(0,k_{z})$, we need the following additional assumptions:

  1. (A.I)

A version of $\mathcal{GP}(0,k_{x})$ takes value on a separable subspace
$\mathbb{B}\subset L_{\infty}(P(dx))$.

  2. (A.II)

The operator $E_{r}$ (Remark B.1) maps $\mathbb{B}$ to a space of continuous,
bounded functions on $\mathcal{Z}$.

Note (A.I) will hold given our Assumption E.2, in which case we can take the
subspace as a power RKHS $\mathcal{H}^{\alpha}$; see Steinwart [66]. As
discussed around that assumption, for some valid choice of $\alpha$,
$\mathcal{H}^{\alpha}$ should match the regularity of the second-stage RKHS
assumed in previous work on kernelized IV, so such a boundedness assumption
matches previous work.

(A.II) will hold if we assume $E_{r}$ maps $f\in\mathcal{H}^{\alpha}$ to
another RKHS over $\mathcal{Z}$, with a continuous, bounded reproducing
kernel. Such an RKHS is often assumed in previous work; note that it does not
have have the optimal regularity. Alternatively, the assumption can also be
fulfilled by the assumption that $P(dx\times dz)$ have a continuous Lebesgue
density and the marginal density $p(z)$ does not vanish.

We now establish the following lemma. It shows the $\mathcal{I}$ defined in
Sec. 3 fulfills the conditions in Asm. 4.1. It also shows that by defining
$k_{z}$ with $E_{r}$ as below, we can remove the null set indeterminancies in
Sec. 3: all possible $k_{z}$’s have the same integral operator (i.) and are
thus equivalent up to null sets (Claim B.2 (i)), yet they are shown to be
continuous (iii.).

######  Lemma B.5 (bounded kernel and GP draws).

Let $f$ be a Gaussian measure with marginal distributions matching
$\mathcal{GP}(0,k_{x})$. Let $E_{r}$ be defined in Remark B.1. Then under
(A.I),

  1. i.

There exists a kernel $k_{z}$, s.t. $E_{r}f\sim\mathcal{GP}(0,k_{z})$, with
integral operator equaling $ET_{x}E^{\top}$.

  2. ii.

$k_{z}$ is bounded, and there exists a version of $g\sim\mathcal{GP}(0,k_{z})$
which always has a finite sup norm.

  3. iii.

If additionally (A.II) holds, $k_{z}$ will be continuous, and its RKHS
$\mathcal{I}$ will also satisfy Lemma 3.2.

###### Proof.

ii.: Observe that by (12), for all $z_{0}\in\mathcal{Z}$, the linear
functional $e_{z_{0}}:f\mapsto(E_{r}f)(z_{0})$ is bounded on $L_{\infty}$. As
$k_{x}$ is a bounded kernel, we have
$\|\cdot\|_{\mathcal{H}}\geq(\sup_{x\in\mathcal{X}}k_{x}(x,x))^{-1}\|\cdot\|_{\infty}$;
thus, $e_{z_{0}}$ is bounded on $\mathcal{H}$, and its Riesz representer
$h_{z_{0}}\in\mathcal{H}$ has norm
$\|h_{z_{0}}\|_{\mathcal{H}}\leq\sup_{x\in\mathcal{X}}k_{x}(x,x)=:\sigma_{x}$.
Moreover, for any $\\{z_{1},z_{1},\ldots,z_{m}\\}\subset\mathcal{Z}$ and
$a\in\mathbb{R}^{m}$, the linear map $f\mapsto\sum_{j=1}^{m}a_{j}e_{z_{j}}(f)$
is also bounded on $L_{\infty}$, and thus $\mathcal{H}$; and its representer
$h_{\\{z_{j},a_{j}\\}}$ has $\mathcal{H}$-norm bounded by
$\|a\|_{2}\sigma_{x}$. By our assumptions on $f\sim\mathcal{GP}(0,k_{x})$, we
can invoke Ghosal and Van der Vaart [71, Definition 11.12-11.13, Lemma 11.14]
which show that, for all $m,\\{z_{i}\\}\in\mathcal{Z}^{m},a\in\mathbb{R}^{m}$
and $f\sim\mathcal{GP}(0,k_{x})$,

$\sum_{j=1}^{m}a_{j}(E_{r}f)(z_{j})\sim\mathcal{N}(0,\|h_{\\{z_{j},a_{j}\\}}\|_{\mathcal{H}}^{2}),\quad\text{where}~{}\|h_{\\{z_{j},a_{j}\\}}\|_{\mathcal{H}}\leq\|a\|_{2}\sigma_{x},$

meaning that $E_{r}f$ distributes as a GP. Its reproducing kernel [71,
Definition 11.12] $k_{z}$ satisfies

$\sup_{z\in\mathcal{Z}}k_{z}(z,z)\leq\sigma_{x}.$ (14)

We also have

$\|E_{r}f\|_{\infty}\overset{\eqref{eq:E-bounded-sup-
norm}}{\leq}\|f\|_{\infty}<\infty.$ (15)

As $E_{r}f$ is a version of $\mathcal{GP}(0,k_{z})$, (14) and (15) prove the
second claim.

i.: Let $k_{z}$ be defined as above, $T_{x}$ denote the integral operator of
$k_{x}$. We claim $k_{z}$ has integral operator $ET_{x}E^{\top}$: this is
because by Claim B.1 applied to $f\sim\mathcal{GP}(0,k_{x})$, we have
$[f]_{\sim}\sim N(0,T_{x})$; moreover, we have
$E_{r}f\sim\mathcal{GP}(0,k_{z})$, and $[E_{r}f]_{\sim}=E[f]_{\sim}\sim
N(0,ET_{x}E^{\top})$ by definition of Gaussian measure, and the boundedness of
$E$. Thus, by Claim B.1 applied to $k_{z}$, its integral operator is
$ET_{x}E^{\top}$.

iii.: The continuity of $k_{z}$ follows from (A.II), and the fact that
continuous GP samples must correspond to an RKHS with a continuous kernel [33,
Example 8.1]. Now it remains to re-establish Lemma 3.2.

Following the proofs for Lemma 3.1, 3.2, let
$\\{(\lambda_{i}^{z},[e_{i}^{z}]_{\sim}):i\in I\\}$ be a set of eigenfunctions
for $ET_{x}E^{\top}$. By Lemma A.1, $\\{[e_{i}^{z}]_{\sim}:i\in I\\}$ then
determine an ONS $\\{\sqrt{\lambda_{i}^{z}}e_{i}^{z}:i\in I\\}$ for
$\mathcal{I}$. It suffices to show this is an ONB, after which we can follow
the proof of Lemma 3.2. But as $k_{z}$ is bounded and continuous, the
inclusion operator $\iota_{z}:\mathcal{I}\to L_{2}(P(dz))$ is now injective
[44, Exercise 4.6]; thus, by Steinwart and Scovel [65, Theorem 3.1],
$\\{\sqrt{\lambda_{i}^{z}}e_{i}^{z}:i\in I\\}$ is an ONB. This completes the
proof. ∎

###  Appendix C Deferred Proofs: Kernel Learning

####  C.1 Notations and Preliminary Observations

Let $m^{\prime}=[m/2]$, $\mathrm{Proj}_{m^{\prime}}{(\cdot)}$ denote the
projection onto the top $m^{\prime}$ Mercer basis
$\psi_{1},\ldots,\psi_{m^{\prime}}$, and the respective Mercer eigenvalues be
$\lambda_{i}$. Then there exists i.i.d. normal rvs $\bar{e}_{ij}$ s.t.

$\begin{pmatrix}\mathrm{Proj}_{m^{\prime}}{g^{(1)}}\\\ \ldots,\\\
\mathrm{Proj}_{m^{\prime}}{g^{(m)}}\end{pmatrix}=\begin{pmatrix}\bar{e}_{11}&\ldots&\bar{e}_{1m^{\prime}}\\\
\ldots&\ldots&\ldots\\\
\bar{e}_{m1}&\ldots&\bar{e}_{mm^{\prime}}\end{pmatrix}\begin{pmatrix}\sqrt{\lambda_{1}}\psi_{1}\\\
\ldots\\\ \sqrt{\lambda_{m^{\prime}}}\psi_{m^{\prime}}\end{pmatrix}.$

Denote the $m\times m^{\prime}$ matrix as $\Xi$. Introduce the notation
$\hat{G}:=(\hat{g}^{(1)}_{n_{1}};\ldots;\hat{g}^{(m)}_{n_{1}})$, and
$\bar{G},\Psi$ so that we can write the above as

$\hat{G}+(\bar{G}-\hat{G})=\Xi\Psi.$

Note our slight abuse of notation: throughout the proof, we will use $\hat{G}$
to refer to both the vector-valued function as in the main text, and a “column
of $m$ functions”, i.e., a linear map from $L_{2}(P(dx))$ (or other suitable
function spaces) to $\mathbb{R}^{m}$. We define the norm

$\|\hat{G}\|_{2}^{2}:=\sum_{i=1}^{m}\|\hat{g}^{(i)}_{n_{1}}\|_{2}^{2}=\int\sum_{i=1}^{m}(\hat{g}^{(i)}_{n_{1}}(z))^{2}P(dz),$

and similarly for $\bar{G},\Psi$; and use notations such as $\hat{G}_{j}$ to
refer to the $j$-th row of this “column vector of functions”, so e.g.,
$\hat{G}_{j}$ refers to $\hat{g}^{(j)}_{n_{1}}$.

As $\Xi$ is a $m\times m^{\prime}$ Gaussian random matrix, where
$m^{\prime}=[m/2]$, we have the high-probability singular value bounds

$c^{\prime\prime}\sqrt{m}\leq s_{min}(\Xi)\leq
s_{max}(\Xi)=s_{max}(\Xi^{\top})\leq c^{\prime}\sqrt{m},$ (16)

where $c^{\prime}>c^{\prime\prime}>0$ are universal constants [76]. Thus,

$s_{max}((\Xi^{\top}\Xi)^{-1}\Xi^{\top})\leq(c^{\prime\prime})^{-2}c^{\prime}m^{-1/2}=:c_{r}m^{-1/2}.$
(17)

And for all $j\leq m^{\prime}$,

$s_{j}((\Xi^{\top}\Xi)^{-1}\Xi^{\top})\overset{\eqref{eq:trivial-sv-
ineq}}{\geq}s_{j}(\Xi^{\top})\|\Xi^{\top}\Xi\|^{-1}\geq
s_{min}(\Xi)\|\Xi^{\top}\Xi\|^{-1}\overset{\eqref{eq:singular-value-
raw}}{\geq}c^{\prime\prime}\sqrt{m}\cdot(c^{\prime}\sqrt{m})^{-2}=:c_{r}^{\prime}m^{-1/2}.$
(18)

(By the other inequality in (11), the remaining singular values are zero.)

We will condition on the event defined in (16) throughout all proofs below. On
this event, we have $\Psi=(\Xi^{\top}\Xi)^{-1}\Xi^{\top}\bar{G}$, and we can
define the transformed feature map

$\hat{\Psi}:=(\Xi^{\top}\Xi)^{-1}\Xi^{\top}\hat{G}.$

####  C.2 Proof for Theroem 4.1

Recall the notations and observations in Appendix C.1.

As we will work with the truncated estimators $\hat{g}^{(j)}_{n_{1}}$, we
first show that the truncation does not affect the error rate. By Borell’s
inequality, we have, for any $B>1$

$\mathbb{P}(\max_{i\in[m]}\|g^{(i)}\|_{\infty}\geq B)\leq
m\,\mathbb{P}(\|g^{(i)}\|_{\infty}\geq B)<C_{1}me^{-C_{2}B^{2}}.$ (19)

Thus the above event has high probability for $B=4C_{2}^{-1/2}\sqrt{\log m}$.
On this event, the truncated estimator will have the same $L_{2}$ error as the
original estimators, leading to

$\mathbb{E}_{\mathcal{D}^{(n_{1})},G}\|\hat{G}-G\|_{2}^{2}:=\mathbb{E}_{\mathcal{D}^{(n_{1})},G}\sum_{j=1}^{m}\|\hat{g}^{(j)}_{n_{1}}-g^{(j)}\|_{2}^{2}=m\xi_{n}^{2}.$

And by Markov’s inequality we have, with high probability

$\|\hat{G}-G\|_{2}^{2}\leq m\xi_{n}^{2}\log n.$ (20)

We further restrict to the event on which the above holds. Now, for any
$g^{*}\in L_{2}(P(dx))$, let

$\bar{e}_{*}:=(\lambda_{1}^{-1/2}\langle
g^{*},\psi_{1}\rangle_{2},\ldots,\lambda_{m^{\prime}}^{-1/2}\langle
g^{*},\psi_{m^{\prime}}\rangle_{2})\in\mathbb{R}^{m^{\prime}}$

so that
$\|\bar{e}_{*}\|_{2}=\|\mathrm{Proj}_{m^{\prime}}{g^{*}}\|_{\mathcal{I}},\bar{e}_{*}^{\top}\Psi=\mathrm{Proj}_{m^{\prime}}{g^{*}}$.
Then for $\tilde{g}^{*}:=\bar{e}_{*}^{\top}\hat{\Psi}$, we have

$\displaystyle\|\tilde{g}^{*}\|_{\tilde{\mathcal{I}}}$
$\displaystyle=m^{1/2}\|\Xi(\Xi^{\top}\Xi)^{-1}\bar{e}_{*}\|_{2}\overset{\eqref{eq:singular-
value}}{\leq}c_{r}\|\bar{e}_{*}\|_{2}=c_{r}\|\mathrm{Proj}_{m^{\prime}}{g^{*}}\|_{\mathcal{I}},$
(21) $\displaystyle\|\tilde{g}^{*}-g^{*}\|_{2}$
$\displaystyle\leq\|\bar{e}_{*}^{\top}(\hat{\Psi}-\Psi)\|_{2}+\|g^{*}-\mathrm{Proj}_{m^{\prime}}{g^{*}}\|_{2}$
$\displaystyle\leq
c_{r}\|\mathrm{Proj}_{m^{\prime}}{g^{*}}\|_{\mathcal{I}}(\xi_{n}+m^{-({\bar{b}}+1)/2})\sqrt{\log
n}+\|g^{*}-\mathrm{Proj}_{m^{\prime}}{g^{*}}\|_{2},$ (22)

where (22) holds on a high-probability event independent of $g^{*}$, by Lemma
C.2 which we prove below. ∎

######  Remark C.1.

The failure probability comes from four sources: a Borell’s inequality in
(19), a Markov’s inequality in (20), the singular value bound (17), and a
Markov’s inequality in the above lemma, on $\bar{G}-G$. The failure
probability of (19) is $O(m^{-3})$, which can be easily improved by increasing
$B$. The failure probability of (17) is exponentially small [77, Theorem
4.6.1]. The Markov’s inequality on $\bar{G}-G$ can be sharpened with another
use of Borell-TIS inequality. Therefore, the main source of failure
probability comes from our lack of further assumptions on the black-box
learner, and can be improved given such assumptions.

######  Lemma C.1.

For any $j\in[m]$,

$\mathbb{E}_{g^{(j)}}\|g^{(j)}-\mathrm{Proj}_{m^{\prime}}{g^{(j)}}\|_{2}^{2}=\sum_{j=m^{\prime}+1}^{\infty}\lambda_{j}\asymp
m^{-{\bar{b}}}.$ (23)

###### Proof.

By the $L_{2}$ convergence of K-L expansion, and the eigendecay assumption. ∎

######  Lemma C.2.

On a high-probability event determined by $G$ and $\mathcal{D}^{(n_{1})}$, we
have, for any $\bar{e}_{*}\in\mathbb{R}^{m}$,

$\|\bar{e}_{*}^{\top}(\hat{\Psi}-\Psi)\|_{2}\leq
c_{r}\|\bar{e}_{*}\|_{2}(\xi_{n}+m^{-({\bar{b}}+1)/2})\sqrt{\log n}.$

###### Proof.

Introduce the notation $\tilde{e}:=\Xi(\Xi^{\top}\Xi)^{-1}\bar{e}_{*}$, so
that we can write the LHS as

$\displaystyle\|\bar{e}_{*}^{\top}(\hat{\Psi}-\Psi)\|_{2}$
$\displaystyle=\|\bar{e}_{*}^{\top}(\Xi^{\top}\Xi)^{-1}\Xi^{\top}(\bar{G}-\hat{G})\|_{2}=\|\tilde{e}(\bar{G}-\hat{G})\|_{2}\leq\|\tilde{e}(\bar{G}-G)\|_{2}+\|\tilde{e}(\hat{G}-G)\|_{2}.$

We consider the two terms in turn.

  1. 1.

For the first term, observe

$\mathbb{E}_{G-\bar{G}}\langle(G-\bar{G})_{i},(G-\bar{G})_{j}\rangle=\mathbf{1}_{\\{i=j\\}}\mathbb{E}\|(G-\bar{G})_{1}\|_{2}^{2},~{}~{}\text{where}~{}~{}\mathbb{E}\|(G-\bar{G})_{1}\|_{2}^{2}\overset{\eqref{eq:kl-
trunc-error}}{\asymp}m^{-{\bar{b}}}.$ (24)

Moreover, $\Xi$ and $G-\bar{G}$ depends on disjoint subsets of the generalized
Fourier coefficients of the GP samples, and are thus independent. Thus,
$\tilde{e}$ and $G-\bar{G}$ are also independent, and

$\displaystyle\mathbb{E}_{G-\bar{G}}\|(G-\bar{G})^{\top}\tilde{e}\|_{2}^{2}$
$\displaystyle=\mathbb{E}_{G-\bar{G}}\>\tilde{e}^{\top}(G-\bar{G})(G-\bar{G})^{\top}\tilde{e}\asymp
m^{-{\bar{b}}}\|\tilde{e}\|_{2}^{2}.$

By the Markov inequality we have, w.h.p. w.r.t. $G-\bar{G}$,

$\|(G-\bar{G})^{\top}\tilde{e}\|_{2}^{2}\leq m^{-{\bar{b}}}\log
n\|\tilde{e}\|_{2}^{2}.$ (25)

  2. 2.

For the second term we have

$\displaystyle\|\tilde{e}^{\top}(G-\hat{G})\|_{2}^{2}$
$\displaystyle=\int(\tilde{e}^{\top}(G-\hat{G})(z))^{2}P(dz)\leq\int\|\tilde{e}\|_{2}^{2}\|(G-\hat{G})(z)\|^{2}P(dz)$
$\displaystyle=\|\tilde{e}\|_{2}^{2}\|G-\hat{G}\|_{2}^{2}\overset{\eqref{eq:G-mean-
bound}}{\leq}\|\tilde{e}\|_{2}^{2}m\xi_{n}^{2}\log n.$ (26)

In the above, the first inequality is the Cauchy-Schwarz inequality in
$\mathbb{R}^{m}$.

Combining these results and the observation that

$\|\tilde{e}\|_{2}=\|\Xi(\Xi^{\top}\Xi)^{-1}\bar{e}_{*}\|_{2}\overset{\eqref{eq:singular-
value}}{\leq}c_{r}m^{-1/2}\|\bar{e}_{*}\|_{2},$ (27)

we have, with high probability w.r.t. $G$ and $\mathcal{D}^{(n_{1})}$,

$\|\bar{e}_{*}^{\top}(\hat{\Psi}-\Psi)\|_{2}^{2}\leq\|\tilde{e}\|_{2}^{2}(m\xi_{n}^{2}+m^{-{\bar{b}}})\log
n\leq c_{r}^{2}\|\bar{e}_{*}\|_{2}^{2}(\xi_{n}^{2}+m^{-({\bar{b}}+1)})\log n.$

∎

####  C.3 Proof for Proposition 4.2

Recall the notations and observations in Appendix C.1. In the following, the
first two lemmas bound the entropy number of $\tilde{\mathcal{I}}_{1}$, and
the final proof uses it to bound the local Rademacher complexity.

######  Lemma C.3.

There exists an RKHS $\bar{\mathcal{I}}$ s.t.
$\tilde{\mathcal{I}}_{1}\subset\bar{\mathcal{I}}_{1}$, and, on the event in
Theorem 4.1, we have, for all $0\leq i_{s}\leq i_{e}$, the Mercer eigenvalue
bound

$\sum_{i=i_{s}}^{i_{e}}\lambda_{i}(\bar{C})\lesssim
i_{s}^{-\bar{b}}+\chi_{n}^{2},$ (28)

where $\chi_{n}:=(m^{-\frac{{\bar{b}}+1}{2}}+\xi_{n})\sqrt{\log n},$ and the
constant hidden in $\lesssim$ is independent of $j,m$ or $n$.

###### Proof.

We first define $\bar{\mathcal{I}}$. Define
$\Delta_{1}G_{n}=\hat{G}_{n}-G_{n},\Delta_{2}G_{n}=G_{n}-\bar{G}_{n}.$ Then

$\displaystyle\tilde{\mathcal{I}}_{1}$
$\displaystyle=\\{\theta^{\top}(\bar{G}_{n}+\Delta_{1}G_{n}+\Delta_{2}G_{n})(\cdot):\|\theta\|_{2}^{2}\leq
m^{-1}\\}$
$\displaystyle\subset\\{\theta_{1}^{\top}\bar{G}_{n}(\cdot)+\theta_{2}^{\top}\Delta_{1}G_{n}(\cdot)+\theta_{3}^{\top}\Delta_{2}G_{n}(\cdot):\|\theta_{1}\|_{2}^{2}+\|\theta_{2}\|_{2}^{2}+\|\theta_{3}\|_{2}^{2}\leq
3m^{-1}\\}.$

The RHS is the unit-norm ball of an RKHS, denoted as $\bar{\mathcal{I}}$. Its
reproducing kernel and integral operator are (recall our notations in Appendix
C.1)

$\displaystyle\bar{k}(z,z^{\prime})$
$\displaystyle=\frac{3}{m}(\bar{G}_{n}(z)^{\top}\bar{G}_{n}(z^{\prime})+\Delta_{1}G_{n}(z)^{\top}\Delta_{1}G_{n}(z^{\prime})+\Delta_{2}G_{n}(z)^{\top}\Delta_{2}G_{n}(z^{\prime})),$
$\displaystyle\bar{C}$
$\displaystyle=\frac{3}{m}(\bar{G}_{n}^{\top}\bar{G}_{n}+(\Delta_{1}G_{n})^{\top}\Delta_{1}G_{n}+(\Delta_{2}G_{n})^{\top}\Delta_{2}G_{n}).$

By Wielandt [78, Theorem 2], we have, for any $0\leq i_{s}\leq i_{e}$,

$\sum_{i=i_{s}}^{i_{e}}\lambda_{i}(\bar{C})\leq\frac{3}{m}\left(\sum_{i=i_{s}}^{i_{e}}\lambda_{i}(\bar{G}_{n}^{\top}\bar{G}_{n})+\operatorname{Tr}((\Delta_{1}G_{n})^{\top}\Delta_{1}G_{n})+\operatorname{Tr}((\Delta_{2}G_{n})^{\top}\Delta_{2}G_{n})\right).$

For the first part above, recall $\bar{G}_{n}=\Xi\Psi$; and on the event
defined in Theorem 4.1, we have, by (16),
$(c^{\prime\prime})^{2}m\leq\lambda_{i}(\Xi^{\top}\Xi)\leq(c^{\prime})^{2}m.$
Thus,

$\lambda_{i}(\bar{G}_{n}^{\top}\bar{G}_{n})=\lambda_{i}(\Psi^{\top}\Xi^{\top}\Xi\Psi)\overset{\eqref{eq:trivial-
sv-ineq}}{\asymp}m\lambda_{i}(\Psi^{\top}\Psi)\lesssim mi^{-({\bar{b}}+1)},$

The second term is bounded as

$\operatorname{Tr}((\Delta_{1}G_{n})^{\top}(\Delta_{1}G_{n}))=\|\Delta_{1}G_{n}\|_{2}^{2}\overset{\eqref{eq:G-mean-
bound}}{\lesssim}m\xi_{n}^{2}\log n.$

For the third, by Markov’s inequality on (24), we have, w.h.p. w.r.t.
$G_{n}-\bar{G}_{n}$,

$\operatorname{Tr}((\Delta_{2}G_{n})^{\top}(\Delta_{2}G_{n}))\lesssim
m^{-{\bar{b}}}\log n.$

Note this event has been included in Theorem 4.1 via Lemma C.2. Combining the
above, we have, for all $0\leq i_{s}\leq i_{e}$,

$\sum_{i=i_{s}}^{i_{e}}\lambda_{i}(\bar{C})\lesssim\left(\sum_{i=i_{s}}^{i_{e}}i^{-({\bar{b}}+1)}\right)+(\xi_{n}^{2}+m^{-({\bar{b}}+1)})\log
n\lesssim i_{s}^{-{\bar{b}}}+(\xi_{n}^{2}+m^{-({\bar{b}}+1)})\log n\leq
i_{s}^{-{\bar{b}}}+\chi_{n}^{2}.$

∎

######  Lemma C.4.

Let $Z^{n_{2}}:=\\{z_{1},\ldots,z_{n_{2}}\\}$ be a sample of $n_{2}$ iid
inputs independent of $\mathcal{D}^{(n_{1})}$, and $L_{2}(Z^{n_{2}})$ denote
the $L_{2}$ space defined by the respective empirical measure. Then, on the
high-probability event defined in Theorem 4.1, we have

$\mathbb{E}_{Z^{n_{2}}}e_{j}(\mathrm{id}:\tilde{\mathcal{I}}\to
L_{2}(Z^{n_{2}}))\lesssim
j^{-({\bar{b}}+1)}(\min\\{j,n_{2}\\}^{\frac{{\bar{b}}+1}{2}}+\min\\{j,n_{2}\\}^{{\bar{b}}+\frac{1}{2}}\chi_{n})),$

where $\chi_{n}:=(m^{-\frac{{\bar{b}}+1}{2}}+\xi_{n})\sqrt{\log n},$ and the
constant hidden in $\lesssim$ is independent of $j,m$ or $n$.

###### Proof.

As $\tilde{\mathcal{I}}_{1}\subset\bar{\mathcal{I}}_{1}$, it suffices to
establish the bound for $\bar{\mathcal{I}}_{1}$. For this, we first invoke
Steinwart and Christmann [44, Theorem 7.30], but with the RHS of the last
display on p. 275 replaced by (28). Following the proof we find, for all
$p\in(0,1)$,

$\mathbb{E}_{Z^{n_{2}}}e_{j}(S_{\bar{k},D}^{*})\lesssim
j^{-1/p}\sum_{i=1}^{\min\\{j,n_{2}\\}}i^{1/p-1}\sqrt{i^{-1}(i^{-{\bar{b}}}+\chi_{n}^{2})}.$

As stated in their Corollary 7.31, the LHS above is our desired empirical
entropy number. Following the proof for that corollary, we set
$p=({\bar{b}}+1)^{-1}$ above, leading to

$\displaystyle\mathbb{E}_{Z^{n_{2}}}e_{j}(\tilde{\mathcal{I}}_{1}\to
L_{2}(Z^{n_{2}}))$ $\displaystyle\lesssim
j^{-({\bar{b}}+1)}\sum_{i=1}^{\min\\{j,n_{2}\\}}i^{{\bar{b}}}(i^{-({\bar{b}}+1)/2}+i^{-1/2}\chi_{n})$
$\displaystyle\lesssim
j^{-({\bar{b}}+1)}(\min\\{j,n_{2}\\}^{\frac{{\bar{b}}+1}{2}}+\min\\{j,n_{2}\\}^{{\bar{b}}+\frac{1}{2}}\chi_{n}).$

∎

We now prove Proposition 4.2, by plugging in our new entropy number bound to
the chaining argument in Steinwart and Christmann [44, Theorem 7.12 – Theorem
7.16].

######  Proposition C.5 (Proposition 4.2, restated).

Let

$\bar{R}_{n_{2}}(\tilde{\mathcal{I}}_{1};\sigma):=\mathbb{E}_{Z^{n_{2}},\varepsilon_{i}}\sup_{g\in\tilde{\mathcal{I}}_{1},\|g\|_{2}\leq\sigma}\left|\frac{1}{n_{2}}\sum_{i=1}^{n_{2}}\varepsilon_{i}g(z_{i})\right|$

be the local Rademacher complexity, and
$B:=\sup_{g\in\tilde{\mathcal{I}}_{1}}\|g\|_{\infty}$. Then its critical
radius, defined as the solution $\delta_{n_{2}}$ to

$\bar{R}_{n_{2}}(\tilde{\mathcal{I}}_{1};\delta)\leq\delta^{2},$ (29)

is bounded as

$\delta_{n_{2}}^{2}\lesssim
B^{\frac{{\bar{b}}}{{\bar{b}}+2}}n_{2}^{-\frac{{\bar{b}}+1}{{\bar{b}}+2}}+\frac{\log^{2}n_{2}}{n_{2}}+\chi_{n}^{2}.$

###### Proof.

Define $I_{\sigma}:=\\{g\in\tilde{\mathcal{I}}_{1},\|g\|_{2}\leq\sigma\\}$. By
Steinwart and Christmann [44, Lemma 7.14 and the last display on p. 254], we
have

$\displaystyle e_{i}(I_{\sigma},L_{2}(Z^{n_{2}}))$ $\displaystyle\leq
e_{i-1}(I_{\sigma},L_{2}(Z^{n_{2}})),$
$\displaystyle\mathbb{E}_{Z^{n_{2}}}e_{1}(I_{\sigma},L_{2}(Z^{n_{2}}))$
$\displaystyle\leq\mathbb{E}_{Z^{n_{2}}}\sup_{g\in
I_{\sigma}}\|g\|_{L_{2}(Z^{n_{2}})}\leq(\sigma^{2}+8B\bar{R}_{n_{2}}(\tilde{\mathcal{I}}_{1};\sigma))^{1/2}=:s_{1}.$

Taking expectation on both sides of Steinwart and Christmann [44, Theorem
7.13], we find

$\displaystyle\bar{R}_{n_{2}}(\tilde{\mathcal{I}};\sigma)$
$\displaystyle\leq\sqrt{\frac{\ln
16}{n_{2}}}\mathbb{E}_{Z^{n_{2}}}\Big{(}\sum_{i=1}^{\infty}2^{i/2}e_{2^{i}}(I_{\sigma},L_{2}(Z^{n_{2}}))+\underbrace{\sup_{g\in
I_{\sigma}}\|g\|_{L_{2}(Z^{n_{2}})}}_{\leq s_{1}}\Big{)}.$

Plugging in our Lemma C.4 to the first term in RHS above, we have, for
$j:=2^{i}$ and $m_{2}:=[\log_{2}n_{2}]$,

$\displaystyle{\phantom{=}}\mathbb{E}_{Z^{n_{2}}}\sum_{i=1}^{\infty}2^{\frac{i}{2}}e_{2^{i}}(I_{\sigma},L_{2}(Z^{n_{2}}))$
$\displaystyle\leq\sum_{i=1}^{\infty}2^{\frac{i}{2}}\min\\{s_{1},j^{-\frac{{\bar{b}}+1}{2}}\\}+\sum_{i=1}^{m_{2}}\cancel{2^{\frac{i}{2}}\cdot
j^{-\frac{1}{2}}}\chi_{n}+\sum_{i=m_{2}+1}^{\infty}2^{\frac{i}{2}}n_{2}^{{\bar{b}}+\frac{1}{2}}j^{-({\bar{b}}+1)}\chi_{n}$
$\displaystyle=:S_{1}+\chi_{n}\log n_{2}+S_{2}.$

By Steinwart and Christmann [44, Lemma 7.15] with
$p\leftarrow({\bar{b}}+1)^{-1}$, we have $S_{1}\lesssim
s_{1}^{\frac{{\bar{b}}}{{\bar{b}}+1}},$ where the constant in $\lesssim$ is
independent of all sample sizes. For $S_{2}$, we have

$\displaystyle S_{2}$ $\displaystyle\leq
n_{2}^{{\bar{b}}+\frac{1}{2}}\chi_{n}\int_{x=m_{2}}^{\infty}2^{\frac{x}{2}-({\bar{b}}+1)x}dx\lesssim
n_{2}^{{\bar{b}}+\frac{1}{2}}\chi_{n}2^{\frac{m_{2}}{2}-({\bar{b}}+1)m_{2}}\leq\chi_{n}.$

Combining, we have

$\displaystyle\mathbb{E}_{Z^{n_{2}}}\sum_{i=1}^{\infty}2^{\frac{i}{2}}e_{2^{i}}(I_{\sigma},L_{2}(Z^{n_{2}}))$
$\displaystyle\lesssim s_{1}^{\frac{{\bar{b}}}{{\bar{b}}+1}}+\chi_{n}\log
n_{2},$ $\displaystyle\bar{R}_{n_{2}}(\tilde{\mathcal{I}};\sigma)$
$\displaystyle\leq
Cn_{2}^{-\frac{1}{2}}\left((\sigma^{2}+B\bar{R}_{n_{2}}(\tilde{\mathcal{I}};\sigma))^{\frac{{\bar{b}}}{2({\bar{b}}+1)}}+\chi_{n}\log
n_{2}\right),$

for some constant $C$ which does not depend on $\sigma$ or any sample sizes.
To find solutions to (29) using the above bound, we will restrict to

$\sigma^{2}\geq 2C\Big{(}\frac{\log^{2}n_{2}}{n_{2}}+\chi_{n}^{2}\Big{)}\geq
4Cn_{2}^{-\frac{1}{2}}\chi_{n}\log n_{2}.$ (30)

For such $\sigma$ we now have

$\bar{R}_{n_{2}}(\tilde{\mathcal{I}};\sigma)\leq
Cn_{2}^{-\frac{1}{2}}(\sigma^{2}+B\bar{R}_{n_{2}}(\tilde{\mathcal{I}};\sigma))^{\frac{{\bar{b}}}{2({\bar{b}}+1)}}+\frac{1}{4}\sigma^{2}\leq\max\Big{\\{}3Cn_{2}^{-\frac{1}{2}}\sigma^{\frac{{\bar{b}}}{{\bar{b}}+1}},3CB^{\frac{{\bar{b}}}{{\bar{b}}+2}}n_{2}^{-\frac{{\bar{b}}+1}{{\bar{b}}+2}},\frac{1}{2}\sigma^{2}\Big{\\}}.$

For the RHS to be $\leq\sigma^{2}$ it suffices to consider the first two
terms, leading to $\sigma^{2}\gtrsim
B^{\frac{{\bar{b}}}{{\bar{b}}+2}}n_{2}^{-\frac{{\bar{b}}+1}{{\bar{b}}+2}}$.
Combining with (30) complete the proof. ∎

####  C.4 Average-Case Analysis

######  Corollary C.6 (average-case approximation error).

In the setting of Theorem 4.1, with $\mathcal{D}^{(n_{1})}$-probability $\to
1$ we have

$\mathbb{E}_{g^{*}\sim\mathcal{GP}}\|\tilde{g}^{*}-g^{*}\|_{2}\leq\mathbb{E}_{g^{*}\sim\mathcal{GP}}\epsilon_{n_{1},m}(g^{*}),$

where $\epsilon_{n_{1},m}$ is defined by replacing the $\xi_{n_{1}}$ in (6)
with $m^{-1/2}\xi_{n_{1}}$, and $\|\tilde{g}^{*}\|_{\tilde{\mathcal{I}}}$
satisfies (5). In particular, for $m=[n_{1}^{1/({\bar{b}}+1)}]$, we have

$\mathbb{E}_{g^{*}\sim\mathcal{GP}}\|\tilde{g}^{*}-g^{*}\|_{2}=\tilde{O}\big{(}\xi_{n_{1}}+n_{1}^{-\frac{{\bar{b}}}{2({\bar{b}}+1)}}\big{)}.$
(31)

While we will not use (31) for IV, it hints at the possibility of improvement,
through more careful analyses and/or additional assumptions. Remark C.2
discusses this in more detail.

###### Proof.

Recall the proof of Theorem 4.1, and the notations therein. We can easily
check that all statements in the proposition still hold true, with the high-
probability events now defined w.r.t. ($G$ and) the combined training samples.
Thus, (21) will still hold, and it remains to prove the two approximation
error bounds. We will improve Lemma C.2 and plug into (22). We can check the
proof of the lemma also remains valid, and (25) is still good enough; thus, it
suffices to improve (26) about the term $\|\tilde{e}^{\top}(G-\hat{G})\|_{2}$.

The proof of Theorem 4.1 only invokes Lemma C.2 with

$\bar{e}_{*}=(\lambda_{1}^{-1/2}\langle
g^{*},\psi_{1}\rangle_{2},\ldots,\lambda_{m^{\prime}}^{-1/2}\langle
g^{*},\psi_{m^{\prime}}\rangle_{2}),$

where $\\{(\lambda_{i},\psi_{i})\\}$ are the Mercer decomposition of $k_{x}$.
For $g^{*}\sim\mathcal{GP}(0,k_{x})$, $\bar{e}_{*}$ will distribute as
$N(0,I)$. We will thus modify Lemma C.2 into:

$\mathbb{E}_{\bar{e}_{*}\sim N(0,I)}\|\bar{e}_{*}(\hat{\Psi}-\Psi)\|_{2}\leq
c_{r}(m^{-1/2}\xi_{n}+m^{-({\bar{b}}+1)/2})\sqrt{\log
n}\cdot\mathbb{E}_{\bar{e}_{*}\sim N(0,I)}\|\bar{e}_{*}\|_{2}.$ (32)

(The statement is still restricted to a high-probability event w.r.t. training
samples and $G$.) To prove the above, note that (25) still holds for any fixed
$\bar{e}_{*}$, so the second term above remains correct.
It remains to deal with the first term. Introduce the notation
$S:=\Xi(\Xi^{\top}\Xi)^{-1}$, so that $\tilde{e}=S\bar{e}_{*}$. We have, for
_any fixed $(\mathcal{D}^{(n_{1})},G)$,_

$\displaystyle\mathbb{E}_{\bar{e}_{*}}(\|\tilde{e}^{\top}(G-\hat{G})\|_{2}^{2})$
$\displaystyle=\mathbb{E}_{\bar{e}_{*}}\operatorname{Tr}(\tilde{e}\tilde{e}^{\top}(G-\hat{G})(G-\hat{G})^{\top})$
$\displaystyle=\mathbb{E}_{\bar{e}_{*}}\operatorname{Tr}(S\bar{e}_{*}\bar{e}_{*}^{\top}S^{\top}(G-\hat{G})(G-\hat{G})^{\top})$
$\displaystyle=\operatorname{Tr}(S\cdot\mathbb{E}_{\bar{e}_{*}}(\bar{e}_{*}\bar{e}_{*}^{\top})\cdot
S^{\top}(G-\hat{G})(G-\hat{G})^{\top})$
$\displaystyle=\operatorname{Tr}(SS^{\top}(G-\hat{G})(G-\hat{G})^{\top})$
$\displaystyle\leq\|SS^{\top}\|\operatorname{Tr}((G-\hat{G})(G-\hat{G})^{\top})=\|SS^{\top}\|\sum_{j=1}^{m}\|g^{(j)}-\hat{g}^{(j)}_{n_{1}}\|_{2}^{2}.$

The inequality above is the von Neumann inequality
$\operatorname{Tr}(AB)\leq\|A\|_{op}\operatorname{Tr}(B)$. Conditioned on the
$(\mathcal{D}^{(n_{1})},G)$-measurable event on which (17) and (20) hold, the
first term above is bounded by $c_{r}^{2}m^{-1}$, and the second by
$m\xi_{n}^{2}\log n$. Thus we have, with
$(\mathcal{D}^{(n_{1})},G)$-probability $\to 1$,

$\mathbb{E}_{g^{*}}\|\tilde{e}^{\top}(G-\hat{G})\|_{2}^{2}\leq
c_{r}^{2}\xi_{n}^{2}\log n\lesssim c_{r}^{2}\xi_{n}^{2}\log
n\frac{(\mathbb{E}_{\bar{e}_{*}}\|\bar{e}_{*}\|_{2})^{2}}{m},$

where the last inequality follows from Vershynin [77, Theorem 3.1.1]. Since
$\mathbb{E}X^{2}\geq(\mathbb{E}X)^{2}$, we complete the new bound for the
first term in (32), and subsequently (32). Following the original proof, we
can see that (22) becomes

$\mathbb{E}_{g^{*}\sim\mathcal{GP}}\|\tilde{g}^{*}-g^{*}\|_{2}\leq\mathbb{E}_{g^{*}\sim\mathcal{GP}}\left[c_{r}\|\mathrm{Proj}_{m^{\prime}}{g^{*}}\|_{\mathcal{I}}(m^{-1/2}\xi_{n}+m^{-({\bar{b}}+1)/2})\sqrt{\log
n}+\|g^{*}-\mathrm{Proj}_{m^{\prime}}{g^{*}}\|_{2}\right],$

which proves the first claim.

For the second claim, from Jensen’s inequality, we know
$\mathbb{E}_{g^{*}\sim\mathcal{GP}}\|\mathrm{Proj}_{m^{\prime}}{g^{*}}\|_{\mathcal{I}}=\mathbb{E}_{\bar{e}_{*}}\|\bar{e}_{*}\|_{2}\leq\sqrt{\mathbb{E}_{\bar{e}_{*}}\|\bar{e}_{*}\|_{2}^{2}}=m^{\prime}$,
and thus the first term in the above display is
$\tilde{O}(\xi_{n}+m^{-\bar{b}/2})$. Similarly, we have

$\mathbb{E}_{g^{*}\sim\mathcal{GP}}\|g^{*}-\mathrm{Proj}_{m^{\prime}}{g^{*}}\|_{2}\leq\left(\mathbb{E}_{g^{*}\sim\mathcal{GP}}\|g^{*}-\mathrm{Proj}_{m^{\prime}}{g^{*}}\|_{2}^{2}\right)^{1/2}\overset{\eqref{eq:kl-
trunc-error}}{\leq}O(m^{-\bar{b}/2}),$

which completes the proof. ∎

######  Remark C.2.

The only change in this proof is a new bound for the estimation error from the
oracle. The Cauchy-Schwarz inequality used in the previous bound (26) was is
typically loose: it requires the vector $\tilde{e}$, determined by the test
function, to be parallel with
$(g^{(1)}-\hat{g}^{(1)}_{n_{1}},\ldots,g^{(m)}-\hat{g}^{(m)}_{n_{1}})$, the
estimation residuals. This is at odds with the intuition that residual
functions for each independently drawn GP samples to be in some sense
uncorrelated.

For an extreme example, suppose our target RKHS
$\mathcal{I}=\mathrm{span}\\{\psi_{1}\\}$ is one-dimensional, so that the
standard GP prior draw can be written as $\epsilon\psi_{1}$ for some
$\epsilon\sim\mathcal{N}(0,1)$. Then, independent draws from
$\mathcal{GP}(0,k_{z})$ should have independent signs, and it seems very
strange if a useful regression algorithm often returns residual functions with
correlated signs determined by $\tilde{e}$, for such $g^{(\cdot)}$. The fact
that such oracles are allowed by our assumption suggest there is room for
improvement, although we do not pursue this path, since the suboptimality is
relatively mild (see the end of Appendix C.5.2). Empirically, our estimator
works well for $n_{1}=n_{2}$.

Note that while similar problems have been studied in multi-task learning,
such works often assume independent inputs for the different tasks, which
corresponds to the different GP prior draws in our setting. Corollary C.6
appears new in this aspect, at least among works that established fast rates.

####  C.5 Deferred Derivations and Additional Discussions

#####  C.5.1 Derivation for Example 4.1

We first derive the expression of $\xi_{n}$. For regression functions of the
form $g=\bar{g}\circ\Phi$, where $\Phi$ is defined as in the text, and
$\bar{g}$ is $\underline{\beta}_{2}$ Hölder regular and have $d_{l}$ inputs,
Schmidt-Hieber [10] establishes the convergence rate

$\xi_{n}\lesssim(n^{-\frac{\beta_{1}}{2\beta_{1}+d_{z}}}+n^{-\frac{\underline{\beta}_{2}}{\underline{\beta_{2}}+d_{l}}})\|\bar{g}\|_{C^{\underline{\beta}_{2}}}\log^{3/2}n+\epsilon_{opt}.$
(33)

(We view $\|\Phi\|_{C^{\beta_{1}}}$ as a constant.) The result holds uniformly
for all $\bar{g}\in C^{\underline{\beta}_{2}}$ with uniformly bounded norm.

Recall Example A.3: for any $\epsilon>0$, there exists a version121212We can
work with any version of the GP prior since the regression oracle only
accesses its evaluation on a finite number of points, which have the same
distribution among all versions. of the Matérn GP $\Pi_{z}$ for $\bar{g}$ s.t.
$\bar{g}\in C^{\beta_{2}-\epsilon}$ a.s. We set
$\underline{\beta}_{2}=\beta_{2}-\epsilon$. By Lemma A.4, the random variable
$\|\bar{g}\|_{C^{\underline{\beta}_{2}}}$ has a subgaussian tail. Thus for any
$p>1$, we can choose some $C_{1}>0$ which depend on $p$ and $\epsilon$,
leading to

$\Pi_{z}(E):=\Pi_{z}(\\{\bar{g}:\|\bar{g}\|_{C^{\underline{\beta}_{2}}}^{2}\leq
C_{1}\log n\\})\geq 1-n^{-p}.$

Let the NN model in Schmidt-Hieber [10] be constructed for
$g=\bar{g}\circ\Phi$ satisfying the above norm bound, and $\hat{g}_{n}$ denote
the resulted estimator. Consider

$\displaystyle\mathbb{E}_{\Pi_{z}}\|\hat{g}_{n}-g\|_{2}^{2}$
$\displaystyle\leq\mathbb{E}_{\Pi_{z}}[\mathbf{1}_{E}\|\hat{g}_{n}-g\|_{2}^{2}]+2\mathbb{E}_{\Pi_{z}}[\mathbf{1}_{E^{c}}(\|\hat{g}\|_{\infty}^{2}+\|g\|_{2}^{2})]$
$\displaystyle\leq\mathbb{E}_{\Pi_{z}}[\mathbf{1}_{E}\|\hat{g}_{n}-g\|_{2}^{2}]+\frac{2\mathbb{E}_{\Pi_{z}}[\|\hat{g}\|_{\infty}^{2}\mid
E^{c}]}{n_{2}^{p}}+2\mathbb{E}_{\Pi_{z}}[\mathbf{1}_{E^{c}}\|g\|_{2}^{2}].$
(34)

The first term clearly has the desired bound, with two extra logarithms. For
the second term, note the estimator in Schmidt-Hieber [10] has sup norm
bounded as $\|\hat{g}_{n}\|_{\infty}\leq C_{1}\log n$, so it also has the
desired bound. For the last term, another application of Lemma A.4 shows that
$\|\bar{g}\|_{2}$ norm also has a subgaussian tail. Let $\Phi_{2}$ denote its
CDF. Then we have

$\mathbb{E}(\mathbf{1}_{E^{c}}\|g\|_{2}^{2})\leq\int_{\Phi_{2}^{-1}(1-n_{2}^{-p})}^{+\infty}x^{2}d\Phi_{2}(x)\lesssim\int_{C_{2}\sqrt{p\log
n_{2}}}^{\infty}x^{2}e^{-C_{3}x^{2}}dx\lesssim n_{2}^{-C_{4}p}\sqrt{p\log
n_{2}}.$

In the above, $C_{2},C_{3},C_{4}$ are determined by $\Pi_{z}$, so we can
choose a sufficiently large $p$ so that the RHS is $\lesssim n^{-1}$. This
completes the derivation for $\xi_{n}$. ∎

We now turn to the regression error using a fixed-form kernel. The Matérn
process $\bar{g}\sim\mathcal{GP}(0,\bar{k}_{z})$ takes value in the Hölder
space $C^{\underline{\beta}_{2}}(\bar{\mathcal{Z}})$ w.p.1, and the order
cannot be made larger (Example A.3). Thus, the random function
$g=\bar{g}\circ\Phi$ can only take value in
$C^{\min\\{\underline{\beta}_{2},\beta_{1}\\}}(\mathcal{Z})$ w.p.1. Our
claimed regression rate follows from [38] (for Matérn kernels), or [79] (for
Gaussian/RBF kernels with an adaptive bandwidth).131313Note that we did not
rule out the possibility of a better rate being attainable: we did not prove a
lower bound, instead only presented the best known upper rate given our Hölder
regularity condition of $g$. However, improvement is most likely impossible:
as discussed in introduction, separability results have been established in
similar feature learning settings for certain fixed-form kernels; in the
setting of this example, [10] established a lower bound for a wavelet model
with a similar order. ∎

#####  C.5.2 Derivations for Example 4.2

By properties of the Matérn RKHS (Example A.1), the latent-space kernel
$\bar{k}_{z}$ has eigendecay $\lambda_{i}\lesssim i^{-(1+2\beta_{2}/d_{l})}$.
Thus we have ${\bar{b}}+1=1+\nicefrac{{2\beta_{2}}}{{d_{l}}}$, and

$\xi_{n_{1}}=\tilde{\mathcal{O}}\Big{(}\epsilon_{fea,n_{1}}+n_{1}^{-\frac{{\bar{b}}-\epsilon}{2({\bar{b}}-\epsilon+1)}}\Big{)}+\epsilon_{opt},$

for all $\epsilon>0$. We first establish the three approximation error bounds.

  1. 1.

The first claim (Case (i) in the text) follows by observing
$m^{-{\bar{b}}}=n_{1}^{-{\bar{b}}/({\bar{b}}+1)}\ll\xi_{n_{1}}$.

  2. 2.

For Case (ii), let $\tilde{g}_{0}\in\tilde{\mathcal{I}}$ be the approximation
returned by Corollary C.6, then we have

$\mathbb{E}_{g}\|\tilde{g}_{0}-g\|_{2}\overset{\eqref{eq:avg-approx-
bound}}{=}\tilde{O}(\xi_{n_{1}}+n_{1}^{-\frac{{\bar{b}}}{2({\bar{b}}+1)}})=\tilde{O}(\xi_{n_{1}}).$
(35)

It remains to provide an RKHS norm bound. By convergence of the Karhunen-Loève
expansion [66, Thm. 3.1] we have
$\mathrm{Proj}_{m}{g}=\sum_{i=1}^{m}\epsilon_{i}(\lambda^{z}_{i})^{1/2}\psi^{z}_{i}$,
where $\\{(\lambda^{z}_{i},\psi^{z}_{i})\\}$ denote the Mercer representation
of $k_{z}$, and $\epsilon_{i}\sim\mathcal{N}(0,1)$ are iid. As
$(\lambda^{z}_{i})^{1/2}\psi^{z}_{i}$ is an ONB of $\mathcal{H}$ (Lem. A.1),
we have
$\|\mathrm{Proj}_{m}{g}\|_{\mathcal{H}}^{2}=\sum_{i=1}^{m}\epsilon_{i}^{2}.$
Thus, $\chi^{2}$-concentration bounds [46, Example 2.11] yields

$\mathbb{P}(\|\mathrm{Proj}_{m}{g}\|_{\mathcal{H}}^{2}>2\sqrt{2}m)\leq
e^{-m}=\exp\big{(}-n^{\frac{1}{{\bar{b}}+1}}\big{)}.$

Thus, we let $\tilde{g}=0$ on the above event, and $\tilde{g}_{0}$ otherwise.
Repeating the argument of (34) we can show that the additional $L_{2}$ error
is negligible, and $\tilde{g}$ satisfies the same bound of (35). In this case,
the average-case $L_{2}$ rate of DNN is $\xi_{n_{1}}$, by definition.

  3. 3.

Additionally, we claim that in Case (ii), for
$m=[n^{\frac{{\bar{b}}}{({\bar{b}}+1)^{2}}}]$, there exists
$\tilde{g}^{*}\in\tilde{\mathcal{I}}$ s.t.

$\mathbb{P}_{g^{*}\sim\mathcal{GP}}\big{(}\|\tilde{g}^{*}\|_{\tilde{\mathcal{I}}}\leq
c_{r}n_{1}^{1/2({\bar{b}}+1)},~{}\|\tilde{g}^{*}-g^{*}\|_{2}=\tilde{\mathcal{O}}(m^{1/2}\xi_{n_{1}}+n_{1}^{-{\bar{b}}^{2}/2({\bar{b}}+1)^{2}})\big{)}\geq
1-e^{-m}\to 1.$

To show this, we shall apply Theorem 4.1 to a different target function
$g^{*}$. Let us denote by $g$ the GP prior draw to be approximated. Then we
have the standard prior mass bound for the “sieve set” [36]: for some $C>0$,

$\forall\bar{n}\in\mathbb{N},\mathbb{P}_{g\sim\mathcal{GP}(0,k_{z})}\big{(}\exists
g^{\dagger}_{\bar{n}},\|g^{\dagger}_{\bar{n}}\|_{\mathcal{I}}\leq
C\bar{n}^{1/2({\bar{b}}+1)},\|g^{\dagger}_{\bar{n}}-g\|_{2}\leq
C\bar{n}^{-{\bar{b}}/2({\bar{b}}+1)}\big{)}\geq
1-\exp\big{(}-\bar{n}^{\frac{1}{{\bar{b}}+1}}\big{)}.$ (36)

(This is proved by combining the Borell inequality [e.g., 71, Prop. 11.17] and
the $L_{2}$ small-ball probability bound [e.g., 66, Cor. 4.9, with
$\beta=1$].) We set $\bar{n}:=[n_{1}^{{\bar{b}}/({\bar{b}}+1)}]$ and restrict
to $g$ in the above event. Invoking Lemma A.3 with $g\leftarrow
g^{\dagger}_{\bar{n}},m\leftarrow[\bar{n}^{\frac{1}{b+1}}]$ leads to

$\|\mathrm{Proj}_{m}{g^{\dagger}_{\bar{n}}}\|_{\mathcal{I}}\leq\|g^{\dagger}_{\bar{n}}\|_{\mathcal{I}}\lesssim\bar{n}^{\frac{1}{2({\bar{b}}+1)}},~{}~{}\|\mathrm{Proj}_{m}{g^{\dagger}_{\bar{n}}}-g^{\dagger}_{\bar{n}}\|_{2}^{2}\lesssim\bar{n}^{-\frac{{\bar{b}}}{2({\bar{b}}+1)}}.$

Invoking Theorem 4.1 with $m\leftarrow
2[\bar{n}^{\frac{1}{{\bar{b}}+1}}],g^{*}\leftarrow\mathrm{Proj}_{m}{g^{\dagger}_{\bar{n}}}$
yields

$\|\tilde{g}\|_{\tilde{\mathcal{I}}}\leq
c_{r}\bar{n}^{\frac{1}{2({\bar{b}}+1)}},~{}~{}\|\tilde{g}-g^{*}\|_{2}\leq
c_{r}\bar{n}^{\frac{1}{2({\bar{b}}+1)}}(\xi_{n_{1}}+\bar{n}^{-1/2})\log
n_{1}+0.$

Two applications of triangle inequality yield

$\|\tilde{g}-g\|_{2}\lesssim\bar{n}^{\nicefrac{{{\bar{b}}}}{{2({\bar{b}}+1)^{2}}}}(\xi_{n_{1}}+\bar{n}^{-\nicefrac{{1}}{{2}}})\log
n+\bar{n}^{-\nicefrac{{{\bar{b}}}}{{2({\bar{b}}+1)}}}=\tilde{\mathcal{O}}\Big{(}\epsilon_{fea,n_{1}}+n_{1}^{-\nicefrac{{({\bar{b}}^{2}-\epsilon)}}{{2({\bar{b}}+1)^{2}}}}\Big{)},$

for all $\epsilon>0$. This result is most useful for functions satisfying the
conditions in (36), which may or may not be random GP prior draws.

In Case (ii) the average-case DNN rate is $\xi_{n_{1}}$, by definition. Case
(iii) can thus be suboptimal when feature learning is easy, as the exponent of
the second term above is worse by a multiplicative factor of
$({\bar{b}}-\epsilon)/({\bar{b}}+1)$. Note the suboptimality can also be
removed if we can take $n_{1}\gtrsim
n_{2}^{\nicefrac{{({\bar{b}}+1)}}{{{\bar{b}}}}}$, at which point estimation
error starts to dominate the combined regression error. The difference in
order diminishes as ${\bar{b}}\to\infty$. The situation in Case (i) is less
clear, since it is unclear if the DNN can take advantage of the improved
regularity of $g^{*}\in\mathcal{I}$, see below.∎

######  Remark C.3 (optimal rate for $g^{*}\in\mathcal{I}$).

We discuss whether the DNN oracle in [10] may estimate $g^{*}\in\mathcal{I}$
with a rate better than $\xi_{n}$.

The condition that $g^{*}\in\mathcal{I}$, or equivalently
$\bar{g}^{*}\in\bar{\mathcal{I}}$, implies that $\bar{g}^{*}$ now possesses a
better Sobolev regularity comparing with typical GP samples (Example A.1-A.3),
but it provides no additional information on Hölder regularity, which is
needed in [10]. If we additionally assume the Hölder regularity is also
improved by the same amount (i.e., $\bar{g}^{*}\in
C^{\underline{\beta}_{2}+(d_{l}/2)}(\bar{\mathcal{Z}})$), the DNN rate will
improve to
$\tilde{\mathcal{O}}(\epsilon_{fea,n_{1}}+n_{1}^{-\nicefrac{{{\bar{b}}+1}}{{2({\bar{b}}+2)}}})$,
by (33); this is smaller than $\xi_{n_{1}}$ if feature learning is
sufficiently easy. Without the extra Hölder regularity, we can only invoke the
Sobolev embedding theorem to obtain $\bar{g}^{*}\in
C^{\underline{\beta}_{2}}(\bar{\mathcal{Z}})$; the derived DNN rate can only
be infinitesimally better than $\xi_{n_{1}}$.

#####  C.5.3 Alternative Kernel Learning Schemes

Let us compare the proposed algorithm with a more obvious alternative. For
simplicity, we assume the feature learning term in $\xi_{n_{1}}$ does not
dominate, and restrict to Case (i) in Example 4.2 ($g^{*}\in\mathcal{I}$).

The established bound requires $O(n_{1}^{{\bar{b}}/({\bar{b}}+1)^{2}})$ calls
to the regression oracle.141414Or equivalently, solving a vector-valued
regression problem with a similar output dimension. In the IV setting, we have
$\bar{b}\geq\max\\{2p-1,b\\}$, and possibly $\bar{b}=b+2p$ (App. B.1). An
alternative instrument learning procedure is to construct a $d$-dimensional
approximation to $\mathcal{H}$, invoke the oracle on its basis functions, and
use the estimates to construct an alternative $\tilde{\mathcal{I}}$. For an
approximation error of $n_{1}^{-(b+1)/2(b+2)}$, we need $d\asymp
n_{1}^{1/(b+2)}$ bases [80, 81], which becomes much worse than
$n_{1}^{1/({\bar{b}}+2)}$ when $p$ is moderately large. Moreover, finding such
an approximation is difficult: the same references note that intuitive
choices, such as uniformly sampled random feature, or Nyström inducing points,
actually require $d\asymp n^{\frac{b}{b+1}}\to n$. It may be also possible
that such an algorithm actually has worse performance, as we have not checked
if its dependency on the oracle estimation error remains unchanged.

Our approximation gracefully adapts to the ill-posedness of the problem, while
also being simpler to implement. Intuitively, this is because it directly
focuses on the approximation of the optimal $\mathcal{I}$ as opposed to
$\mathcal{H}$. In other words, it works with a basis
$\\{\varphi_{i}\\}\subset\mathcal{H}$ that is not necessarily optimal for the
approximation of $\mathcal{H}$, but is designed so that $\\{E\varphi_{i}\\}$
approximate $\mathcal{I}$ well. Interestingly, we do not require the knowledge
of such a basis, instead automatically adapting to it. This is made possible
by the isotropy of the $\mathcal{GP}(0,k_{x})$ prior: the GP can be viewed as
being defined by an aribtrary basis, convergence technicalities
notwithstanding.

###  Appendix D Approximation and Estimation Results for Gaussian Process
Regression

This section includes various technical results for using
$\tilde{\mathcal{I}}$ in the “GP regime”, i.e., when the regression function
does not live in $\tilde{\mathcal{I}}$, but only satisfies an approximability
condition like Asm. 2.2 (iii).b. Regularization becomes weaker in the GP
scheme than the kernel scheme; to see this, consider a regression task with
data $\\{(z_{i},\bar{y}_{i}=g_{0}(z_{i})+\epsilon_{i}):i\in[\bar{n}]\\}$.
Observe that both the KRR estimate and the GP posterior mean estimator can
have the form of

$\hat{g}_{n}:=\arg\min_{g\in\tilde{\mathcal{I}}}\frac{1}{\bar{n}}\sum_{i=1}^{\bar{n}}(g(z_{i})-\bar{y}_{i})^{2}+\bar{\nu}\|g\|_{\tilde{\mathcal{I}}}^{2},$

but the kernel scheme uses $\bar{\nu}\gg\bar{n}^{-1}$ (unless the RKHS has
fixed finite dimensionality), to dominate the critical radius of the RKHS
[46], whereas the GP scheme uses $\bar{\nu}\asymp\bar{n}^{-1}$. See [38, 67]
for more discussions. It is thus understandable that the GP scheme requires
additional assumptions.

We first impose the following assumption on the _true RKHS_ $\mathcal{I}$:

######  Assumption D.1.

Let the Mercer eigenvalues of $\mathcal{I}$ be $\lambda_{i}\asymp
i^{-({\bar{b}}+1)}$. Then (EMB) holds for
$\gamma=\frac{{\bar{b}}-1}{{\bar{b}}+1}$.

Recall the “top-level RKHS” $\bar{\mathcal{I}}$ and $\mathcal{I}$ satisfy
embedding properties of the same order. As discussed in Example A.1, this
assumption holds when $\bar{\mathcal{I}}$ is a Matérn RKHS with ${\bar{b}}>2$,
and its kernel is supported on a bounded, open subset of $\mathbb{R}^{d_{l}}$
with smooth boundaries.

The above assumption controls estimation error using the true RKHS
$\mathcal{I}$. To see this, observe Lemma A.2 now implies

$\|g\|_{\infty}\lesssim\|g\|_{\mathcal{I}}^{\frac{{\bar{b}}-1}{{\bar{b}}+1}}\|g\|_{2}^{\frac{2}{{\bar{b}}+1}},\quad\forall
g\in\mathcal{I}.$ (37)

Fast-rate convergence analyses typically require sup norm bounds on a
“relevant subset” of the hypothesis space, such as
$\\{g-g_{0}:g\in\mathcal{I},\|g-g_{0}\|_{2}\text{ is small}\\}$. The
interpolation inequality (37) provides such a bound if we were using the true
model $\mathcal{I}$. Note that such assumptions are not needed in the kernel
regularization scheme, because as discussed above, in that case the ridge
regularizer has a greater magnitude, thus providing more control for the RKHS
norm on the relevant subset.

However, we will use the approximate model $\tilde{\mathcal{I}}$ for
estimation, so to control similar sup norm quantities, we also need some crude
bound on the sup-norm approximation error of $\tilde{\mathcal{I}}$. This, in
turn, requires the following crude sup-norm error bound for the regression
oracle.

######  Assumption D.2.

On a high-probability event determined by $\mathcal{D}^{(n_{1})}$, we have,
for $m\asymp n^{\nicefrac{{1}}{{{\bar{b}}+1}}}$,
$\sum_{j=1}^{m}\|g^{(j)}-\hat{g}^{(j)}_{u,n_{1}}\|_{\infty}^{2}\lesssim 1.$

The assumption requires sup norm error to scale at
$n^{-\nicefrac{{1}}{{2({\bar{b}}+1)}}}$. We expect it to be mild, at least
when ${\bar{b}}$ is reasonably large (and feature selection is not too hard).
This is because in the classical kernel literature, the sup norm error rate
for Sobolev kernels is $O(n^{-\nicefrac{{({\bar{b}}-1)}}{{2({\bar{b}}+1)}}})$
[48]. Moreover, if the GP samples are
$\nicefrac{{({\bar{b}}-\epsilon)d}}{{2}}$-Hölder regular as in the Matérn
example, error rates in $L_{2}$ and sup norm will both be
$O(n^{-\nicefrac{{({\bar{b}}-\epsilon)}}{{2({\bar{b}}+1-\epsilon)}}})$.

####  D.1 Sup Norm Approximation

All results in this subsection assume Assumption D.1.

######  Lemma D.1.

Let $\Pi_{g}$ denote the standard GP prior defined by $k_{z}$, and
$\tau_{n}=n^{-\frac{{\bar{b}}/2}{{\bar{b}}+1}}$. There exist constants
$C,C^{\prime}>0$ s.t. (i) for all $n\in\mathbb{N}$,

$\Pi_{g}(\Theta_{n}):=\Pi_{g}(\\{g=g_{h}+g_{e}:\|g_{h}\|_{\mathcal{I}}^{2}\leq
Cn\tau_{n}^{2},\|g_{e}\|_{2}^{2}\leq C\tau_{n}^{2},\|g_{e}\|_{\infty}^{2}\leq
Cn^{-\frac{1}{{\bar{b}}+1}}\\})\geq 1-e^{-C^{\prime}n\tau_{n}^{2}}.$ (38)

(ii) In the above display, we can always choose $g_{h}=\mathrm{Proj}_{m}{g}$,
for $m=[n^{\frac{1}{{\bar{b}}+1}}]$.

###### Proof.

(i) is a consequence of Steinwart [66] and is proved in Wang et al. [26,
Corollary 13, with $b\leftarrow{\bar{b}},b^{\prime}\leftarrow b-1$]. To show
(ii), it suffices to verify that

$\tilde{g}_{h}:=\mathrm{Proj}_{m}{g_{h}}+\mathrm{Proj}_{m}{g_{e}},~{}~{}\tilde{g}_{e}:=\mathrm{Proj}_{>m}{g_{h}}+\mathrm{Proj}_{>m}{g_{e}}$

still satisfy the above display. For $\tilde{g}_{h}$, let $\lambda_{m}\asymp
m^{-({\bar{b}}+1)}$ be the $m$-th Mercer eigenvalue for $k_{z}$. Then

$\displaystyle\|\tilde{g}_{h}\|_{\mathcal{I}}$
$\displaystyle\leq\|g_{h}\|_{\mathcal{I}}+\|\mathrm{Proj}_{m}{g_{e}}\|_{\mathcal{I}}\leq\sqrt{Cn\tau_{n}^{2}}+\lambda_{m}^{-1/2}\|\mathrm{Proj}_{m}{g_{e}}\|_{2}\lesssim
n^{\frac{1/2}{{\bar{b}}+1}}=\sqrt{n\tau_{n}^{2}}.$
$\displaystyle\|\mathrm{Proj}_{>m}{g_{h}}\|_{2}$
$\displaystyle\leq\lambda_{m}^{1/2}\|g_{h}\|_{\mathcal{I}}\lesssim\tau_{n}.$
(39) $\displaystyle\|\tilde{g}_{e}\|_{2}$
$\displaystyle\leq\|g_{e}\|_{2}+\|\mathrm{Proj}_{>m}{g_{h}}\|_{2}\leq
C\tau_{n}+\lambda_{m}^{1/2}\|g_{h}\|_{\mathcal{I}}\lesssim\tau_{n}.$
$\displaystyle\|\tilde{g}_{e}\|_{\infty}$
$\displaystyle\leq\|g_{e}\|_{\infty}+\|\mathrm{Proj}_{>m}{g_{h}}\|_{\infty}\overset{\eqref{eq:interpolation}}{\leq}Cn^{-\frac{1/2}{{\bar{b}}+1}}+(\|\mathrm{Proj}_{>m}{g_{h}}\|_{\mathcal{I}})^{\frac{{\bar{b}}-1}{{\bar{b}}+1}}(\|\mathrm{Proj}_{>m}{g_{h}}\|_{2})^{\frac{2}{{\bar{b}}+1}}$
$\displaystyle\leq
Cn^{-\frac{1/2}{{\bar{b}}+1}}+n^{\frac{1/2}{{\bar{b}}+1}\cdot\frac{{\bar{b}}-1}{{\bar{b}}+1}}(\|\mathrm{Proj}_{>m}{g_{h}}\|_{2})^{\frac{2}{{\bar{b}}+1}}\overset{\eqref{eq:sieve-