= K_{*x} (\lambda  + LK_{xx} )^{-1}L Y.
\]

As $\tilde Y-Y, g_0$ and $f_0$ are independent, the covariance operator of $f^*$ is

\begin{align*}
\cC^\prime &= {\lambda^{-1}\cC}(
\empCxz\invEmpCz({n\lambda\empCz} + \lambda \nu I)
\invEmpCz\empCzx + \lambda^2 I
)\lambda^{-1}\cC  \\
&= \lambda^{-1}\cC(
\lambda n\empCxz\invEmpCz\empCzx + \lambda^2 I
\overset{\eqref{eq:equiv-cov-form}}{=}
\cC. %
\end{align*}

In view of (<ref>), we know

\[ \begin{aligned}
S_*\cC^\prime S^*_*
&=  S_*S_*^* - S_*S_x^*(\lambda L^{-1} + S_x S_x^*)^{-1}S_xS_*^*
=  K_{**} - K_{*x}(\lambda L^{-1} + K_{xx})^{-1}K_{x*},
\end{aligned} \]

which matches the posterior covariance matrix (<ref>) with replaced kernels.

For the discussion in Section <ref>,
observe that (<ref>) implies the posterior mean estimator $\hat f_n$ satisfies

\begin{equation}\label{eq:posterior-mean-functional}
\hat f_n =
(n\empCxz\invEmpCz\empCzx+\lambda I)^{-1} n\empCxz\invEmpCz
\Biggl(\empCzx f_0 + \frac{\empSz^*(Y-f_0(X))}{n}\Biggr),
\end{equation}

The derivations do not rely on the RKHSes being finite-dimensional.

§.§ Assumptions in Proposition <ref>

The subsequent analysis will rely
on the following assumptions on the random feature expansion. We only state them for $x$ for conciseness; the requirements for $z$ are similar.

The following assumption holds for, e.g., random Fourier features [Rahimi et al., 2007].

As $m\to\infty$, it holds that
\sup_{x,x'\in\cX} \left | k_x(x,x') - \tilde{k}_{x,m}(x,x') \right |  \overset{p}{\rightarrow} 0.

There exists a constant $\tilde{\kappa}>0$ such that
$\max_{m\in\mb{N}} \sup_{x\in\cX} \tilde{k}_{x,m}(x,x)\le \tilde{\kappa}$.

§.§ Analysis of Random Feature Approximation

We recall the following facts: for $A,B\in\mb{R}^{n\times n}$,

$$\|A\| \le \|A\|_F \le \sqrt{n}\|A\|, ~~ A^{-1}-B^{-1} = A^{-1}(B-A)B^{-1}.$$

For all $m\in\mb{N}$, let $k_{x,m}$ be a random feature approximation to $k_x$ such that Assumption <ref> holds, and
let $\tilde{k}_{z,m}$ be an approximation to $k_z$ satisfying a similar requirement as above. Then the random feature-approximated posterior
\Pi_m(f(x_*) \givendata) = \cN(\tilde{\mu}, \tilde{S})

\lim_{m\rightarrow\infty} \sup_{x^*\in\cX^l} \|\mu - \tilde{\mu}\|_2 = 0, ~~
\lim_{m\rightarrow\infty} \sup_{x^*\in\cX^l} \|\tilde{S} - S\|_F= 0,

for any fixed training data $(X,Y,Z)$, $l\in\mb{N}$, and $\lambda,\nu>0$. In the above, $\tilde \mu$ and $\tilde S$ are defined as in (<ref>)-(<ref>), but with
the Gram matrices redefined
using $\tilde{k}_{x,m}$ and $\tilde{k}_{z,m}$.

\epsilon_m := \max\Big(
\sup_{x,x'\in\cX} \left | k(x,x') - \tilde{k}_{x,m}(x,x') \right | ,
\sup_{z,z'\in\cZ} \left | k(z,z') - \tilde{k}_{z,m}(z,z') \right | \Big).

By assumption $\epsilon_m\overset{p}{\rightarrow} 0$.
For $\tilde{S}$ we consider the decomposition

\begin{align*}
\|\tilde{S} - S\| %
&\le \|\tilde{K}_{**} - K_{**}\|   \\
&\phantom{{}={}}+\|\tilde{K}_{*x} - K_{*x}\|  \|
\tilde{L}\|\|(\lambda I + \tilde{K}_{xx}\tilde{L})^{-1}\tilde{K}_{x*} \| \\
&\phantom{{}={}} +\|K_{*x}\| \| \tilde{L} - L\| \|(\lambda I+\tilde{K}_{xx}\tilde{L})^{-1}\tilde{K}_{x*}\|  \\
&\phantom{{}={}} +
\|K_{*x} L\|\|(\lambda I+\tilde{K}_{xx}\tilde{L})^{-1}
- (\lambda I+ K_{xx}L)^{-1}\| \|\tilde{K}_{x*}\| \\
&\phantom{{}={}} +
\|K_{*x}L(\lambda I+K_{xx}L)^{-1}\|\|\tilde{K}_{x*}-K_{x*}\| \\
&=: (\text{I}) + (\text{II}) + (\text{III}) + (\text{IV}) + (\text{V}).
\end{align*}

In the following, we use $O(\cdot)$ and $O_p(\cdot)$ to represent the asymptotic behaviour when $m \to \infty$.
Since $n$ and $l$ are fixed, the operator norms of the matrices $K_{*x},L,K_{xx}$ are $O(1)$.
Observe that $\|K_{zz}-\tilde{K}_{zz}\|\le \sqrt{n}\epsilon_m$.
By the triangle inequality, the inequality $\|\cdot\|\le \|\cdot\|_F$ and the boundedness of $\tilde{k}_{x,m}$ and $\tilde{k}_{z,m}$, we have $\|\tilde{K}_{*x}\| = O(1)$. Both $O(\cdot)$ terms above are independent of $x^*$. Finally, recall that $\|L\| = \|K_{zz}(K_{zz}+\nu I)^{-1}\|\le 1$ and similarly $\|\tilde{L}\|\le 1$. Using these facts, we have

\begin{align*}
(\text{I}) &\le \|\tilde{K}_{**}-K_{**}\|_F \le l \epsilon_m  \rightarrow 0. \\
(\text{II}) &\le \sqrt{ln}\epsilon_m \cdot 1 \cdot \lambda^{-1}\cdot O(1) \rightarrow 0. \\
\|\tilde{L}-L\| &= \|K_{zz}(K_{zz}+\nu I)^{-1} - \tilde{K}_{zz}(\tilde{K}_{zz}+\nu I)^{-1}\| \rightarrow 0. \\
\|K_{zz}-\tilde{K}_{zz}\|\cdot \nu^{-1} + \|\tilde{K}_{zz}(\tilde{K}_{zz}+\nu I)^{-1}\|\|(K_{zz}-\tilde{K}_{zz})(K_{zz}+\nu I)^{-1}\| \\
&\le 2\sqrt{n}\epsilon_m \cdot \nu^{-1} \rightarrow 0.
\\
(\text{III}) &\le O(1)\cdot \|\tilde{L}-L\|\cdot \lambda^{-1}O(1)
\rightarrow 0 \\
(\text{IV}) &= O(1)\cdot \|(\lambda I+\tilde{K}_{xx}\tilde{L})^{-1}\|\|\tilde{K}_{xx}\tilde{L}-K_{xx}L\|\|(\lambda I+ K_{xx}L)^{-1}\| \\
&\le  O(1)\cdot\lambda^{-2}\cdot (\|\tilde{K}_{xx}-K_{xx}\|\|\tilde{L}\| + \|K_{xx}\|\|\tilde{L}-L\|)
\rightarrow 0. \\
(\text{V}) &= O(1) \cdot \sqrt{ln}\epsilon_m \rightarrow 0.
\end{align*}

Moreover, the converges above are all independent of the choice of $x^*$.
Thus we have
\sup_{x^*\in\cX^l} \|\tilde{S}-S\|_F\le l \sup_{x^*\in\cX^l} \|\tilde{S}-S\|\rightarrow 0.

Using a similar argument we have
\sup_{x^*\in\cX^l} \|\tilde{\mu}-\mu\|_2 \rightarrow 0.

§.§ Analysis of the Optimization Algorithm

Hyperparameters $\nu,\lambda\in\RR$. Random feature models
$\theta\mapsto f(\cdot;\theta)$, $\varphi\mapsto g(\cdot;\varphi)$.

A single sample from the approximate posterior
Initialize: draw $\theta_0\sim \cN(0,I), \varphi_0\sim \cN(0,\lambda\nu^{-1}I),\tilde{Y}\sim\cN(Y,\lambda I)$
()$\ell\gets 1,\ldots,L-1$
$\hat{\theta}_\ell \gets \theta_{\ell-1} - \eta_\ell \hat \nabla_{\theta} \cL_{\mrm{rf}}(\theta_{\ell-1},\varphi_{\ell-1},\theta_0,\varphi_0)$
$\hat{\varphi}_\ell \gets \varphi_{\ell-1} + \eta_\ell \hat \nabla_{\varphi} \cL_{\mrm{rf}}(\theta_{\ell-1},\varphi_{\ell-1},\theta_0,\varphi_0)$
$\theta_{\ell+1} \gets \mrm{Proj}_{B_f}(\hat{\theta}_{\ell})$
$\varphi_{\ell+1} \gets \mrm{Proj}_{B_g}(\hat{\varphi}_{\ell})$

Modified randomized prior algorithm for approximate inference.

For the purpose of the analysis we consider the standard SGDA algorithm as outlined in Algorithm <ref>.
In the algorithm $\cL_{\mrm{rf}}$ denotes the objective in (<ref>),
and $\mrm{Proj}_B$ denotes the projection into the $\ell_2$-norm ball with radius $B$,
and $\hat \nabla \cL_{\mrm{rf}}$ represents a stochastic (unbiased) approximation of the gradient $\nabla \cL_{\mrm{rf}}$.
In the following, we will suppress the dependency of $\cL_{\mrm{rf}}$ on $\theta_0, \varphi_0$ for simplicity.

we introduce the notations

\Phi_f := \frac{1}{\sqrt{m}}\begin{bmatrix}
\phi_{x,m}(x_1)^\top \\
\vdots \\
\phi_{x,m}(x_n)^\top
\end{bmatrix}\in\mb{R}^{n\times m},
\quad
\Phi_g := \frac{1}{\sqrt{m}}\begin{bmatrix}
\phi_{z,m}(z_1)^\top \\
\vdots \\
\phi_{z,m}(z_n)^\top
\end{bmatrix}\in\mb{R}^{n\times m},

where we recall $X := (x_1,\ldots,x_n)$ and $Z := (z_1,\ldots,z_n)$ are the training data.

Observe that $\Phi_f \theta = f(X;\theta),\Phi_g \varphi = g(Z;\varphi)$, we can rewrite the objective (<ref>) as

\begin{align}
\cL_{\mrm{rf}}(\theta, \varphi) &=
\theta^\top \Phi_f^\top \Phi_g\varphi -\tilde{Y}^\top \Phi_g \varphi - \frac{1}{2} \varphi^\top \Phi_g^\top \Phi_g \varphi - \frac{\nu}{2}\|\varphi - \varphi_0\|_2^2 + \frac{\lambda}{2}\|\theta - \theta_0\|_2^2.
\end{align}

We additionally define

\begin{align*}
\cL_i(\theta, \varphi) = n\left(\theta^\top \Phi_f^\top E_i \Phi_g\varphi -\tilde{Y}^\top E_i \Phi_g \varphi - \frac{1}{2} \varphi^\top \Phi_g^\top E_i \Phi_g \varphi \right)- \frac{\nu}{2}\|\varphi - \varphi_0\|_2^2 + \frac{\lambda}{2}\|\theta - \theta_0\|_2^2,
\end{align*}

where $E_i := e_ie_i^\top$ and $\{e_i\}_{i\in[n]}$ is the standard orthogonal basis of $\mathbb{R}^n$.
We can see that
\cL_{\mrm{rf}}(\theta, \varphi) = \frac{1}{n} \sum_{i\in[n]} \cL_{i}(\theta, \varphi).
Therefore, the stochastic gradient in Algorithm <ref> can be defined as

\begin{equation}
\label{eqn:stochastic-gradient-rf}
\hat \nabla \cL_{\mrm{rf}}(\theta, \varphi) :=
\nabla \cL_\cI (\theta, \varphi)
=\sum_{i\in[n]} \nabla \cL_i(\theta, \varphi) \mathbf{1}_{i = \cI},
\end{equation}

where $\cI$ is a random variable sampled from the uniform distribution of the set $[n]$.

In practice we run the algorithm concurrently on $J$ sets of parameters, starting from independent draws of initial conditions $\{\theta^{(j)}_0,\phi^{(j)}_0\}$; moreover, the projection is not implemented, and there are various other modifications to further improve stability, as described in Appendix <ref>.

The following lemma is a convergence theorem of Algorithm <ref> under the choice of stochastic gradient defined in (<ref>).

Fix an $m\in\mb{N}$.
Denote by $\theta^*$ the optima of (<ref>) and take $\eta_\ell := \frac{1}{\mu(\ell + 1)}$ with $\mu = \min\{ \lambda, \nu \}$.
Then for any $\epsilon,B_1,B_2,B_3>0$, there exist $B_f, B_g >0$
such that when $L = \Omega( \delta^{-1} \epsilon^{-2})$,
the approximate optima $\theta_{L}$ returned by Algorithm <ref> satisfies

\PP\left (\{\|\theta_L - \theta^*\|_2 > \epsilon\}
\cap E_n\right ) \leq \delta,

E_n :=
\left \{\|\theta_0\|_2 +  \|\varphi_0\|_2\le B_1, \|\tilde{Y}\|_2\le B_2, \sup_{z\in \cZ} \tilde{k}_{z,m}(z,z) +   \sup_{x\in\cX} \tilde{k}_{x,m}(x,x)\le B_3
\right \},

and $\tilde{k}_{\cdot,m}$ denotes the random feature-approximated kernel. The randomness in the statement above is from the sampling of the initial values $\theta_0,\varphi_0$, the gradient noise.

Recall from (<ref>) that
$\theta^*$ is a sum of bounded linear transforms of $\theta_0,\varphi_0$ and $\tilde{Y}_0$. Thus on the event $E_n$,
$\|\theta^*\|_2$ is bounded.
Similarly, $\|\varphi^*\|_2$ is also bounded on $E_n$ by (<ref>).
We choose $B_f$ and $B_g$ to be their maximum values on the event $E_n$.

$\cL_{\mrm{rf}}$ is strongly-convex in $\theta$, and strongly-concave in $\varphi$, it has the unique stationary point $(\theta^*, \varphi^*)$.
We will then bound
\|\theta_\ell - \theta^{*}\|_2^2 + \|\varphi_\ell - \varphi^*\|_2^2
Let $\sigma_f,\sigma_g$ be the minimal constants such that $\|\nabla_{\theta} \cL_{i}(\theta, \varphi)\|_2^2 \leq \sigma_f^2$, $\|\nabla_{\varphi} \cL_{i}(\theta, \varphi)\|_2^2 \leq \sigma_g^2$ for all $i \in [n], \| \theta \|_2 \leq B_f$ and $\| \varphi \|_2 \leq B_g$. Denote
$B:=\max\{B_f,B_g\}$, so we have $\|\theta\|_2, \|\varphi\|_2 \leq B$.

\begin{align*}
r_\ell = \mathbb{E} \left[\|\theta_\ell - \theta^*\|_2^2 + \|\varphi_\ell - \varphi^*\|_2^2\right].
\end{align*}

We want to know how $r_\ell$ contracts.
We first make a stochastic gradient step on $\theta_\ell$ with step size $\eta_\ell$, i.e., $\hat \theta_{\ell + 1} := \theta_\ell - \eta_\ell \hat \nabla_\theta \cL_{\mrm{rf}}(\theta_\ell, \varphi_\ell)$ with $\hat\nabla \cL_{\mrm{rf}}$ defined in (<ref>).

\begin{align*}
\mathbb{E}[\|\hat{\theta}_{\ell+1} - \theta^*\|_2^2 ~|~ \theta_{\ell}, \varphi_{\ell}]
\leq  \|\theta_{\ell} - \theta^*\|_2^2 - 2\eta_\ell\langle \theta_{\ell} - \theta^*, \nabla_{\theta} \cL(\theta_{\ell}, \varphi_{\ell})\rangle + \eta_\ell^2 \sigma_f^2,
\end{align*}

where the expectation is taken with respect to the randomness of the gradient.
For the above inner product term, we have that

\begin{align*}
\langle \theta_{\ell} - \theta^*, \nabla_{\theta} \cL_{\mrm{rf}}(\theta_{\ell}, \varphi_{\ell})\rangle
&= \langle \theta_{\ell} - \theta^*, \nabla_{\theta} \cL_{\mrm{rf}}(\theta_{\ell}, \varphi_{\ell}) - \nabla_{\theta} \cL_{\mrm{rf}}(\theta^*, \varphi^*)\rangle\\
&= \lambda\|\theta_{\ell} - \theta^*\|_2^2 + \langle \theta_{\ell} - \theta^*, \Phi_f^\top \Phi_g (\varphi_{\ell} - \varphi^*) \rangle.
\end{align*}

Next, we consider the gradient step on $\varphi_\ell$ with step size $\eta_\ell$, i.e.,
$\hat \varphi_{\ell + 1} := \varphi_\ell + \eta_\ell \hat \nabla_\varphi \cL_{\mrm{rf}}(\theta_\ell, \varphi_\ell)$.
Then, we have that

\begin{align*}
\mathbb{E}[\|\hat{\varphi}_{\ell + 1} - \varphi^*\|_2^2 ~|~ \theta_{\ell}, \varphi_{\ell}]
\leq  \|\varphi_{\ell} - \varphi^*\|_2^2 + 2\eta_\ell \langle \varphi_{\ell} - \varphi^*, \nabla_{\varphi} \cL_{\mrm{rf}}(\theta_{\ell}, \varphi_{\ell}) \rangle + \eta_\ell^2\sigma_g^2.
\end{align*}

We similarly deal with the inner product term:

\begin{align*}
\langle \varphi_{\ell} - \varphi^*, \nabla_{\varphi} \cL_{\mrm{rf}}({\theta}_{\ell}, \varphi_{\ell}) \rangle
&= \langle \varphi_{\ell} - \varphi^*, \nabla_{\varphi} \cL_{\mrm{rf}}({\theta}_{\ell}, \varphi_{\ell}) - \nabla_{\varphi} \cL_{\mrm{rf}}(\theta^*, \varphi^*)\rangle\\
&= -\langle \varphi_{\ell} - \varphi^*, (\Phi_g^\top \Phi_g + \nu I)(\varphi_{\ell} - \varphi^*)\rangle + \langle \varphi_{\ell} - \varphi^*, \Phi_g^\top \Phi_f({\theta}_{\ell} - \theta^*) \rangle \\
&\leq -\nu\| \varphi_{\ell} - \varphi^* \|^2_2 + \langle \varphi_{\ell} - \varphi^*, \Phi_g^\top \Phi_f({\theta}_{\ell} - \theta^*) \rangle,
\end{align*}

Combining the above results, we have

\begin{align*}
r_{\ell + 1} \leq \EE [\|\hat{\theta}_{\ell+1}-\theta^*\|_2^2+\|\hat{\varphi}_{\ell+1}-\varphi^*\|_2^2\mid \theta_\ell,\varphi_\ell]
\le (1 - 2 \mu \eta_\ell) r_\ell
+ \eta_\ell^2(\sigma_f^2 + \sigma_g^2),
\end{align*}

where we have set $\mu := \min\{ \nu, \lambda \}$, and the first inequality follows from the fact that the projection onto a convex set is a contraction map, i.e., $\| \mrm{Proj}_B(x) - \mrm{Proj}_B(y) \| \leq \|x - y \|$.

Let $\sigma^2 = \sigma_f^2 + \sigma_g^2$ and $\eta_\ell = \frac{\xi}{\ell+1}$ for some $ \xi > \frac{1}{2\mu}$, by induction we have

\begin{align*}
r_{\ell} \leq \frac{c_{\xi}}{\ell + 1}, \quad \text{where }
c_{\xi} = \max\left \{ r_0, \frac{2\xi^2 \sigma^2}{2\mu\xi - 1}\right \}.
\end{align*}

Specifically, taking $\xi = \mu^{-1}$, we have

\begin{equation}
\label{eqn:r_ell-contraction}
r_\ell \leq \frac{1}{\ell + 1} \max\left \{ r_0, \frac{2\sigma^2}{\mu^2} \right \}.
\end{equation}

We now track the constants we have used in (<ref>).
Note that on the event $E_n$,

\[ r_0
\leq 2 \left(\|\theta_0\|_2^2 + \|\theta^*\|_2^2 + \|\varphi_0\|_2^2 + \|\varphi^*\|_2^2\right)
\leq 4(B_1^2 + B^2). \]

Recall the definition of $\sigma^2$:

\begin{align*}
\sigma^2= \max_{i\in [n], \|\theta\|_2, \|\varphi\|_2\leq B} \|\nabla_{\theta}\cL_i(\theta, \varphi)\|_2^2 + \max_{i\in [n], \|\theta\|_2, \|\varphi\|_2\leq B} \|\nabla_{\varphi} \cL_i(\theta, \varphi)\|_2^2 =: (\mrm I) + (\mrm{II}).
\end{align*}

For the first term, we have

\begin{align*}
(\mrm I)
& =  \max_{i\in [n], \|\theta\|_2, \|\varphi\|_2\leq B} \|\lambda (\theta - \theta_0) + n \Phi_f^\top E_i \Phi_g\varphi \|_2^2 \\
& \leq \max_{i\in[n], \|\theta\|_2, \|\varphi\|_2\leq B} \left(2 \lambda^2 \|\theta - \theta_0\|_2^2 + 2n^2 \|\Phi_f^\top E_i \Phi_g\varphi\|_2^2\right)
\leq  4\lambda^2 (B^2 + B_1^2) + 2n^2B_3^2B^2.
\end{align*}

Similarly, for the second term, we have

\begin{align*}
(\mrm{II})  %
& =  \max_{i\in [n], \|\theta\|_2, \|\varphi\|_2\leq B}\|n(\theta^\top \Phi_f^\top E_i \Phi_g - \tilde{Y}^\top E_i\Phi_g - \Phi_g^\top E_i \Phi_g\varphi) - \nu(\varphi - \varphi_0)\|_2^2 \\
&\leq 4n^2 B_3^2B^2 + 2n^2 B_2^2 B^2 + 4\nu^2(B^2 + B_1^2).
\end{align*}

Thus, we know that

\begin{align*}
\sigma^2 \leq 8(\lambda^2 + \nu^2) (B^2 + B_1^2) + 6 n^2 B_3^2B^2 + 2n^2 B_2^2B^2 =: \tilde C.
\end{align*}

Taking $L_\delta = \delta^{-1} \epsilon^{-2} \max\{4B_1^2 + 4B^2, \tilde C \mu^{-1} \}$ and $\eta_{ \ell} = \frac{1}{\mu(\ell + 1)}$, by (<ref>), we know that

\begin{align*}
\mathbb{P}(\|\theta_{L} - \theta^*\|_2 > \epsilon)
\leq \epsilon^{-2} \EE \| \theta_L - \theta^* \|_2^2
\leq \epsilon^{-2} r_\ell
\leq \delta.
\end{align*}

§.§ Proof of Proposition <ref>

By Lemma <ref>, for any $\epsilon_1>0$ we have

\begin{equation}\label{eq:approxinf-0}
\lim_{m\rightarrow\infty}
\PP\left (
\left \{\sup_{x^*\in\cX^l}\|\tilde{\mu}-\mu\|_2 > \epsilon_1 \right \}\cup
\left \{\sup_{x^*\in\cX^l}\|\tilde{S}-S\|_F > \epsilon_1     \right \}
\right ) = 0,
\end{equation}

where the randomness is from the sampling of random feature bases.

Fix an arbitrary set of $\epsilon_1>0,\delta_0>0$. Then we can find
$m\in\mb{N}$ such that the event in (<ref>) has probability smaller than $\delta_0$.
Combining Assumption <ref> with the fact that $\theta_0,\phi_0,\tilde{Y}_0$ are now Gaussian random variables with fixed dimensionality,
for any $\delta_1>0$,
we can choose $B_1,B_2,B_3$ such that
the event $E_n$ defined in Lemma <ref> has probability $1-\delta_1$.
Thus for any $\epsilon_2>0$, when the number of iteration steps exceeds $\Omega( \delta_1^{-1} \epsilon^{-2}_2)$, we have

\begin{equation}\label{eq:approxinf-2}
\PP(\|\hat{\theta}_m-\theta^*_m\|_2 > \epsilon_2)
\leq \PP( \{ \|\hat{\theta}_m-\theta^*_m\|_2 > \epsilon_2 \} \cap E_n) + \PP(E_n^c)
\leq 2\delta_1,
\end{equation}

where $\hat{\theta}_m$ denotes the approximate optima returned by Algorithm <ref> after $\Omega( \delta_1^{-1} \epsilon^{-2}_2)$ iterations, $\theta^*_m$ denotes the exact optima of the minimax objective, and the randomness is from the gradient noise as well as the perturbations $f_0,g_0,\tilde{Y}$.
Thus we have

\EE \|\hat{\theta}_m-\theta^*_m\|_2
\leq \epsilon_2 + 2\delta_1 ( \EE \| \hat \theta_m \|_2 + \EE \| \theta^*_m \|_2)
\leq \epsilon_2 + 4\delta_1 B.

From the choice of $B$ in Lemma <ref>,
we can see that $
\delta_1 B \le \EE(\|\theta^*_m\| \cdot (1-\mbf{1}_{E_n}))$, and thus converges to $0$ as $\delta_1\rightarrow 0$.
$\EE \|\hat{\theta}_m-\theta^*_m\|_2 $ converges to 0, and for any $x^*\in\cX^l$,

\[ \begin{aligned}
\EE \sup_{x^*\in\cX^l}\|f(x^*;\hat{\theta}_m)-f(x^*;\theta^*_m)\|_2
&= \EE \sup_{x^*\in\cX^l}\|\phi_{x,m}(x^*)^\top (\hat{\theta}_m-\theta^*_m)\|_2 \\
&\le l \sqrt{\tilde{\kappa}} \cdot \EE \|\hat{\theta}_m-\theta^*_m\|_2
\rightarrow 0,
\end{aligned} \]

where the expectation is taken with respect to the gradient noise, perturbations, and random feature draws.
Hence, the mean and covariance of $f(x^*;\hat{\theta}_m)$ converges to that of $f(x^*;\theta^*_m)$ as intended,
and we know that  the following holds with probability at least $1 - \delta_0$

\sup_{x^*\in\cX^l} \max\left\{
\|\EE(f(x^*;\hat{\theta}_m)) - \EE(f(x^*;\theta_m)) \|_2,
\|\mrm{Cov}(f(x^*;\hat{\theta}_m)) - \mrm{Cov}(f(x^*;\theta_m)) \|_F
\right\} \leq \epsilon_1

Combining this with (<ref>) completes the proof.

§ IMPLEMENTATION DETAILS, EXPERIMENT SETUP AND ADDITIONAL RESULTS

§.§ Hyperparameter Selection

We follow the strategy in previous work <cit.> and select hyperparameters by minimizing the observable first or second stage loss, depending on which part they directly correspond to.

For the first stage, the loss is

\cL_{v1} = \mrm{Tr}(K_{xx} - 2K_{x\tilde{x}} L + K_{\tilde x\tilde x}L^\top L) = \EE_{f\sim\mc{GP}(0,k)} \|f(X) - L f(\tilde{X})\|_2^2

where $L := K_{z \tilde z}(K_{\tilde{z}\tilde{z}}+\nu I)^{-1}$, and tilde indicates the held-out data.
From the above equality we can see that a Monte-Carlo estimator for $L_1$ can be constructed with the following procedure:

* Draw $f\sim\mc{GP}(0,k_x)$.

* Perform kernel ridge regression on the dataset $\{(\tilde z_i, f(\tilde x_i))\}$.

* Return the mean squared error on the dataset $(X,Z)$.

This procedure can also be implemented for the NN-based models.

For the second stage, the loss $\sum_{i=1}^n \hat{d}_n(\hat{E}_n f, \hat{b})$ can be computed directly, for both the closed-form quasi-posterior and the random feature approximation.
For the approximate inference algorithm, as we can see from (<ref>) that the dual functions $\{g(\cdot;\varphi^{(k)})\}$ are samples from Gaussian process posteriors centered at the needed point estimates $\hat{E}_n f(\cdot;\theta^{(k)})$, instead of the point estimates themselves,
we train separate validator models to approximate the latter. The validator models have the architecture to the dual functions used for training, and follow the same learning rate schedule. The validator models are trained before estimating the validation statistics, and we run SGD until convergence to ensure an accurate estimate.

§.§ Details in the Approximate Inference Algorithm

To draw multiple samples from the quasi-posterior efficiently, our algorithm runs $J$ SGDA chains in parallel, with different perturbations $
\{(\tilde{Y}^{(j)}, f_0^{(j)}, g_0^{(j)}): j\in [J]\}
While the convergence analysis works with the extremely simple Algorithm <ref>, in practice we extend it to improve stability and accelerate convergence:

* we employ early stopping based on the validation statistics;

* before the main optimization loop we initialize the dual parameters at the approximate optima $\arg\min_\varphi L_{\mrm{rf}}(f^{(j)}, g(\cdot;\varphi))$, by running SGD until convergence;

* in each SGDA iteration, we use $K_1>1$ GD steps on $g$ and one GA step for $f$;

* after every $K_2$ epochs, we fix $\theta^{(j)}$ and train the dual parameters $\varphi^{(j)}$ for one epoch.

All the above choices are
shown to improve the observable validation statistics. We fix $K_1=3,K_2=2$ which are determined on the 1D datasets using the validation statistics.

§.§ Nyström Approximation for the Closed-form Quasi-Posterior

Recall the definition of the Nyström-approximated quasi-posterior: we change the definition of the quasi-likelihood to Assumption <ref>, where $\cI_n := \mrm{span}\{k(z_{i_j},\cdot): j\in [m]\}$, and
$z_{i_1},\ldots,z_{i_m}$ are inducing points sampled from the training data. The new closed-form expressions are then

\begin{align}
\Pi(f(x_*)\givendata) &= \cN(
K_{*x} \tilde\Lambda Y,
K_{**} - K_{*x} \tilde\Lambda K_{x*}),  \\
\text{where}~~ \tilde \Lambda &:=
(\lambda I + K_{xx}\tilde L)^{-1}\tilde L, \\
\tilde L &:= K_{z\tilde z}(\nu K_{\tilde z\tilde z} + K_{\tilde z z}K_{z\tilde z})^{-1}K_{\tilde z z},
\end{align}

and $\tilde z$ denotes the inducing points.
In the implementation, we compute a low-rank factorization $A^\top A := \tilde L$, defined through the Cholesky factorization of $\effnu K_{\tilde z\tilde z}+K_{\tilde zz}K_{z\tilde z}$, and use the following equivalent expression of $\tilde\Lambda$ obtained from the Woodbury identity:

\tilde \Lambda = \lambda^{-1} A^\top (I - (\lambda I + A K_{xx}A^\top)^{-1} A K_{xx}A^\top) A =
A^\top U \mrm{diag}\!\left(
(\lambda + \gamma_1)^{-1}, \ldots, (\lambda+\gamma_n)^{-1}
\right) U^\top A,

where $U \mrm{diag}(\gamma_1,\ldots,\gamma_n) U^\top = AK_{xx}A^\top$ denotes the eigendecomposition. This ensures that matrix inversion is only performed on $m\times m$ matrices.

We use uniform sampling for inducing points, which is equivalent to the leverage score sampling scheme studied in [Rudi et al., 2016] as we work on the symmetric $\cU[\mb{S}^1]$.

§.§ Experiment in Section <ref>: Additional Results

Additional results are plotted in Figure <ref>. For the contraction rate of the derivative, observe the power space of $W^{s,2}(\mb{T})$ matches lower-order Sobolev spaces, so the contraction rate follows that of the power space, as in Section <ref>.

Additional results.
The left three plots correspond to the estimation of $f$; the rightmost plot corresponds to the estimation of $f'$. Dashed line indicates the theoretical rate.

§.§ 1D Simulation: Experiment Setup Details

In constructing the datasets, let $\tilde f_0$ denote the sine, step, abs or linear ($\tilde f_0(x)=x$) function; we then set $f_0 = \tilde f_0(4\cdot(2x-1))$ if $\tilde f_0$ is sine, abs or linear, $\mbf{1}_{\{2x-1<0\}}+2.5\cdot \mbf{1}_{\{2x-1\ge 0\}}$ otherwise.
These choices are made to maintain consistency with previous work [Lewis and Syrgkanis, 2018, Bennett et al., 2019], which used the same transformed step function and defined $\bx$ so that it has a range of approximately $[-4,4]$.

For 2SLS and the kernelized IV methods,
we determine $\lambda$ and $\nu$ following <ref>. To improve stability, we repeat the procedures on $50$ random partitions of the combined training and validation set, and choose the hyperparameters that minimize the average loss. The hyperparameters are chosen from a log-uniform grid consisting of $10$ values in the range of $[0.1, 30]$. The occasional instability of hyperparameter selection is also reported in [Muandet et al., 2020].
For BayesIV, we run the MCMC sampler for 25000 iterations, discard the first
5000 iterations for burn in, and take one sample out of every 80 consecutive iterations to construct the approximate posterior.
For bootstrap we use 20 samples. In both cases, further increasing the computational budget will not improve performance.

We normalize the dataset to have zero mean and unit variance. For all kernel methods we set the kernel bandwidth using the median trick.

§.§ 1D Simulation: Full Results

Normalized MSE, average CI coverage and CI width on test set are provided in Table <ref>-<ref>.
We can see that the quasi-posterior consistently provides the most reliable coverage, and
is most advantageous in the small-sample setting, or when instrument strength is weak.
All methods suffer from lower coverage on the  design, due to its discontinuity, and in the weak IV setting, due to instabilities in the hyperparameter selection process.

We provide the following visualizations: (i) Figure <ref> plots the quasi-posterior on the  design, for the Matérn and RBF kernels; we also include the approximate inference algorithm in the plot.
(ii) Figure <ref> compares the quasi-posterior and bootstrap in the weak instrument setting.
The hyperparameters $\lambda,\nu$ are set to the median of all values selected in the repeated experiments, and hyperparameters in the approximate inference algorithm are set as in the demand experiment below.
Additional visualizations are provided in the conference version of this work.

Matérn-$3/2$ kernel

Matérn-$5/2$ kernel

RBF, closed form

RBF, RF approx

1D datasets: visualization of the quasi-posterior on the  design using various models. Dots indicate training data. We fix $N=1000,\alpha=0.5$.

1D datasets: visualization of the weak instrument setting, on the  design with the Matérn$-3/2$ kernel. Results for other kernels are qualitatively similar.

§.§ Demand Simulation: Experiment Setup Details

All variables in the dataset are normalized to have zero mean and unit variance.
For BayesIV, we run the MCMC sampler for 50000 iterations, discard the first 10000 samples as burn in, and take every 80th sample for inference.
For the kernelized methods, hyperparamater selection follows the 1D experiments.
For the NN-based methods, implementation details are discussed in Appendix <ref>;
for both our method and bootstrap, we draw $J=10$ samples from the predictive distribution.

We select hyperparameters by applying the procedure in Appendix <ref> to a fixed train / validation split, since on this dataset we observe little variation in its results. Hyperparameters include $\lambda,\nu$, and the learning rate schedule (initial learning rate $\eta_0$ and period of learning rate decay $\tau$). The learning rate is adjusted by multiplying it by a factor of $0.8$ every $\tau$ iterations.
We fix the optimizer to Adam, and train until validation statistics no longer improves.

For the lower-dimensional setup, we select
$\lambda$ and $\nu$ from a log-linearly scaled grid of 10 values, with the range of $[5\times 10^{-3}, 5]$ and $[0.05, 1]$, respectively. The ranges are chosen based on preliminary experiments using the range of $[0.1, 30]$. We determine $\eta$ from $\{5\times 10^{-4}, 10^{-3}, 5\times 10^{-3}, 1\times 10^{-2}, 5\times 10^{-2}\}$, and $\tau$ from $\{80,160,320,640\}$. We fix the batch size at $256$. The NN architecture consists of two fully-connected layers, with 50 hidden units and the tanh activation. We also experimented with NNs with 3 hidden layers or with ReLU activation, and made this choice based on the validation statistics.

For the image-based setup, the range of $\lambda$ and $\eta$ follows the above. For $\nu,\tau$ we consider $
\nu\in [1,100], \tau\in\{640,1280,2560,5120\},
$ based on preliminary experiments. We fix the batch size at $80$.
The network architecture is adapted from [Hartford et al., 2017], and consists of two $3\times 3$ convolutional layers with $64$ filters, followed by max pooling, dropout, and three fully-connected layers with $64,32$ and $1$ units.

Following the setup in all previous work, we use a uniform grid on $[5,30]\times [0,10]\times \{0,\ldots,6\}$ as the test set.

§.§ Demand Simulation: Full Results and Visualizations

Results in the large-sample settings are presented in Table <ref>,
which are consistent with the discussion in the main text.

We plot the predictive distributions for all methods in Figure <ref>, on the same cross-section as in the main text, for $N=1000$.
(We omit the plot for $N=10^4$ and the image experiment, since in those settings bootstrap and the quasi-posterior have similar behaviors.) As we can see,
all non-NN baselines except BayesIV produce overly smooth predictions, presumably due to the lack of flexibility in these models.
Note that the visualizations only correspond to an intersection of the true function $f(x_0,t_0,s)$, with $x_0,t_0$ fixed; the complete function has the form of $x\cdot s\cdot \psi(t)$, ignoring the less significant terms, and thus may incur a large norm penalty in the less flexible RKHSes.
The issue is further exacerbated by the discrepancy between the training and test distributions: the former is non-uniform due to confounding. As we can see from Figure <ref>, in the region where $t$ is close to $5$, the data is scarce for most values of $x$, which may explain the reason that the closed-form kernels fail to provide good coverage around $t=5$ (and $s=3,x=17.5$, as used in the visualizations), and the reason that both NN-based methods assign higher uncertainty around this location.

BayesIV has a different failure mode: as it employs additive regression models for both stages $p(\bx\mid\bz), p(\by\mid\bx)$, it approximates this cross-section relatively well. However, as the true structural function does not have an additive decomposition, its prediction in other regions can be grossly inaccurate;
we plot one such cross-section in Figure <ref>(a).

When implemented with the NN model, bootstrap CIs are sharper in regions with more training data, although the difference is often insignificant. The difference in out-of-distribution regions is more significant, where bootstrap is often less robust. An example is shown in Figure <ref>.

Setting     3cLow-dimensional, $N=10^4$     3cImage, $N=5\times 10^4$

Method     BS-2SLS     BS-NN     QB-NN      BS-2SLS     BS-NN     QB-NN

CI Cvg.

CI Wid.
$.014\pm .002$
$.136\pm .015$
$.203\pm .013$
$.132\pm .039$
$.636\pm .027$
$.597\pm .024$

Deferred results on the demand design. Results are averaged over 20 trials for the low-dimensional experiment, and 10 trials for the image experiment.

Demand experiment: scatter plot of $10^4$ samples from the training data distribution $p(x,t\mid s=4)$. The dashed line indicates the cross-section used in Figure <ref>.

Demand experiment: visualizations of the predictive distributions for $N=1000$, on the same cross-section as in Figure <ref>.

Demand experiment: visualizations of the predictive distributions for $N=1000$ on a out-of-distribution cross-section, obtained by fixing $t=9,s=6$ and varying $x$.

§.§.§ Computational cost
We report the typical training time for a single set of hyperparameters, excluding JIT compilation time, on a GeForce GTX 1080Ti GPU.
In the lower-dimensional experiments, training takes around 25 minutes for a single set of hyperparameters when $N=10^3$, or around 30 minutes when $N=10^4$; in both cases 6 experiments can be carried out in parallel on a single GPU.
In the image experiment, training takes around 7.5 hours.

The time cost above is for the optimal hyperparameter configuration; experiments using suboptimal hyperparameters usually take a shorter period of time due to early stopping.
It can also be improved by switching to low-precision numerical operations, or with various heuristics in the hyperparameter search (e.g., using a smaller $J$ in an initial search).

[Adams and Fournier, 2003]
Robert A Adams and John JF Fournier.
Sobolev Spaces.
Elsevier, 2003.

[Agapiou et al., 2013]
Sergios Agapiou, Stig Larsson, and Andrew M Stuart.
Posterior contraction rates for the bayesian approach to linear
ill-posed inverse problems.
Stochastic Processes and their Applications, 1230
(10):0 3828–3860, 2013.

[Angrist and Pischke, 2008]
Joshua D Angrist and Jörn-Steffen Pischke.
Mostly Harmless Econometrics: An Empiricist's Companion.
Princeton University Press, 2008.

[Bennett and Kallus, 2020]
Andrew Bennett and Nathan Kallus.
The variational method of moments.
arXiv preprint arXiv:2012.09422, 2020.

[Bennett et al., 2019]
Andrew Bennett, Nathan Kallus, and Tobias Schnabel.
Deep generalized method of moments for instrumental variable
In Advances in Neural Information Processing Systems, 2019.

[Blundell et al., 2007]
Richard Blundell, Xiaohong Chen, and Dennis Kristensen.
Semi-nonparametric IV estimation of shape-invariant Engel curves.
Econometrica, 750 (6):0 1613–1669, 2007.

[Borovitskiy et al., 2020]
Viacheslav Borovitskiy, Alexander Terenin, Peter Mostowsky, and Marc Peter
Matérn Gaussian processes on Riemannian manifolds.
In Advances in Neural Information Processing Systems, 2020.

[Cavalier, 2008]
Laurent Cavalier.
Nonparametric statistical inverse problems.
Inverse Problems, 240 (3):0 034004, 2008.

[Chen and Christensen, 2015]
Xiaohong Chen and Timothy M. Christensen.
Optimal sup-norm rates and uniform inference on nonlinear functionals
of nonparametric IV regression.
arXiv preprint arXiv:1508.03365, 2015.

[Chen and Pouzo, 2012]
Xiaohong Chen and Demian Pouzo.
Estimation of nonparametric conditional moment models with possibly
nonsmooth generalized residuals.
Econometrica, 800 (1):0 277–321, 2012.

[Chen and Reiss, 2011]
Xiaohong Chen and Markus Reiss.
On rate optimality for ill-posed inverse problems in econometrics.
Econometric Theory, pages 497–521, 2011.

[Chernozhukov and Hong, 2003]
Victor Chernozhukov and Han Hong.
An MCMC approach to classical estimation.
Journal of Econometrics, 1150 (2):0 293–346,

[Conley et al., 2008]
Timothy G Conley, Christian B Hansen, Robert E McCulloch, and Peter E Rossi.
A semi-parametric Bayesian approach to the instrumental variable
Journal of Econometrics, 1440 (1):0 276–305,

[Cuzick et al., 1997]
Jack Cuzick, Robert Edwards, and Nereo Segnan.
Adjusting for non-compliance and contamination in randomized clinical
Statistics in Medicine, 160 (9):0 1017–1029,

[Dai et al., 2016]
Bo Dai, Niao He, Yunpeng Pan, Byron Boots, and Le Song.
Learning from conditional distributions via dual embeddings.
arXiv preprint arXiv:1607.04579, 2016.

[Davidson and Mackinnon, 2010]
Russell Davidson and James G. Mackinnon.
Wild bootstrap tests for IV regression.
Journal of Business & Economic Statistics, 280
(1):0 128–144, 2010.

[Dikkala et al., 2020]
Nishanth Dikkala, Greg Lewis, Lester Mackey, and Vasilis Syrgkanis.
Minimax estimation of conditional moment models.
In Advances in Neural Information Processing Systems,
volume 33, pages 12248–12262, 2020.

[Dudík et al., 2007]
Miroslav Dudík, Steven J Phillips, and Robert E Schapire.
Maximum entropy density estimation with generalized regularization
and an application to species distribution modeling.
Journal of Machine Learning Research, 2007.

[Eldredge, 2016]
Nathaniel Eldredge.
Analysis and probability on infinite-dimensional spaces.
arXiv preprint arXiv:1607.03591, 2016.

[Fischer and Steinwart, 2020]
Simon Fischer and Ingo Steinwart.
Sobolev norm learning rates for regularized least-squares algorithms.
Journal of Machine Learning Research, 21:0 205–1,

[Florens and Simoni, 2012]
Jean-Pierre Florens and Anna Simoni.
Nonparametric estimation of an instrumental regression: A
quasi-Bayesian approach based on regularized posterior.
Journal of Econometrics, 1700 (2):0 458–475,

[Florens and Simoni, 2021]
Jean-Pierre Florens and Anna Simoni.
Gaussian processes and Bayesian moment estimation.
Journal of Business & Economic Statistics, 390
(2):0 482–492, 2021.

[Flores-Lagunes, 2007]
Alfonso Flores-Lagunes.
Finite sample evidence of IV estimators under weak instruments.
Journal of Applied Econometrics, 220 (3):0
677–694, 2007.

[Ghosal and Van der Vaart, 2017]
Subhashis Ghosal and Aad Van der Vaart.
Fundamentals of Nonparametric Bayesian Inference, volume 44.
Cambridge University Press, 2017.

[Ghosal et al., 2000]
Subhashis Ghosal, Jayanta K Ghosh, and Aad W Van Der Vaart.
Convergence rates of posterior distributions.
Annals of Statistics, pages 500–531, 2000.

[Giné and Nickl, 2021]
Evarist Giné and Richard Nickl.
Mathematical Foundations of Infinite-Dimensional Statistical
Cambridge University Press, 2021.

[Greenland, 2000]
Sander Greenland.
An introduction to instrumental variables for epidemiologists.
International Journal of Epidemiology, 290
(4):0 722–729, 2000.

[Hartford et al., 2017]
Jason Hartford, Greg Lewis, Kevin Leyton-Brown, and Matt Taddy.
Deep IV: A flexible approach for counterfactual prediction.
In International Conference on Machine Learning, pages
1414–1423, 2017.

[He et al., 2020]
Bobby He, Balaji Lakshminarayanan, and Yee Whye Teh.
Bayesian deep ensembles via the neural tangent kernel.
arXiv preprint arXiv:2007.05864, 2020.

[Horowitz, 2011]
Joel L. Horowitz.
Applied nonparametric instrumental variables estimation.
Econometrica, 790 (2):0 347–394, 2011.

[Horowitz, 2012]
Joel L. Horowitz.
Specification testing in nonparametric instrumental variable
Journal of Econometrics, 1670 (2):0 383–396,

[Jaakkola et al., 1999]
Tommi Jaakkola, Marina Meila, and Tony Jebara.
Maximum entropy discrimination.
In Proceedings of the 12th International Conference on Neural
Information Processing Systems, pages 470–476, 1999.

[Jacot et al., 2018]
Arthur Jacot, Franck Gabriel, and Clement Hongler.
Neural tangent kernel: Convergence and generalization in neural
In Advances in Neural Information Processing Systems,

[Kanagawa et al., 2018]
Motonobu Kanagawa, Philipp Hennig, Dino Sejdinovic, and Bharath K
Gaussian processes and kernel methods: A review on connections and
arXiv preprint arXiv:1807.02582, 2018.

[Kato, 2013]
Kengo Kato.
Quasi-Bayesian analysis of nonparametric instrumental variables
The Annals of Statistics, 410 (5):0
2359–2390, 2013.

[Kim, 2002]
Jae-Young Kim.
Limited information likelihood and Bayesian analysis.
Journal of Econometrics, page 19, 2002.

[Kleibergen and Van Dijk, 1998]
Frank Kleibergen and Herman K Van Dijk.
Bayesian simultaneous equations analysis using reduced rank
Econometric Theory, pages 701–743, 1998.

[Kleibergen and Zivot, 2003]
Frank Kleibergen and Eric Zivot.
Bayesian and classical approaches to instrumental variable
Journal of Econometrics, 1140 (1):0 29–72,

[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, 390 (5):0
2626–2657, 2011.

[Knapik and Salomond, 2018]
Bartek Knapik and Jean-Bernard Salomond.
A general approach to posterior contraction in nonparametric inverse
Bernoulli, 240 (3), 2018.

[Lee et al., 2019]
Jaehoon Lee, Lechao Xiao, Samuel Schoenholz, Yasaman Bahri, Roman Novak, Jascha
Sohl-Dickstein, and Jeffrey Pennington.
Wide neural networks of any depth evolve as linear models under
gradient descent.
In Advances in Neural Information Processing Systems,
volume 32, 2019.

[Lewis and Syrgkanis, 2018]
Greg Lewis and Vasilis Syrgkanis.
Adversarial generalized method of moments.
arXiv preprint arXiv:1803.07164, 2018.

[Liao et al., 2020]
Luofeng Liao, You-Lin Chen, Zhuoran Yang, Bo Dai, Zhaoran Wang, and Mladen
Provably efficient neural estimation of structural equation model: An
adversarial approach.
arXiv preprint arXiv:2007.01290, 2020.

[Liao and Jiang, 2011]
Yuan Liao and Wenxin Jiang.
Posterior consistency of nonparametric conditional moment restricted
The Annals of Statistics, 390 (6):0
3003–3031, 2011.

[Lopes and Polson, 2014]
Hedibert F Lopes and Nicholas G Polson.
Bayesian instrumental variables: priors and likelihoods.
Econometric Reviews, 330 (1-4):0 100–121,

[Mastouri et al., 2021]
Afsaneh Mastouri, Yuchen Zhu, Limor Gultchin, Anna Korba, Ricardo Silva, Matt J
Kusner, Arthur Gretton, and Krikamol Muandet.
Proximal causal learning with kernels: Two-stage estimation and
moment restriction.
arXiv preprint arXiv:2105.04544, 2021.

[Miao et al., 2018]
Wang Miao, Zhi Geng, and Eric J Tchetgen Tchetgen.
Identifying causal effects with proxy variables of an unmeasured
Biometrika, 1050 (4):0 987–993, 2018.

[Moreira et al., 2004]
Marcelo Moreira, Jack R Porter, and Gustavo A Suarez.
Bootstrap and higher-order expansion validity when instruments may be
weak, 2004.

[Muandet et al., 2020]
Krikamol Muandet, Arash Mehrjou, Si Kai Lee, and Anant Raj.
Dual instrumental variable regression.
arXiv preprint arXiv:1910.12358, 2020.

[Neal, 2012]
Radford M Neal.
Bayesian Learning for Neural Networks, volume 118.
Springer Science & Business Media, 2012.

[Newey and Powell, 2003]
Whitney K Newey and James L Powell.
Instrumental variable estimation of nonparametric models.
Econometrica, 710 (5):0 1565–1578, 2003.

[Osband et al., 2018]
Ian Osband, John Aslanides, and Albin Cassirer.
Randomized prior functions for deep reinforcement learning.
In Proceedings of the 32nd International Conference on Neural
Information Processing Systems, pages 8626–8638, 2018.

[Pearce et al., 2020]
Tim Pearce, Felix Leibfried, and Alexandra Brintrup.
Uncertainty in neural networks: Approximately Bayesian ensembling.
In Proceedings of the 23rd International Conference on
Artificial Intelligence and Statistics, volume 108, pages 234–244, 2020.

[Rahimi et al., 2007]
Ali Rahimi, Benjamin Recht, et al.
Random features for large-scale kernel machines.
In Advances in Neural Information Processing Systems.
Citeseer, 2007.

[Rousseau and Szabo, 2017]
Judith Rousseau and Botond Szabo.
Asymptotic behaviour of the empirical Bayes posteriors associated
to maximum marginal likelihood estimator.
The Annals of Statistics, 450 (2), 2017.

[Rudi et al., 2016]
Alessandro Rudi, Raffaello Camoriano, and Lorenzo Rosasco.
Less is more: Nyström computational regularization.
arXiv preprint arXiv:1507.04717, 2016.

[Singh et al., 2019]
Rahul Singh, Maneesh Sahani, and Arthur Gretton.
Kernel Instrumental variable regression.
arXiv preprint arXiv:1906.00232, 2019.

[Srinivas et al., 2012]
Niranjan Srinivas, Andreas Krause, Sham M. Kakade, and Matthias Seeger.
Gaussian process optimization in the bandit setting: No regret and
experimental design.
IEEE Transactions on Information Theory, 580
(5):0 3250–3265, 2012.

[Sriperumbudur et al., 2011]
Bharath K Sriperumbudur, Kenji Fukumizu, and Gert RG Lanckriet.
Universality, characteristic kernels and RKHS embedding of
Journal of Machine Learning Research, 120 (7), 2011.

[Steinwart, 2019]
Ingo Steinwart.
Convergence types and rates in generic Karhunen-Loève expansions
with applications to sample path properties.
Potential Analysis, 510 (3):0 361–395, 2019.

[Steinwart and Christmann, 2008]
Ingo Steinwart and Andreas Christmann.
Support Vector Machines.
Springer Science & Business Media, 2008.

[Steinwart and Scovel, 2012]
Ingo Steinwart and Clint Scovel.
Mercer’s theorem on general domains: On the interaction between
measures, kernels, and RKHSs.
Constructive Approximation, 350 (3):0
363–417, 2012.

[Steinwart et al., 2009]
Ingo Steinwart, Don R Hush, Clint Scovel, et al.
Optimal rates for regularized least squares regression.
In COLT, pages 79–93, 2009.

[Stock et al., 2002]
James H Stock, Jonathan H Wright, and Motohiro Yogo.
A survey of weak instruments and weak identification in generalized
method of moments.
Journal of Business & Economic Statistics, 200
(4):0 518–529, 2002.

[Strichartz, 1983]
Robert S Strichartz.
Analysis of the Laplacian on the complete Riemannian manifold.
Journal of Functional Analysis, 520 (1):0
48–79, 1983.

[Stuart, 2010]
Andrew M Stuart.
Inverse problems: a Bayesian perspective.
Acta Numerica, 19:0 451–559, 2010.

[Szabó et al., 2015]
Botond Szabó, A. W. van der Vaart, and J. H. van Zanten.
Frequentist coverage of adaptive nonparametric Bayesian credible
The Annals of Statistics, 430 (4):0
1391–1428, 2015.

[Triebel, 1983]
Hans Triebel.
Theory of Function Spaces.
Springer Basel, Basel, 1983.
ISBN 978-3-0346-0415-4 978-3-0346-0416-1.

[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, 360 (3), 2008a.

[van der Vaart and van Zanten, 2008]
A. W. van der Vaart and J. H. van Zanten.
Reproducing kernel Hilbert spaces of Gaussian priors.
arXiv preprint arXiv:0805.3252, pages 200–222,

[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, 12:0 2095–2119,

[Wainwright, 2019]
Martin J Wainwright.
High-Dimensional Statistics: A Non-Asymptotic Viewpoint,
volume 48.
Cambridge University Press, 2019.

[Wang et al., 2021]
Ziyu Wang, Yuhao Zhou, Tongzheng Ren, and Jun Zhu.
Scalable quasi-bayesian inference for instrumental variable
In Advances in Neural Information Processing Systems, 2021.

[Wiesenfarth et al., 2014]
Manuel Wiesenfarth, Carlos Matías Hisgen, Thomas Kneib, and Carmen
Bayesian nonparametric instrumental variables regression based on
penalized splines and dirichlet process mixtures.
Journal of Business & Economic Statistics, 320
(3):0 468–482, 2014.

[Xu et al., 2020]
Liyuan Xu, Yutian Chen, Siddarth Srinivasan, Nando de Freitas, Arnaud Doucet,
and Arthur Gretton.
Learning deep features in instrumental variable regression.
arXiv preprint arXiv:2010.07154, 2020.

[Yang et al., 2017]
Yun Yang, Anirban Bhattacharya, and Debdeep Pati.
Frequentist coverage and sup-norm convergence rate in Gaussian
process regression.
arXiv preprint arXiv:1708.04753, 2017.

[Zellner, 1995]
Arnold Zellner.
Bayesian method of moments (BMOM) analysis of mean and regression
arXiv preprint arXiv:bayes-an/9511001, 1995.

[Zhang et al., 2020]
Rui Zhang, Masaaki Imaizumi, Bernhard Schölkopf, and Krikamol Muandet.
Maximum moment restriction for instrumental variable regression.
arXiv preprint arXiv:2010.07684, 2020.

[Zhang, 2006]
Tong Zhang.
From $\epsilon$-entropy to KL-entropy: Analysis of minimum
information complexity density estimation.
The Annals of Statistics, 340 (5):0
2180–2210, 2006.

[Zhu and Xing, 2009]
Jun Zhu and Eric P Xing.
Maximum entropy discrimination Markov networks.
Journal of Machine Learning Research, 100 (11), 2009.