text:
\begin{center}
    {\bf {\large Supplementary Material:\\
    Minimax Estimation of Conditional Moment Models}}
\end{center}

\etocdepthtag.toc{mtappendix}
\etocsettagdepth{mtchapter}{none}
\etocsettagdepth{mtappendix}{subsection}
\tableofcontents

\section{Further Discussion on Related Work}
\label{app:related}

The non-parametric IV problem has a long history in econometrics \citep{newey2003instrumental,blundell2007semi,chen2012estimation,chen2018optimal,hall2005nonparametric,horowitz2007asymptotic,horowitz2011applied,darolles2011nonparametric,chen2009efficient}. Arguably the closest to our work is that of \cite{chen2012estimation} (in particular their Theorem~4.1), who consider estimation of non-parametric function classes and estimation via the method of sieves and a penalized minimum distance estimator of the form: $\min_{h \in \mcH} \E[\E[y-h(x)\mid z]^2] + \lambda R(h)$, where $R(h)$ is a regularizer. The authors approximate the function class $\mcH$ by linear functions in a growing feature space. Subsequently, they also estimate the function $m(z) = \E[y-h(x)\mid z]$ based on another growing sieve.

Though it may seem at first that the approach in that paper and ours are quite distinct, the population limit of our objective function coincides with theirs.
To see this, consider the simplified version of our estimator presented in \Cref{eqn:minimax-no-norms}, where the function classes are already norm-constrained and no norm based regularization is imposed. Moreover, for a moment consider the population version of this estimator, i.e.
\begin{equation}
    \min_{h\in \mcH} \max_{f\in \mcF} \Psi(h, f) - \|f\|_2^2 = \min_{h\in \mcH} \max_{f\in \mcF} \E[(y - h(x)) f(z) - f(z)^2]
\end{equation}
Observe that if $\mcF$ is expressive enough (if $T(h_0-h)\in \mcF$), then the maximizing test function is $\frac{1}{2}\E[y-h(x)\mid z]=\frac{1}{2}\E[h_0(x)-h(x)\mid z]$.
Then by the law of iterated expectations, the population criterion becomes:
\begin{equation}
   \min_{h\in \mcH} \E\left[(y - h(x)) \frac{1}{2}\E[y-h(x)\mid z] - \frac{1}{4}\E[y-h(x)\mid z]^2\right] = \min_{h\in \mcH} \frac{1}{4} \E\left[\E[y-h(x)\mid z]^2\right]
\end{equation}
Thus in the population limit and without norm regularization on the test function $f$, our criterion is equivalent to the minimum distance criterion analyzed in \cite{chen2012estimation}.
Another point of similarity is that we prove convergence of the estimator in terms of the pseudo-metric, the projected MSE defined in Section 4 of \cite{chen2012estimation} - and like that paper we require additional conditions to relate the pseudo-metric to the true MSE.

The present paper differs in a number of ways: (i) the finite sample criterion is different; (ii) we prove our results using localized Rademacher analysis which allows for weaker assumptions; (iii) we consider for a broader range of estimation approaches than linear sieves, necessitating more of a focus on optimization.

Digging into the second point, \cite{chen2012estimation} take a more traditional parameter recovery approach which requires several minimum eigenvalue conditions and several regularity conditions to be satisfied for their estimation rate to hold (see e.g. their Assumptions 3.1, 3.2, 3.3, 4.1 and C.1). This is like a mean squared error proof in an exogenous linear regression setting, that requires a minimum eigenvalue of the feature co-variance to be bounded.
Moreover, such parameter recovery methods seem limited to the growing sieve approach, since only then one has a clear finite dimensional parameter vector to work on for each fixed $n$.

In contrast we work with infinite dimensional parameter spaces directly and our analysis makes no further assumptions other than boundedness of the random variables and the conditional moment restriction in order to provide a projected MSE rate.
We do not require that the hypothesis space be a convex set, nor that the moment is path-wise differentiable with respect to $h$.
Relaxing these assumptions is important, since they are violated in three of our leading examples: linear hypothesis spaces with hard sparsity constraints or for neural network spaces or for tree based regressors.
Another benefit of the localized Rademacher analysis is that we do not require a preliminary proof of consistency, which is typical of more classical approaches to MSE rates.
Such proofs typically require that $n$ be larger than some constant before the convergence rate kicks in, so that the estimator is within some small ball around the truth. This constant can sometimes be prohibitively large. Our convergence rate is global and holds without any lower bound condition on $n$.
The sieve method is most closely related to our RKHS section (and the expository sieve \Cref{sec:sieves}), where essentially we consider infinite dimensional linear function spaces.
However, unlike the sieve method, we do not clip the eigenfunctions to a finite set that is growing, but rather impose an RKHS penalty.
We show that this approach has advantages in auto-tuning to the ill-posedness of the problem.
Finally, we do not require a bound on the ill-posedness of the problem in order to prove convergence rates in terms of the pseudo-metric - this bound is only needed in post-processing to relate the pseudo-metric to the MSE.
By contrast \cite{chen2012estimation} use the bounded ill-posedness condition (Assumption~4.1) to prove convergence in the pseudo-metric.

As a concrete example of the differences in the analysis, we apply our main
\Cref{thm:reg-main-error} for the case where $\mcH$ and $\mcF$ are growing sieves, equipped with the parameter $\ell_2$ norms, i.e. $\mcH=\left\{\ldot{\theta}{\phi_n(\cdot)}: \theta\in \R^{k_n}\right\}$, $\mcF=\left\{\ldot{\beta}{\psi_n(\cdot)}: \beta\in \R^{m_n}\right\}$, $\|\ldot{\theta}{\phi_n(\cdot)}\|_{\mcH}=\|\theta\|_2$, $\|\ldot{\beta}{\psi_n(\cdot)}\|_{\mcF}=\|\beta\|_2$, for some fixed and growing feature maps $\phi_n(\cdot)$, $\psi_n(\cdot)$. In that case $\eta_n$ will correspond to the approximation error of the sieve $\psi_n$ that is used for the test function space and, if we choose $h_*=\argmin_{h\in \mcH} \|h_*-h_0\|_2$, then $\|T(h_* -h_0)\|_2\leq \|h_*-h_0\|_2 =: \epsilon_n$, will correspond to the approximation error of the sieve $\phi_n$ that is used for approximating the model $h_0$. In that case, \Cref{thm:reg-main-error} gives a bound of $O\left(\delta_n \|\theta_*\|_2 + \eta_n + \epsilon_n\right)$, where $\theta_*$ is the $\ell_2$ norm of the parameter of the projection of $h_0$ on the sieve space for the model, i.e $\argmin_{\theta\in \R^{k_n}} \|\ldot{\theta}{\phi(\cdot)} - h_0\|_2^2$. Moreover, $\delta$ is a bound on the critical radius of $\mcF_U$ and $\mcG_{B, U}$. Since both are finite dimensional linear functions, via standard covering arguments (see \Cref{lem:covering}), we can bound $\delta = O\left(\sqrt{\frac{\max\{k_n, m_n\}\, \log(n)}{n}}\right)$.\footnote{The $\log(n)$ factor can also be saved with a more careful analysis of the critical radius for finite dimensional linear function spaces (see Section~\ref{sec:sieves}).} Combined with ill-posedness conditions provided in \citep{chen2012estimation}, our results can thus give an alternative proof to the results in \citep{chen2012estimation} that i) do not make minimum eigenvalue conditions, ii) provide adaptivity to $\|\theta_*\|_2$, without knowledge of it, thereby justifying theoretically the use of the regularization term $R(h)$, that was mostly proposed for experimental improvement in \citep{chen2012estimation}. We provide a more thorough exposition of how our main theorem applies to the case of growing sieves in \Cref{sec:sieves}.

The localized Rademacher analysis also allows us to consider hypothesis spaces that are not linear sieves, such as neural nets and random forests.
This introduces some new optimization difficulties, as the estimator cannot be written in closed form (as it can for linear sieves).
Our work gives several solutions for these difficulties, via iterative first order algorithms.
Intuitively, our optimization algorithms gradually and iteratively make gradient steps towards solving both optimization problems (of regressing $y-h(x)$ on $z$ and minimizing $\E[\E[y-h(x)\mid z]^2]$ over $\mcH$), as opposed to calculating full solutions of either problem. This formulation allows us to work with arbitrary hypothesis spaces and not just linear sieves.

There is also a growing body of work in the machine learning literature on the non-parametric instrumental variable regression problem \citep{deepiv,bennett2019deep,singh2019kernel,muandet2019dual,muandet2020kernel}. The seminal work of \cite{deepiv} provided a methodology for training neural networks that solve the instrumental variable problem by taking a non-parametric analogue of the two stage least squares method.

\cite{bennett2019deep} also consider a minimax criterion with a variance penalty. Albeit the variance penalty they impose is not the second moment of the test functions and depends on a preliminary estimate of the true model. Moreover, they only show asymptotic consistency of their estimate and not finite sample rates and primarily focus on neural network applications (see Section~\ref{sec:neural-networks} for more details). \cite{singh2019kernel} consider a RKHS analogue of \cite{deepiv}, where the hypothesis space $h$ fall in an RKHS and the conditional distribution of $X$ conditional on $Z$ is represented via a conditional kernel mean embedding. They offer very strong finite RKHS-norm rates on the estimated $h$, which typically imply sup-norm rates of the recovered function. Albeit, we focus on projected MSE and MSE rates and achieve faster rates as a function of the eigendecay of the kernel and the degree of ill-posedness. Moreover, the work of \cite{singh2019kernel} makes several stronger \emph{prior} assumptions, that control the smoothness of the function within the kernel, assumptions that are typical of RKHS norm guarantees in kernel ridge regression \citep{caponnetto2007optimal}, but which are not required for the weaker MSE metric. \cite{muandet2019dual} also propose a method that is very related to the second moment penalized method that we propose, albeit the motivation stems from a different dual formulation of the two-stage-least-squares problem presented in \citep{deepiv} and similar to \citep{bennett2019deep} only offer asymptotic consistency of the estimator and only focus on RKHS function spaces. Finally, \cite{muandet2020kernel} consider the version of the minimax criterion that does not impose the second moment penalty on $f$, and make the important observation that for RKHS function spaces, the internal maximization takes a closed form, leading to a pairwise sample criterion (see Equation~\eqref{eqn:kloss-rkhs} and Equation~\eqref{eqn:kloss-mmd}). Moreover, they focus primarily on hypothesis testing as opposed to estimation. The un-penalized criterion can have sub-optimal convergence guarantees, as it does not posses the property that as the hypothesis of the learner gets close to the truth, then the adversary is testing smaller functions in terms of variance. The inability to achieve the fast rates attained via the critical radius was the main reason why we introduced the second moment penalty. The suboptimality of the un-penalized kernel based criterion was also proven in the context of hypothesis testing by \cite{balasubramanian2017optimality}, who also show that a form of second moment penalization can yield hypothesis tests with optimal power, when the alternative is very close to the null. Moreover, for RKHS, we show that the penalized method still admits a closed form solution, albeit now the closed form depends on the inverse of a kernel matrix, which makes it less amenable to gradient training as we discuss in \ref{sec:neural-networks}.

\section{Beyond the IV Moments}\label{app:beyond-iv}

Our results easily extend to arbitrary moments that are linear in $h$, which can capture several other problems in econometrics and causal inference, but for simplicity of exposition we focus on the case of moments of the form $y-h(x)$. Moreover, our results can also be extended to non-linear and non-smooth moments $\psi(y;h(x))$, albeit in that case our convergence rates will be with respect to the distance metric: $d(\hat{h}, h)=\sqrt{\E[\E[\psi(y; \hat{h}(x)) - \psi(y; h(x))\mid z]^2]}$ as opposed to the projected MSE distance. For instance, in the case of $\alpha$-quantile IV regression: $\psi(y; h(x)) = a - 1\{y\leq h(x)\}$ and the distance metric corresponds to: $d(\hat{h}, h_0)=\sqrt{\E[\E[1\{y\leq h(x)\} - \alpha\mid z]^2]}$.

\section{Supplementary Discussion of Main Theorems}\label{app:further-main}

\subsection{Adaptivity of Regularized Estimator}\label{app:adaptivity}
Suppose that we know that for $B, U=1$, we have that functions in $\mcH_B, \mcF_U$ have ranges in $[-1, 1]$ as their inputs range in $\mcX$ and $\mcZ$ correspondingly. Then our Theorem requires that we set: $\lambda \geq \delta^2$ and $\mu\geq 2\lambda (4L^2 + 27)$, where $\delta^2$ depends on the critical radius of the function class $\mcF_1$ and $\mcG_1$. Observe that none of these values depend on the norm of the benchmark hypothesis $\|h_*\|_{\mcH}$, which can be arbitrary and not constrained by our theorem. For instance, if we knew that the true model $h_0\in \mcH$ and $T(h-h_0)\in \mcF_{L^2\|h-h_0\|_{\mcH}^2}$, then we can apply the latter theorem to get rates of the form:
\begin{equation}
    O\left(\delta \max\left\{1, \|h_0\|_{\mcH}^2\right\}\right)
\end{equation}
with $\lambda = \delta^2$ and $\mu = 2\delta^2 (4L^2 + 27)$. This hyperparameter tuning only requires knowledge of the critical radius of the function classes $\mcF_1$ adn $\mcH_1$ and the Lipschitz constant of the operator $T$, but does not require knowledge of the norm of the true model $\|h_0\|_{\mcH}$, nor upper bounds on it. If the true model does not fall in the hypothesis $\mcH$, then observe that we also require knowledge of the unconstrained approximation error, i.e. if we knew that:
\begin{equation}
    \inf_{h\in \mcH} \|h-h_0\|_2 \leq \epsilon_n
\end{equation}
and that $T(h-h_0)\in \mcF_{L\|h-h_0\|_{\mcH}}$, then we can choose $\delta\geq \epsilon_n$ to get rates of the form:
\begin{equation}
    O\left(\delta \max\left\{1, \|h_*\|_{\mcH}^2\right\} + \epsilon_n\right)
\end{equation}
where $h_*=\arginf_{h\in \mcH} \|h-h_0\|_2$. Again we do not require knowledge of the norm of the unconstrained projection, $\|h_*\|_{\mcH}$, just bounds on the approximation error of the unconstrained function space. Then the regularized estimator adapts to the norm of the projection of the true model on $\mcH$. These results are inline with recent work on statistical learning theory \citep{lecue2017regularization,Lecue2018} for square losses and extend these qualitative insights to the minimax objectives that we deal with.

\subsection{Critical Radius and Rademacher Complexity via Covering}

The critical radius of a function class is characterized to within a constant factor by it's empirical localized Rademacher critical radius, which subsequently is chracterized by the empirical entropy integral. The empirical Rademacher complexity of a function class $\mcG:\mcV\to [-1,1]$, for a given set of samples $S=\{v_i\}_{i=1}^n$ is defined as:
\begin{equation}
    R_S(\delta; \mcG) = \E_{\{\epsilon_i\}_{i=1}^n}\left[\sup_{g\in \mcG: \|g\|_{2,n}\leq \delta} \frac{1}{n}\sum_i \epsilon_i g(v_i)\right]
\end{equation}
The empirical critical radius is defined as any solution $\hat{\delta}_n$ to:
\begin{equation}
    R_S(\delta; \mcG) \leq \delta^2
\end{equation}
Proposition 14.1 of \cite{wainwright2019high} shows that w.p. $1-\zeta$, \begin{equation}\label{eqn:emp-v-pop-critical}
\delta_n=O\left(\hat{\delta}_n + \sqrt{\frac{\log(1/\zeta)}{n}}\right).
\end{equation}
Thus we can choose $\delta$ in our main theorems based on the empirical critical radius $\hat{\delta}_n$.

Moreover, an upper bound on the empirical critical radius can be obtained via the empirical covering integral defined as follows. An empirical $\epsilon$-cover of $\mcG$, is any function class $\mcG_{\epsilon}$, such that for all $g\in \mcG$, $\inf_{g_\epsilon\in \mcG_{\epsilon}} \|g_\epsilon-g\|_{2,n}\leq \epsilon$. We denote with $N(\epsilon, \mcG, S)$ as the size of the smallest empirical $\epsilon$-cover of $\mcG$. The empirical metric entropy of $\mcG$ is defined as $H(\epsilon, \mcG, S)=\log(N(\epsilon, \mcG, S))$. An empirical $\delta$-slice of $\mcG$ is defined as $\mcG_{S,\delta}=\{g\in \mcG: \|g\|_{2,n}\leq \delta\}$. Then the empirical critical radius of $\mcG$ is upper bounded by any solution to the inequality:
\begin{equation}\label{eqn:empirical-covint}
    \int_{\delta^2/8}^{\delta} \sqrt{\frac{H(\epsilon, \mcG_{S, \delta}, S)}{n}} d\epsilon \leq \frac{\delta^2}{20}
\end{equation}
Observe that a conservative upper bound on $\hat{\delta}_n$ comes from replacing $\mcG_{S,\delta}$ inside the integral with $\mcG$, i.e. when we do not restrict the function class to be in an empirical $\delta$-slice, when calculating it's empirical metric entropy. For many function classes (e.g. parametric $\ell_2$-balls, RKHS, high-dimensional sparse parametric spaces, VC-subgraph classes) this still yields tight results. For some other cases, such as $\ell_1$-balls centered around a sparse parameter, this can be loose.

When we make this relaxation, then observe that we can derive an upper bound on the critical radius of $\mcG_{B,U}$, as a function of the empirical metric entropy of $\mcH$ and $\mcF$. Observe that if $\mcH_{\epsilon}$ is an empirical $\epsilon$-cover of $\mcH$ and $\mcF_{\epsilon}$ is an empirical $\epsilon$-cover of $\mcF_U$, then since $\mcH$ contains functions uniformly bounded in $[-1,1]$, we have that:
\begin{equation}
    \inf_{h\in \mcH_{\epsilon}, f_{\epsilon}\in \mcF_{\epsilon}} \|(h_{\epsilon}-h) f_{\epsilon} - (h-h_*) f_h^U\|_{2,n} \leq 2 \|h_\epsilon - h_*\|_{2,n} + 2\|f_{\epsilon} - f_h^U\|_{2,n}\leq 4\epsilon
\end{equation}
Thus, the product of these two spaces is an $\epsilon$-cover of the function class $\mcG$ defined in Equation~\cref{eqn:reg-class-G}. Hence, the empirical metric entropy of $\mcG$ satisfies:
\begin{equation}
    H(\epsilon, \mcG_{B,U}, S) \leq H(\epsilon/4, \mcH_B, S) + H(\epsilon/4, \mcF_U, S)
\end{equation}

Thus by applying Proposition 14.1 of \cite{wainwright2019high} we get the following corollary.
\begin{corollary}\label{lem:covering}
Suppose that $\hat{\delta}_n$ satisfies the inequality:
\begin{equation}
    \int_{\delta^2/8}^{\delta} \sqrt{\frac{H(\epsilon/4, \mcH_{2B}, S) + H(\epsilon/4, \mcF_{3U}, S)}{n}} d\epsilon \leq \frac{\delta^2}{20}
\end{equation}
Then w.p. $1-\zeta$, $\delta_n\leq  O\left(\hat{\delta}_n + \sqrt{\frac{\log(1/\zeta)}{n}}\right)$, where $\delta_n$ is the maximum of the critical radii of $\mcF_{3U}$, $\mcG_{B, U}$ and $\hat{\mcG}_{B,U}$.
\end{corollary}
For instance, if $\mcH$ and $\mcF$ is assumed to be a VC-subgraph class with constant VC dimension, then the above is satisfied for $\hat{\delta}_n = O\left(\sqrt{\frac{\log(n)}{n}}\right)$.

\subsection{Solving the Min-Max Optimization Problem}

In this section we outline some strategies for addressing the empirical min-max problem required by the estimators described in Equations~\eqref{eqn:reg-estimator} and \eqref{eqn:reg-estimator-2}. In subsequent sections, we will present instances of these optimization approaches for each of the function classes that we consider.

First observe that if the hypothesis space can be parameterized as $h(x; \theta)$, such that the moment $\psi(y; h(x;\theta))$ is convex in $\theta$ and the inner optimization problem is solvable in closed form
then we can solve the empirical problem via subgradient descent: i.e. letting
\begin{align}
    f_*(\cdot;h) :=~& \argsup_{f\in \mcF} \Psi_n(h,f) - \lambda \Phi(f), \\
    \theta_{t+1} :=~& \theta_t - \eta\, \left(\E_n\left[ f_*(z; h(x;\theta_t)) \nabla_\theta h(x;\theta_t)\right] + \mu\nabla_{\theta} R(h(\cdot; \theta_t))\right)
\end{align}
where $\Phi, R$ are the regularizers on $f$ and $h$ correspondingly. After $T$ iterations, the average parameter $\bar{\theta}=\frac{1}{T}\sum_{t=1}^T \theta_t$, will correspond to an $O\left(T^{-1/2}\right)$ approximate solution to the min-max problem. This approximate solution will satisfy the same guarantees as $\hat{h}$ presented in \Cref{thm:reg-main-error} and \Cref{thm:reg-main-error-2}, augmented by an extra $O\left(T^{-1/2}\right)$ additive factor.

Many times, even if the hypothesis space is not parameterizable by a finite dimensional parameter vector $\theta$, universally, we can invoke characterizations (typically referred to as \emph{representer theorems}), that prove that the empirical solution can always be expressed in terms of a finite set of parameters (many times of the order of the number of samples). This is for instance the case when $\mcF$ and $\mcH$ belong to a Reproducing Kernel Hilbert space, as we will see in Section~\ref{sec:rkhs}. In such settings, we will see that even the overall min-max optimization problem can be expressed in closed form, involving only matrix inversions and mutliplications, with matrices of size of the order of $n^2$.

Since the min-max problem does not have a smooth gradient, one can also benefit by invoking algorithms that are tailored to saddle point problems. These improvements typically assume some structure on the inner optimization problem. For instance, if the function $f$ can be parameterizedd as $f(\cdot; w)$ such that the inner maximization problem is concave in $w$ then faster than $T^{-1/2}$ optimization rates can be achieved. We will see examples of such settings in the high-dimensional linear function class setting in Section~\ref{sec:high-dim-linear}. The following set of papers provide examples of algorithms that achieve $T^{-1}$ approximation rates (see e.g. \cite{nesterov2005smooth,Nemirovski2004,Rakhlin2013,mokhtari2019unified}).

One simple such algorithm is the simultaneous optimistic mirror descent algorithm proposed in \cite{Rakhlin2013} and also recently analyzed by several papers, both theoretically and empirically, in the context of non-convex optimization problems (see e.g. \cite{Daskalakis2017,Mertikopoulos2018}). In this algorithm, instead of fully solving the internal optimization problem, we only take gradient steps. However, it modifies the gradient descent algorithm to incorporate a notion of \emph{optimism} (i.e. that the next gradient will look similar to the last gradient). In particular, if we use the short-hand notation $\Psi_n(\theta, w):=\Psi_n(h(\cdot;\theta), f(\cdot;w))$, then in the simplified setting where we have no regularization on $\theta, w$, the algorithm is described via the following update dynamics:
\begin{align}
    \theta_{t+1} =~& \theta_t - 2 \eta \nabla_\theta \Psi_n(\theta_t, w_t) + \eta \nabla_\theta \Psi_n(\theta_{t-1}, w_{t-1})\\
    w_{t+1} =~& w_t + 2 \eta \nabla_{w} \Psi_n(\theta_t, w_t) - \eta \nabla_w \Psi_n(\theta_{t-1}, w_{t-1})
\end{align}
Convex constraints on $\theta$ and $w$ can be easily incorporated via projection steps and we defer to \cite{Rakhlin2013} for the formal definition of the algorithm in that setting. Similarly, for the regularized versions one would simply replace $\Psi_n$ with its regularized counterparts.

Unlike the sub-gradient descent approach, the simultaneous optimistic gradient dynamics, with the regularized version of our estimator, can also be implemented in a stochastic gradient manner, where a mini-batch of samples are drawn at each step (with replacement), from the empirical set of samples and $\Psi_n$ is replaced with the empirical expectation over that sub-sample. This can enable applications where storing all the dataset in-memory is prohibitive. Moreover, this algorithm has variants that have been proven beneficial for neural nets (see, e.g. the Optimistic Adam algorithm of \cite{Daskalakis2017}, also used in the related work of \cite{bennett2019deep} in a generalized method of moments setup). Properties of simultaneous gradient dynamics in non-convex/non-concave settings have also been a topic of recent interest in the machine learning community and recent techinques from this line of work can be invoked to empirically solve the optimization problem (see e.g. \cite{Jin2019,Nouiehed2019,Thekumparampil2019,yang2020global,lin2020near}).

\subsection{From Projected MSE to MSE: Measure of Ill-Posedness}

If we want to get a bound on the RMSE of $\hat{h}$, i.e. $\|h-h_0\|_2$, then we need to bound the quantity:
\begin{equation}
    \tau^*(\delta) = \sup_{h\in \mcH_B: \|T(h-h_*)\|_2 \leq \delta} \|h-h_*\|_2
\end{equation}
In fact, it suffices to bound the measure of ill-posedness of the operator $T$ with respect to the function class $\mcH_B$, defined as:
\begin{equation}
\label{eqn:illposedness}
    \tau := \sup_{h\in \mcH_B} \frac{\|h-h_*\|_2}{\|T(h-h_*)\|}.
\end{equation}
Both of these measures have been used in the literature on conditional moment models. For instance, \cite{chen2012estimation} defines both of these measures for the case where $\mcH_B$ is a space of growing linear sieves. In that case, the second measure $\tau$ is typically referred to as the \emph{sieve measure of ill-posedness}. Then observe that \Cref{thm:reg-main-error} implies that:
\begin{equation}
    \|\hat{h}-h_*\|_2 \leq \tau \|T(\hat{h}-h_*)\|_2 \leq O\left( \tau\, \delta_n + \tau\|T(h_*-h_0)\|_2\right) \leq O\left( \tau\, \delta_n + \tau\,\|h_*-h_0\|_2\right)
\end{equation}
which by a triangle inequality also implies that:
\begin{equation}
    \|\hat{h}-h_0\|_2 \leq O\left( \tau\, \delta_n + (\tau + 1) \|h_*-h_0\|_2\right)
\end{equation}
Choosing $h_*=\argmin_{h\in \mcH: \|h\|_\mcH\leq B} \|h_*-h_0\|$, yields the bound:
\begin{equation}
    \|\hat{h}-h_0\|_2 \leq O\left( \tau\, \delta_n + (\tau + 1) \inf_{h\in \mcH: \|h\|_{\mcH}} \|h-h_0\|_2\right)
\end{equation}
Subsequently one can appropriately choose $\mcH$ and $B$ so as to trade-off the ill-posedness constant and the bias term.

Moreover, we show that when we have a bounded ill-posedness measure, then we can prove a more convenient version of \Cref{thm:reg-main-error}, that only requires bounds on the critical radius of the centered function classes $\shull(\mcH_B-h_*)=\{r(h-h_*): h\in \mcH_B, r\in [0, 1]\}$ and $\shull(T(\mcH-h_*))=\{T(h-h_*): h\in \mcH_B, r\in [0, 1]\}$, as opposed to the space $\mcG$ that contains products of these functions.
\begin{theorem}\label{thm:main-error-ill}
Let $\mcF$ be a symmetric and star-convex set of test functions and
consider the estimator in Equation~\eqref{eqn:reg-estimator}.
Let $h_0$ be any hypothesis (not necessarily in $\mcH$) that satisfies the Conditional Moment~\eqref{eqn:cond-moments} and suppose that $\mcH$ satisfies that:
\begin{equation}
    \inf_{h\in \mcH} \|h-h_0\|_2 \leq \epsilon_n
\end{equation}
and let $h_*=\arginf_{h\in \mcH} \|h-h_0\|_2$. Moreover, suppose that:
\begin{equation}
\forall h\in \mcH: \min_{f\in \mcF_{L^2\|h-h_*\|_{\mcH}^2}} \|f - T(h-h_*)\|_2 \leq \eta_n
\end{equation}
Assume that functions in $\mcH_B$ and $\mcF_{3U}$ have uniformly bounded ranges in $[-1, 1]$ and that:
\begin{equation}
\delta:=\delta_n + \eta_n + \epsilon_n + c_0 \sqrt{\frac{\log(c_1/\zeta)}{n}}
\end{equation}
for universal constants $c_0, c_1$, and $\delta_n$ an upper bound on the critical radii of the classes $\mcF_{3U}$ and
\begin{align}
\shull(\mcH_B-h_*) :=~& \{r\, (h-h_*): h - h_*\in \mcH_B, r\in [0, 1]\}\\
\shull(T(\mcH_B-h_*)) :=~& \{r\, f_h: h - h_*\in \mcH_B, r\in [0, 1]\}
\end{align}
where $f_h=\argmin_{f\in \mcF_U} \|f-T(h-h_*)\|_2$.
If $O(\delta^2)\geq \lambda \geq \delta^2/U$ and $O(\delta^2)\geq \mu\geq 2\lambda (4 L^2 + 27U/B)$, then $\hat{h}$ satisfies w.p. $1-3\,\zeta$:
\begin{equation}
     \|h-h_0\|_2 = O\left(\tau^2 \delta \max\{1, \|h_*\|_{\mcH}\}\right)
\end{equation}
\end{theorem}

\subsection{Minimax Optimality of Estimation Rate}

In this section we take the viewpoint of establishing minimax optimal rates for the estimation problem of interest and discuss under which circumstances the upper bound we provide will typically be tight (i.e. achieving the statistically best possible projected RMSE). Suppose that the only prior assumptions we are willing to make about our data generating process is that it satisfies the moment condition, that $h_0\in \mcH$ and that $T_0 \in \mcT$ for some function class $\mcH$ and linear operator class $\mcT$. Moreover, let $\mcF:=\{Th: T\in \mcT, h\in \mcH\}$. What is the minimax estimation rate, with respect to the projected MSE norm, achievable in this setting? More concretely, let $D(h, T)$ be any distribution consistent with function $h$, linear operator $T$ and conditional moment condition $Th=E[y\mid z]$. Then for any estimator $\hat{h}$, that takes as input a training sample $S$ of size $n$, drawn i.i.d. from $D(h, T)$, and returns a function $\hat{h}_S$, we want to lower bound the minimax optimal rate:
\begin{equation}
    \min_{\hat{h}} \max_{h_0\in \mcH, T_0 \in \mcT}  \E_{S\sim D(h_0,T_0)^n}\left[ \|T_0(\hat{h}_S-h_0)\|_2^2 \right]
\end{equation}

If the space $\mcT$ contains the identity, then this is lower bounded by the RMSE rates of a non-parametric regression problem over hypothesis space $\mcH$. Thus by standard results on regression problems, the critical radius of $\mcH$ is insurmountable for many classes $\mcH$ of interest (see e.g. \cite{massart2000,bartlett2005local,rakhlin2017minimax}.

Moreover, suppose that there exists a $T\in \mcT$ such that: for all $f$ there exists $h\in \mcH$, such that $Th=f$, i.e. $T$ is the worst mapping that allows one to span all of $\mcF$. Then even if we knew $T=T_0$, we could not bypass the critical radius of $\mcF$ for many classes $\mcF$ of interest (see e.g. \cite{bartlett2005local,rakhlin2017minimax}). More generally, we can lower bound the minimax risk as:
\begin{equation}
    \max_{T_0\in \mcT} \min_{\hat{h}} \max_{h_0\in \mcH}  \E_{S\sim D(h_0,T_0)^n}\left[ \|T_0(\hat{h}_S-h_0)\|_2^2 \right]
\end{equation}
Let $\mcF_T=\{Th: h\in \mcH\}$. Then the above can be re-written:
\begin{equation}
    \max_{T_0\in \mcT} \min_{\hat{f} \in \mcF_T} \max_{f_0\in \mcF}  \E_{S\sim D(f_0)^n}\left[ \|\hat{f}_S-f_0\|_2^2 \right]
\end{equation}
where $D(f_0)$ is any distribution that satisfies $\E[y\mid z]=f_0$. This is the minimax lower bound for the regression problem of predicting $y$ from $z$, assuming that $\E[y\mid z]\in \mcF_T$. Thus we have that the minimax rate is at least $\max_{T} \delta(\mcF_T)$. If we knew that there was a finite set of $k$ representative linear operators $T_1, \ldots, T_k$ in $\mcT$, such that $\mcF = \mcF_{T_1} \cup \ldots \cup \mcF_{T_k}$, then observe that the critical radius of $\mcF$ is at most $O(\log(k))$ more than the maximum critical radius of each of the $\mcF_{T_i}$. Thus the only case that remains open where our upper bound might not be providing tight results is when there is not such finite small set of representative operators in $\mcT$. In many of our settings, we will have that $\delta(\mcF) \sim \delta(\mcH)$, which is achieved for the single identity operator $\mcT=I$. The case where our upper bound is loose, is essentially the case when knowing the operator, or some equivalence class of the operator, can significantly reduce the sample complexity of the problem. Potentially in such settings fitting a first stage model of $T$ to identify the equivalence class or a finite number of viable equivalence classes and focus only on a remaining set of $k$ candidate $\mcF_{T_1}\cup\ldots\cup \mcF_{T_i}$ in a second stage can be beneficial. However, in most of our applications this setting does not arise. One for instance can follow techniques similar to aggregation algorithms \cite{rakhlin2017minimax}, that applies our minimax estimator on an $\epsilon$ partition of the original hypothesis $\mcH$ and then aggregates the resulting winning hypothesis from each partition. However, this would typically be a computationally inefficient algorithm.

\section{Application: Growing Linear Sieves}\label{sec:sieves}

Consider the case where $\mcH$ and $\mcF$ are growing linear sieves, i.e.
\begin{align}
\mcH=\mcH_n:=\left\{\ldot{\theta}{\phi_n(\cdot)}: \theta\in \R^{k_n}\right\},\\
\mcF=\mcF_n:=\left\{\ldot{\beta}{\psi_n(\cdot)}: \beta\in \R^{m_n}\right\},
\end{align}
equipped with norms $\|\ldot{\theta}{\phi_n(\cdot)}\|_{\mcH}=\|\theta\|_2$, $\|\ldot{\beta}{\psi_n(\cdot)}\|_{\mcF}=\|\beta\|_2$, for some known and growing feature maps $\phi_n(\cdot)$, $\psi_n(\cdot)$.

Moreover, we denote with $\eta_n$ the approximation error of the sieve $\psi_n$ that is used for the test function space, i.e. for all $h, h_*\in \mcH$:
\begin{equation}
    \inf_{f\in \mcF} \|f - T(h-h_*)\|_{2} \leq \eta_n
\end{equation}
and, let $\epsilon_n$ the approximation error of the sieve $\phi_n$ used for the model, i.e.:
\begin{equation}
\inf_{h \in \mcH} \|h - h_0\|_2 \leq \epsilon_n
\end{equation}
In that case, applying \Cref{thm:reg-main-error} with $h_* = \arginf_{h\in \mcH} \|h - h_0\|_2$, gives a bound w.p. $1-\zeta$ of:
\begin{equation}
 \|T(\hat{h}-h_0)\|_2 \leq O\left(\left(\delta_n +\epsilon_n + \sqrt{\frac{\log(1/\zeta)}{n}}\right)\max\{1, \|\theta_*\|_2^2\} + \eta_n \right)
\end{equation}
where $\theta_*$ is the $\ell_2$ norm of the parameter that corresponds to $h_*$.

Moreover, $\delta_n$ is a bound on the critical radius of $\mcF_U$ and $\mcG_{B, U}$. Since both are finite dimensional linear functions, via standard covering arguments (see \Cref{lem:covering}), we can bound $\delta_n = O\left(\sqrt{\frac{\max\{k_n, m_n\} \log(n)}{n}}\right)$. We also now provide a more intricate argument that removes the $\log(n)$ from this rate. Observe that $\mcF_U$ is a simple linear model space and therefore existing results directly apply to show that the critical radius of $\mcF_U$ is at most $\sqrt{\frac{m_n}{n}}$ (see e.g. Example 13.5 of \cite{wainwright2019high}). The function space $\mcG_{B,U}$ is a bit more subtle. We will in fact bound the critical radius of the following larger class:
\begin{equation}
    \tilde{\mcG}_{B, U} = \{ (x, z) \to \ldot{\theta - \theta_*}{\phi_n(x)} \ldot{\beta}{\psi_n(z)}: \theta \in \R^{k_n}, \beta \in \R^{m_n}, \|\theta-\theta_*\|_2\leq B, \|\beta\|_{2}\leq U\}
\end{equation}
We will use the empirical covering integral bound on the critical radius, presented in Equation~\eqref{eqn:empirical-covint}. Thus we need to bound the metric entropy of the function class $\tilde{\mcG}_{B,U}(\delta)=\{g\in \tilde{\mcG}_{B,U}: \|g\|_{2,n}\leq \delta\}$. Let $\Psi_n$ denote the $n\times k_n$ matrix whose $i$-th row corresponds to the vector $\psi_n(x_i)$ and similarly $\Phi_n$. Observe that the norm empirical $\ell_{2,n}$ norm can then be written as:
\begin{equation}
    \|\ldot{\theta-\theta_*}{\phi_n(\cdot)}\ldot{\beta}{\psi_n(\cdot)}\|_{2,n} = \frac{\|\Psi_n (\theta-\theta_*)\|_2 \|\Phi_n \beta\|_2}{\sqrt{n}}
\end{equation}
Thus $\ell_{2,n}$ defines a norm on the space defined by the Hadamard (coordinate-wise) product $v_1 \circ v_2$ of two vectors $v_1, v_2$ in $\text{range}(\Psi_n)$ and $\text{range}(\Phi_n)$, correspondingly, i.e. $\|v_1\circ v_2\| = \frac{\|v_1\|_2\, \|v_2\|_2}{\sqrt{n}}$. Moreover, $\tilde{\mcG}_{B,U}(\delta)$ is isomorphic to a $\delta$-ball in this space. Moreover, observe that the dimension of the space $\{v_1\circ v_2: v_1\in \text{range}(\Psi_n), v_2\in \text{range}(\Phi_n)\}$ is at most $\text{rank}(\Psi_n)\, \text{rank}(\Phi_n) \leq k_n\cdot m_n$. Therefore by the volumetric argument presented in Example~5.4 of \cite{wainwright2019high}, we get that for any set of samples $S$ of size $n$, $\log(H(\epsilon, \tilde{\mcG}_{B,U}(\delta), S)\leq k_n\, m_n\, \log\left(1 + \frac{2 \delta}{\epsilon}\right)$. Moreover, observe that:
\begin{align}
    \int_{0}^{\delta} \log(H(\epsilon, \tilde{\mcG}_{B,U}(\delta), S) d\epsilon \leq~& \sqrt{\frac{k_n\, m_n}{n}}\int_{0}^{\delta} \sqrt{\log\left(1 + \frac{2 \delta}{\epsilon}\right)} d\epsilon\\
    \leq~& \delta  \sqrt{\frac{k_n\, m_n}{n}}\int_{0}^{1} \sqrt{\log\left(1 + \frac{2}{u}\right)} du = c\, \delta  \sqrt{\frac{k_n\, m_n}{n}}
\end{align}
for some constant $c$. Thus Equation~\eqref{eqn:empirical-covint} is satisfied for $\delta = O\left(\sqrt{\frac{k_n\, m_n}{n}}\right)$. Combining all these we get a projected MSE rate w.p. $1-\zeta$ of:
\begin{equation}
    \|T(\hat{h}-h_0)\|_2 = O\left(\left(\sqrt{\frac{k_n\, m_n}{n}} + \eta_n + \epsilon_n + \sqrt{\frac{\log(1/\zeta)}{n}}\right)\max\{1, \|\theta_*\|_2^2\}\right)
\end{equation}
Invoking standard bounds on the approximation error of classical sieves (e.g. wavelets) and optimally balancing $k_n, m_n$, yields concrete rates (see e.g. \cite{chen2012estimation} for particular approximation rates of known sieves).

Combined with ill-posedness conditions provided in \citep{chen2012estimation}, our results can thus give an alternative proof to the results in \citep{chen2012estimation} that i) do not make minimum eigenvalue conditions, ii) provide adaptivity to $\|\theta_*\|_2$, without knowledge of it, thereby justifying theoretically the use of the regularization term $R(h)$, that was mostly proposed for experimental improvement in \citep{chen2012estimation}. For instance, one concrete ill-posedness condition is that $\lambda_{\min}\left(\E\left[\E[\phi_n(x)\mid z]\E[\phi_n(x)\mid z]^{\top}\right]\right)\geq \gamma_n$ and $\lambda_{\max}\left(\E\left[\psi_n(x)\psi_n(x)^{\top}\right]\right)\leq \sigma_n$. Then the ill-posedness constant is upper bounded by $\tau_n=\sigma_n/\gamma_n$. Moreover, if one assumes a bound on ill-posedness, then \Cref{thm:main-error-ill} requires $\delta$ to be an upper bound of simpler function spaces, that all correspond to simple linear function spaces in finite dimensions. Thus a smaller bound of $O\left(\sqrt{\frac{\max\{k_n,m_n\}}{n}}\right)$, suffices, leading to an error w.p. $1-\zeta$ of the form:
\begin{equation}
    \|\hat{h}-h_0\|_2 = O\left(\tau_n^2 \left( \sqrt{\frac{\max\{k_n,m_n\}}{n}} + \eta_n + \epsilon_n + \sqrt{\frac{\log(1/\zeta)}{n}}\right)\max\{1, \|\theta_*\|_2^2\}\right)
\end{equation}

\section{Application: Reproducing Kernel Hilbert Spaces}\label{app:rkhs}

In this section we deal with the case where $h_0$ lies in a Reproducing Kernel Hilbert space (RKHS) with kernel $K_\mcH:\mcX\times \mcX \to \R$, denoted with $\bbH_K$ and $T h_0$ lies in another RKHS with kernel $K_{\mcF}: \mcZ \times \mcZ \to \R$. We present the three components required to apply our general theory.

First we characterize the set of test functions that are sufficient to satisfy the requirement that $T(h-h_0)\in \mcF_U$; \emph{under non-parametric assumptions on the conditional density $p(x\mid z)$ then we can have $K_{\mcH}=K_{\mcF}$}. Second, by recent results in statistical learning theory, the critical radius of the function classes $\mcF$ and $\mcG$ can be characterized as a function of the \emph{eigendecay of the kernel $K$ and the product kernel $K_{\times}((x, z), (x',z'))=K(x,x')\cdot K(z,z')$ and in the worst-case is of the order of $n^{-1/4}$}. Combining these two facts, we can then apply \Cref{thm:reg-main-error}, to get a bound on the estimation error of the minimax or regularized minimax estimator. Finally, we show that for this set of test functions and hypothesis spaces, \emph{the empirical min-max optimization problem can be solved in closed form}; in particular the inner maximization problem can be shown to correspond roughly to a regularized version of a pairwise metric of the form: $\sum_{i,j} \psi_i K(z_i, z_j) \psi_j$, where $\psi_i = \psi(y_i; h(x_i))$.

\subsection{Characterization of Sufficient Test Functions}

In general, it suffices to assume that the linear operator $T$ is regular enough that it satisfies that for any $h\in \mcH$, we have that $Th \in \bbH_{\KF}$ for some known kernel $\KF$ and that it is an $L$-Lipschitz operator with respect to the pair of RKHS norms $\Hnorm{\cdot}$, $\KFnorm{\cdot}$. Then observe that we  satisfy the requirement that $T(h-h_*)\in \mcF_{L^2\,\Hnorm{h-h_*}^2}$, if we take $\mcF=\bbH_{\KF}$. We now present two complementary sets of sufficient conditions for which the aforementioned property holds.

The first set of conditions applies to a generic function class $\mcH$ and asks principally that $p(x|\cdot)$ belongs to a common RKHS for each $x$.

\begin{lemma}\label{lem:rkhs2}
Suppose that, for each $x$, $p(x|\cdot)$ is an element of an RKHS $\bbH_{\KF}$ and $h \in \mcH$ satisfies $|h(x)| \leq \kappa(x) \Hnorm{h}$ for some $\kappa: \mcX \to \R$.
If $L \defeq \int \kappa(x) \KFnorm{p(x|\cdot)} dx < \infty$,
then $Th \in \bbH_{\KF}$ with
$\KFnorm{Th} \leq L \Hnorm{h}$.
\end{lemma}
\begin{proof}
For any nonnegative $h$, Jensen's inequality implies that
\balign
\label{eq:Th_integrated_norm_bound}
\textstyle
\KFnorm{Th}
    = \Knorm{\int h(x) p(x|\cdot) dx}
    \leq \int |h(x)| \KFnorm{p(x|\cdot)} dx.
\ealign
The same result \cref{eq:Th_integrated_norm_bound} holds for arbitrary signed $h$ due to the decomposition $h = h_+ - h_-$ for $h_+(x) = \max(h(x),0)$ and $h_-(x) = \max(-h(x),0)$, the identity
$|h(x)| = |h_+(x)| + |h_-(x)|$,
and the triangle inequality $\KFnorm{Th} \leq \KFnorm{Th_+} + \KFnorm{Th_-}$.

Now consider any $h \in \mcH$ satisfying $|h(x)| \leq \kappa(x) \Hnorm{h}$ for some $\kappa: \mcX \to \R$
By our inequality \cref{eq:Th_integrated_norm_bound}, we have
\begin{align}
\KFnorm{Th}
    \leq \Hnorm{h} \int \kappa(x) \KFnorm{p(x|\cdot)} dx
    = L \Hnorm{h}.
\end{align}
\end{proof}

The second set of conditions applies when $h$ belongs to a translation-invariant RKHS and ensures that $Th$ belongs to the same RKHS.
Suppose that the kernel $\KH(x, y)=k(x-y)$. Moreover, suppose that $p(x\mid z)=\rho(x-z)$. Then the following lemma states that $Th\in \bbH_{\KH}$ and hence also $T(h-h_*)\in \bbH_{\KH}$ for any $h, h_*\in \bbH_{\KH}$.
\begin{lemma}\label{lem:rkhs}
Suppose the conditional distribution of $X$ given $Z=z$ has continuous density $p(x|z) = \rho(x-z)$ and that $\KH(x,y) = k(x-y)$ for $k$ positive definite and continuous.
If the generalized Fourier transform of $k$ is continuous on $\R^d \backslash \{0\}$,
then $Th \in \bbH_{\KH}$ for all $h \in \bbH_{\KH}$ with $\KHnorm{Th} \leq L \KHnorm{h}$ for $L = \norm{\hat{\rho}}_{\infty}$.
\end{lemma}
\begin{proof}
Fix any $h\in H_K$.
By \cite[Thm.~10.21]{wendland2004scattered}, $\KHnorm{h} = \norm{\hat{h}/\sqrt{\hat{k}}}_2 < \infty$.
Moreover, since $\rho$ is in $L^1$, the Hausdorff-Young inequality implies that $\hat{\rho}\in L^\infty$.
Hence, since $Th = h * \rho$,
\[
\KHnorm{Th}^2
= \int \widehat{Th}(\omega)^2/\hat{k}(\omega) d\omega
= \int \hat{h}(\omega)^2\hat{\rho}(\omega)^2/\hat{k}(\omega) d\omega \leq \norm{\hat{\rho}}_\infty^2 \norm{\hat{h}/\sqrt{\hat{k}}}_2^2
= L^2 \KHnorm{h}^2
<\infty,
\]
so that $Th \in \bbH_{\KH}$ by \cite[Thm.~10.21]{wendland2004scattered}.
\end{proof}

Thus in \Cref{thm:reg-main-error} we can use $\mcH=\mcF=\bbH_K$ for $K=K_{\mcH}$. Moreover, we can set $B$ to be an upper bound on the squared RKSH norm of $h_0$, i.e. $\|h_0\|_{\mcH}^2\leq B$ so that we can take $h_*=h_0$ and have $\|T(h_*-h_0)\|_2=0$, i.e. zero bias. Moreover, by Lemma \ref{lem:rkhs} we also know that $\|Th_0\|_{\mcF}^2\leq L B$ for some constant $L$. Thus we can set $U=2LB$ in Theorem~\ref{thm:reg-main-error} and have that Equation~\eqref{cond:reg-unnormalized} holds with $\eta_n=0$. Thus by \Cref{thm:reg-main-error}, we can get that the estimator in Equation~\eqref{eqn:reg-estimator} satisfies w.p. $1-3\zeta$:
\begin{equation}
\|T(\hat{h}-h_0)\|_2 \leq \delta_n + c_0 \sqrt{\frac{\log(c_1/\delta)}{n}}
\end{equation}
where $\delta_n$ is an upper bound on the critical radii of $\mcF_{6LB}$ and $\mcG_B$, which simplify to:
\begin{align}
    \mcF_{3U} :=~& \left\{f\in \bbH_K: \Knorm{f}^2 \leq  6LB\right\} \\
    \mcG_B :=~& \left\{ (x, z) \to (h(x) - h_0(x))\, T(h-h_0) (z): h\in \bbH_K, \Knorm{h-h_0}^2 \leq B\right\}
\end{align}
Similar rates can also be established for the regularized estimator analogue in \Cref{thm:reg-main-error}, without explicit knowledge of $B$.

\subsection{Critical Radius of $\mcF_{3U}$ and $\mcG_B$}

We now turn to analyze the critical radii of $\mcF_{3U}$ and $\mcG_B$. We first show that these function spaces are also RKHS with appropriate kernels and have bounded RKHS norms. This is trivial for $\mcF_U$. Moreover, observe that the space $\mcG_B$, contains the product of two functions $h f$, where $h:\mcX\to [-1, 1]$ and $f:\mcZ \to [-1, 1]$ and such that $h\in \mcH$ and $f=Th\in \mcF$. Thus the space $\mcG$, with inner product $\ldot{hf}{h'f'}_{\mcG}=\ldot{h}{h'}_{\mcH}\, \ldot{f}{f'}_\mcF$, also admits a reproducing kernel, defined as (see Proposition 12.2 of \cite{wainwright2019high}):
\begin{equation}
K_\mcG((x;z), (x';z'))=\KH(x,x')\, \KF(z,z')
\end{equation}
Moreover, $\|hf\|_{\mcG} = \|h\|_{\mcH}\, \|f\|_{\mcF}$. Thus if $h$, satisfies $\|h\|_{\mcH}^2\leq B$, then by Lemma \ref{lem:rkhs}, $\|Th\|_\mcF^2 \leq L \Knorm{h}^2\leq L B$ for some constant $L$ and $\|hf\|_{\mcG}^2\leq LB^2$.

Assuming that the RKHS spaces $\mcF$ and $\mcG$, also have a sufficiently fast eigendecay then existing results in statistical learning theory also bound the generalization error \cite{wainwright2019high}. In particular, Corollary 14.2 of \cite{wainwright2019high}, shows that for any RKHS $\bbH_K$, if we let
\[
\bbH_K^B:=\{h\in \bbH_K: \Knorm{h}\leq B\},
\]
then we can bound the localized Rademacher and empirical Rademacher complexity as:
\begin{align}
    \mcR(\delta; \bbH_K^B) \leq~& B\sqrt{\frac{2}{n}}\sqrt{\sum_{j=1}^\infty \min\{\lambda_j, \delta^2\}} &
    \mcR_S(\delta; \bbH_K^B) \leq~& B\sqrt{\frac{2}{n}}\sqrt{\sum_{j=1}^{n} \min\{\lambda_j^S, \delta^2\}}
\end{align}
where $\lambda_j$ are the eigenvalues of the kernel and $\lambda_j^S$ are the empirical eigenvalues of the empirical kernel matrix ${\bf K}$ defined as ${\bf K}_{ij}=K(x_i, x_j)/n$. Moreover, the unrestricted Rademacher complexity is upper bounded as (see Lemma 26.10 of \cite{shalev2014understanding}):
\begin{align}
    \mcR(\bbH_K^B) \leq~& O\left(B\sqrt{\frac{\max_{x\in \mcX} K(x, x)}{n}}\right)
\end{align}
Thus in the worst case we can take $\delta_n=O\left(\sqrt{B}\left(\frac{\max_{x\in \mcX} K(x, x)}{n}\right)^{1/4}\right)$, to get a non-parametric rate of convergence.\footnote{Observe that: $K(x, x)=\sum_{j=1}^\infty \lambda_j e_j(x)^2$ and therefore: $\sum_{j=1}^\infty \lambda_j = \sum_{j=1}^{\infty} \lambda_j \E_x[e_j(x)^2] = \E_x[\sum_{j=1}^{\infty}\lambda_j e_j(x)^2] = \E_x[K(x,x)]\leq \max_{x\in \mcX} K(x,x)$. Thus in the worst case, when $\lambda_j \geq \delta^2$ for most $j$, we still recover the non-localized from the localized bounds.} However, for many kernels, the eigendecay will be sufficiently fast, that $\delta^2$ will not be binding in the minimum. For instance, for the Gaussian kernel in one dimension on the domain $[0, 1]$, with bandwidth of $1$, i.e. $K(x,x' )=e^{-\frac{(x-x')^2}{2}}$, we have that $\delta_n=O\left(B\sqrt{\frac{\log(n+1)}{n}}\right)$ (see Example 14.4 of \cite{wainwright2019high}).

\paragraph{Data-adaptive estimation} Moreover, by Equation~\eqref{eqn:emp-v-pop-critical}, we can choose $\delta$ in \Cref{thm:reg-main-error} based on the empirical critical radius. Observe that the empirical eigenvalues are directly computable from the data and hence, we can calculate a data-adaptive quantity $\hat{\delta}_n$ and choose $\delta$ in \Cref{thm:reg-main-error}, based on this data-adaptive quantity plus an $O\left(\sqrt{\frac{\log(1/\zeta)}{n}}\right)$ term. Moreover, if we use the regularized estimator, then we also do not require knowledge of $B$, which leads to a very data-adaptive estimation scheme. The only thing required is knowledge of an upper bound on the Lipschitz constant $L$ of the operator $T$ with respect to the RKHS norm.

\subsection{Closed-Form Solution to Optimization Problem}\label{app:rkhs-optimization}
Finally, we show that the optimization problem that defines the estimator in Equation~\eqref{eqn:reg-estimator} can be computed in closed form. We present the results for the constrained estimator, but exact analogues also hold for the regularized version.
The proof can be found in \cref{sec:proof-rkhs-closed-form-max}.

\begin{proposition}[Closed-form maximization]\label{rkhs-closed-form-max}
Suppose $\mcF$ is an RKHS with kernel $K$ equipped with the canonical RKHS norm $\Fnorm{\cdot} = \Knorm{\cdot}$. Then for any $h$
\begin{align}\label{eqn:maximization}
    \sup_{f\in \mcF}\,\, \Psi_n(h, f)^2 - \lambda\left(\KFnorm{f}^2 + \frac{U}{n\delta^2}\|f\|_{2,n}^2\right)
    &= \frac{1}{4\lambda} \psi_n^\top  K_n^{1/2}(\textfrac{U}{n\delta^2}K_n + I)^{-1} K_n^{1/2} \psi_n\\
    &= \frac{1}{4\lambda} \psi_n^\top  K_n\, (\textfrac{U}{n\delta^2}K_n + I)^{-1} \psi_n
\end{align}
where $K_n = (K(z_i,z_j))_{i,j=1}^n$ is the empirical kernel matrix and $\psi_n = (\frac{1}{n}\psi(y_i\midsem{} h(x_i)))_{i=1}^n$.
\end{proposition}

We note that if we did not enforce the extra $\ell_{2,n}$ norm constraint on $f$ (i.e. $\delta\to\infty$, then the above inner optimization problem simplifies to:
\begin{equation}\label{eqn:kloss-rkhs}
    \sup_{f\in \mcF}\,\, \Psi_n(h, f)^2  - \lambda\, \KFnorm{f}^2 = \frac{1}{4\lambda} \psi_n^\top  K_n \psi_n = \frac{1}{4\lambda\, n^2}\sum_{i, j} \psi(y_i;h(x_i)) K(z_i, z_j) \psi(y_j; h(x_i))
\end{equation}
i.e. we get a pair-wise residual loss, weighted by a kernel matrix that is only a function of the conditioning set $z$.

Thus the solution $\hat{h}$ of the estimator in Equation~\eqref{eqn:reg-estimator} is equivalent to:
\begin{equation}
    \hat{h} = \argmin_{h\in \mcH} \frac{1}{4\lambda} \psi_n^\top  M \psi_n + \mu \|h\|_{\mcH}^2 = \argmin_{h\in \mcH} \psi_n^\top  M \psi_n + 4\mu\, \lambda \|h\|_{\mcH}^2
\end{equation}
where $M:= K_n^{1/2}(\textfrac{U}{n\delta^2}K_n + I)^{-1} K_n^{1/2}$. Finally, we show that this outer maximization also has a closed form solution. See \cref{sec:proof-closed-form-min} for the proof.

\begin{proposition}[Closed-form minimization]\label{prop:closed-form-min}\label{prop:closed-form-min-reg}
Suppose that $\mcH$ and $\mcF$ are the RKHSes of the kernels $K_\mcH$ and $K_\mcF$, equipped with the canonical RKHS norms $\Hnorm{\cdot} = \KHnorm{\cdot}$ and $\Fnorm{\cdot} = \KFnorm{\cdot}$.
Define the empirical kernel matrices $K_{\mcH,n} = (K_{\mcH}(x_i,x_j))_{i,j=1}^n$ and $K_{\mcF,n} = (K_{\mcF}(z_i,z_j))_{i,j=1}^n$.
Then the following estimator is an optimizer of Equation~\eqref{eqn:reg-estimator}:
\begin{align}\label{eqn:optimal-min-vector-rkhs-reg}
\hat{h} =~& \sum_{i=1}^n \alpha_{\lambda_*, i} \KH(x_i, \cdot)
&
\alpha_{\lambda}
:=~& (\KHn M\KHn  + 4\,\lambda\, \mu \KHn)^{\dagger}\KHn M y
\end{align}
for $M = K_{\mcF,n}^{1/2}(\textfrac{U}{n\delta^2}K_{\mcF,n} + I)^{-1} K_{\mcF,n}^{1/2} \equiv K_{\mcF,n}(\textfrac{U}{n\delta^2}K_{\mcF,n} + I)^{-1}$ and $A^\dagger$ is the Moore-Penrose pseudoinverse of a matrix $A$.

\end{proposition}

\paragraph{Hyper-parameter tuning} Observe that \Cref{thm:reg-main-error} states that as long as the regularization strength satisfies that $\lambda \mu =\Theta(\delta^4 L^2)$, then this estimator will provide results that automatically scale with the RKHS norm of true hypothesis $h_0$. Moreover, the regularization hyperparameter $\lambda\cdot \mu$ can also be tuned in practice by evaluating the loss function $\psi_n^\top  M \psi_n$ on a left-out sample, with parameters $n, \delta$ set to the appropriate ones for the size of that sample.

\begin{figure}
    \centering
    \begin{subfigure}[b]{0.3\textwidth}
        \includegraphics[width=\textwidth]{figures/2dpoly_fn_rkhs.png}
        \caption{$-1.5\cdot x + .9\cdot x^2$}
        \label{fig:2dpoly}
    \end{subfigure}
    ~

    \begin{subfigure}[b]{0.3\textwidth}
        \includegraphics[width=\textwidth]{figures/step_fn_rkhs.png}
        \caption{$1+1.5\cdot 1\{x>0\}$}
        \label{fig:step}
    \end{subfigure}
    ~

    \begin{subfigure}[b]{0.3\textwidth}
        \includegraphics[width=\textwidth]{figures/sin_fn_rkhs.png}
        \caption{$\sin(x)$}
        \label{fig:sin}
    \end{subfigure}
    \begin{subfigure}[b]{0.3\textwidth}
        \includegraphics[width=\textwidth]{figures/abs_fn_rkhs.png}
        \caption{$|x|$}
        \label{fig:abs}
    \end{subfigure}
    \begin{subfigure}[b]{0.3\textwidth}
        \includegraphics[width=\textwidth]{figures/sigmoid_fn_rkhs.png}
        \caption{$\frac{2}{1+e^{-2x}}$}
        \label{fig:sigmoid}
    \end{subfigure}
    \begin{subfigure}[b]{0.3\textwidth}
        \includegraphics[width=\textwidth]{figures/3dpoly_fn_rkhs.png}
        \caption{$-1.5\cdot x + .9\cdot x^2 + x^3$}
        \label{fig:3dpoly}
    \end{subfigure}
    \caption{Estimated functions based on our minimax estimator for different true functions. We use an rbf kernel with parameter $\gamma=.1$ and $1000$ samples. We chose critical radius parameter $\delta=5/n^{.4}$ and the regularization hyper-parameter $\tau$ is chosen via k-fold cross-validation. The data generating process was: $x=.6\, z + .4\, u + \delta$ and $y=h_0(x) + u + \epsilon$ and $z,u\sim N(0, 2)$ and $\epsilon, \delta\sim N(0, .1)$.}\label{fig:rkhs}
\end{figure}

\paragraph{Low-Rank Approximation and Nystrom's Method} The solution to the empirical optimization problem requires inverting an $n\times n$ kernel matrix, which takes time $O(n^3)$. This can be prohibitive for moderate sample sizes of the order of tens of thousands. We note here that one can construct very good approximations to the solution in \Cref{prop:closed-form-min} by considering low-rank approximations of the kernel matrix $K$. We present here one such low-rank approximation, based on Nystrom's method, but we note that the plethora of recent literature on low-rank kernel approximation methods are applicable to our problem too (see e.g. \cite{kumar2012sampling,bach2005predictive,musco2017recursive,oglic2017nystrom}).

Suppose that we can express our kernel matrices as $\KHn$ and $\KFn$ as $\KFn = D D^\top $ and $\KHn= V V^\top $, where $D$ and $V$ are of dimensions $n\times r$ and such that we can express the \emph{kernel row} of any new test sample as:
\begin{align}
(\KH(x_1, x), \ldots, \KH(x_n, x)) =  V \phi(x)
\end{align}
for some $r$-dimensional vector $\phi(x)$. Then we can express $h(x) = \phi(x)^\top  V^\top  a_\lambda$. If we then define $\gamma = V^\top  \alpha_\lambda$. Then we can re-write the closed form solutions to the min and max problems as follows:
\begin{align}
    \sup_{f\in \mcF}\,\, \Psi_n(h, f)^2 - \lambda\left(\KFnorm{f}^2 + \frac{U}{n\delta^2}\|f\|_{2,n}^2\right) = \frac{1}{4\lambda} \psi_n^\top  D \left(\frac{U}{n\delta^2}D^\top D + I \right)^{-1} D^\top  \psi_n
\end{align}
and if we let $Q:=\left(\frac{U}{n\delta^2}D^\top D + I \right)^{-1}$ and $A=V^\top D$, then:
\begin{align}
    \gamma :=~& \left( A Q A^\top  + 4\,\lambda\, \mu\, I\right)^{-1} A Q D^\top y\\
    \hat{h}(x) :=~& \phi(x)^\top \gamma
\end{align}
Observe that every matrix calculation in the above expressions requires time at most $O(n\, r^2)$ to be computed. Thus if $r\ll n$, we have massively reduced the computation time from $\Theta(n^3)$ to $O(n\, r^2)$, making the method practical even very large data regimes.

Even though $r$ in the worst-case can be of size $n$, we can typically well-approximate the kernel matrices with $r \ll  n$. One popular approach for achieving this is Nystrom's method, which essentially sub-samples a set of $r$ points and uses the normalized kernel distances with respect to this subset of points as $D$ and $V$, respectively. In particular, let $S$ denote an $n\times r$ matrix whose $i$-th column contains a $1$ in position $j$ for some randomly sampled index $j$. Then $KS$ is an $n\times r$ sub-matrix of $K$, where a subset $S$ of the columns of $K$ are chosen at random.\footnote{Several sampling strategies have been proposed in the literature to improve upon pure uniform sampling (see e.g. \cite{kumar2012sampling,musco2017recursive,oglic2017nystrom}). One popular practical and simple method is to perform some version of unsupervised clustering of the samples, such as kmeans clustering, and choosing the points as the cluster centroids.} Then we can approximate $K$ via $VV^\top $, where $V = KS M^{1/2}$ and $M = (S^\top KS)^+$ (i.e. $V$ is contains normalized kernel-based similarities to the subset $S$ of $r$ randomly chosen points). Moreover, for any new test point, we can set $\phi(x)=M^{1/2} (\KH(x_i, x))_{i\in S}$.

\begin{figure}
    \centering
    \begin{subfigure}[b]{0.3\textwidth}
        \includegraphics[width=\textwidth]{figures/nystrom_sin_fn_rkhs.png}
        \caption{$\sin(x)$}
        \label{fig:nystrom_sin}
    \end{subfigure}
    ~

    \begin{subfigure}[b]{0.3\textwidth}
        \includegraphics[width=\textwidth]{figures/nystrom_step_fn_rkhs.png}
        \caption{$1+1.5\cdot 1\{x>0\}$}
        \label{fig:nystromg_step}
    \end{subfigure}
    ~

    \begin{subfigure}[b]{0.3\textwidth}
        \includegraphics[width=\textwidth]{figures/nystrom_3dpoly_fn_rkhs.png}
        \caption{$-1.5\cdot x + .9\cdot x^2 + x^3$}
        \label{fig:nystrom_3dpoly}
    \end{subfigure}
    \caption{Estimates based on Nystrom approximation, with $50$ nystrom samples, for the same dgp and parameter setup as in Figure~\ref{fig:rkhs}.}\label{fig:nystrom}
\end{figure}

\subsection{Bounds on Ill-Posedness Measure}\label{app:rkhs-illposedness}
The results so far in the section provide bounds on the projected RMSE. In this last section, we show that under further assumptions on the strength of the instrument (i.e. the correlation of $x$ and $z$), then the projected RMSE rates also imply rates for the RMSE. We give an example such set of conditions, mostly as an example of a sufficient set of assumptions that lead to RMSE rates and in order to provide qualitative insights on what RMSE rates one can expect in different regimes of the instrument strength and the eigendecay of the kernel. In this section we will assume that the space $\mcH$ is also augmented with a hard constraint on the RKHS norm, i.e. $\mcH = \bbH_K^B = \{h\in \bbH_K: \|h\|_{K} \leq B\}$. Assuming $\|h_0\|_{K}\leq B$ this does not change the statistical guarantees and moreover the closed form optimization theorems, can easily be amended to incorporate a hard constraint on top of the regularization (due to the equivalent between hard constraints and regularization). Imposing this hard constraint will simplify the analysis of this section.\footnote{We note that the proof of \Cref{thm:reg-main-error} implies that even without a hard constraint, with high probability $\|\hat{h}\|_{K}^2 \leq \|h_0\|_K^2 + \frac{\delta^2 + \lambda U}{\mu}$. Thus the results of this section hold for $B=\|h_0\|_K^2 + \frac{\delta^2 + \lambda U}{\mu}$ even without the extra hard constraint.}

By Mercer's theorem we can express any function in the RKHS $\bbH_K^B$, in terms of the eigenfunctions of the kernel:
\begin{equation}
    h = \sum_{j\in J} a_j e_j
\end{equation}
with $e_j: \mcX \to \R$, such that $\E[e_j(x)^2]=1$ and $\E[e_i(x)\, e_j(x)]=0$ and $J$ a countable set. Moreover, we have $\|h\|_2^2 = \sum_{j\in J} a_j^2$ and $\Knorm{h}=\sum_{j\in J} \frac{a_j^2}{\lambda_j}\leq B$. Thus we have that $\|h\|_{\mcH}^2 \leq B$ implies that for all $m\in \mathbb{N}_+$: $\sum_{j\geq m} a_j^2 \leq \lambda_m B$. Moreover, we have:
\[
\|Th\|_2^2 = \sum_{i, j\in J} a_i a_j \E[\E[e_i(x)\mid z]\E[e_{j}(x)\mid z]].
\]

For any $m\in \mathbb{N}_+$, let $I:=\{1, \ldots, m\}$, $e_I=(e_1, \ldots, e_m)$, $a_I=(a_1, \ldots, a_m)$ and:
\begin{equation}
    V_m := \E[\E[e_I(x)\mid z]\, \E[e_I(x) \mid z]^\top ]
\end{equation}
and suppose that $\lambda_{\min}(V_m) \geq \tau_{m}$, i.e. that these finite eigenfunctions maintain some fraction of their independent components, even when they are smoothened through the conditional expectation $p(x\mid z)$. Furthermore suppose that for all $i \leq m < j$: $\left|\E[\E[e_i(x)\mid z]\E[e_{j}(x)\mid z]]\right| \leq \gamma_m \leq c\, \tau_m$ (for some constant $c$), i.e. the smoothening performed by the conditional expectation does not ruin a lot the orthogonality of the first $m$ eigenfunctions with eigenfunctions for indices larger than $m$. Observe that if we had a perfect instrument, i.e. $z$ was perfectly correlated with $x$, then $V_m=I_m$ and $\E[\E[e_i(x)\mid z]\E[e_{j}(x)\mid z]]=\E[e_i(x) e_j(x)]=0$. Thus for a perfect instrument $\tau_m=1$ and $\gamma_m=0$. Therefore the latter requirements are implicit assumptions on the strength of the instrument.\footnote{Potentially the strongest assumption of these is that $\gamma_m\leq \tau_m$. This could be avoided by restricting the hypothesis space $\mcH_B$ to only be supported on the first $m$ eigenfunctions. However, this would require being able to diagonalize the kernel and also to tune the estimator to the unknown parameters $\tau_m$.} We show that under these assumptions, we can bound the measure of ill-posedness as follows.

\begin{lemma}\label{lem:ill-posedness-rkhs}
Suppose that $\lambda_{\min}(V_m) \geq \tau_m$ and for some constant $c>0$, for all $i\leq m < j$,
\[
\left|\E[\E[e_i(x)\mid z]\E[e_{j}(x)\mid z]]\right| \leq c\, \tau_m
\]
Then:
\begin{equation}
    \tau^*(\delta)^2 := \max_{h\in \bbH_K^B: \|Th\|_2\leq \delta} \|h\|_2^2 \leq \min_{m\in \mathbb{N}_+} \left(\frac{4\delta^2}{\tau_m} + (4c^2+1) B \lambda_{m+1} \right)
\end{equation}
\end{lemma}

The optimal choice of $m_*$ roughly solves the equation: $\tau_m \lambda_{m+1}= \delta^2/B$. If for instance $\lambda_m \leq m^{-b}$ for $b>1$, and $\tau_{m} \geq m^{-a}$ for $a>0$, then: $m_* \sim \delta^{2/(a+b)}$, leading to a rate of:
\begin{equation}
    \|\hat{h}-h_*\|_2 = O\left(\delta^{b/(a+b)}\right)
\end{equation}
We see that the RMSE rate is of a slower order than the projected MSE rate. If $\lambda_m$ has an exponential eigendecay, i.e. $\lambda_m \sim 2^{-m}$ (e.g. such as in the case of a Gaussian kernel), and $\tau_m \geq m^{-a}$, then $m_* \sim \log(1/\delta^2)$ and we get:
\begin{equation}
    \|\hat{h}-h_*\|_2 = O\left(\delta\, (\log(1/\delta))^{a/2}\right)
\end{equation}
Thus we only get a logarithmic increase in the RMSE rate as compared to the Projected RMSE rate. However, we note that if also $\tau_m \sim 2^{-a\, m}$ and $\lambda_m \sim 2^{-b\, m}$, then we get rates of $O\left(\delta^{b/(a+b)}\right)$, by settings $m_* \sim \log(1/\delta^{2/(a+b)})$. Finally, in the severely ill-posed setup, where $\tau_m \sim 2^{-m}$ and $\lambda_m \sim m^{-b}$, then we have $m_* \sim \log(1/\delta^2)$ and:
\begin{equation}
    \|\hat{h}-h_*\|_2 = O\left( \frac{1}{\log(1/\delta)^{b}} \right)
\end{equation}
leading to a very slow rate of convergence that will typically be of the order of $1/\log(n)$.

Observe that we achieve the rate for the optimal choice of $m$, without the need to tune our algorithm. The RKHS norm penalty implicitly clips the weight that our functions can put on eigenfunctions with large index and hence controls the measure of ill-posedness for whatever is the decay rates of the eigenvalues $\lambda_m$ and $\tau_m$.

\section{Application: High-Dimensional Sparse Linear Function Spaces}\label{app:high-dim-linear}

In this section we deal with high-dimensional linear function classes, i.e. the case when $\mcX, \mcZ\subseteq \R^p$ for $p\gg n$ and $h_0(x) = \ldot{\theta_0}{x}$. We will address the case when the function $\theta_0$ is assumed to be sparse, i.e. $\|\theta_0\|_{0}:=\{j\in [p]: |\theta_j|>0\}\leq s$. We will be denoting with $S$ the subset of coordinates of $\theta_0$ that are non-zero and with $S^c$ its complement. For simplicity of exposition we will also assume that $\E[x_i\mid z]=\ldot{\beta}{z}$, though most of the results of this section also extend to the case where $\E[x_i\mid z]\in \mcF_i$ for some $\mcF_i$ with small Rademacher complexity. We provide two sets of results, dependent on whether we make further minimum eigenvalue assumptions on the covariance matrix of the random variables $\E[x_i\mid z]$.

\subsection{Hard Sparsity Constraints without Minimum Eigenvalue}\label{app:high-dim-linear-hard}

In the first result, we apply \Cref{thm:reg-main-error} to show that even without any further assumptions on the eigenvalues of the covariance matrix
\[
V:=\E[\E[x\mid z]\E[x\mid z]^\top ],
\]
we can attain fast rates of the order of $n^{-1/2}$ that are logarithmic in $p$ and only linear in the sparsity $s$ of $h_0$ and the sparsity $r$ of the conditional expectation functions $\E[x_i\mid z]$. Albeit the optimization problem we need to solve to get these rates is non-convex and has running time that is exponential in $r, s$. This setting covers and extends the linear moment case of the setting analyzed in \citep{fan2014endogeneity}; albeit we only provide RMSE and projected RMSE rates.

\begin{corollary}
Suppose that $h_0(x)=\ldot{\theta_0}{x}$ with $\|\theta_0\|_0\leq s$ and $\E[x_i\mid z]=\ldot{\beta_{0}^i}{z}$ with $\|\beta_{0}^{i}\|_0\leq r$. Then let $\mcH$ consist of all $s$-sparse linear functions of $x$ and $\mcF$ consist of all $(s\cdot r)$-sparse linear functions of $z$ with coefficients in $[-1,1]$. in $p$ dimensions with only $s$ non-zero coefficients and $\mcF$ consists of linear functions in $q$ dimensions with $r$ non-zero coefficients. Then the estimator presented in Equation~\eqref{eqn:reg-estimator}, satisfies that w.p. $1-\zeta$:
\begin{equation}
    \|T(\hat{h}-h_0)\|_2 \leq O\left( \sqrt{\frac{r\, s \log(p\, n)}{n}} + \sqrt{\frac{\log(1/\zeta)}{n}}\right)
\end{equation}
\end{corollary}

The proof follows immediately from the fact that the metric entropy of $r\, s$-sparse linear functions in $p$-dimensions, with coefficients in $[-1,1]$ is of the order of $O\left(r\, s \log(p/\epsilon)\right)$. Thus we can invoke \Cref{lem:covering} to get a bound of $O\left(\sqrt{\frac{r\, s\, \log(p\, n)}{n}}\right)$ on the critical radii of classes $\mcF_{3U}$ and $\mcG_{B, U}$ and apply \Cref{thm:reg-main-error}.

\subsection{$\ell_1$-Relaxation under Minimum Eigenvalue Condition}\label{app:high-dim-linear-ell1}

In the second set of results we assume a restricted minimum eigenvalue of $\gamma$ on the matrix $V$ and apply \Cref{thm:reg-main-error-2} to get fast rates of the order of $n^{-1/2}$, that also scale logarithmically in $p$, linearly in $r, s$ and $\gamma^{-1}$. Moreover, the optimization problem required is now a convex problem as we replace the hard sparsity constraint with an $\ell_1$ constraint. This dichotomy of computationally efficient vs computationally hard estimation dependent on whether we make minimum eigenvalue assumptions is a well established result in exogenous regression problems \citep{zhang2014lower} and hence we provide here analogous positive results for the endogenous regression setup. We also note that without the minimum eigenvalue condition, our \Cref{thm:reg-main-error} still provides slow rates of the order of $n^{-1/4}$, for computationally efficient estimators that replace the hard sparsity constraint with an $\ell_1$-norm constraint. Our results based on the $\ell_1$-constraint are also closely related to the work of \cite{gautier2011high}, who analyzes an endogenous analogue of the Dantzig selector. Our work proposes an alternative to the Dantzig selector that enjoys similar estimation rate guarantees.

\begin{repcorollary}{cor:sparse-linear-reg-ell1}
Suppose that $h_0(x)=\ldot{\theta_0}{x}$ with $\|\theta_0\|_0\leq s$ and $\|\theta_0\|_1\leq B$ and $\|\theta_0\|_{\infty}\leq 1$. Moreover, suppose that $\E[x_i\mid z] = \ldot{\beta_0^i}{z}$, with $\beta_0^i\in \R^p$ and $\|\beta_0^i\|_1\leq U$ and that the co-variance matrix $V$ satisfies the following restricted eigenvalue condition:
\begin{equation}
    \forall \nu\in \R^p \text{ s.t. } \|\nu_{S^c}\|_1 \leq \|\nu_S\|_1 + 2\,\delta_{n,\zeta}: \nu^\top  V\nu \geq \gamma \|\nu\|_2^2
\end{equation}
Then let $\mcH = \{x\to \ldot{\theta}{x}: \theta \in \R^p\}$, $\|\ldot{\theta}{\cdot}\|_{\mcH}=\|\theta\|_1$, $\mcF_U=\{z \to \ldot{\beta}{z}: \beta\in \R^p, \|\beta\|_1\leq U\}$ and $\|\ldot{\beta}{\cdot}\|_{\mcF}=\|\beta\|_1$. Then the estimator presented in Equation~\eqref{eqn:reg-estimator-2} with $\lambda\leq \frac{\gamma}{8s}$, satisfies that w.p. $1-\zeta$:
\begin{equation}
    \|T(\hat{h}-h_0)\|_2 \leq O\left( \max\left\{1, \frac{1}{\lambda}\frac{\gamma}{s}\right\} \sqrt{\frac{s}{\gamma}} \left((B + U + 1)\sqrt{\frac{\log(p)}{n}} + \sqrt{\frac{\log(p/\zeta)}{n}}\right)\right)
\end{equation}
If instead we assume that $\|\beta_0^i\|_2\leq U$ and $\sup_{z\in \mcZ} \|z\|_2\leq R$ then by setting $\mcF_U=\{z\to \ldot{\beta}{z}: \|\beta\|_2\leq U\}$ and $\|\ldot{\beta}{\cdot}\|_\mcF= \|\beta\|_2$, we have:
\begin{equation}
    \|T(\hat{h}-h_0)\|_2 \leq O\left( \max\left\{1, \frac{1}{\lambda}\frac{\gamma}{s}\right\} \sqrt{\frac{s}{\gamma}} \left((B + 1)\sqrt{\frac{\log(p)}{n}} + \frac{U\, R}{\sqrt{n}} + \sqrt{\frac{\log(p/\zeta)}{n}}\right)\right)
\end{equation}
\end{repcorollary}

\paragraph{Second order influence from $\E[x_i\mid z]$ model complexity} Notably, observe that in the case of $\|\beta_0^i\|_2\leq U$, we note that if one wants to learn the true $\beta$ with respect to the $\ell_2$ norm or the functions $\E[x_i\mid z]$ with respect to the RMSE, then the best rate one can achieve (by standard results for statistical learning with the square loss), even when one assumes that $\sup_{z\in \mcZ} \|z\|_2\leq R$ and that $\E[zz^{\top}]$ has minimum eigenvalue of at least $\gamma$, is: $\min\left\{\sqrt{\frac{p}{n}}, \left(\frac{U\, R}{n}\right)^{1/4}\right\}$. For large $p\gg n$ the first rate is vacuous. Thus we see that even though we cannot accurately learn the conditional expectation functions at a $1/\sqrt{n}$ rate, we can still estimate $h_0$ at a $1/\sqrt{n}$ rate, assuming that $h_0$ is sparse. Therefore, the minimax approach offers some form of robustness to nuisance parameters, reminiscent of the type of robustness of Neyman orthogonal methods (see e.g. \citep{Chernozhukov2018double}).

\subsection{Solving the $\ell_1$-Relaxation Optimization Problem via First-Order Methods}\label{app:high-dim-linear-opt}

The estimator presented in \Cref{cor:sparse-linear-reg-ell1} require solving optimization problems of the form:
\begin{align}\label{eqn:minimax-ell1}
\min_{\theta: \|\theta\|_1\leq B} & \max_{\beta: \|\beta\| \leq U} \ldot{\E_n\left[(y - \ldot{\theta}{x})z\right]}{\beta} + \mu \|\theta\|_1
\end{align}
for some $R, \mu$ and for norm $\|\cdot\|$ either $\|\cdot\|_1$ or $\|\cdot\|_2$ (in the constrained estimator $\mu=0$; while in the regularized $R=\infty$ - though in practice we can set it to some large value for stability of the optimization process). Observe that inner optimization simplifies to:
\begin{align}
   \min_{\theta: \|\theta\|_1\leq B} & \left\|\E_n\left[(y - \ldot{\theta}{x})z\right]\right\|_{*} +  \frac{\mu}{U} \|\theta\|_1
\end{align}
where $\|\cdot\|_*$ is the dual norm of $\|\cdot\|$ (i.e. the $\ell_{\infty}$ norm in the case where $\|\cdot\|$ is the $\ell_1$ norm and the $\ell_2$ norm in the case where $\|\cdot\|$ is the $\ell_2$ norm). One approach to solving these optimization problems is using projected sub-gradient descent:
\begin{align}
    \beta_t =~& \argmax_{\beta: \|\beta\| \leq U} \ldot{\E_n\left[(y - \ldot{\theta_t}{x})z\right]}{\beta}\\
    \theta_{t+1} =~& \Pi\left(\theta_t + \eta\, \E_n\left[x\, z^\top \right]\beta_t\, - \frac{\mu}{U} \sign(\theta_t)\right)\\
    \Pi(\theta) =~& \argmin_{\theta' : \|\theta' \|_1\leq B} \|\theta - \theta' \|_2
\end{align}
Moreover, for both $\ell_1$ and $\ell_2$ norm, the solution to $\beta_t$ can be easily found in closed form.\footnote{For the case of the $\ell_1$ norm: $\beta_t=U e_{i_t} \sign(\E_n[y-\ldot{\theta_t}{x})z_{i_t}])$, with $i_t=\argmax_{i} |\E_n[(y-\ldot{\theta_t}{x})z_i]|$. For the case of the $\ell_2$ norm: $\beta_t = \E_n[(y-\ldot{\theta_t}{x})z]\cdot U/\|\E_n[(y-\ldot{\theta_t}{x})z]\|_2$} After $O(1/\epsilon^2)$ iterations and for $\eta=\Theta(\epsilon)$, we will have that $\bar{\theta}=\frac{1}{T}\sum_{t=1}^T \theta_t$, is an $\epsilon$-approximate solution to the optimization problem.

\paragraph{Improved Iteration Complexity with Optimistic FTRL Dynamics} The sub-gradient descent approach has two caveats: i) the rate of $1/\epsilon^2$ is considerably slow and would require a large number of iterations to converge to a reasonable solution, ii) the gradient does not admit an unbiased stochastic version (due to the non-linearity introduced by the $\argmax$ operation that defines $\beta_t$), and therefore the algorithm does not admit a stochastic variant, which is useful for large samples. We can improve the error rate by invoking algorithms that address non-smooth optimization problems that take the form of a min-max objective of some underlying smooth loss.

First, we show that we can remove the non-smoothness of the $\ell_1$-regularization by lifting the parameter $\theta$ to a $2p$-dimensional positive orthant. Consider two vectors $\rho^+, \rho^-\geq 0$ and then setting $\theta=\rho^+ - \rho^-$, with $\rho = (\rho^+; \rho^-)$ and $\|\rho\|_1\leq B$. Observe that for any feasible $\theta$, the solution $\rho_i^+=\theta_i 1\{\theta_i>0\}$ and $\rho_i^-=\theta_i 1\{\theta_i \leq 0\}$ is still feasible and achieves the same objective. Moreover, any solution $\rho$, maps to a feasible solution $\theta$ (since $\|\theta\|_1 \leq \|\rho_+ - \rho_-\|_1\leq \|\rho^+\|_1 +\|\rho^-\|_1\leq B$) and thus the two optimization programs have the same optimal solutions. Then, if we define with $v=(x; -x)$, then the optimization problem can be re-stated as:
\begin{equation}
    \min_{\rho\geq 0: \|\rho\|_1\leq B}  \max_{\beta: \|\beta\|\leq U} \ell(\rho, \beta)
\end{equation}
where:
\begin{align}
    \ell(\rho, \beta) :=~& \beta^\top \E_n[zy] - \beta^\top  \E_n[zv^\top ] \rho + \mu \sum_{i=1}^{2p} \rho_i
\end{align}
This falls exactly into the class of problems analyzed in a line of work on bi-linear minimax optimization, starting from the seminal work of \cite{nesterov2005smooth}. For instance, we can view the problem as a two-player bi-linear zero-sum game and invoke the Optimistic Follow-the-Regularized-Leader (OFTRL) or Optimistic Mirror Descent (OMD) paradigm of \cite{Rakhlin2013,syrgkanis2015fast}, to find an $\epsilon$-approximate solution for $\rho$ in $O(1/\epsilon)$ iterations. The algorithm repeats for $T$ iterations the updates:
\begin{align}
    \rho_{t+1} =~& \argmin_{\rho \geq 0: \|\rho\|_1\leq B} \sum_{\tau\leq t} \ell(\rho, \beta_{\tau}) + \ell(\rho, \beta_t) + \frac{1}{\eta} R_{\min}(\rho)\\
    \beta_{t+1} =~& \argmax_{\beta: \|\beta\|_1\leq U} \sum_{\tau\leq t} \ell(\rho_{\tau}, \beta) + \ell(\rho_{t}, \beta) - \frac{1}{\eta} R_{\max}(\beta)
\end{align}
and returns $\bar{\rho}=\frac{1}{T} \sum_{t=1}^T \rho_t$, $\bar{\beta}=\frac{1}{T}\sum_{t=1}^T \beta_t$.\footnote{Finally, if we want to compare with $s$-sparse solutions and we want to enhance sparsity of the returned solution, then we can always truncate to zero at the end of training any coordinate of $\bar{\theta}=\bar{\rho}^+-\bar{\rho}^-$ that was smaller than $1/(s\, n^{1/2+\epsilon})$. This can introduce an extra lower order approximation error of at most $1/n^{1/2+\epsilon}$ in our projected MSE theorem, since by this shrinkage procedure, the error with respect to a sparse solution $\theta_0$ can only increase on the non-zero entries of $\theta_0$ and it can only increase by at most $1/(sn^{1/2+\epsilon})$ on every such entry.} We note that if we did not double count the last period's loss and we used $R_{\min}(x)=R_{\max}(x)=\frac{1}{2}\|x\|_2^2$, then this would correspond to running simultaneous gradient descent dynamics for both parameters $\rho, \beta$. Moreover, the parameters $\bar{\rho}, \bar{\beta}$ can be thought as primal and dual solutions and we can use the duality gap as a certificate for convergence of the algorithm.\footnote{In particular, $\bar{\rho}$ and $\bar{\beta}$ are an $\epsilon$-equilibrium of the zero-sum game.}
\begin{equation}
    \text{tol} = \max_{\beta: \|\beta\| \leq U} \ell(\bar{\rho}, \beta) - \min_{\rho: \|\rho\|_1\leq B} \ell(\rho, \bar{\beta})
\end{equation}

This approach addresses both problems with projected sub-gradient descent: i) as we will show below, the iteration complexity is $O\left((B+U^2)\log(B\,p)/\epsilon\right)$, instead of $1/\epsilon^2$, ii) the per-iteration losses $\ell(\rho, \beta_t)$, $\ell(\rho_t, \beta)$ in the FTRL formulation can be replaced with unbiased estimates, while still maintaining theoretical guarantees and therefore the algorithm admits a stochastic analogue which makes it scalable to very large data sets.\footnote{We note that the fast rate of $1/\epsilon$ will deteriorate with the size of the mini-batch, but a $1/\epsilon^2$ rate is always achievable and the step-size $\eta$ should be appropriately tuned to account for the mini-batch sampling noise.}

To instantiate this paradigm we need to find appropriate regularizers for the strategy spaces of the two players. Below we outline two concrete such algorithms for the two cases of the norm of $\beta$ and provide worst-case convergence rates.

\paragraph{$\ell_1$-ball adversary} For the case when $\|\beta\|=\|\beta\|_1$, we can further simplify the problem by showing that the inner optimization can be performed over a $2p$-dimensional simplex. If we let $u = (z; -z)$, then we can re-write the optimization problem as:
\begin{align}
    \ell(\rho, w) :=~& w^\top \E_n[uy] - w^\top  \E_n[uv^\top ] \rho + \frac{\mu}{U} \sum_{i=1}^{2p} \rho_i
\end{align}
\begin{equation}
    \min_{\rho\geq 0: \|\rho\|_1\leq B}  \max_{w: \|w\|= 1} \ell(\rho, w)
\end{equation}
Since both player strategies $\rho$, $w$ are constrained to be in an $\ell_1$-ball, we can get iteration complexity that only grows logarithmically with the dimension $p$, if for each player we use OFTRL with an entropic regularizer: i.e. $R_{\min}(x) = R_{\max}(x) = \sum_{i=1}^{2p} x_i \log(x_i)$, denotes the negative entropy.

\begin{proposition}\label{prop:sparse-optimization-ell1}
Consider the algorithm that for $t=1, \ldots, T$, sets:
\begin{align}

    \tilde{\rho}_{i, t+1} =~& \tilde{\rho}_{i, t} e^{- 2\frac{\eta}{B}\, \left(- \E_n[v_i u^\top w_{t}] + \frac{\mu}{U}\right) + \frac{\eta}{B}\, \left(- \E_n[v_i u^\top w_{t-1}] + \frac{\mu}{U}\right)} &
    \rho_{t+1} =~& \tilde{\rho}_{t+1}\, \min\left\{1, \frac{B}{\|\tilde{\rho}_{t+1}\|_1}\right\}\\
    \tilde{w}_{i, t+1} =~& w_{i, t} e^{2\, \eta\, \E_n[(y - \rho_t^\top v)\, u_i] - \eta\, \E_n[(y - \rho_{t-1}^\top v)\, u_i]} & w_{t+1} =~& \frac{\tilde{w}_{t+1}}{\|\tilde{w}_{t+1}\|_1}
\end{align}
with $\tilde{\rho}_{i,-1}=\tilde{\rho}_{i,0}=1/e$ and $\tilde{w}_{i,-1}=\tilde{w}_{i,0}=1/(2p)$ and returns $\bar{\rho}=\frac{1}{T} \sum_{t=1}^T \rho_t$. Then for $\eta=\frac{1}{4\|\E_n[vu^\top]\|_{\infty}}$,\footnote{For a matrix $A$, we denote with $\|A\|_{\infty}=\max_{i, j} |A_{ij}|$} after
\begin{equation}
T=16\|\E_n[vu^\top]\|_{\infty} \frac{4B^2 \log(B\vee 1) + (B+1) \log(2p)}{\epsilon}
\end{equation}
iterations, the parameter $\bar{\theta}=\bar{\rho}^{+} - \bar{\rho}^-$ is an $\epsilon$-approximate solution to the minimax problem in Equation~\eqref{eqn:minimax-ell1}.
\end{proposition}

Moreover, every update step requires computation time $O(\min\{n\, p, p^2\})$.\footnote{If $p\geq n$, then at every iteration we can calculate $m^{(j)}=v^{(j)}\cdot w_t$, for each sample $v^{(j)}$; which takes $O(n\cdot p)$ time; and then update each $\tilde{\rho}_{i, t+1}$ based on the quantity $\E_n[v_i u^\top w_{t}]=\frac{1}{n} \sum_j v_i^{(j)} m^{(j)}$. If $p<n$, then we can calculate $\Sigma_n = \E_n[vu^\top ]$ ahead of time and at each period calculate $\E_n[v_i u^\top w_{t}]=(\Sigma w_t)_i$; which would require $O(p^2)$ time.} Using techniques for sparse gradient updates, one could also potentially improve the iteration complexity to not depend linearly on the dimension $p$ (see e.g. \cite{langford2009sparse,Duchi2008l1,Duchi2009forward,mcmahan2011follow}), but we defer such approaches to future work.

\paragraph{$\ell_2$-ball adversary} For the case when $\|\beta\|=\|\beta\|_2$, then we can use $R_{\max}(\beta)=\frac{1}{2}\|\beta\|_2^2$, which leads to an alternative update rule for the maximizing player.
In this case, the update of the maximizing player is essentially optimistic gradient descent, modulo the normalization so as to respect the $\ell_2$-norm constraint.

\begin{proposition}\label{prop:sparse-optimization-ell2}
Consider the algorithm that for $t=1, \ldots, T$, sets:
\begin{align}

    \tilde{\rho}_{i, t+1} =~& \tilde{\rho}_{i, t} e^{- 2\frac{\eta}{B}\, \left(- \E_n[v_i z^\top \beta_{t}] + \frac{\mu}{U}\right) + \frac{\eta}{B}\, \left(- \E_n[v_i z^\top \beta_{t-1}] + \frac{\mu}{U}\right)} &
    \rho_{t+1} =~& \tilde{\rho}_{t+1} \min\left\{1, \frac{B}{\|\tilde{\rho}_{t+1}\|_1}\right\}\\

    \tilde{\beta}_{t+1} =~& \tilde{\beta}_{t+1} + 2\eta \E_n[(y - \rho_{t}^\top v)\, z] - \eta \E_n[(y - \rho_{t-1}^\top v)\, z]  & \beta_{t+1} =~& \tilde{\beta}_{t+1} \min\left\{1, \frac{U}{\|\tilde{\beta}_{t+1}\|_2}\right\}
\end{align}
with $\tilde{\rho}_{i,-1}=\tilde{\rho}_{i,0}=1/e$ and $\tilde{\beta}_{-1}=\tilde{\beta}_0 = 0$. Then for $\eta=\frac{1}{4\|\E_n[zv^\top\|_{2,\infty}}$,\footnote{For a matrix $A$, we denote with $\|A\|_{2,\infty}=\sqrt{\sum_i \max_{j} A_{ij}^2}$} after
\begin{equation}
T=16\|\E_n[zv^\top]\|_{2,\infty} \frac{4B^2 \log(B\vee 1) + B \log(2p) + U^2/2}{\epsilon}.
\end{equation}
iterations, the parameter $\bar{\theta}=\bar{\rho}^{+} - \bar{\rho}^-$ is an $\epsilon$-approximate solution to the minimax problem in Equation~\eqref{eqn:minimax-ell1}.
\end{proposition}
Observe that if $v_j\in [-H, H]$ then the quantity $\|\E_n[zv^\top]\|_{2,\infty}\|$ can be upper bounded by $H\, \sqrt{\E_n[\|z\|_2^2]}$, which under the assumptions of \Cref{cor:sparse-linear-reg-ell1} is at most a constant.

\begin{figure}
    \centering
    \begin{subfigure}[b]{0.3\textwidth}
        \includegraphics[width=\textwidth]{figures/n600_true_v_est_sparse_linear.png}
        \caption{true vs. est. $\theta$ ($n=600$)}
        \label{fig:true_v_est}
    \end{subfigure}
    \begin{subfigure}[b]{0.3\textwidth}
        \includegraphics[width=\textwidth]{figures/true_v_est_sparse_linear.png}
        \caption{true vs. est. $\theta$ ($n=1000$)}
        \label{fig:true_v_est2}
    \end{subfigure}
    ~

    \begin{subfigure}[b]{0.3\textwidth}
        \includegraphics[width=\textwidth]{figures/duals.png}
        \caption{dual variables $w^+ - w^-$}
        \label{fig:duals}
    \end{subfigure}
    \caption{Estimates based on minimax estimator proposed in \Cref{prop:sparse-optimization-ell1}.  The left figure depicts the $p=2000$ estimated coefficients compared to the true coefficients; we also include the coefficients of i) a direct lasso regression to portray the importance of dealing with the endogeneity problem (Lasso), ii) a two-stage lasso regression where we regress each $x_i$ on $z$ and then regress $y$ on $\E[x\mid z]$, all regressions performed with lasso where the first stage regularization was fixed to $0.01$ and the final stage was chosen via cross-validation (2SLasso), iii) the algorithm in \Cref{prop:sparse-optimization-ell1} (SparseIV), iv) a stochastic variant of the algorithm in \Cref{prop:sparse-optimization-ell1} where a mini-batch of $10$ samples is used at each iteration (StochasticSparseIV). The right pictures depicts the coefficients of the dual test function learned by the adversary at equilibrium, which is of the form: $f(z)=\sum_{i=1}^p (w_i^+ - w_i^-) z_i$. The data generating process was: $x, z, u\in \R^p$, $x = z + u$, $y = \ldot{x + u}{\theta}$, $z, u\sim N(0, I_d)$, $\theta=(1, -1, 0, \ldots, 0)$, $p=2000$.}\label{fig:sparse_linear}
\end{figure}

\subsection{Bounds on Ill-Posedness Measure}\label{app:high-dim-linear-ill}

Let $h(x)=\ldot{\theta}{x}$, $h_0(x)=\ldot{\theta_0}{x}$ and $\nu=\theta-\theta_0$. Then observe that we have:
\begin{equation}
    \|T(h-h_0)\|_2^2 = \nu^{\top} \E\left[\E[x\mid z]\E[x\mid z]^\top\right] \nu = \nu^\top V \nu \geq \lambda_{\min}(V) \|\nu\|_2^2
\end{equation}
where we remind that $V:= \E\left[\E[x\mid z]\E[x\mid z]^\top\right]$ and $\lambda_{\min}(V)$ denotes the minimum eigenvalue of $V$. Moreover, if we let $\Sigma=\E\left[xx^\top\right]$ then:
\begin{equation}
    \|h-h_0\|_2^2 = \nu^\top \E\left[x x^\top\right] \nu \leq \lambda_{\max}(\Sigma) \|\nu\|_2^2
\end{equation}
Thus we see that the measure of ill-posedness can be upper bounded as:
\begin{equation}
\tau \leq \sqrt{\frac{\lambda_{\max}(\Sigma)}{\lambda_{\min}(V)}}
\end{equation}
Thus assuming that these eigenvalues are upper and lower bounded correspondingly, then the results of this section extend also to RMSE guarantees for the recovered $\hat{h}$ and not just projected RMSE guarantees, at the cost of an extra multiplicative factor of $\tau$.

Moreover, we note that in both our hard sparsity and $\ell_1$-relaxed estimators we have further constraints on the vector $\nu$ and thus we only require the minimum and maximum eigenvalue to be bounded subject to these constraints. For instance, in the case of hard sparsity, we know that $\nu$ is a $2s$-sparse vector. Thus it suffices to require the minimum eigenvalue of $V$ and the maximum eigenvalue of $\Sigma$ to be bounded only for such $2s$-sparse vectors (i.e. they should hold for all $2s\times 2s$ square sub-matrices of $\Sigma$ and $V$). Similarly, for the $\ell_1$ based estimators we know that the vector $\nu$ falls in a restricted cone, such that most of the $\ell_1$ norm of $\nu$ is concentrated on the $s$ coordinates of the true coefficient $\theta_0$. Thus we solely need the $\lambda_{\min}$ and $\lambda_{\max}$ constraints to be valid only in this restricted cone of vectors.

\section{Application: Shape Constrained Functions}\label{sec:shape}

In this section, we consider the case when $x\in [0, 1]$ and we make shape constraints on $h_0$. We look at both monotonicity/total variation bound constraints and convexity constraints.

\subsection{Monotone functions and functions with small total variation}\label{sec:tviv}
Consider the case when $h_0$ is a function with range in $[0, 1]$ and of bounded total variation, $BV(h_0)\leq 1$.\footnote{Our results easily extend to arbitrary intervals $x\in [a,b]$ and ranges $[-H, H]$, though we restrict to $[0, 1]$ for simplicity of exposition.} We let $\mcH:=BV(1)$ denote the latter class of functions. Moreover, we assume that the operator $T$ satisfies that $Th$ is a monotone non-decreasing (or non-increasing) function of $z$ for any monotone non-decreasing (or non-increasing) function $h$ of $x$. Total variation function classes in linear inverse problems with a known linear operator have also been recently analyzed by \cite{del2019total} and a minimax loss based estimator was also considered, similar in spirit to our general framework.

Observe that any function $h$ with range in $[0, 1]$ and total variation at most $1$ can be written as the difference of two non-decreasing functions $h_+, h_-$ with ranges in $[0, 1]$, i.e. $h = h_+ - h_-$. Thus we note that our assumption on $T$ implies that if $h\in BV(1)$, then $Th = Th_+ - Th_- = f_+ - f_-$, where $f_+$ and $f_-$ are monotone non-decreasing functions in $[0, 1]$. Thus $Th\in BV(1)$ and $T(h-h_0) \in BV(2)$. Thus in order to apply our main theorems, it suffices to take $\mcF=BV(2)$, i.e. the class of functions that can be expressed as the difference of two monotone non-decreasing functions with range in $[0, 2]$. Alternatively, we could also define the norm of a function in the function classes $\mcF$ and $\mcH$ as the total variation, which would enable the regularized estimator to adapt to the total variation of the true hypothesis. For simplicity, we assume a known upper bound.

Furthermore, we note that by standard results in statistical learning theory (see e.g. exercise 18, p.153 of \cite{VanDerVaartWe96} or excercise 3.6.7 of \cite{Gine2015}), that the class of monotone functions with range in $[0, 2]$ have metric entropy of the order of $O(1/\epsilon)$. Thus the same holds for the class $BV(2)$, leading to a critical radius of $\delta_n = O\left( n^{-1/3} \right)$, by invoking \Cref{lem:covering}. Thus by applying our \Cref{thm:reg-main-error}, we get that the corresponding estimators presented in these sections, when $\mcH=BV(1)$ and $\mcF=BV(2)$ (and no norm constraints, which can be emulated by setting $B=U=\infty$), satisfy w.p. $1-\zeta$:
\begin{equation}
    \|T(\hat{h}-h_0)\|_2 = O\left( \frac{1}{n^{1/3}} + \sqrt{\frac{\log(1/\zeta)}{n}}\right)
\end{equation}
The latter rate matches known lower bounds on the achievable RMSE for monotone functions even in the case of exogenous regression problems \cite{chatterjee2015}.

\paragraph{Efficiently solving the optimization problem} We can solve the empirical optimization problem by using piece-wise constant monotone functions (or piece-wise linear), i.e. when running the estimator on $n$ samples, we can describe the function $h$ via a $2n$-dimensional vector $\theta=(\theta^+; \theta^-)$, such that $1\geq \theta_1^+\geq \ldots \geq \theta_n^+\geq 0$ and $1\geq \theta_1^-\geq \ldots \geq \theta_n^-\geq 0$.\footnote{If we want to enforce a monotone non-decreasing $h$, then we can set $\theta^-=0$ and similarly, for a monotone non-increasing algorithm $\theta^+=0$.} Let $\Theta$ describe the set of $\theta$ that satisfy these constraints. Similarly,
we can describe $f$ via a vector $w=(w^+; w^-)$, such that $2\geq w_1^+\geq \ldots \geq w_n^+\geq 0$ and $2\geq w_1^-\geq \ldots \geq w_n^-\geq 0$. Let $W$ describe the set of $w$ that satisfy these constraints.

Then for every sample $i$, if we let $q_x(i)$ be the rank of sample $i$ (i.e. sample $i$ has the $q_x(i)$ highest $x$), when we order all samples based on $x$, we can set $h(x_i)=\theta_{q_x(i)}^+ - \theta_{q_x(i)}^-$. Similarly, if we let $q_z(i)$ be the rank of sample $i$, when we order all samples based on $z$, we can set $f(z_i)=w_{q_z(i)}^+ - w_{q_z(i)}^-$. For simplicity of exposition and w.l.o.g. we will assume that samples are ordered in terms of $x$, i.e. $q_x(i) = i$. Thus we can simplify the optimization problem in \Cref{thm:reg-main-error} as:
\begin{equation}
    \min_{\theta\in \Theta} \max_{w\in W} \sum_i (y_i - (\theta_{i}^+-\theta_{i}^-)) (w_{q_z(i)}^+ - w_{q_z(i)}^-) - \lambda \sum_{i=1}^n (w_i^+ - w_i^-)^2
\end{equation}
where the conclusions of the theorem hold if $\lambda \geq 1$. Since the loss:
\begin{equation}
    \ell(\theta, w) = \sum_i (y_i - (\theta_{i}^+-\theta_{i}^-)) (w_{q_z(i)}^+ - w_{q_z(i)}^-) - \lambda \sum_{i=1}^n (w_i^+ - w_i^-)^2
\end{equation}
is convex in $\theta$ and concave in $w$ and the spaces $\Theta, W$ are convex sets, we can solve this problem by running simultaneous projected gradient descent for $\theta$ and $w$ separately and returning the average solutions, i.e.: for $t=1,\ldots, T$:
\begin{align}
    \theta_t =~& \Pi_{\Theta}(\theta_{t-1} - \eta \nabla_{\theta} \ell(\theta_{t-1}, w_{t-1}))\\
    w_t =~& \Pi_W(w_{t-1} + \eta \nabla_w \ell(\theta_{t-1}, w_{t-1}))
\end{align}
and return $\bar{\theta}=\frac{1}{T} \sum_{t=1}^T \theta_t$. After $O(n/\epsilon^2)$ iterations this would return an $\epsilon$-approximate solution to the minimax problem. Each iteration step would require running a projection on the spaces $\Theta, W$. If we let $\tilde{\theta}\in \R^{2n}$, then we need to find a solution to the problem:
\begin{equation}
    \min_{\theta \in \Theta} \frac{1}{2n} \sum_i (\tilde{\theta}_i^+ -\theta_i^+)^2 + (\tilde{\theta}_i^- - \theta_i^-)^2
\end{equation}
Since the objective and the constraints decompose for the two parts of the vector, this corresponds to running two isotonic regressions for $\theta_i^+$ and $\theta_i^-$ with observations $\tilde{\theta}_i^+$ and $\tilde{\theta}_i^-$. Thus each problem can be solved via the well-known Pool-Adjacent-Violator (PAV) algorithm, which requires $O(n)$ computation time. Similarly, we can deal with the projection of $w$. Thus each iteration of the simultaneous projected gradient descent algorithm requires four calls to the PAV algorithm. If we further want to impose Lipschitzness constraints on our estimates, then we can instead use the Lipschitz-PAV algorithm (see \cite{Yeganova2009LPAV,Kakade2011Isotonic}) to project onto spaces $\Theta$ and $W$ that are augmented with lipschitzness constraints, e.g. $0\leq \theta_i^+ - \theta_j^+ \leq L(x_i - x_j)$ for all $i\leq j$. Albeit the LPAV algorithm requires computation of $O(n^2)$.

\begin{figure}
    \centering
    \begin{subfigure}[b]{0.3\textwidth}
        \includegraphics[width=\textwidth]{figures/iso_direct.png}
        \caption{Isotonic Regression $y\sim x$}
        \label{fig:iso_direct}
    \end{subfigure}
    ~

    \begin{subfigure}[b]{0.3\textwidth}
        \includegraphics[width=\textwidth]{figures/iso_iv.png}
        \caption{Isotonic IV}
        \label{fig:iso_iv}
    \end{subfigure}
    ~

    \begin{subfigure}[b]{0.3\textwidth}
        \includegraphics[width=\textwidth]{figures/iso_lip_iv.png}
        \caption{Lipschitz Isotonic IV}
        \label{fig:iso_lip_iv}
    \end{subfigure}
    \caption{Estimated functions based on our minimax estimator under monotonicity constraints. The first figure depicts a direct isotonic regression that ignores endogeneity. The second figure depics our isotonic IV regression, without any lipschitz constraints and the final figure depicts our isotonic IV regression with Lipschitzness constraints. The data generating process was: $h_0(x)=x^2\, 1\{x>0\}$, $x=.6\, z + .4\, u + \delta$ and $y=h_0(x) + u + \epsilon$ and $z,u\sim N(0, 2)$ and $\epsilon, \delta\sim N(0, .1)$. ($n=1000$)}\label{fig:iso}
\end{figure}

\paragraph{Generality of computational approach} We note that the above approach of solving the endogenous regression problem with shape constraints via our minimax estimator essentially applies to any type of shape constraints and reduces the minimax problem to a standard square loss problem subject to the same shape constraints (assuming that both $\mcH$ and $\mcF$ satisfy the same shape constraints; i.e. that these constraints are invariant to the application of the operator $T$). Thus to solve the minimax problem we simply require an oracle for the square loss problem. In the the setting described in this section we used the PAV and LPAV algorithm as such oracles. In the next section we will be using a quadratic optimization subject to linear constraints solver as our oracle.

\paragraph{Ill-posedness} We note that the recent work of \cite{Chetverikov2017}, shows that when $x, z\in [0,1]$ and the distributions of $x$ and $z$ have full support and lower-bounded density, then for any function $h$, that is $\alpha$-approximately monotone and continuously differentiable, then $\|Th\|_2 \geq \frac{1}{\tau} \|h\|_{2,t}$, where $\|h\|_{2,t}=\int_{x_1}^{x_2} h(x)^2 dx$, for some $0<x_1 < x_2 <1$. The result requires several more regularity conditions on the operator $T$ and the constant $\tau$ depends on constants in these regularity conditions (e.g. the lower bound on the density, the quantities $x_1$ and $1-x_2$, the constant $\alpha$, etc). Thus under these further regularity conditions, we have that for any $h_*$ that is $\alpha$-approximately constant and for $h$ being a monotone function $\|T(h-h_*)\|_2 \geq \frac{1}{\tau} \|h\|_{2,t}$. Thus our bound on $\|T(h-h_*)\|_2$ also implies a bound on $\|h-h_*\|_{2,t}$. This claim, roughly recovers the main estimation rate result of \cite{Chetverikov2017}.

\subsection{Convex functions}\label{sec:convexiv}

In this section we consider the case when $h_0$ is assumed to be a convex function in $[0,1]$, $\Gamma$-Lipschitz and with range in $[0, 1]$. Moreover, we asusme that the linear operator $T$ satisfies that for any convex $\Gamma$-Lipschitz function $h$, $Th$ is also convex and $\Gamma$-Lipschitz. Observe that if $T$ is a symmetric density, i.e. $Th=h\star \rho$ (where $\star$ denotes the convolution operator), for some conditional density function $\rho$, then we have $(Th)''(z) = (h'') \star \rho \geq 0$, since $h''(x)\geq 0$ and $\rho(x)\geq 0$ for all $x$. Thus any such symmetric density satisfies our constraints.

The work of \cite{bronshtein1976varepsilon} shows that the metric entropy this function class, even in the $d$-dimensional hypercube, with respect to the $\ell_{\infty}$ norm, and therefore also with respect to the $\ell_{2,n}$ norm, is of the order of $\epsilon^{-d/2}$ (see also the recent work of \cite{guntuboyina2012covering}). Thus we get that by invoking \Cref{lem:covering}, for $d=1$, we can choose $\delta_n$ in \Cref{thm:reg-main-error} in the order of $O(n^{2/5})$, leading to the corollary that the estimator in \Cref{thm:reg-main-error}, for the case when $\mcH$ is the space of convex, $\Gamma$-Lipscthiz functions with range in $[0,1]$ and $\mcF$ is the space of differences of two convex functions, each $\Gamma$-Lipschitz and with range in $[0, 1$, then w.p. $1-\zeta$:
\begin{equation}
    \|T(\hat{h}-h_0)\|_2 = O\left( \frac{1}{n^{2/5}} + \sqrt{\frac{\log(1/\zeta)}{n}}\right)
\end{equation}

\paragraph{Solving the optimization problem} Moreover, we can address the optimization problem in manner similar to the previous section. We can choose estimators that optimize over piece-wise linear functions and hence can be uniquely determined by their values on the $n$ samples, i.e. we can describe $h$ by a $n$-dimensional vector $\theta$, such that $h(x_i) = \theta_{q_x(i)}$ (where $q_x(i)$ as defined in the previous section). Similarly, we can descirbe $f\in \mcF$ via a $2n$-dimensional vector $w=(w^+; w^-)$, such that $f(z_i)=w_{q_z(i)}^+ - w_{q_z(i)}^-$. Subsequently, we can apply the simultaneous projected gradient descent approach, which reduces the minimax optimization problem to solving the projection problem. Observe that we can describe the constraints that describe the vectors $\theta$ and $w$ as linear constraints. Using the same idea as the one described in Example 13.4 of \cite{wainwright2019high}, we can express the convexity constraint as the existence of a subgradient, i.e. there must exist sub-gradients $u, \mu^+, \mu^- \in \R^{n}$ such that for all $i,j\in [n]$:
\begin{align}
    \theta_j \geq~& \theta_i + \ldot{u_i}{x_{q_x^{-1}(j)} - x_{q_x^{-1}(i)}}\\
    w_j^+ \geq~& w_i^+ + \ldot{\mu_i^+}{z_{q_z^{-1}(j)} - z_{q_z^{-1}(i)}}\\
    w_j^- \geq~& w_i^- + \ldot{\mu_i^-}{z_{q_z^{-1}(j)} - z_{q_z^{-1}(i)}}
\end{align}
This is a set of linear constraints of $\theta, w^+, w^-, u, \mu^+, \mu^-$. Moreover, the lipschitz constraints corresponds to another set of linear constraints, for all $i\in [n]$:
\begin{align}
    -\Gamma (x_{q_x^{-1}(i+1)} - x_{q_x^{-1}(i)}) \leq \theta_{i+1} - \theta_i \leq \Gamma (x_{q_x^{-1}(i+1)} - x_{q_x^{-1}(i)})
\end{align}
and similarly for $w^+, w^-$. Thus projecting onto onto $\Theta$ or $W$, corresponds to a convex quadratic optimization problem with $2n$ variables and $O(n^2)$ linear constraints. Therefore, we can compute such projections in polynomial time at every iteration of the simultaneous projected gradient descent algorithm. In practice, one can achieve substantial speedup by subsampling a set of $s \ll n$ points and restricting the curve to a piece-wise linear function in between these points. This would reduce the number of variables and constraints to $2s$ and $O(s^2)$, correspondingly.

\begin{figure}
    \centering
    \begin{subfigure}[b]{0.3\textwidth}
        \includegraphics[width=\textwidth]{figures/bounded_tv.png}
        \caption{Bounded TV}
        \label{fig:bv}
    \end{subfigure}
    ~

    \begin{subfigure}[b]{0.3\textwidth}
        \includegraphics[width=\textwidth]{figures/bounde_tv_lipschitz.png}
        \caption{Bounded TV and $1$-Lipschitz}
        \label{fig:bv_lip}
    \end{subfigure}
    ~

    \begin{subfigure}[b]{0.3\textwidth}
        \includegraphics[width=\textwidth]{figures/convex_lipschitz.png}
        \caption{Convex and $1$-Lipschitz}
        \label{fig:convex_lip}
    \end{subfigure}
    \caption{Estimated functions based on our minimax estimator for different sets of shape constraints. In the last figure we also depict the direct regression estimate subject to the same constraints, i.e. if we regressed $y$ on $x$, ignoring endogeneity. The data generating process was: $h_0(x)=|x|$ and $x=.5\, z + .5\, u + \delta$ and $y=h_0(x) + u + \epsilon$ and $z,u\sim N(0, 2)$ and $\epsilon, \delta\sim N(0, .1)$. ($n=1000$)}\label{fig:shape}
\end{figure}

\section{Neural Networks}\label{app:neural-networks}

In this section we describe how one can apply the theoretical findings from the previous sections to understand how to train neural networks that solve the conditional moment problem. We will consider the case when our true function $h_0$ can be represented (or well-approximated) by a deep neural network function of $x$, for some given domain specific network architecture, and we will represent it as $h_0(x)=h_{\theta_0}(x)$, where $\theta_0$ are the weights of the neural net. Moreover, we will assume that the linear operator $T$, satisfies that for any set of weights $\theta$, we have that $T h_{\theta}$ belongs to a set of functions that can be represented (or well-approximated) as another deep neural network architecture, and we will denote these functions as $f_w(z)$, where $w$ are the weights of the neural net.
\paragraph{Adversarial GMM Networks (AGMM)} Thus we can apply our general approach presented in \Cref{thm:reg-main-error} and consider the estimator:
\begin{equation}\label{eqn:agmm-nn}
    \hat{\theta} = \argmin_{\theta} \sup_{w}  \E_n[\psi(y_i;h_{\theta}(x_i)) f_{w}(z)] - \lambda \left(\|f_w\|_{\mcF}^2 + \frac{U}{n \delta^2}\sum_i f_w(z_i)^2\right) + \mu \|h_{\theta}\|_{\mcH}^2
\end{equation}
where $\lambda, \mu, U, \delta$ are hyperparameters that need to satisfy the conditions of the theorem. In particular, if we know that the neural nets $h_{\theta}, f_{w}$ output functions in $[0,1]$, then we can choose $U=B=1$, $\lambda = \delta^2$, $\mu = 2\delta^2 (4L^2 + 27)$, where $L$ is a bound on the lipschitzness of the operator $T$ with respect to the two function space norms and $\delta$ is a bound on the critical radius of the function spaces $\mcF_{3}$ and $\hat{\mcG}_{1,L^2}$. Then problem takes the form:
\begin{equation}
    \hat{\theta} = \argmin_{\theta} \sup_{w}  \E_n[\psi(y_i;h_{\theta}(x_i)) f_{w}(z)] - \delta^2 \|f_w\|_{\mcF}^2 - \frac{1}{n}\sum_i f_w(z_i)^2 + c\, \delta^2 \|h_{\theta}\|_{\mcH}^2
\end{equation}
for some constant $c>1$ that depends on the lipschitzness of the operator $T$. Moreover, theoretically we can set the critical radius $\delta$ by invoking \Cref{lem:covering}, and using existing results on the pseudo-dimension of the neural network architecture, for which there exist known bounds \cite{anthony2009neural} that scale with the number of nodes and edges of the neural net. Moreover, one can also use the recent work of \cite{bartlett2017spectrally,Golowich2018}, to provide size independent bounds on the critical radius of these classes, that only depend on spectral properties of the learned weight matrices of the neural nets.

The work of \cite{bennett2019deep} also proposed the use of second moment penalization of the test function, albeit from a different perspective. In particular, their approach stems from a reasoning based on the optimally weighted GMM estimator. In this work we show that second moment penalization arises also when one wants to achieve fast rates of convergence in terms of mean squared error of the learned function. Moreover, the regularization presented in \cite{bennett2019deep} is not a simple second moment penalization, but the second moment of each sample is re-weighted based on the moment evaluated at a preliminary estimate of $\theta$, i.e. $\sum_{i} f_w(z_i)^2 \psi(y_i; h_{\tilde{\theta}}(x_i))^2$. The preliminary estimate of $\tilde{\theta}$ is an extra burden and typically requires sample splitting and first stage estimation. Here we show that such re-weighting is not required if one simply wants fast projected MSE rates. Moreover, this alternative penalty has the property that as the model $h$ becomes very accurate, then $\psi(y_i; h(x_i))\approx 0$ and hence the penalty vanishes as the model becomes accurate. This is a big qualitative difference of the two penalties and it is not clear that the penalty that rescales with the moment enjoys the same theoretical guarantees in terms of projected MSE as the simpler second moment penalty.

In the remainder of the section, we will mostly focus on the practical aspect of training neural networks, such as what would be appropriate architectures for the test function space, based on the intuition developed in the prior theoretical developments of the paper and what would be appropriate optimization algorithms for solving the optimization problem.

\subsection{MMD-GMM: A Neural Network Architecture for Adversarial GMM}
\label{app:neural-networks-mmd}

\paragraph{Maximum Mean Discrepancy GMM Networks (MMD-GMM).} Our results for RKHS function spaces, suggest that one class of test functions are functions that fall in an RKHS. Observe that \Cref{lem:rkhs2} shows that, even when $h$ is an arbitrary function represented by a neural network, as long as $p(x\mid \cdot)$ is a function that belongs to an RKHS $\bbH_K$, with some kernel $K$, then $Th\in \bbH_K$. Thus we can choose test functions in $\bbH_K$.

In many neural network applications, we might have that $p(x\mid \cdot)$ is not in an RKHS (or might have very large RKHS norm), when we use the raw instrument $z$, as $z$ might be very high-dimensional and structured (e.g. an image). However, it might be natural to assume that there is some latent representation $g(z)$ of the instrument $z$, such that: $p(x\mid z) = \rho(x\mid g(z))$ and such that $\rho(x\mid \cdot)$ is in an RKHS.

Thus we will generalize our RKHS approach to augment the adversary with the ability to simultaneously learn the representation $g_w$ (represented as a neural network with weights $w$), and also choose the best function in the RKHS of the implied kernel $K_w(z,z') := K(g_w(z), g_w(z'))$. With this generalization, we are still guaranteeing that $T(h-h_0)\in \mcF$, whenever $p(x\mid \cdot)=\rho(x\mid g(\cdot))$ and $\rho(x\mid \cdot)$ is in $\bbH_K$.

Using the variational characterization of the best function in the RKHS presented in Equation~\eqref{eqn:maximization} we get that the optimization of the adversary can be rephrased as optimizing over test functions of the form $f(z)=\frac{1}{n}\sum_{i=1}^n \beta_i K_w(z_i, z)$, leading to an objective for the adversary of the form:
\begin{align}
    \sup_{\beta, w}~&  \frac{1}{n^2}\sum_{i, j} \left(\psi(y_i;h_{\theta}(x_i)) K_w(z_i, z_j) \beta_j
    - \delta^2 \beta_i\,  K_w(z_i, z_j)\, \beta_j\right) - \frac{1}{n}\sum_i \left(\sum_{j} \frac{\beta_j}{n} K_w(z_i, z_j)\right)^2
\end{align}
which can be written as an average over triplets of samples:
\begin{align}
    \frac{1}{n^3}\sum_{i, j, k} \left( \psi(y_i;h_\theta(x_i)) K_w(z_i, z_j) \beta_j  - \beta_i \left(\delta^2 K_w(z_i, z_j) + K_w(z_i, z_k)K_w(z_k, z_j)  \right) \beta_j \right)
\end{align}
Kernels applied to learned representations have been applied in the context of distribution learning (see e.g. the work on MMD-GANs \cite{li2017mmd,binkowski2018demystifying}) and distribution testing (see the recent work of \cite{liu2020learning}).

\paragraph{Unregularized MMD-GMM.} When we omit the $\ell_{2,n}$ regularization then the optimal solution for $\beta$ can be found in closed form (see \Cref{rkhs-closed-form-max}) and the MMD-GMM simplifies to:
\begin{equation}\label{eqn:kloss-mmd}
    \argmin_{\theta} \sup_{w} \frac{1}{n^2}\sum_{i, j} \psi(y_i;h_{\theta}(x_i)) K_w(z_i, z_j) \psi(y_j; h_{\theta}(x_j)) + c \delta^4 \|h_\theta\|_{\mcH}^2
\end{equation}
This version (without fixed kernel parameters $w$) was also independently analyzed from the perspective of testing by \cite{muandet2020kernel}. However, the $\ell_{2,n}$ penalty is crucial for obtaining fast rates (e.g. rates that adapt to the eigendecay in the case of RKHS spaces). On the other hand, the unregularized MMD-GMM admits a much easier implementation as we do not need to deal with the $n$ parameters $\beta$ and in the case where we use fixed kernel parameters $w$ we don't even need adversarial training.

\paragraph{Kernel Approximation} Moreover, as we saw in the RKHS section, it can be beneficial from a computational perspective to approximate the kernel function by sampling a set of training points (either at random or more cleverly based on either leverage scores or k-means clustering) and restrict the space of functions to be supported only on this subset of the points, i.e. $f(z) = \frac{1}{s}\sum_{i=1}^s \beta_i K(g_w(z_i^*), g_w(z))$, where $z_i^*$ is a set of representative samples and approximating the RKHS norm penalty with $\sum_{i,j\in S} \beta_i K_w(z_i^*, z_j^*) \beta_j$. This has the benefit of only depending on an $|S|$-dimensional vector $\beta$, that the adversary needs to optimize over, as opposed to $n$-dimensional. Moreover, in practice, instead of constraining the centers to be of the form $g_w(z_i^*)$, we could instead consider arbitrary centers $c_i$ in the space of the output of $g_w$ and consider test functions of the form: $f(z) = \frac{1}{s}\sum_{i=1}^s \beta_i K(c_i, g_w(z))$, where $c_i$ are parameters that could also be trained via gradient descent. The latter essentially corresponds to adding what is known as an RBF layer at the end of the adversary neural net. This simplified architecture seems the most appealing from a practical point of view (as it does not require any pre-selection of representative samples $z_i^*$) and is depicted in Figure~\ref{fig:mmdgmm-architecture-simple}.

\begin{figure}
    \centering
        \includegraphics[width=.7\textwidth]{figures/architecture.png}
    \caption{MMD-GMM architecture of adversary's test function.}\label{fig:mmdgmm-architecture}
\end{figure}

\paragraph{Multi-Kernel MMD-GMM.} The case of sparse linear representations portrays that it might be important to test many different classes of functions, each potentially trained on a separate part of the input space, since different instruments might be correlated with different treatments and many of these treatments can be irrelevant.
\begin{equation}
    \sup_{w_1, \ldots, w_m, t\in [m]}  \E_n[\psi(y_i;h(x_i)) f_{w_t}(z_{S_t})] - \delta^2 \|f_{w_t}\|_{\mcF}^2 -  \frac{1}{n}\sum_i f_{w_t}(z_{S_t, i})^2
\end{equation}
where $S_t$ are pre-defined subsets of the instruments and $z_{S_t}$ corresponds to the sub-vector of instruments. Each of these functions $f_{w_t}$ corresponding to a neural net.

One can also combine the above approaches and set $f_{w_t}(z_{S_t}) = \frac{1}{n}\sum_{j} \beta_{tj} K_{w_t}(z_{S_t, j}, z_{S_t})$, i.e. allow for the test function that takes as input the subset of the instruments $S_t$ to be in an RKHS of a learned kernel $w_t$. This leads to taking a supremum over a set of kernels in the MMD-GMM objective, where each kernel calculates similarity based on a subset of the input instruments, i.e.:
\begin{align}
    \sup_{\beta, w, t}~&  \frac{1}{n^2}\sum_{i, j} \left(\psi(y_i;h_{\theta}(x_i)) K_{w}^t(z_{i}, z_{j}) \beta_{tj}
    - \delta^2 \beta_{ti}\,  K_w^t(z_i, z_j)\, \beta_{tj}\right) - \frac{1}{n}\sum_i \left(\sum_{j} \frac{\beta_{tj}}{n} K_w^t(z_i, z_j)\right)^2
\end{align}
where $K_w^t(z_i, z_j)$ is shorthand notation for $K_{w_t}(z_{S_t, i}, z_{S_t, j})$. The adverary's objective can also be written as choosing a distribution $p_t$ over the $t$ kernels, leading to an adversary objective of:
\begin{align}
    \frac{1}{n^3}\sum_{i, j, k} \left( \psi(y_i;h_\theta(x_i)) \sum_{t} p_t K_w^t(z_i, z_j) \beta_j  - \beta_i\, \beta_j \sum_{t} p_t (\delta^2 K_w^t(z_i, z_j) + K_w^t(z_i, z_k)K_w^t(z_k, z_j))  \right)
\end{align}
We can again reduce the complexity of the optimization problem by restricting to a subset of samples to represent the test functions.

This combined method targets settings where different instruments are correlated with different latent ``treatment factors'', treatment factors are high-dimensional but only a small subset of them having a large and additively separable effect on the outcome and the relationship between the treatment factor and the instrument is non-linear. Thus it tackles several sources of high-dimensionality in the instrumental variable regression problem.

\begin{figure}
    \centering
        \includegraphics[width=.7\textwidth]{figures/architecture_simple.png}
    \caption{Simplified MMD-GMM architecture of adversary's test function with kernel final activation layer.}\label{fig:mmdgmm-architecture-simple}
\end{figure}

\subsection{Adversarial Training: Simultaneous Optimistic First-Order Stochastic Optimization}
\label{app:neural-networks-opt}

The optimization problem that we are facing is similar to the optimization problem that is encountered in training Generative Adversarial Networks, i.e. we need
to solve a non-convex, non-concave zero-sum game, where the strategy of each of the two players are the parameters of a neural net. This is obviously a computationally intractable problem from a worst-case perspective. However, typical instances are far from worst-case and there has been a surge of recent work proposing iterative optimization algorithms inspired by the convex-concave zero-sum game theory (see, e.g. the Optimistic Adam algorithm of \cite{Daskalakis2017}). For instance, one can expect that in practice most early layers of a neural net will change very slowly or will not have a face transition in their non-linearities. In that case, the main parameters that matter are the parameters of the final layers of the two neural nets. However, the zero-sum game is convex-concave in these parameters. Hence, assuming that the features constructed in the final layer of the two neural nets, change slowly, then one should expect convex-concave zero-sum game optimization theory to apply. Such arguments have been recently exploited in the case of square loss minimization with deep over-parameterized neural networks (see e.g. \cite{AllenZhu2018,du2018gradient,Soltanolkotabi2019}). It is highly plausible and an interesting question for future research, whether such guarantees extend to the minimax problem that we are facing here. For instance, recent work of \cite{Lei2019}, provides an instance of a minimax objective, related to training Wasserstein GANs, where stochastic iterative optimization of neural nets provably converges to an optimal solution.

In our implementation and experiments we used the optimistic Adam algorithm as was also proposed in \cite{bennett2019deep}. Other algorithms that could prove useful for our problem are the extra-gradient or stochastic extra-gradient algorithm (see e.g. \cite{Hsieh2019,Mishchenko2019}).

\begin{figure}
    \centering
    \begin{subfigure}[b]{0.3\textwidth}
        \centering
        \captionsetup{justification=centering}
        \includegraphics[width=\textwidth]{figures/neural_net_test_n1.png}
        \caption{AGMM\\ ($p=1$, $n=4000$)}
    \end{subfigure}
    \begin{subfigure}[b]{0.3\textwidth}
        \centering
        \captionsetup{justification=centering}
        \includegraphics[width=\textwidth]{figures/kernel_gmm_fixed_kernel_n1.png}
        \caption{MMD-GMM\\ ($p=1$, $n=4000$)}
    \end{subfigure}
    \begin{subfigure}[b]{0.3\textwidth}
        \centering
        \captionsetup{justification=centering}
        \includegraphics[width=\textwidth]{figures/kernel_gmm_trained_kernel_n1.png}
        \caption{Learned Kernel MMD-GMM\\ ($p=1$, $n=4000$)}
    \end{subfigure}
    \begin{subfigure}[b]{0.3\textwidth}
        \centering
        \captionsetup{justification=centering}
        \includegraphics[width=\textwidth]{figures/neural_net_test.png}
        \caption{AGMM\\ ($p=50$, $n=4000$)}
    \end{subfigure}
    \begin{subfigure}[b]{0.3\textwidth}
        \centering
        \captionsetup{justification=centering}
        \includegraphics[width=\textwidth]{figures/kernel_gmm_fixed_kernel.png}
        \caption{MMD-GMM\\ ($p=50$, $n=4000$)}
    \end{subfigure}
    \begin{subfigure}[b]{0.3\textwidth}
        \centering
        \captionsetup{justification=centering}
        \includegraphics[width=\textwidth]{figures/kernel_gmm_trained_kernel.png}
        \caption{Learned Kernel MMD-GMM\\ ($p=50$, $n=4000$)}
    \end{subfigure}
    \caption{Estimated function based on our minimax estimator with neural networks as a function of the relevant treatment. The $h_{\theta}$ function was a two layer neural net with $100$ hidden units. In the first figure an two-layer neural net was used as a test function $f_w$. In the second and third, we used the MMD-GMM test functions with a low rank approximation. In the second we used test functions of the form: $f_{\beta}(z) = \sum_{i=1}^s \beta_i K_{\gamma}(c_i, z)$, with $c_i$ a fixed grid of test points in $[-3,3]^p$ and $K$ is the rbf kernel with parameter $\gamma=.2$, i.e. $K(z,z')=\exp(-\gamma \|z-z'\|_2^2)$. In the third we learned the kernel, i.e. we used test functions of the form: $f_{w,\beta}(z)=\sum_{i=1}^s \beta_i K_{\gamma}(c_i, g_w(z))$ and $g_w(z)=\text{relu}(Az+b)$ (all the parameters $A, b, \beta, c_i, \gamma$ where trained). The networks were trained via the simultaneous Optimistic Adam algorithm. The data generating process was: $h_0(x)=|x[0]|$ and $x=.6\, z + .4\, u + \delta$ and $y=h_0(x) + u + \epsilon$ and $z\sim N(0, 2 I_p)$, $u\sim N(0, 2)$ and $\epsilon, \delta\sim N(0, .1)$.}\label{fig:deepgmm}
\end{figure}

\begin{figure}
    \centering
    \begin{subfigure}[b]{0.3\textwidth}
        \centering
        \captionsetup{justification=centering}
        \includegraphics[width=\textwidth]{figures/weak_inst.png}
        \caption{Weak Instruments\\ ($p=2$, $n=4000$)}
    \end{subfigure}
    \caption{Estimated function based on our minimax estimator with neural networks as a function of the relevant treatment. The setup is the same as in Figure~\ref{fig:deepgmm}, but we now made the instrument very weak. The data generating process was: $h_0(x)=|x[0]|\, 1\{x[0]>0\}$ and $x=.05\, z + .95\, u + \delta$ and $y=h_0(x) + u + \epsilon$ and $z\sim N(0, 2 I_p)$, $u\sim N(0, 2)$ and $\epsilon, \delta\sim N(0, .1)$.}\label{fig:deepgmm-weak}
\end{figure}

\section{Random Forests via a Reduction Approach}\label{app:ensemble}

In this section we deal with the problem of training random forests that solve the non-parametric IV problem. In particular, we aim to develop a learning procedure that learns a hypothesis $h$ that solves the Conditional Moment~\eqref{eqn:cond-moments}, that is represented as an ensemble of regression trees. Prior work on random forests for causal inference problems has primarily focused on learning forests that capture the heterogeneity of the treatment effect of a treatment, but did not account for non-linear relationships between the treatment and the outcome variable. We will provide a theoretical foundation of the proposed method by taking a reductions approach to the minimax problem defined by our estimator.

For simplicity, throughout this section we will assume that the hypothesis spaces $\mcH$ and $\mcF$ are bounded and have bound critical radius and will make no further norm constraints. Thus the estimator proposed in \Cref{thm:reg-main-error}\footnote{By setting $\lambda=\delta^2/U$, $\mu=2\lambda \left(4L^2 + 27U/B\right)$ using an $\ell_{\infty}$ norm in both function spaces and taking $U, B\to \infty$. Observe that we can also take $L=1$, since $\|Th\|_{\infty} \leq \|h\|_{\infty}$ for any $T$.} takes the simple form of:
\begin{equation}
    \hat{h} = \argmin_{h\in \mcH} \sup_{f\in \mcF}  \E_n[\psi(y_i;h(x_i)) f(z_i)] - \E_n[f(z_i)^2]
\end{equation}

Since the statistical properties of random forests is an active area of investigation, we will solely focus on the optimization problem and leave the statistical properties (e.g. bounding the critical radius or bias of Random Forest methods) to future work. Our goal is to reduce the aforementioned optimization problem to classification and regression oracles over arbitrary hypothesis spaces. Subsequently in practice we can use random forests as oracles.

\paragraph{Reducing the Optimization to Regression and Classification Oracles} To achieve this reduction we will make the assumption that the space $\mcF$ defines a convex image set on the samples, i.e. the set $A=\{(f(z_1),\ldots, f(z_n)): f\in \mcF\}$ is a convex set. This can potentially be violated for tree based methods, but in practice will be alleviated when training a forest with a large set of trees.

We will show that we can reduce the problem to a regression oracle over the function space $\mcF$ and a classification oracle over the function space ${\cal B}$. We will assume that we have a regression oracle that solves the square loss problem over $\mcF$: for any set of labels and features $z_{1:n}, u_{1:n}$ it returns
\begin{equation}
    \text{Oracle}_{\mcF}(z_{1:n}, u_{1:n}) = \argmin_{f\in \mcF} \frac{1}{n}\sum_{i=1}^n \left(u_i -  f(z_i)\right)^2
\end{equation}
Moreover, we will assume that we have a classification oracle that solves the weighted binary classification problem over ${\cal B}$: for any set of sample weights $w_{1:n}$, binary labels $v_{1:n}$ in $\{0, 1\}$ and features $x_{1:n}$:
\begin{equation}
    \text{Oracle}_{\mcH}(x_{1:n}, v_{1:n}, w_{1:n}) = \argmax_{h\in {\cal \mcH}} \frac{1}{n}\sum_{i=1}^n w_i\, \Pr_{z_i\sim \text{Bernoulli}\left(\frac{1+h(x_i)}{2}\right)}\left[v_i = z_i\right]
\end{equation}
Observe that the objective in the equation above is equivalent to a classification accuracy objective, assuming that $h$ outputs values in $[-1, 1]$ and it corresponds to an expected accuracy objective if one interprets $(h(x)+1)/2$ as the probability of label $1$ conditional on $x$. Having access to these oracles we can then show the following computational result:
\begin{reptheorem}{thm:reduction-algorithm}
Consider the algorithm where for $t=1, \ldots, T$: let
\begin{align}
    u_i^t =~& \frac{1}{2} \left( y_i - \frac{1}{t-1} \sum_{\tau=1}^{t-1} h_{\tau}(x_i)\right), &
    f_t =~& \text{Oracle}_{\mcF}\left(z_{1:n}, u_{1:n}^t\right)\\
    v_i^t =~& 1\{f_t(z_i) > 0\},
    w_i^t = |f_t(z_i)| &
    h_t =~& \text{Oracle}_{\mcH}\left(x_{1:n}, v_{1:n}^t, w_{1:n}^t\right)
\end{align}
Then the ensemble hypothesis: $\bar{h}=\frac{1}{T} \sum_{t=1}^T h_t$, is a $\frac{8\, (\log(T)+1)}{T}$-approximate solution to the minimax problem in Equation~\eqref{eqn:minimax-no-norms}.
\end{reptheorem}

In practice, we will consider a random forest regression method as the oracle over $\mcF$ and a binary decision tree classification method as the oracle for $\mcH$.

Moreover, we observe that if the hypothesis space $\mcH$ can be expressed as linear span of base hypothesis, i.e. $\mcH=\{\sum_i w_i b_i: b_i \in B\}$, then observe that because the best-response problem of the learner is linear in the output of the hypothesis, it suffices to optimize only over the space of base hypothesis. Then the algorithm will return a linear span, supported on $T$ base hypothesis that solves the minimax problem over the whole linear span. This improvement can also lead to statistical rate improvements. For instance, if the base hypothesis $B$ is a VC class with VC dimension $d$ (e.g. a binary decision tree with small depth, see e.g. \citep{Mansour2000}), then the algorithm returns a convex combination of $T$ base hypothesis, which has VC dimension at most $d\, T$ \citep{shalev2014understanding}. Thus the entropy integral of $\mcH$ is of the order of $\sqrt{\frac{T\, d \log(n)}{n}}$. If we further have that the entropy integral of $\mcF$ is at most $\kappa(\mcF)$, then we get a final rate of the order of:
\begin{equation}
    \sqrt{\frac{T\, d \log(n)}{n}} + \kappa(\mcF) + \frac{\log(T)}{T}
\end{equation}
Setting, $T=O(n^{1/4})$, one can achieve rates of the order of $n^{-1/4} + \kappa(\mcF)$.

In practice, we will leverage the above observation and train a single binary classification tree at each period of the algorithm, as our $\text{Oracle}_{\mcH}$. In the end the final prediction will be the prediction of the random forest represented by the ensemble of the $T$ trees trained at each period. We refer to this algorithm as Random Forest IV (RFIV).

\newpage

\section{Experimental Analysis}
\label{app:funcform}

We consider the following data generating processes: for $n_x=1$ and $n_z\geq 1$

\begin{align}
    y =~& h_0(x[0]) + e + \delta, & \delta \sim N(0, .1)\\
    x =~& \gamma\, z[0] + (1-\gamma)\, e + \gamma, &
    z \sim N(0, 2\, I_{n_z}), e\sim N(0, 2), \gamma\sim N(0, .1)
\end{align}
While, when $n_x=n_z>1$, then we consider the following modified treatment equation:
\begin{align}
x =~& \gamma\, z + (1-\gamma)\, e + \gamma, &
\end{align}
We consider several ranges of the number of samples $n$, number of treatments $n_x$, number of instruments $n_z$ and instrument strength $\gamma$ and the following functional forms for $h_0$:
\begin{enumerate}
    \item abs: $h_0(x) = |x|$
    \item 2dpoly: $h_0(x)= -1.5\, x + .9\, x^2$
    \item sigmoid: $h_0(x) = \frac{2}{1 + e^{-2 x}}$
    \item sin: $h_0(x)=\sin(x)$
    \item frequentsin: $h_0(x)=\sin(3\,x)$
    \item abssqrt: $h_0(x)=\sqrt{|x|}$
    \item step: $h_0(x)=1\{x < 0\} + 2.5\, 1\{x\geq 0\}$
    \item 3dpoly: $h_0(x)=-1.5\, x + .9\, x^2 + x^3$
    \item linear: $h_0(x)=x$
    \item randpw: piece wise linear function drawn at random
    \item abspos: $h_0(x)= x\, 1\{x\geq 0\}$
    \item sqrpos: $h_0(x)= x^2\, 1\{x\geq 0\}$
    \item band: $h_0(x)=1\{-.75\leq x \leq .75\}$
    \item invband: $h_0(x) = 1-1\{-.75\leq x \leq .75\}$
    \item steplinear: $h_0(x)=2\, 1\{x\geq 0\} - x$
    \item pwlinear: $h_0(x)=(x+1)\,1\{x\leq -1\} + (x-1)\, 1\{x>=1\}$
\end{enumerate}

We consider as classic benchmarks 2SLS with a polynomial features of degree $3$ (2SLS) and a regularized version of 2SLS where ElasticNetCV is used in both stages (Reg2SLS). We have implemented several of the algorithms described in the paper:
\begin{enumerate}
    \item NystromRKHS: The method described in \Cref{app:rkhs}, with the Nystrom approximation described in \Cref{app:rkhs-optimization}. We used $100$ Nystrom samples for the approximation.
    \item ConvexIV: The variant of the method described in \Cref{sec:convexiv} with both lipscthiz and convexity constraints (lipschitz bound of $L=2$).
    \item TVIV: The variant of the method described in \Cref{sec:tviv} without a lipschitz constraint and only total variation constraint.
    \item LipTVIV: The variant of the method described in \Cref{sec:tviv} with lipscthiz constraint and total variation constraint (lipscthiz bound of $L=2$)
    \item RFIV: The method described in \Cref{app:ensemble}, where a Random Forest Regressor is used as an oracle for the adversary (with $40$ trees, max depth $2$, bootstrap sub-sampling enabled, and minimum leaf size of $40$) and Random Forest Classifier (with $5$ trees, max depth $2$, minimum leaf size of $40$ and bootstrap subsampling disabled) was used as an oracle for the learner. The optimization was run for $T=200$ iterations.
    \item SpLin: The method described in \Cref{app:high-dim-linear-ell1} with the specific optimization method described in \Cref{prop:sparse-optimization-ell1}.
    \item StSpLin: A stochastic gradient descent variant of SpLin, where a mini-batch of $100$ samples is used at every step to calculate the co-variance matrices.
    \item AGMM: The method described in Equation~\eqref{eqn:agmm-nn}. A two-layer neural net with $100$ hidden units at each layer and leaky ReLU units was used for both the learner and the adversary architecture. Optimization was done via the Optimistic Adam.
    \item KLayerFixed: The variant of the method described in \Cref{app:neural-networks-mmd}, where an RBF layer is attached at the end of the adversary's architecture with fixed centers, i.e. testing functions of the form: $f(z) = \sum_{j=1}^{n_{\text{centers}}} K(c_j, g_w(z)) \beta_j$, with $n_{\text{centers}}=100$. The centers $c_j$ are placed in a $100$ dimensional feature space and the function $g_w$ is a two-layer neural net with $100$ hidden units in each layer.
    \item KLayerTrained: The same as KLayerFixed, but the centers of the RBF layer are trained.
    \item CentroidMMD: The version of the MMMD-GMM in \Cref{app:neural-networks-mmd}, where we select a subset of the data points to use as centers in the Kernel approximation, i.e. testing functions of the form: $f(z) = \sum_{j=1}^{n_{\text{centers}}} K(g_w(z_j^*, g_w(z)) \beta_j$. $z_j^*$ are chosen as the centroids of a KMeans clustering and $n_{\text{centers}}=100$. $g_w$ is the same architecture as in KLayerFixed.
    \item KLossMMD: The method described in Equation~\eqref{eqn:kloss-mmd}, where no $\ell_{2,n}$ penalty is imposed on the adversary test function. $g_w$ is the same architecture as in KLayerFixed.
\end{enumerate}

In addition to these regimes, we consider high-dimensional experiments with images, following the scenarios proposed in \cite{bennett2019deep} where either the instrument $z$ or treatment $x$ or both are images from the MNIST dataset consisting of grayscale images of $28\times28$ pixels. We compare the performance of our approaches to that of \cite{bennett2019deep}, using their code.
A full description of the DGP is given in \Cref{app:mnist}.

\paragraph{Results.} The main findings are: i) for small number of treatments, the RKHS method with a Nystrom approximation (NystromRKHS), outperforms all methods (Figure~\ref{fig:table1}) with only exception being functions that are highly non-smooth or non-continuous, in which case the methods that are based on shape constraints (ConvexIV, TVIV, LipTVIV) are better, ii) for moderate number of instruments and treatments, Random Forest IV (RFIV) significantly outperforms most methods, with second best being neural networks (AGMM, KLayerTrained) (Figure~\ref{fig:table2}), iii) the estimator for sparse linear hypotheses can handle an ultra-high dimensional regime (Figure~\ref{fig:table3}), iv) neural network methods (AGMM, KLayerTrained) outperform the state of the art in prior work \citep{bennett2019deep} for tasks that involve images (Figure~\ref{fig:mnist_dgps}). The figures below present the average MSE across $100$ experiments ($10$ experiments for Figure~\ref{fig:mnist_dgps}) and two times the standard error of the average MSE.

\begin{figure}[H]
    \centering
    \scriptsize
    \begin{tabular}{c||c|c|c|c|c|c|c}
    & NystromRKHS & 2SLS & Reg2SLS & ConvexIV & TVIV & LipTVIV & RFIV\\
    \hline
    abs & {\bf 0.045 $\pm$ 0.010 } & 0.100 $\pm$ 0.035 & 1.733 $\pm$ 2.981 & 0.054 $\pm$ 0.005 & 0.089 $\pm$ 0.005 & 0.047 $\pm$ 0.004 & 0.084 $\pm$ 0.007 \\
2dpoly & 0.121 $\pm$ 0.014 & {\bf 0.036 $\pm$ 0.022 } & 9.068 $\pm$ 16.071 & 0.060 $\pm$ 0.007 & 0.090 $\pm$ 0.009 & 0.069 $\pm$ 0.009 & 0.379 $\pm$ 0.022 \\
sigmoid & {\bf 0.016 $\pm$ 0.003 } & 0.071 $\pm$ 0.037 & 0.429 $\pm$ 0.244 & 0.029 $\pm$ 0.005 & 0.067 $\pm$ 0.004 & 0.034 $\pm$ 0.003 & 0.044 $\pm$ 0.006 \\
sin & {\bf 0.023 $\pm$ 0.003 } & 0.090 $\pm$ 0.042 & 0.801 $\pm$ 0.420 & 0.055 $\pm$ 0.005 & 0.074 $\pm$ 0.004 & 0.036 $\pm$ 0.003 & 0.057 $\pm$ 0.007 \\
frequentsin & 0.129 $\pm$ 0.005 & 0.193 $\pm$ 0.040 & 0.145 $\pm$ 0.017 & 0.143 $\pm$ 0.008 & 0.115 $\pm$ 0.005 & {\bf 0.106 $\pm$ 0.005 } & 0.126 $\pm$ 0.010 \\
abssqrt & {\bf 0.033 $\pm$ 0.004 } & 0.099 $\pm$ 0.039 & 0.117 $\pm$ 0.046 & 0.045 $\pm$ 0.007 & 0.096 $\pm$ 0.006 & 0.047 $\pm$ 0.004 & 0.064 $\pm$ 0.008 \\
step & {\bf 0.035 $\pm$ 0.003 } & 0.103 $\pm$ 0.043 & 0.497 $\pm$ 0.276 & 0.054 $\pm$ 0.005 & 0.073 $\pm$ 0.004 & 0.044 $\pm$ 0.003 & 0.056 $\pm$ 0.007 \\
3dpoly & 0.220 $\pm$ 0.037 & {\bf 0.004 $\pm$ 0.003 } & 0.066 $\pm$ 0.014 & 0.396 $\pm$ 0.051 & 0.138 $\pm$ 0.028 & 0.190 $\pm$ 0.036 & 0.687 $\pm$ 0.069 \\
linear & 0.019 $\pm$ 0.003 & 0.038 $\pm$ 0.021 & 0.355 $\pm$ 0.189 & {\bf 0.017 $\pm$ 0.005 } & 0.042 $\pm$ 0.002 & 0.027 $\pm$ 0.002 & 0.048 $\pm$ 0.005 \\
randpw & 0.067 $\pm$ 0.012 & 0.092 $\pm$ 0.024 & 3.810 $\pm$ 5.878 & 0.162 $\pm$ 0.032 & 0.073 $\pm$ 0.009 & {\bf 0.046 $\pm$ 0.006 } & 0.121 $\pm$ 0.015 \\
abspos & {\bf 0.022 $\pm$ 0.003 } & 0.060 $\pm$ 0.027 & 0.299 $\pm$ 0.157 & 0.022 $\pm$ 0.004 & 0.062 $\pm$ 0.004 & 0.033 $\pm$ 0.003 & 0.055 $\pm$ 0.006 \\
sqrpos & 0.064 $\pm$ 0.013 & {\bf 0.026 $\pm$ 0.015 } & 0.490 $\pm$ 0.494 & 0.030 $\pm$ 0.006 & 0.034 $\pm$ 0.003 & 0.033 $\pm$ 0.005 & 0.181 $\pm$ 0.013 \\
band & {\bf 0.059 $\pm$ 0.003 } & 0.125 $\pm$ 0.051 & 0.085 $\pm$ 0.017 & 0.086 $\pm$ 0.008 & 0.102 $\pm$ 0.006 & 0.059 $\pm$ 0.004 & 0.071 $\pm$ 0.008 \\
invband & {\bf 0.056 $\pm$ 0.003 } & 0.130 $\pm$ 0.041 & 0.138 $\pm$ 0.051 & 0.075 $\pm$ 0.008 & 0.102 $\pm$ 0.006 & 0.059 $\pm$ 0.004 & 0.073 $\pm$ 0.008 \\
steplinear & 0.141 $\pm$ 0.009 & 0.231 $\pm$ 0.085 & 0.203 $\pm$ 0.063 & 0.138 $\pm$ 0.008 & 0.156 $\pm$ 0.009 & {\bf 0.100 $\pm$ 0.006 } & 0.141 $\pm$ 0.011 \\
pwlinear & {\bf 0.032 $\pm$ 0.004 } & 0.051 $\pm$ 0.024 & 0.058 $\pm$ 0.025 & 0.037 $\pm$ 0.006 & 0.061 $\pm$ 0.003 & 0.035 $\pm$ 0.003 & 0.068 $\pm$ 0.006
    \end{tabular}
    \caption{$n=300$, $n_z=1$, $n_x=1$, $\gamma=.6$}
\end{figure}

\begin{figure}[H]
    \centering
    \scriptsize
    \begin{tabular}{c||c|c|c|c|c|c|c}
    & NystromRKHS & 2SLS & Reg2SLS & ConvexIV & TVIV & LipTVIV & RFIV\\
    \hline
    abs & {\bf 0.010 $\pm$ 0.001 } & 0.025 $\pm$ 0.001 & 0.025 $\pm$ 0.002 & 0.031 $\pm$ 0.001 & 0.031 $\pm$ 0.001 & 0.021 $\pm$ 0.001 & 0.026 $\pm$ 0.002 \\
2dpoly & 0.022 $\pm$ 0.005 & {\bf 0.002 $\pm$ 0.000 } & 0.043 $\pm$ 0.039 & 0.052 $\pm$ 0.004 & 0.034 $\pm$ 0.004 & 0.037 $\pm$ 0.004 & 0.286 $\pm$ 0.013 \\
sigmoid & {\bf 0.005 $\pm$ 0.001 } & 0.007 $\pm$ 0.001 & 0.021 $\pm$ 0.017 & 0.011 $\pm$ 0.000 & 0.018 $\pm$ 0.001 & 0.008 $\pm$ 0.001 & 0.015 $\pm$ 0.001 \\
sin & {\bf 0.005 $\pm$ 0.001 } & 0.013 $\pm$ 0.002 & 0.033 $\pm$ 0.025 & 0.035 $\pm$ 0.001 & 0.020 $\pm$ 0.001 & 0.009 $\pm$ 0.001 & 0.017 $\pm$ 0.001 \\
frequentsin & 0.118 $\pm$ 0.001 & 0.117 $\pm$ 0.001 & 0.115 $\pm$ 0.001 & 0.116 $\pm$ 0.001 & 0.089 $\pm$ 0.002 & 0.105 $\pm$ 0.002 & {\bf 0.087 $\pm$ 0.004 } \\
abssqrt & {\bf 0.011 $\pm$ 0.001 } & 0.018 $\pm$ 0.001 & 0.018 $\pm$ 0.001 & 0.020 $\pm$ 0.001 & 0.028 $\pm$ 0.001 & 0.016 $\pm$ 0.001 & 0.022 $\pm$ 0.002 \\
step & 0.022 $\pm$ 0.001 & 0.029 $\pm$ 0.001 & 0.043 $\pm$ 0.017 & 0.034 $\pm$ 0.001 & 0.026 $\pm$ 0.001 & {\bf 0.020 $\pm$ 0.001 } & 0.026 $\pm$ 0.002 \\
3dpoly & 0.028 $\pm$ 0.012 & {\bf 0.000 $\pm$ 0.000 } & 0.010 $\pm$ 0.003 & 0.325 $\pm$ 0.026 & 0.086 $\pm$ 0.019 & 0.121 $\pm$ 0.020 & 0.375 $\pm$ 0.036 \\
linear & 0.004 $\pm$ 0.001 & {\bf 0.002 $\pm$ 0.000 } & 0.022 $\pm$ 0.022 & 0.002 $\pm$ 0.000 & 0.013 $\pm$ 0.001 & 0.007 $\pm$ 0.001 & 0.012 $\pm$ 0.001 \\
randpw & 0.031 $\pm$ 0.006 & 0.057 $\pm$ 0.010 & 0.131 $\pm$ 0.111 & 0.150 $\pm$ 0.032 & 0.032 $\pm$ 0.004 & {\bf 0.029 $\pm$ 0.004 } & 0.054 $\pm$ 0.010 \\
abspos & 0.006 $\pm$ 0.001 & 0.007 $\pm$ 0.001 & 0.015 $\pm$ 0.009 & {\bf 0.005 $\pm$ 0.000 } & 0.016 $\pm$ 0.001 & 0.008 $\pm$ 0.001 & 0.016 $\pm$ 0.001 \\
sqrpos & 0.011 $\pm$ 0.003 & {\bf 0.004 $\pm$ 0.000 } & 0.010 $\pm$ 0.006 & 0.011 $\pm$ 0.002 & 0.011 $\pm$ 0.001 & 0.012 $\pm$ 0.002 & 0.091 $\pm$ 0.007 \\
band & 0.031 $\pm$ 0.001 & 0.046 $\pm$ 0.001 & 0.046 $\pm$ 0.001 & 0.059 $\pm$ 0.001 & 0.039 $\pm$ 0.002 & {\bf 0.031 $\pm$ 0.002 } & 0.032 $\pm$ 0.002 \\
invband & {\bf 0.031 $\pm$ 0.001 } & 0.046 $\pm$ 0.001 & 0.046 $\pm$ 0.001 & 0.049 $\pm$ 0.001 & 0.039 $\pm$ 0.002 & 0.031 $\pm$ 0.001 & 0.032 $\pm$ 0.002 \\
steplinear & 0.066 $\pm$ 0.002 & 0.085 $\pm$ 0.003 & 0.089 $\pm$ 0.005 & 0.104 $\pm$ 0.001 & 0.074 $\pm$ 0.002 & {\bf 0.064 $\pm$ 0.002 } & 0.066 $\pm$ 0.003 \\
pwlinear & {\bf 0.007 $\pm$ 0.001 } & 0.009 $\pm$ 0.000 & 0.012 $\pm$ 0.001 & 0.017 $\pm$ 0.001 & 0.016 $\pm$ 0.001 & 0.009 $\pm$ 0.001 & 0.016 $\pm$ 0.001
    \end{tabular}
    \caption{$n=2000$, $n_z=1$, $n_x=1$, $\gamma=.6$}
\end{figure}

\begin{figure}[H]
    \centering
    \scriptsize
    \begin{tabular}{c||c|c|c|c|c|c|c}
    & NystromRKHS & 2SLS & Reg2SLS & ConvexIV & TVIV & LipTVIV & RFIV\\
    \hline
    abs & {\bf 0.008 $\pm$ 0.001 } & 0.027 $\pm$ 0.001 & 0.027 $\pm$ 0.001 & 0.024 $\pm$ 0.000 & 0.016 $\pm$ 0.001 & 0.012 $\pm$ 0.001 & 0.017 $\pm$ 0.001 \\
2dpoly & 0.009 $\pm$ 0.002 & {\bf 0.001 $\pm$ 0.000 } & 0.016 $\pm$ 0.007 & 0.036 $\pm$ 0.003 & 0.018 $\pm$ 0.002 & 0.022 $\pm$ 0.003 & 0.151 $\pm$ 0.010 \\
sigmoid & {\bf 0.004 $\pm$ 0.000 } & 0.007 $\pm$ 0.000 & 0.017 $\pm$ 0.005 & 0.013 $\pm$ 0.000 & 0.011 $\pm$ 0.001 & 0.007 $\pm$ 0.000 & 0.012 $\pm$ 0.001 \\
sin & {\bf 0.003 $\pm$ 0.000 } & 0.023 $\pm$ 0.002 & 0.033 $\pm$ 0.006 & 0.055 $\pm$ 0.001 & 0.013 $\pm$ 0.001 & 0.009 $\pm$ 0.001 & 0.014 $\pm$ 0.001 \\
frequentsin & 0.114 $\pm$ 0.001 & 0.114 $\pm$ 0.001 & 0.113 $\pm$ 0.001 & 0.114 $\pm$ 0.001 & 0.048 $\pm$ 0.001 & 0.051 $\pm$ 0.001 & {\bf 0.024 $\pm$ 0.001 } \\
abssqrt & {\bf 0.008 $\pm$ 0.000 } & 0.017 $\pm$ 0.001 & 0.017 $\pm$ 0.001 & 0.017 $\pm$ 0.000 & 0.015 $\pm$ 0.001 & 0.011 $\pm$ 0.001 & 0.015 $\pm$ 0.001 \\
step & 0.021 $\pm$ 0.000 & 0.031 $\pm$ 0.001 & 0.039 $\pm$ 0.004 & 0.038 $\pm$ 0.000 & 0.015 $\pm$ 0.001 & {\bf 0.012 $\pm$ 0.001 } & 0.018 $\pm$ 0.001 \\
3dpoly & 0.030 $\pm$ 0.006 & {\bf 0.000 $\pm$ 0.000 } & 0.001 $\pm$ 0.000 & 0.344 $\pm$ 0.025 & 0.081 $\pm$ 0.015 & 0.114 $\pm$ 0.016 & 0.366 $\pm$ 0.031 \\
linear & 0.003 $\pm$ 0.000 & {\bf 0.001 $\pm$ 0.000 } & 0.016 $\pm$ 0.008 & 0.002 $\pm$ 0.000 & 0.009 $\pm$ 0.000 & 0.008 $\pm$ 0.000 & 0.010 $\pm$ 0.001 \\
randpw & 0.021 $\pm$ 0.004 & 0.055 $\pm$ 0.009 & 0.069 $\pm$ 0.010 & 0.157 $\pm$ 0.032 & 0.015 $\pm$ 0.002 & {\bf 0.013 $\pm$ 0.002 } & 0.028 $\pm$ 0.004 \\
abspos & 0.004 $\pm$ 0.000 & 0.007 $\pm$ 0.000 & 0.013 $\pm$ 0.003 & {\bf 0.003 $\pm$ 0.000 } & 0.010 $\pm$ 0.001 & 0.007 $\pm$ 0.000 & 0.013 $\pm$ 0.001 \\
sqrpos & 0.008 $\pm$ 0.002 & {\bf 0.004 $\pm$ 0.000 } & 0.008 $\pm$ 0.003 & 0.025 $\pm$ 0.003 & 0.013 $\pm$ 0.002 & 0.018 $\pm$ 0.002 & 0.109 $\pm$ 0.008 \\
band & 0.026 $\pm$ 0.001 & 0.044 $\pm$ 0.001 & 0.044 $\pm$ 0.001 & 0.056 $\pm$ 0.001 & 0.018 $\pm$ 0.001 & {\bf 0.014 $\pm$ 0.001 } & 0.020 $\pm$ 0.001 \\
invband & 0.026 $\pm$ 0.001 & 0.044 $\pm$ 0.001 & 0.044 $\pm$ 0.001 & 0.046 $\pm$ 0.001 & 0.018 $\pm$ 0.001 & {\bf 0.015 $\pm$ 0.001 } & 0.020 $\pm$ 0.001 \\
steplinear & 0.042 $\pm$ 0.001 & 0.064 $\pm$ 0.001 & 0.066 $\pm$ 0.002 & 0.079 $\pm$ 0.001 & 0.036 $\pm$ 0.001 & {\bf 0.032 $\pm$ 0.001 } & 0.032 $\pm$ 0.001 \\
pwlinear & {\bf 0.005 $\pm$ 0.000 } & 0.010 $\pm$ 0.000 & 0.013 $\pm$ 0.002 & 0.019 $\pm$ 0.000 & 0.011 $\pm$ 0.001 & 0.008 $\pm$ 0.001 & 0.014 $\pm$ 0.001
    \end{tabular}
    \caption{$n=2000$, $n_z=1$, $n_x=1$, $\gamma=.8$}
\end{figure}

\begin{figure}[H]
    \centering
    \scriptsize
    \begin{tabular}{c||c|c|c|c}
    & NystromRKHS & 2SLS & Reg2SLS & RFIV\\
    \hline
    abs & 0.026 $\pm$ 0.010 & 0.025 $\pm$ 0.001 & 0.054 $\pm$ 0.007 & {\bf 0.023 $\pm$ 0.001 } \\
2dpoly & 0.033 $\pm$ 0.006 & {\bf 0.002 $\pm$ 0.000 } & 0.361 $\pm$ 0.059 & 0.292 $\pm$ 0.012 \\
sigmoid & 0.015 $\pm$ 0.006 & {\bf 0.006 $\pm$ 0.000 } & 0.096 $\pm$ 0.016 & 0.014 $\pm$ 0.001 \\
sin & 0.019 $\pm$ 0.007 & {\bf 0.012 $\pm$ 0.001 } & 0.142 $\pm$ 0.024 & 0.016 $\pm$ 0.001 \\
frequentsin & 0.131 $\pm$ 0.007 & 0.117 $\pm$ 0.001 & 0.116 $\pm$ 0.003 & {\bf 0.069 $\pm$ 0.003 } \\
abssqrt & 0.027 $\pm$ 0.010 & {\bf 0.018 $\pm$ 0.001 } & 0.026 $\pm$ 0.004 & 0.019 $\pm$ 0.001 \\
step & 0.036 $\pm$ 0.006 & 0.028 $\pm$ 0.001 & 0.116 $\pm$ 0.017 & {\bf 0.021 $\pm$ 0.001 } \\
3dpoly & 0.018 $\pm$ 0.008 & {\bf 0.000 $\pm$ 0.000 } & 0.021 $\pm$ 0.003 & 0.416 $\pm$ 0.041 \\
linear & 0.015 $\pm$ 0.005 & {\bf 0.002 $\pm$ 0.000 } & 0.120 $\pm$ 0.019 & 0.012 $\pm$ 0.001 \\
randpw & {\bf 0.047 $\pm$ 0.010 } & 0.057 $\pm$ 0.011 & 0.448 $\pm$ 0.185 & 0.050 $\pm$ 0.009 \\
abspos & 0.019 $\pm$ 0.007 & {\bf 0.007 $\pm$ 0.001 } & 0.060 $\pm$ 0.010 & 0.014 $\pm$ 0.001 \\
sqrpos & 0.025 $\pm$ 0.005 & {\bf 0.004 $\pm$ 0.001 } & 0.065 $\pm$ 0.010 & 0.092 $\pm$ 0.007 \\
band & 0.056 $\pm$ 0.012 & 0.046 $\pm$ 0.001 & 0.053 $\pm$ 0.003 & {\bf 0.027 $\pm$ 0.002 } \\
invband & 0.051 $\pm$ 0.012 & 0.046 $\pm$ 0.001 & 0.052 $\pm$ 0.004 & {\bf 0.027 $\pm$ 0.002 } \\
steplinear & 0.087 $\pm$ 0.006 & 0.084 $\pm$ 0.001 & 0.103 $\pm$ 0.005 & {\bf 0.059 $\pm$ 0.002 } \\
pwlinear & 0.023 $\pm$ 0.008 & {\bf 0.010 $\pm$ 0.001 } & 0.026 $\pm$ 0.004 & 0.014 $\pm$ 0.001
    \end{tabular}
    \caption{$n=2000$, $n_z=5$, $n_x=1$, $\gamma=.6$}
\end{figure}

\begin{figure}[H]
    \centering
    \scriptsize
    \begin{tabular}{c||c|c|c|c}
    & NystromRKHS & 2SLS & Reg2SLS & RFIV\\
    \hline
    abs & 0.027 $\pm$ 0.011 & 0.035 $\pm$ 0.002 & 0.107 $\pm$ 0.016 & {\bf 0.021 $\pm$ 0.001 } \\
2dpoly & 0.050 $\pm$ 0.019 & {\bf 0.006 $\pm$ 0.000 } & 0.545 $\pm$ 0.080 & 0.282 $\pm$ 0.014 \\
sigmoid & 0.017 $\pm$ 0.009 & 0.014 $\pm$ 0.001 & 0.115 $\pm$ 0.023 & {\bf 0.013 $\pm$ 0.001 } \\
sin & 0.023 $\pm$ 0.009 & 0.020 $\pm$ 0.001 & 0.181 $\pm$ 0.045 & {\bf 0.017 $\pm$ 0.001 } \\
frequentsin & 0.136 $\pm$ 0.012 & 0.126 $\pm$ 0.001 & 0.117 $\pm$ 0.003 & {\bf 0.065 $\pm$ 0.003 } \\
abssqrt & 0.026 $\pm$ 0.008 & 0.030 $\pm$ 0.002 & 0.038 $\pm$ 0.006 & {\bf 0.018 $\pm$ 0.002 } \\
step & 0.035 $\pm$ 0.008 & 0.036 $\pm$ 0.001 & 0.135 $\pm$ 0.025 & {\bf 0.021 $\pm$ 0.002 } \\
3dpoly & 0.022 $\pm$ 0.018 & {\bf 0.001 $\pm$ 0.000 } & 0.035 $\pm$ 0.005 & 0.402 $\pm$ 0.045 \\
linear & 0.022 $\pm$ 0.008 & {\bf 0.007 $\pm$ 0.001 } & 0.123 $\pm$ 0.020 & 0.011 $\pm$ 0.001 \\
randpw & {\bf 0.047 $\pm$ 0.009 } & 0.061 $\pm$ 0.010 & 0.457 $\pm$ 0.165 & 0.051 $\pm$ 0.011 \\
abspos & 0.022 $\pm$ 0.008 & 0.015 $\pm$ 0.001 & 0.082 $\pm$ 0.015 & {\bf 0.013 $\pm$ 0.001 } \\
sqrpos & 0.042 $\pm$ 0.017 & {\bf 0.008 $\pm$ 0.001 } & 0.129 $\pm$ 0.020 & 0.086 $\pm$ 0.006 \\
band & 0.056 $\pm$ 0.013 & 0.056 $\pm$ 0.001 & 0.062 $\pm$ 0.007 & {\bf 0.027 $\pm$ 0.002 } \\
invband & 0.052 $\pm$ 0.012 & 0.058 $\pm$ 0.002 & 0.060 $\pm$ 0.006 & {\bf 0.026 $\pm$ 0.002 } \\
steplinear & 0.102 $\pm$ 0.013 & 0.097 $\pm$ 0.002 & 0.099 $\pm$ 0.005 & {\bf 0.059 $\pm$ 0.003 } \\
pwlinear & 0.031 $\pm$ 0.008 & 0.017 $\pm$ 0.002 & 0.033 $\pm$ 0.006 & {\bf 0.014 $\pm$ 0.001 }
    \end{tabular}
    \caption{$n=2000$, $n_z=10$, $n_x=1$, $\gamma=.6$}
\end{figure}

\begin{figure}[H]
    \centering
    \scriptsize
    \begin{tabular}{c||c|c|c|c}
    & NystromRKHS & 2SLS & Reg2SLS & RFIV\\
    \hline
    abs & 0.051 $\pm$ 0.002 & 0.262 $\pm$ 0.076 & {\bf 0.031 $\pm$ 0.002 } & 0.038 $\pm$ 0.001 \\
2dpoly & 0.226 $\pm$ 0.012 & 0.106 $\pm$ 0.033 & {\bf 0.105 $\pm$ 0.027 } & 0.316 $\pm$ 0.013 \\
sigmoid & 0.025 $\pm$ 0.002 & 0.198 $\pm$ 0.060 & 0.056 $\pm$ 0.002 & {\bf 0.015 $\pm$ 0.001 } \\
sin & 0.035 $\pm$ 0.002 & 0.222 $\pm$ 0.066 & 0.077 $\pm$ 0.006 & {\bf 0.022 $\pm$ 0.001 } \\
frequentsin & 0.140 $\pm$ 0.002 & 0.386 $\pm$ 0.084 & 0.114 $\pm$ 0.001 & {\bf 0.108 $\pm$ 0.002 } \\
abssqrt & 0.037 $\pm$ 0.002 & 0.288 $\pm$ 0.087 & 0.025 $\pm$ 0.001 & {\bf 0.025 $\pm$ 0.001 } \\
step & 0.045 $\pm$ 0.002 & 0.234 $\pm$ 0.064 & 0.076 $\pm$ 0.002 & {\bf 0.025 $\pm$ 0.001 } \\
3dpoly & 0.308 $\pm$ 0.030 & {\bf 0.009 $\pm$ 0.003 } & 0.027 $\pm$ 0.004 & 0.414 $\pm$ 0.034 \\
linear & 0.040 $\pm$ 0.002 & 0.124 $\pm$ 0.039 & 0.058 $\pm$ 0.006 & {\bf 0.014 $\pm$ 0.001 } \\
randpw & 0.131 $\pm$ 0.015 & 0.266 $\pm$ 0.163 & 0.161 $\pm$ 0.028 & {\bf 0.077 $\pm$ 0.011 } \\
abspos & 0.034 $\pm$ 0.002 & 0.185 $\pm$ 0.057 & 0.043 $\pm$ 0.002 & {\bf 0.017 $\pm$ 0.001 } \\
sqrpos & 0.111 $\pm$ 0.008 & 0.088 $\pm$ 0.028 & {\bf 0.029 $\pm$ 0.002 } & 0.097 $\pm$ 0.006 \\
band & 0.060 $\pm$ 0.002 & 0.327 $\pm$ 0.085 & 0.055 $\pm$ 0.001 & {\bf 0.038 $\pm$ 0.001 } \\
invband & 0.060 $\pm$ 0.002 & 0.311 $\pm$ 0.089 & 0.054 $\pm$ 0.001 & {\bf 0.039 $\pm$ 0.001 } \\
steplinear & 0.161 $\pm$ 0.004 & 0.457 $\pm$ 0.115 & 0.100 $\pm$ 0.003 & {\bf 0.090 $\pm$ 0.002 } \\
pwlinear & 0.052 $\pm$ 0.003 & 0.187 $\pm$ 0.058 & {\bf 0.017 $\pm$ 0.001 } & 0.018 $\pm$ 0.001
    \end{tabular}
    \caption{$n=2000$, $n_z=5$, $n_x=5$, $\gamma=.6$}
\end{figure}

\begin{figure}[H]
    \centering
    \scriptsize
    \begin{tabular}{c||c|c|c|c}
    & NystromRKHS & 2SLS & Reg2SLS & RFIV\\
    \hline
    abs & 0.143 $\pm$ 0.005 & 10050.672 $\pm$ 13267.141 & 0.122 $\pm$ 0.011 & {\bf 0.049 $\pm$ 0.001 } \\
2dpoly & 0.595 $\pm$ 0.025 & 5890.128 $\pm$ 8261.553 & 4.510 $\pm$ 1.245 & {\bf 0.346 $\pm$ 0.014 } \\
sigmoid & 0.045 $\pm$ 0.003 & 11712.144 $\pm$ 16799.716 & 0.091 $\pm$ 0.005 & {\bf 0.017 $\pm$ 0.001 } \\
sin & 0.058 $\pm$ 0.003 & 13769.428 $\pm$ 20805.861 & 0.114 $\pm$ 0.006 & {\bf 0.029 $\pm$ 0.001 } \\
frequentsin & 0.136 $\pm$ 0.004 & 12928.749 $\pm$ 19554.361 & 0.144 $\pm$ 0.004 & {\bf 0.120 $\pm$ 0.002 } \\
abssqrt & 0.062 $\pm$ 0.004 & 12764.707 $\pm$ 17195.564 & 0.079 $\pm$ 0.005 & {\bf 0.034 $\pm$ 0.001 } \\
step & 0.064 $\pm$ 0.003 & 12187.342 $\pm$ 17814.756 & 0.109 $\pm$ 0.004 & {\bf 0.027 $\pm$ 0.001 } \\
3dpoly & 0.648 $\pm$ 0.039 & 432.572 $\pm$ 596.731 & {\bf 0.061 $\pm$ 0.005 } & 0.444 $\pm$ 0.029 \\
linear & 0.080 $\pm$ 0.002 & 6964.376 $\pm$ 9566.774 & 0.107 $\pm$ 0.006 & {\bf 0.016 $\pm$ 0.001 } \\
randpw & 0.272 $\pm$ 0.029 & 1882.000 $\pm$ 1998.862 & 0.682 $\pm$ 0.539 & {\bf 0.093 $\pm$ 0.013 } \\
abspos & 0.067 $\pm$ 0.003 & 8841.523 $\pm$ 11921.282 & 0.095 $\pm$ 0.005 & {\bf 0.020 $\pm$ 0.001 } \\
sqrpos & 0.243 $\pm$ 0.010 & 4250.312 $\pm$ 5449.534 & 0.126 $\pm$ 0.014 & {\bf 0.105 $\pm$ 0.006 } \\
band & 0.078 $\pm$ 0.004 & 20401.368 $\pm$ 29655.000 & 0.090 $\pm$ 0.004 & {\bf 0.049 $\pm$ 0.002 } \\
invband & 0.079 $\pm$ 0.004 & 11210.315 $\pm$ 14271.847 & 0.090 $\pm$ 0.005 & {\bf 0.048 $\pm$ 0.002 } \\
steplinear & 0.212 $\pm$ 0.005 & 22217.181 $\pm$ 33274.806 & 0.141 $\pm$ 0.005 & {\bf 0.110 $\pm$ 0.002 } \\
pwlinear & 0.075 $\pm$ 0.003 & 9280.655 $\pm$ 12159.776 & 0.041 $\pm$ 0.004 & {\bf 0.021 $\pm$ 0.001 }
    \end{tabular}
    \caption{$n=2000$, $n_z=10$, $n_x=10$, $\gamma=.6$}
\end{figure}

\begin{figure}[H]
    \centering
    \scriptsize
    \begin{tabular}{c||c|c|c|c|c}
    & AGMM & KLayerFixed & KLayerTrained & CentroidMMD & KLossMMD\\
    \hline
    abs & {\bf 0.062 $\pm$ 0.003 } & 0.190 $\pm$ 0.006 & 0.127 $\pm$ 0.007 & 0.114 $\pm$ 0.007 & 0.193 $\pm$ 0.007 \\
2dpoly & {\bf 0.099 $\pm$ 0.006 } & 0.971 $\pm$ 0.040 & 0.240 $\pm$ 0.014 & 0.204 $\pm$ 0.022 & 0.467 $\pm$ 0.023 \\
sigmoid & 0.040 $\pm$ 0.001 & 0.063 $\pm$ 0.002 & {\bf 0.024 $\pm$ 0.001 } & 0.058 $\pm$ 0.003 & 0.043 $\pm$ 0.003 \\
sin & 0.074 $\pm$ 0.002 & 0.076 $\pm$ 0.002 & {\bf 0.057 $\pm$ 0.002 } & 0.098 $\pm$ 0.003 & 0.083 $\pm$ 0.004 \\
frequentsin & 0.158 $\pm$ 0.002 & {\bf 0.120 $\pm$ 0.002 } & 0.128 $\pm$ 0.002 & 0.181 $\pm$ 0.004 & 0.160 $\pm$ 0.007 \\
abssqrt & 0.060 $\pm$ 0.003 & {\bf 0.058 $\pm$ 0.004 } & 0.060 $\pm$ 0.003 & 0.093 $\pm$ 0.004 & 0.090 $\pm$ 0.007 \\
step & 0.066 $\pm$ 0.002 & 0.076 $\pm$ 0.002 & {\bf 0.050 $\pm$ 0.001 } & 0.088 $\pm$ 0.003 & 0.069 $\pm$ 0.003 \\
3dpoly & {\bf 0.426 $\pm$ 0.027 } & 0.716 $\pm$ 0.037 & 0.491 $\pm$ 0.029 & 0.496 $\pm$ 0.030 & 0.526 $\pm$ 0.032 \\
linear & 0.020 $\pm$ 0.001 & 0.142 $\pm$ 0.003 & {\bf 0.013 $\pm$ 0.001 } & 0.029 $\pm$ 0.002 & 0.027 $\pm$ 0.001 \\
randpw & {\bf 0.127 $\pm$ 0.020 } & 0.449 $\pm$ 0.051 & 0.165 $\pm$ 0.024 & 0.169 $\pm$ 0.025 & 0.218 $\pm$ 0.030 \\
abspos & {\bf 0.034 $\pm$ 0.002 } & 0.090 $\pm$ 0.003 & 0.039 $\pm$ 0.002 & 0.057 $\pm$ 0.003 & 0.060 $\pm$ 0.003 \\
sqrpos & {\bf 0.059 $\pm$ 0.003 } & 0.347 $\pm$ 0.013 & 0.131 $\pm$ 0.007 & 0.113 $\pm$ 0.009 & 0.178 $\pm$ 0.009 \\
band & 0.088 $\pm$ 0.003 & {\bf 0.068 $\pm$ 0.002 } & 0.074 $\pm$ 0.003 & 0.117 $\pm$ 0.004 & 0.130 $\pm$ 0.037 \\
invband & 0.088 $\pm$ 0.003 & {\bf 0.073 $\pm$ 0.005 } & 0.077 $\pm$ 0.003 & 0.114 $\pm$ 0.004 & 0.120 $\pm$ 0.026 \\
steplinear & 0.176 $\pm$ 0.003 & 0.197 $\pm$ 0.004 & {\bf 0.133 $\pm$ 0.003 } & 0.218 $\pm$ 0.005 & 0.170 $\pm$ 0.010 \\
pwlinear & 0.049 $\pm$ 0.001 & 0.074 $\pm$ 0.002 & {\bf 0.033 $\pm$ 0.001 } & 0.063 $\pm$ 0.002 & 0.049 $\pm$ 0.002
    \end{tabular}
    \caption{$n=2000$, $n_z=10$, $n_x=10$, $\gamma=.6$}
\end{figure}

\begin{figure}[H]
    \centering
    \scriptsize
    \begin{tabular}{c||c|c|c|c}
     $p=$ & 1000 & 10000 & 100000 & 1000000\\
    \hline
    SpLin & 0.020 $\pm$ 0.003 & 0.021 $\pm$ 0.003 & - & -\\
    StSpLin & 0.020 $\pm$ 0.002  & 0.023 $\pm$ 0.002 & 0.033 $\pm$ 0.002 & 0.050 $\pm$ 0.004
    \end{tabular}
    \caption{$n=400$, $n_z=n_x:=p$, $\gamma=.6$, $h_0(x[0])=x[0]$}
\end{figure}

\begin{figure}[H]
    \centering
    \scriptsize
    \begin{tabular}{c||c|c|c|c}
      & DeepGMM (\cite{bennett2019deep}) & AGMM  & KLayerTrained \\
    \hline
    $\text{MNIST}_z$ & 0.12 $\pm$ 0.07 & 0.04 $\pm$ 0.03 & 0.05 $\pm$ 0.02 \\
    $\text{MNIST}_x$ & 0.34 $\pm$ 0.21  & 0.24 $\pm$ 0.08 & 0.36 $\pm$ 0.20  \\
    $\text{MNIST}_{xz}$ & 0.26 $\pm$ 0.16 & 0.21 $\pm$ 0.07 & 0.26 $\pm$ 0.11
    \end{tabular}
    \caption{MSE on the high-dimensional DGPs}
\end{figure}

\subsection{Experiments with Image Data}\label{app:mnist}
In this section, we describe the experimental setup for our experiments with high-dimensional data using the MNIST dataset. We replicate the data-generating process of \cite{bennett2019deep}. We present a full description here for completeness.

\paragraph{The Data-Generating Process}
We begin by describing a low-dimensional DGP which will define a mapping for $x$ or $z$ or both to be MNIST images. The data-generating process is:
\begin{align}
    &y = g_0(x^{\text{low}}) + e + \delta \\
    &z^{\text{low}} \sim \text{Uniform}([-3,3]^2) \\
    &x^{\text{low}} = z^{\text{low}}_1 + e + \gamma \\
    &e \sim \mathcal{N}(0,1) , \delta,\gamma \sim \mathcal{N}(0,0.1).
\end{align}
Let $\pi(x) = \text{round}\left(\min(\max(1.5x+5,0),9) \right)$. $\pi$ is a transformation function that maps inputs to an integer between 0 and 9. Let $\text{RandomImage}(d)$ be a function which selects a random MNIST image from the class of images corresponding to digit $d$. The three high-dimensional scenarios are:
\begin{align}
    &\text{MNIST}_Z: x = x^{\text{low}}, z = \text{RandomImage}(\pi(z^{\text{low}}_1)) \\
    &\text{MNIST}_X: x = \text{RandomImage}(\pi(x^{\text{low}})), z = z^{\text{low}} \\
    &\text{MNIST}_{XZ}: x = \text{RandomImage}(\pi(x^{\text{low}})), z = \text{RandomImage}(\pi(z^{\text{low}}_1)).
\end{align}

We use the function $g_0(x)=|x|$ to compare with \cite{bennett2019deep} but in general, the other functional forms described above can also be used. Similar to \cite{bennett2019deep} we normalize the data so that $y$ has zero mean and unit standard deviation.

We evaluate the performance of our AGMM and KLayerTrained estimators on these 3 data-generating processes with 20,000 train samples and 2,000 test samples and compare their performance to that achieved when we evaluate \cite{bennett2019deep}'s code (performance is measured by the average mean squared error of the predictions on test data).

\paragraph{Setup}
We describe more details about our experimental setup for the MNIST experiments here. We run 10 Monte-Carlo runs of each experiment and report the average MSE and the standard deviation in the MSE achieved.

\paragraph{Architectures}
We use a 4-layer convolutional architecture in all cases where the input to the network is an image. This consists of 2 convolutional layers with a 3x3 kernel followed by two fully connected layers with 9216 and 512 hidden units respectively. A ReLU activation is applied after each layer. Along with that, a max-pooling operation is applied after the first two convolutional layers and a dropout operation (with dropout probability 0.1) is applied before each fully connected layer. When the instrument or treatment is low-dimensional we use a 2 layer fully connected neural network with 200 neurons in the hidden layer along with the dropout function as before.
All networks use ReLU as the activation function.

\paragraph{Early Stopping}
We utilize the early stopping procedure proposed in \cite{bennett2019deep} which works as follows. In addition to the 20,000 training samples, 10,000 samples are used for preparing a set of candidate adversary functions prior to training. During training at each epoch, the maximum error incurred by the learner against the candidates in this pre-computed list is recorded. The early stopping selects the model whose maximum error as computer above is the smallest.

\paragraph{Hyper-Parameters}
We use a batch size of 100 samples, and run for 200 epochs where an epoch is defined as one full pass over the train set. We have as hyper-parameters learning rates for the learner and adversary networks, the regularization terms for the weights of the learner and the adversary, and a regularization term on the norm of the output of the adversary network.
For the $\text{MNIST}_x$ experiment, we saw best results when the weight penalizations on both the learner and the adversary were set to very small values as compared to the other two experiments.

\newpage

\section{Proofs from Section~\ref{sec:estimation} and \Cref{app:further-main}}

\subsection{Preliminary Lemmas}

\begin{lemma}\label{lem:population-lower-bound}
Let $f_h$, be any test function that satisfies: $\|f_h - T(h_*-h)\|_2\leq \epsilon$ and let
\[
\Psi(h, f) := \E\brk*{\psi(y\midsem{}h(x))\, f(z)}.
\]
Then:
\begin{equation}
    \frac{1}{\|f_h\|_2} \left(\Psi\left(h, f_h\right) - \Psi(h_*, f_h)\right) \geq \|T(h-h_*)\|_2 - 2\epsilon_n
\end{equation}
\end{lemma}
\begin{proof}
Let $f_h^*=T(h_*-h)$ and observe that by the tower law of expectations:
\begin{equation}
    \frac{1}{\|f_h\|_2}(\Psi(h, f_h) - \Psi(h_*, f_{h}))
    = \frac{\E[(h_*-h)(x)\, f_h(z)]}{\|f_h\|_2} = \frac{\E[f_h^*(z)\, f_h(z)]}{\|f_h\|_2}
\end{equation}
However, observe that by the Cauchy-Schwarz inequality we have:
\begin{align}
    \E[f_h^*(Z)\, f_h(Z)] =~& \E[f_h(Z)^2] + \E[f_h(Z) (f_h^*(Z) - f_h(Z))] \geq \|f_h\|_2^2 -\left| \E[f_h(Z) (f_h^*(Z) - f_h(Z))] \right|\\
    \geq~& \|f_h\|_2^2 -\sqrt{\E[f_h(Z)^2]} \sqrt{\E[(f_h^*(Z) - f_h(Z))^2]}\\
    \geq~& \|f_h\|_2^2 -\|f_h\|_2 \|f_h^* - f_h\|_2\\
    \geq~& \|f_h\|_2^2 - \epsilon_n \|f_h\|_2
\end{align}
Thus we have:
\begin{equation}
    \frac{1}{\|f_h\|_2}(\Psi(h, f_{h}) - \Psi(h_*, f_h))
    \geq \|f_h\|_2-\epsilon_n
\end{equation}
Finally, by a triangle inequality,
\begin{equation}
\|f_h\|_2 \geq \|f_h^*\|_2 - \|f_h^* - f_h\|_2 \geq \|f_h^*\|_2 - \epsilon_n.
\end{equation}
Hence, we can conclude that:
\begin{equation}
    \frac{1}{\|f_h\|_2}(\Psi(h, f_{h}) - \Psi(h_*, f_h)) \geq \|f_h^*\|_2 - 2\epsilon_n = \|T(h-h_*)\|_2 - 2\epsilon_n
\end{equation}
\end{proof}

\subsection{Proof of \Cref{thm:reg-main-error}}

\begin{proof} For convenience let:
\begin{align}
    \Psi(h, f) :=~& \E\brk*{\psi(y\midsem{}h(x))\, f(z)} = \E[T(h_0-h)(z)\, f(z)] \tag{by conditional moment restriction}\\
    \Psi_n(h, f) :=~& \frac{1}{n}\sum_{i=1}^{n} \psi(y_i\midsem{} h(x_i))\,f(z_i)
\end{align}
Moreover, for our choice of $\delta$ as described in the statement of the theorem, let:
\begin{align}
    \mcH_{B} :=~& \left\{ h \in \mcH: \|h\|_\mcH^2 \leq B\right\}\\
    \mcF_{U} :=~& \left\{ f \in \mcF: \|f\|_\mcF^2 \leq U\right\}
\end{align}
Moreover, let:
\begin{align}
    \Psi_n^\lambda(h, f) =~& \Psi_n(h, f) - \lambda \left(\|f\|_{\mcF}^2 + \frac{U}{\delta^2} \|f\|_{2,n}^2\right)\\
    \Psi^{\lambda}(h, f) =~& \Psi(h, f) - \lambda \left(\frac{2}{3}\|f\|_{\mcF}^2 + \frac{U}{\delta^2} \|f\|_{2}^2\right)
\end{align}
Thus our estimate can be written as:
\begin{equation}
    \hat{h} := \argmin_{h\in \mcH} \sup_{f\in \mcF} \Psi_n^{\lambda}(h, f) + \mu \|h\|_{\mcH}^2
\end{equation}

\paragraph{Relating empirical and population regularization.} As a preliminary observation, we have that by Theorem 14.1 of \cite{wainwright2019high}, w.p. $1-\zeta$:
\begin{equation}
    \forall f\in \mcF_{3U}: \left|\|f\|_{n,2}^2 - \|f\|_2^2 \right| \leq \frac{1}{2} \|f\|_2^2 + \delta^2
\end{equation}
for our choice of $\delta:=\delta_n + c_0 \sqrt{\frac{\log(c_1/\zeta)}{n}}$, where $\delta_n$ upper bounds the critical radius of $\mcF_{3U}$ and $c_0, c_1$ are universal constants. Moreover, for any $f$, with $\|f\|_{\mcF}^2\geq 3U$, we can consider the function $f \sqrt{3U}/\|f\|_{\mcF}$, which also belongs to $\mcF_{3U}$, since $\mcF$ is star-convex. Thus we can apply the above lemma to this re-scaled function and multiply both sides by $\|f\|_{\mcF}^2/(3U)$, leading to:
\begin{equation}
    \forall f\in \mcF \text{ s.t. } \|f\|_{\mcF}^2 \geq 3U: \left|\|f\|_{n,2}^2 - \|f\|_2^2 \right| \leq \frac{1}{2} \|f\|_2^2 + \delta^2\frac{\|f\|_{\mcF}^2}{3U}
\end{equation}
Thus overall, we have:
\begin{equation}\label{eqn:reg-pop-emp}
    \forall f\in \mcF: \left|\|f\|_{n,2}^2 - \|f\|_2^2 \right| \leq \frac{1}{2} \|f\|_2^2 + \delta^2 \max\left\{1, \frac{\|f\|_{\mcF}^2}{3U}\right\}
\end{equation}
Thus we have that w.p. $1-\zeta$:
\begin{align}
    \forall f\in \mcF: \|f\|_{\mcF}^2 + \frac{U}{\delta^2} \|f\|_{2,n}^2 \geq~& \|f\|_{\mcF}^2 + \frac{U}{\delta^2} \left(\|f\|_{2}^2 - \delta^2 \max\left\{1, \frac{\|f\|_{\mcF}^2}{3U}\right\}\right)\\
    \geq~& \|f\|_{\mcF}^2 + \frac{U}{\delta^2} \|f\|_{2}^2 - \max\left\{U, \frac{1}{3} \|f\|_{\mcF}^2\right\}\\
    \geq~& \frac{2}{3} \|f\|_{\mcF}^2 + \frac{U}{\delta^2} \|f\|_{2}^2 - U \label{eqn:reg-empirical-vs-pop}
\end{align}

\paragraph{Upper bounding centered empirical sup-loss.} We now argue that the centered empirical sup-loss: $\sup_{f \in \mcF} (\Psi_n(\hat{h}, f) - \Psi_n(h_*, f))$ is small. By the definition of $\hat{h}$:
\begin{equation}\label{eqn:reg-optimality}
\sup_{f\in \mcF} \Psi_n^{\lambda}(\hat{h}, f) \leq \sup_{f\in \mcF} \Psi_n^{\lambda}(h_*, f) + \mu \left(\|h_*\|_{\mcH}^2 - \|\hat{h}\|_{\mcH}^2\right)
\end{equation}

By Lemma~7 of \cite{Foster2019}, the fact that $\phi(y; h_*(x)) f(z)$ is $2$-Lipschitz with respect to $f(z)$ (since $y\in [-1, 1]$ and $\|h_*\|_{\infty}\in [-1, 1]$) and by our choice of $\delta:=\delta_n + c_0 \sqrt{\frac{\log(c_1/\zeta)}{n}}$, where $\delta_n$ is an upper bound on the critical radius of $\mcF_{3U}$, w.p. $1-\zeta$:
\begin{align}
 \forall f\in \mcF_{3U}: \left|\Psi_n(h_*, f) - \Psi(h_*, f)\right| \leq 36\delta \|f\|_2 + 36\delta^2
\end{align}
Thus, if $\|f\|_{\mcF} \geq \sqrt{3U}$, we can apply the latter inequality for the function $f \sqrt{3U}/\|f\|_{\mcF}$, which falls in $\mcF_{3U}$, and then multiply both sides by $\|f\|_{\mcF}/\sqrt{3U}$ to get:
\begin{align}\label{eqn:reg-concentration}
    \forall f \in \mcF: \left|\Psi_n(h_*, f) - \Psi(h_*, f)\right| \leq 36\delta \|f\|_2 + 36\delta^2 \max\left\{1, \frac{\|f\|_{\mcF}}{\sqrt{3U}}\right\}
\end{align}

By Equations~\eqref{eqn:reg-empirical-vs-pop} and \eqref{eqn:reg-concentration}, we have that w.p. $1-2\zeta$:
\begin{align}
   \sup_{f\in \mcF} \Psi_n^{\lambda}(h_*, f)
    =~& \sup_{f\in \mcF} \left(\Psi_n(h_*, f) - \lambda \left(\|f\|_{\mcF}^2 + \frac{U}{\delta^2} \|f\|_{2,n}^2\right)\right)\\
    \leq~& \sup_{f\in \mcF} \left(\Psi(h_*, f) + 36 \delta^2 + \frac{36\delta^2}{\sqrt{3U}}\|f\|_{\mcF} + 36 \delta \|f\|_2 - \lambda \left(\|f\|_{\mcF}^2 + \frac{U}{\delta^2} \|f\|_{2,n}^2\right)\right)\\
    \leq~& \sup_{f\in \mcF} \left(\Psi(h_*, f) + 36 \delta^2 + \frac{36\delta^2}{\sqrt{3U}}\|f\|_{\mcF}  + 36 \delta \|f\|_2 - \lambda \left(\frac{2}{3}\|f\|_{\mcF}^2 + \frac{U}{\delta^2} \|f\|_{2}^2\right) + \lambda U\right)\\
    \leq~& \sup_{f\in \mcF} \Psi^{\lambda/2}(h_*, f) + 36 \delta^2 + \lambda U \\
    &+ \sup_{f\in \mcF} \left( \frac{36\delta^2}{\sqrt{3U}}\|f\|_{\mcF} - \frac{\lambda}{2}\frac{2}{3}\|f\|_{\mcF}^2\right) + \sup_{f\in \mcF} \left(36 \delta \|f\|_2 - \frac{\lambda}{2} \frac{U}{\delta^2} \|f\|_{2}^2\right)\\
\end{align}
Moreover, observe that for any norm $\|\cdot\|$ and any constants $a,b>0$:
\begin{align}
    \sup_{f\in \mcF} \left(a\|f\| - b \|f\|^2\right) \leq \frac{a^2}{4b}
\end{align}
Thus if we assume that $\lambda\geq \delta^2/U$, we have:
\begin{align}
\sup_{f\in \mcF} \left( \delta^2\frac{36}{\sqrt{3U}}\|f\|_{\mcF} - \frac{\lambda}{2}\frac{2}{3}\|f\|_{\mcF}^2\right) \leq~& \frac{36^2}{4}\frac{\delta^4}{U\lambda} \leq 324 \delta^2\\
    \sup_{f\in \mcF} \left(36 \delta \|f\|_2 - \frac{\lambda}{2} \frac{U}{\delta^2} \|f\|_{2}^2\right)\leq~& \frac{36^2 \delta^4}{2\lambda U} \leq 648 \delta^2
\end{align}
Thus we have:
\begin{align}
    \sup_{f\in \mcF} \Psi_n^{\lambda}(h_*, f)\leq \sup_{f\in \mcF} \Psi^{\lambda/2}(h_*, f) + \lambda U + O(\delta^2)
\end{align}
Moreover:
\begin{align}
\sup_{f\in \mcF} \Psi_n^{\lambda}(\hat{h}, f) =~& \sup_{f \in \mcF} \left(\Psi_n(\hat{h}, f) - \Psi_n(h_*, f) + \Psi_n(h_*, f) - \lambda\left(\|f\|_{\mcF}^2 + \frac{U}{\delta^2} \|f\|_{2,n}^2\right)\right)\\ \geq~& \sup_{f \in \mcF} \left(\Psi_n(\hat{h}, f) - \Psi_n(h_*, f) - 2\lambda \left(\|f\|_{\mcF}^2 + \frac{U}{\delta^2} \|f\|_{2,n}^2\right)\right) \\
& + \inf_{f\in \mcF} \left(\Psi_n(h_*, f) + \lambda \left(\|f\|_{\mcF}^2 + \frac{U}{\delta^2} \|f\|_{2,n}^2\right)\right)\\
    =~& \sup_{f \in \mcF} \left(\Psi_n(\hat{h}, f) - \Psi_n(h_*, f)- 2\lambda \left(\|f\|_{\mcF}^2 + \frac{U}{\delta^2} \|f\|_{2,n}^2\right)\right) - \sup_{f\in \mcF} \Psi_n^{\lambda}(h_*, f)
\end{align}
Combining this with Equation~\eqref{eqn:reg-optimality} yields:
\begin{align}
    \sup_{f \in \mcF} \left(\Psi_n(\hat{h}, f) - \Psi_n(h_*, f)- 2\lambda \left(\|f\|_{\mcF}^2 + \frac{U}{\delta^2} \|f\|_{2,n}^2\right)\right) \leq~& 2\, \sup_{f\in \mcF} \Psi_n^{\lambda}(h_*, f) + \mu \left(\|h_*\|_{\mcH}^2 - \|\hat{h}\|_{\mcH}^2\right)\\
    \leq~& O(\delta^2) + \lambda U + 2\sup_{f\in \mcF} \Psi^{\lambda/2}(h_*, f)\\
    & + \mu \left(\|h_*\|_{\mcH}^2 - \|\hat{h}\|_{\mcH}^2\right)
\end{align}

\paragraph{Lower bounding centered empirical sup-loss.} For any $h$, let
\begin{equation}
f_h := \arginf_{f\in \mcF_{L^2\|h-h_*\|_\mcH^2}} \|f-T(h_* - h)\|_2.
\end{equation}
and observe that by our assumption, for any $h\in \mcH$: $\|f_h - T(h_*-h)\|_2\leq \eta_n$.

Suppose that $\|f_{\hat{h}}\|_2\geq \delta$ and let $r = \frac{\delta}{2\|f_{\hat{h}}\|_2}\in [0, 1/2]$. Then observe that since $f_{\hat{h}}\in \mcF_{L\|h-h_*\|_\mcH}$ and $\mcF$ is star-convex, we also have that $r f_h\in \mcF_{L\|h-h_*\|_\mcH}$. Thus we can lower bound the supremum by its evaluation at $r\, f_h$:
\begin{align}
    \sup_{f \in \mcF}\left(\Psi_n(\hat{h}, f) - \Psi_n(h_*, f)- 2\lambda \left(\|f\|_{\mcF}^2 + \frac{U}{\delta^2} \|f\|_{2,n}^2\right)\right)\geq~& r (\Psi_n(\hat{h}, f_{\hat{h}}) - \Psi_n(h_*, f_{\hat{h}}))\\
    &- 2\lambda r^2 \left(\|f_{\hat{h}}\|_{\mcF}^2 + \frac{U}{\delta^2} \|f_{\hat{h}}\|_{2,n}^2\right)
\end{align}

Moreover, since $\delta_n$ upper bounds the critical radius of $\mcF_{3U}$, $\|f_{\hat{h}}\|_{\mcF}\leq L\|\hat{h}-h_*\|_{\mcH}$ and by Equation~\eqref{eqn:reg-pop-emp}:
\begin{align}
    r^2 \left(\|f_{\hat{h}}\|_\mcF^2 + \frac{U}{\delta^2} \|f_{\hat{h}}\|_{2,n}^2 \right)
    \leq~& \|f_{\hat{h}}\|_\mcF^2 + \frac{U}{\delta^2} r^2 \|f_{\hat{h}}\|_{2,n}^2\\
    \leq~& \|f_{\hat{h}}\|_{\mcF}^2 + \frac{U}{\delta^2} r^2 \left(2\|f_{\hat{h}}\|_2^2 + \delta^2 + \delta^2 \frac{\|f_{\hat{h}}\|_{\mcF}^2}{3U}\right)\\
    \leq~& \frac{4}{3} L^2 \|h-h_*\|_{\mcH}^2  + \frac{U}{2} + \frac{U}{4}
    \leq 2L^2 \|h-h_*\|_{\mcH}^2 + U
\end{align}
Thus we get:
\begin{align}
    \sup_{f \in \mcF}\left(\Psi_n(\hat{h}, f) - \Psi_n(h_*, f)- 2\lambda \left(\|f\|_{\mcF}^2 + \frac{U}{\delta^2} \|f\|_{2,n}^2\right)\right)\geq~& r (\Psi_n(\hat{h}, f_{\hat{h}}) - \Psi_n(h_*, f_{\hat{h}}))\\
    & - 4 \lambda L^2 \|h-h_*\|_{\mcH}^2 - 2\lambda U
\end{align}

Observe that:
\begin{align}
    \Psi_n(h, f_{h}) - \Psi_n(h_*, f_{h}) =~& \frac{1}{n} \sum_{i=1}^n  (h_*(x_i) - h(x_i))\, f_h(h_*-h)(z_i)\\
    \Psi(h, f_h) - \Psi(h_*, f_h) =~& \E[(h_*(x_i) - h(x_i))\, f_h(z_i)]
\end{align}
By Lemma~7 of \cite{Foster2019}, and by our choice of $\delta:=\delta_n + c_0 \sqrt{\frac{\log(c_1/\zeta)}{n}}$, where $\delta_n$ upper bounds the critical radius of $\mcG$, we have that w.p. $1-\zeta$: $\forall h$, such that $h-h_*\in \mcH_B$
\begin{align}
     \left|(\Psi_n(h, f_{h}) - \Psi_n(h_*, f_{h})) - (\Psi(h, f_h) - \Psi(h_*, f_h))\right| \leq~& 18\delta \sqrt{\E[(h_*(X) - h(X))^2\, f_h(Z)^2]} + 18\delta^2\\
     \leq~& 18 \delta  \sqrt{\E[f_h(Z)^2]} + 18\delta^2\\
     =~& 18 \delta \|f_h\|_2 + 18\delta^2 \label{eqn:diff-crit-radii}
\end{align}
where in the second inequality we used the fact that $h-h_*$ has range in $[-1,1]$, when $\|h-h_*\|_{\mcH}\leq B$. If $h-h_*$ has $\|h-h_*\|_{\mcH}^2\geq B$, we can apply the latter for $(h-h_*)\sqrt{B}/\|h-h_*\|_{\mcH}$ and multiply both sides by $\|h-h_*\|_{\mcH}^2/B$:
\begin{align}
     \left|(\Psi_n(h, f_{h}) - \Psi_n(h_*, f_{h})) - (\Psi(h, f_h) - \Psi(h_*, f_h))\right| \leq~& 18 \delta \|f_h\|_2 \frac{\|h-h_*\|_{\mcH}}{\sqrt{B}} + 18\delta^2 \frac{\|h-h_*\|_{\mcH}^2}{B}
\end{align}
Thus we have that for all $h\in \mcH$:
\begin{align}
     \left|(\Psi_n(h, f_{h}) - \Psi_n(h_*, f_{h})) - (\Psi(h, f_h) - \Psi(h_*, f_h))\right| \leq~& \left( 18 \delta \|f_h\|_2 + 18\delta^2\right) \max\left\{1, \frac{\|h-h_*\|_{\mcH}^2}{B}\right\}
\end{align}

Applying the latter bound for $h:=\hat{h}$ and multiplying by $r:=\frac{\delta}{2\|f_{\hat{h}}\|_2}\in [0, 1/2]$, yields:
\begin{align}
    r(\Psi_n(\hat{h}, f_{\hat{h}}) - \Psi_n(h_*, f_{\hat{h}}))
    \geq~& r(\Psi(\hat{h}, f_{\hat{h}}) - \Psi(h_*, f_{\hat{h}})) - 18\delta^2\max\left\{1, \frac{\|h-h_*\|_{\mcH}^2}{B}\right\}
\end{align}
Moreover, observe that by Lemma~\ref{lem:population-lower-bound} and the fact that $\|f_{\hat{h}} - T(h_*-\hat{h})\|_2\leq \eta_n$, we have:
\begin{align}
    r(\Psi(\hat{h}, f_{\hat{h}}) - \Psi(h_*, f_{\hat{h}})) \geq \frac{\delta}{2}\|T(h_*-\hat{h})\|_2 - \delta \eta_n
\end{align}
Thus we have:
\begin{align}
     \sup_{f \in \mcF}\left(\Psi_n(\hat{h}, f) - \Psi_n(h_*, f)- 2\lambda \left(\|f\|_{\mcF}^2 + \frac{U}{\delta^2} \|f\|_{2,n}^2\right)\right) \geq~& \frac{\delta}{2}\|T(h_*-h)\|_2 - \delta \eta_n \\
     & - 27\delta^2\max\left\{1, \frac{\|h-h_*\|_{\mcH}^2}{B}\right\} \\
     & - 4\lambda L^2 \|h-h_*\|_{\mcH}^2 - 2\lambda U
\end{align}

\paragraph{Combining upper and lower bound.} Combining the upper and lower bound on the centered population sup-loss we get that w.p. $1-3\zeta$: either $\|f_{\hat{h}}\|_2\leq \delta$ or:
\begin{align}
    \frac{\delta}{2} \|T(\hat{h}-h_*)\|_2 \leq~& O(\delta^2 + \delta \eta_n + \lambda U) + 2\sup_{f\in \mcF} \Psi^{\lambda/2}(h_*, f)\\
    & + 27\delta^2 \frac{\|\hat{h}-h_*\|_{\mcH}^2}{B} + 4\lambda L^2 \|\hat{h}-h_*\|_{\mcH}^2  + \mu \left(\|h_*\|_{\mcH}^2 - \|\hat{h}\|_{\mcH}^2\right)
\end{align}
We now control the last part. Since $\lambda \geq \delta^2/U$, the latter is upper bounded by:
\begin{align}
    \lambda \left(\frac{27\, U}{B} + 4L^2\right) \|\hat{h}-h_*\|_{\mcH}^2  + \mu \left(\|h_*\|_{\mcH}^2 - \|\hat{h}\|_{\mcH}^2\right)\leq~&  2\lambda \left(\frac{27U}{B} + 4L^2\right) \left(\|\hat{h}\|_{\mcH}^2 + \|h_*\|_{\mcH}^2\right)\\
    & + \mu \left(\|h_*\|_{\mcH}^2 - \|\hat{h}\|_{\mcH}^2\right)
\end{align}
Since $\mu \geq 2\lambda \left(\frac{27U}{B} + 4L^2\right)$, the latter is upper bounded by:
\begin{align}
    \left(2\lambda\left( \frac{27U}{B} + 4L^2 \right) + \mu\right) \|h_*\|_{\mcH}^2
\end{align}
Thus as long as $\mu \geq 2\lambda \left(\frac{27U}{B} + 4L^2\right)$ and $\lambda \geq \delta^2/U$, we have:
\begin{equation}
    \frac{\delta}{2} \|T(\hat{h}-h_*)\|_2 \leq O(\delta^2 + \delta \eta_n + \lambda U) + 2\sup_{f\in \mcF} \Psi^{\lambda/2}(h_*, f) +   \left(2\lambda\left( \frac{27U}{B} + 4L^2 \right) + \mu\right) \|h_*\|_{\mcH}^2
\end{equation}
Dividing over by $\delta$ and treating $L, U, B$ as constants, we get:
\begin{equation}
     \|T(\hat{h}-h_*)\|_2 \leq O(\delta + \eta_n + \|h_*\|_{\mcH}^2 \left(\lambda/\delta + \mu/\delta\right)) + \frac{2}{\delta}\sup_{f\in \mcF} \Psi^{\lambda/2}(h_*, f)
\end{equation}
Thus either $\|f_{\hat{h}}\|\leq \delta$ or the latter inequality holds. However, in the case when $\|f_{\hat{h}}\|\leq \delta$, we have by a triangle inequality that: $\|T(\hat{h}-h_*)\|_2\leq \delta + \eta_n$. Thus in any case the latter inequality holds.

\paragraph{Upper bounding population sup-loss at minimum.} Let $f_0 = T(h_0 - h_*)$ and observe that:
\begin{equation}
    \sup_{f\in \mcF} \Psi^{\lambda/2}(h_*, f) = \sup_{f\in \mcF} \E[ f_0(z)\, f(z) ] - \frac{\lambda}{2} \left(\frac{2}{3}\|f\|_{\mcF}^2 + \frac{U}{\delta^2} \|f\|_{2}^2\right) \leq \sup_{f\in \mcF} \E[ f_0(z)\, f(z) ] - \frac{\lambda}{2} \frac{U}{\delta^2} \|f\|_{2}^2
\end{equation}
Then by the Cauchy-Schwarz inequality and since $\lambda \geq \delta^2/U$:
\begin{equation}
    \sup_{f\in \mcF} \E[f_0(z)\, f(z)] - \lambda  \frac{U}{\delta^2} \|f\|_{2}^2 \leq \sup_{f\in \mcF} \|f_0\|_2 \|f\|_2 - \frac{\lambda}{2}  \frac{U}{\delta^2} \|f\|_{2}^2 \leq \frac{\|f_0\|_2^2}{2\lambda U} \delta^2 \leq \frac{\|f_0\|^2}{2}
\end{equation}

\paragraph{Concluding.} Concluding we get that w.p. $1-3\zeta$:
\begin{equation}
    \|T(\hat{h}-h_*)\|_2 \leq O(\delta + \eta_n + \|h_*\|_{\mcH}^2 \left(\lambda/\delta + \mu/\delta\right)) + \frac{\|T(h_*-h_0)\|_2^2}{\delta}
\end{equation}
By a triangle inequality:
\begin{align}
\|T(\hat{h}-h_0)\|_2 \leq~& \|T(\hat{h}-h_*)\|_2 + \|T(h_*-h_0)\|_2\\ \leq~& O(\delta + \eta_n + \|h_*\|_{\mcH}^2 \left(\lambda/\delta + \mu/\delta\right)) + \frac{\|T(h_*-h_0)\|_2^2}{\delta} + \|T(h_*-h_0)\|_2
\end{align}
\end{proof}

\subsection{Proof of \Cref{thm:reg-main-error-2}}

\begin{proof}
By the definition of $\hat{h}$:
\begin{equation}
    0\leq \sup_{f} \Psi_n(\hat{h}, f) \leq \sup_{f} \Psi_n(h_0, f) + \lambda \left(\|h_0\|_{\mcH} - \|\hat{h}\|_{\mcH}\right)
\end{equation}
Let $\mcF_{U}^i=\{f\in \mcF_i: \|f\|_{\mcF}\leq U\}$ and $\delta_{n,\zeta}=\max_{i=1}^d 2\,\mcR(\mcF_U^i) + c_0 \sqrt{\frac{\log(c_1/\zeta)}{n}}$ for some universal constants $c_0,c_1$. By Theorem~26.5 and 26.9 of \cite{shalev2014understanding}, and since $\mcF_U^i$ is a symmetric class and $\sup_{y\in \mcY, x\in \mcX}|y-h_0(x)|\leq 2$, w.p. $1-\zeta$:
\begin{equation}
    f\in \mcF_U^i \left|\Psi_n(h_0, f) - \Psi(h_0, f)\right| \leq \delta_{n,\zeta}
\end{equation}
Since $\Psi(h_0, f)=0$ for all $f$, we have that, w.p. $1-\zeta$:
\begin{equation}
    \|\hat{h}\|_{\mcH} \leq \|h_0\|_{\mcH} + \delta_{n,\zeta}/\lambda
\end{equation}
Let $B_{n,\lambda,\zeta} = (\|h_0\|_{\mcH} + \delta_{n,\zeta}/\lambda)^2$. Then if we let $\epsilon_{n,\lambda, \zeta}=\max_{i}\mcR(\mcH_{B_{n,\lambda,\zeta}}\cdot \mcF_U^i) + c_0 \sqrt{\frac{\log(c_1/\zeta)}{n}}$ for some universal constants $c_0,c_1$.
\begin{equation}
    \forall h\in \mcH_{B_{n,\lambda,\zeta}}, f\in \mcF_U^i \left|\Psi_n(h, f) - \Psi(h, f)\right| \leq \delta_{n,\zeta}
\end{equation}
By a union bound over the $d$ function classes composing $\mcF$, we have that w.p. $1-2\zeta$:
\begin{equation}
    \sup_{f\in \mcF_U} \Psi_n(h_0, f) \leq \sup_{f\in \mcF_U} \Psi(h_0, f) + \delta_{n,\zeta/d} = \delta_{n,\zeta/d}
\end{equation}
and
\begin{equation}
    \sup_{f\in \mcF_U} \Psi_n(\hat{h}, f) \geq \sup_{f\in \mcF_U} \Psi(\hat{h}, f) - \epsilon_{n,\zeta/d}
\end{equation}
Since, by assumption, for any $h\in \mcH_{B_{n,\lambda,\zeta}}$, $\frac{T(h_0-h)}{\|T(h_0-h)\|_2}\in \spanF_{R}(\mcF_U)$, we have $\frac{T(h_0-h)}{\|T(h_0-h)\|_2}=\sum_{i=1}^p w_i f_i$, with $p<\infty$, $\|w\|_1\leq \kappa$ and $f_i\in \mcF_U$. Thus we have:
\begin{align}
    \sup_{f\in \mcF_U} \Psi(\hat{h}, f) \geq~& \frac{1}{\kappa} \sum_{i=1}^p w_i \Psi(\hat{h}, f_i) = \frac{1}{\kappa} \Psi\left(\hat{h}, \sum_i w_i f_i\right)\\
    =~& \frac{1}{\kappa} \frac{1}{\|T(h_0-\hat{h})\|_2}\Psi(\hat{h}, T(h_0 - \hat{h}))\\
    =~& \frac{1}{\kappa} \frac{1}{\|T(h_0-\hat{h})\|_2}\E[T(h_0-\hat{h})(z)^2]\\
    =~& \frac{1}{\kappa} \|T(h_0-\hat{h})\|_2
\end{align}
Combining all the above we have:
\begin{equation}
    \|T(h_0-\hat{h})\|_2 \leq \kappa\, \left(\epsilon_{n, \lambda, \zeta/d} + \delta_{n,\zeta/d} + \lambda \left(\|h_0\|_{\mcH} - \|\hat{h}\|_{\mcH}\right)\right)
\end{equation}
Moreover, since functions in $\mcH$ and $\mcF$ are bounded in $[-1,1]$, we have that the function $h\cdot f$ is $1$-Lipschitz with respect to the vector of functions $(h, f)$. Thus we can apply a vector version of the contraction inequality \cite{maurer2016vector} to get that:
\begin{equation}
\mcR(\mcH_{B_{n,\lambda, z}}\cdot \mcF_U^i) \leq 2\, \left(\mcR(\mcH_{B_{n,\lambda, z}}) + \mcR(\mcF_U^i)\right)
\end{equation}
Finally, we have that since $\mcH$ is star-convex:
\begin{equation}
    \mcR(\mcH_{B_{n,\lambda, z}}) \leq \sqrt{B_{n,\lambda, z}}\,\mcR(\mcH_1)
\end{equation}
Leading the final bound of:
\begin{equation}
    \|T(h_0-\hat{h})\|_2 \leq \kappa \left( 2\left(\|h_0\|_{\mcH} + \delta_{n,\zeta}/\lambda\right) \mcR(\mcH_1) + 2\, \max_{i=1}^d \mcR(\mcF_U^i) + c_0\sqrt{\frac{\log(c_1\, d/\zeta)}{n}} + \lambda \left(\|h_0\|_{\mcH}-\|\hat{h}\|_{\mcH}\right)\right)
\end{equation}
Since $\|h_0\|_{\mcH}\leq R$ and $\lambda \geq \delta_{n,\zeta}$, we get the result.
\end{proof}

\subsection{Proof of Theorem~\ref{thm:main-error-ill}}

The proof is identical to that of Theorem~\ref{thm:reg-main-error} with small modifications. Hence we solely mention these modifications and omit the full proof.

The only part that we change is instead of the set of Equations~\eqref{eqn:diff-crit-radii}, we instead view $\psi(y;h(x))\, f_h(z)$ as a function of the vector valued function $(x, z) \to (h(x), f_h(z))$. Then we note that since $h, f$ take values in $[-1, 1]$ and $y\in [-1, 1]$, we note that this function $2$-Lipschitz with respect to this vector. Then we can apply Lemma~7 of \cite{Foster2019}, and by our choice of $\delta:=\delta_n + c_0 \sqrt{\frac{\log(c_1/\zeta)}{n}}$, where $\delta_n$ upper bounds the critical radius of $\shull(\mcH_B-h_*)$ and $\shull(T(\mcH_B-h_*))$, we have that w.p. $1-\zeta$: $\forall h\in \mcH_B$:
\begin{align}
     \left|(\Psi_n(h, f_{h}) - \Psi_n(h_*, f_{h})) - (\Psi(h, f_h) - \Psi(h_*, f_h))\right| \leq~& 36\delta \left(\|h-h_*\|_2 + \|f_h\|_2\right) + 18\delta^2
\end{align}

Subsequently, we can follow identical steps to conclude that w.p. $1-3\zeta$, either $\|f_{\hat{h}}\|_2\leq \delta$ or:
\begin{equation}
\|T(\hat{h}-h_0)\|_2 \leq O\left(\delta + \delta \frac{\|\hat{h}-h_*\|_2}{\|f_{h}\|_2} + \eta_n + \|h_*\|_{\mcH}^2 \left(\lambda/\delta + \mu/\delta\right) + \frac{\|T(h_*-h_0)\|_2^2}{\delta}\right)
\end{equation}
Subsequently, by the measure of ill-posedness we have:
\begin{equation}
    \|\hat{h}-h_*\|_2 \leq \tau \|T(\hat{h}-h_*)\|_2
\end{equation}
Moreover, observe that when $\|f_{\hat{h}}\|_2 \geq \delta \geq 3\eta_n$, then we have by a triangle inequality that:
\[
\|T(h-h_*)\|_2 \geq \|f_{\hat{h}}\|_2-\eta_n \geq 2\eta_n
\]
and:
\begin{equation}
    \|f_{\hat{h}}\|_2 \geq \|T(h-h_*)\|_2 - \eta_n \geq \frac{1}{2} \|T(h-h_*)\|_2
\end{equation}
Thus we get that:
\begin{equation}
    \frac{\|\hat{h}-h_*\|_2}{\|f_{h}\|_2} \leq \tau \frac{\|T(h-h_*)\|_2}{\|f_{\hat{h}}\|_2} \leq 2\tau
\end{equation}
Thus overall we have that either $\|f_{\hat{h}}\|_2\leq \delta$ or:
\begin{align}
\|\hat{h}-h_*\|_2\leq~&  O\left(\tau\left(\tau\delta + \eta_n + \|h_*\|_{\mcH}^2 \left(\lambda/\delta + \mu/\delta\right) + \frac{\|T(h_*-h_0)\|_2^2}{\delta}\right)\right) \\
\leq~&   O\left(\tau\left(\tau\delta + \eta_n + \|h_*\|_{\mcH}^2 \left(\lambda/\delta + \mu/\delta\right) + \frac{\|h_*-h_0\|_2^2}{\delta}\right)\right) \label{eqn:final-bound-ill-posed}
\end{align}
where the last inequality follows by that fact that Jensen's inequality implies that $\|T(h_*-h_0)\|_2\leq \|h_*-h_0\|_2$. Moreover, if $\|f_{\hat{h}}\|_2\leq \delta$, then by a triangle inequality that $\|T(\hat{h}-h_*)\|_2 \leq \delta + \eta_n$, which, subsquently implies by invoking the bound on the ill-posedness measure that: $\|\hat{h}_*-h\|\leq \tau (\delta + \eta_n)$. Thus in any case the bound in Equation~\eqref{eqn:final-bound-ill-posed} holds.
Choosing $h_*:= \arginf_{h \in \mcH_B} \|h-h_0\|_2$, yields the result.

\section{Proofs from Section~\ref{sec:rkhs} and \Cref{app:rkhs}}

\subsection{Proof of \Cref{rkhs-closed-form-max}}\label{sec:proof-rkhs-closed-form-max}
\begin{proof}
Since $\norm{f}_{2,n}$ depends on $f$ only through the values $f(z_1), \dots, f(z_n)$, and the maximization over $f$ in \cref{eqn:maximization} is the penalized problem
\[
\sup_{f\in \mcF}\,\, \frac{1}{n}\sum_{i=1}^{n} \psi(y_i\midsem{} h(x_i))\,f(z_i) - \lambda (\textfrac{U}{\delta^2}\norm{f}_{2,n}^2 + \Knorm{f}^2)
\]
for some choice of $\lambda \geq 0$, the generalized representer theorem of \citep[Thm. 1]{scholkopf2001generalized} implies that an optimal solution of the constrained problem in \cref{eqn:maximization} takes the form
\[
f^*(z) = \sum_{i=1}^n \alpha_i^* K(z_i, z)
\]
for some weight vector $\alpha^* \in \R^n$.
Now consider a function
\[
f(z) = \sum_{i=1}^n \alpha_i K(z_i, z)
\]
for any $\alpha \in \R^n$.
We have $\Knorm{f}^2 = \alpha^\top  K_n\alpha$,
$f(z_i) = e_i^\top  K_n\alpha$, and
\[
\norm{f}_{2,n}^2
    = \frac{1}{n}\sum_{i=1}^n f(z_i)^2
    = \frac{1}{n}\sum_{i=1}^n \alpha^\top  K_ne_ie_i^\top  K_n\alpha
    = \frac{1}{n} \alpha^\top  K_n^2\alpha.
\]
Thus the penalized problem is equivalent to the finite dimensional maximization problem:
\begin{align}
    \sup_{\alpha \in\R^n}\,\, \psi_n^\top  K_n\alpha - \lambda \alpha^\top \left( \frac{U}{n\delta^2} K_n + I\right) K_n \alpha
\end{align}
by taking the first order condition, the latter has a closed form optimizer of:
\begin{equation}
    \alpha^* = \frac{1}{2\lambda} \left( \frac{U}{n\delta^2} K_n + I\right)^{-1} \psi_n
\end{equation}
and optimal value of:
\begin{align}
    \frac{1}{4\lambda} \psi_n^\top  K_n \left( \frac{U}{n\delta^2} K_n + I\right)^{-1} \psi_n =  \frac{1}{4\lambda} \psi_n^\top  K_n^{1/2} \left( \frac{U}{n\delta^2} K_n + I\right)^{-1} K_n^{1/2}\psi_n
\end{align}
where in the last equality we used a classic matrix inverse identity for kernel matrices.\footnote{The fact that for any matrix $X$: $X (X^\top X + \lambda I)^{-1} = XX^\top (X^\top X + \lambda I)$, and that $K_n = K_n^{1/2} K_n^{1/2}$ and $K_n^{1/2}$ is symmetric.}
\end{proof}

\subsection{Proof of \Cref{prop:closed-form-min}}
\label{sec:proof-closed-form-min}
\begin{proof}
By \cref{rkhs-closed-form-max},
\begin{align}\label{eqn:constrained-rkhs-min-closed-form-max}
    \hat{h} =
    \argmin_{h\in \mcH} \frac{1}{4\lambda} \psi_n^\top  M \psi_n + \mu \KHnorm{h}^2 = \argmin_{h\in \mcH} \psi_n^\top  M \psi_n + 4\lambda\, \mu \KHnorm{h}^2
\end{align}
where $\psi_n = (\frac{1}{n}\psi(y_i\midsem{} h(x_i)))_{i=1}^n$.
Since the objective of  \cref{eqn:constrained-rkhs-min-closed-form-max} depends only on $h$ only through the values $h(x_1), \dots, h(x_n)$, and the problem, the generalized representer theorem of \citep[Thm. 1]{scholkopf2001generalized} implies that an optimal solution of the  problem \cref{eqn:constrained-rkhs-min-closed-form-max} takes the form
\[
h^*(x) = \sum_{i=1}^n \alpha_i^* \KH(x_i, x)
\]
for some weight vector $\alpha^* \in \R^n$.
Now consider a function
\[
h(z) = \sum_{i=1}^n \alpha_i \KH(z_i, z)
\]
for any $\alpha \in \R^n$.
We have $\KHnorm{h}^2 = \alpha^\top  \KHn \alpha$,
$h(z_i) = e_i^\top  \KHn\alpha$, and
$\psi_n = y - \KHn\alpha$.
The problem \cref{eqn:constrained-rkhs-min-closed-form-max} is therefore equivalent to
\begin{align}
    \min_{\alpha \in\R^n}\,\, \alpha^\top  \KHn M \KHn \alpha - 2 y^\top  M \KHn \alpha + 4\lambda\, \mu\, \alpha^\top \KHn \alpha.
\end{align}
By \citep[Ex. 4.22]{boyd2004convex}, this problem is solved by:
\begin{align}
    \alpha^* := \left(\KHn\, M\, \KHn + 4\lambda\, \mu\, \KHn\right)^{\dagger} \KHn M y
\end{align}

\end{proof}

\subsection{Proof of \Cref{lem:ill-posedness-rkhs}}

\begin{proof}
Under these assumptions we have:
\begin{align}
    \|Th\|_2^2 =~& a_I^\top  V_m a_I - 2\sum_{i\leq m < j} a_i a_j \E[\E[e_i(x)\mid z]\E[e_j(x)\mid z]] + \E\left[\left(\sum_{j>m} a_j \E[e_j(x)\mid z]\right)^2\right]\\
    \geq~& a_I^\top  V_m a_I - 2\sum_{i\leq m < j} |a_i a_j| \left|\E[\E[e_i(x)\mid z]\E[e_j(x)\mid z]]\right|\\
    \geq~& a_I^\top  V_m a_I - 2\sum_{i\leq m < j} |a_i a_j| c\,\tau_m\\
    \geq~& \tau_m \|a_I\|_2^2 - 2c\, \tau_m\sum_{i\leq m} |a_i| \sqrt{\sum_{j > m} a_j^2}\\
    \geq~& \tau_m \|a_I\|_2^2 - 2c\, \tau_m \sqrt{\lambda_{m+1} B} \sum_{i\leq m} |a_i| \\
    \geq~& \tau_m \|a_I\|_2^2 - 2c\, \tau_m \sqrt{\lambda_{m+1} B} \sqrt{\sum_{i\leq m} a_i^2}\\
    \geq~& \tau_m \|a_I\|_2^2 - 2c\, \tau_m \sqrt{\lambda_{m+1} B} \|a_I\|_2
\end{align}
Thus if $\|Th\|_2\leq \delta$, then by solving the above quadratic inequality and using the fact that $(a+b)^2\leq 2a^2+2b^2$, we have for all $m$:
\begin{equation}
    \|a_I\|_2^2  \leq \frac{4\delta^2}{\tau_m} + 4 c^2 \lambda_{m+1} B
\end{equation}
Moreover, observe that by the RKHS norm bound:
\begin{equation}
    \|h\|_2^2 = \sum_{j\in J} a_j^2 \leq \|a_I\|_2^2 + \lambda_m B
\end{equation}
Thus we can bound:
\begin{equation}
    \tau^*(\delta)^2 = \min_{h: \|Th\|_2^2\leq \delta^2} \|h\|_2 \leq \min_{m\in \mathbb{N}_+} \frac{4\delta^2}{\tau_m} + (4c^2+1)\lambda_{m+1} B
\end{equation}
\end{proof}

\section{Proofs from Section~\ref{sec:high-dim-linear} and \Cref{app:high-dim-linear}}

\subsection{Proof of \Cref{cor:sparse-linear-reg-ell1}}

\begin{proof}
Let $\mcH=\{\ldot{\theta}{x}: \theta\in \R^d\}$ and $\|h\|_{\mcH}=\|\theta\|_1$. Moreover, suppose that $h_0$ is $s$-sparse. Then if $h\in H_{B_{n,\lambda,\zeta}}$, then:
\begin{align}
    \delta_{n,\zeta}/\lambda + \|\theta_0\|_1 \geq \|\hat{\theta}\|_1 =\|\theta_0 + \nu\|_1 = \|\theta_0 + \nu_S\|_1+\|\nu_{S^c}\|_1 \geq \|\theta_0\|_1 - \|\nu_S\|_1 + \|\nu_{S^c}\|_1
\end{align}
Thus:
\begin{equation}
    \|\nu\|_1 \leq 2 \|\nu_S\|_1 + \delta_{n,\zeta}/\lambda \leq 2\sqrt{s} \|\nu_S\|_2  + \delta_{n,\zeta}/\lambda \leq 2\sqrt{s}\|\nu\|_2  + \delta_{n,\zeta}/\lambda \leq 2 \sqrt{\frac{s}{\gamma} \nu^\top  V \nu}  + \delta_{n,\zeta}/\lambda
\end{equation}
Moreover, observe that:
\begin{equation}
    \|T(h-h_0)\|_2 = \sqrt{\E[ \ldot{\nu}{\E[x\mid z]}^2 ]} = \sqrt{\nu^\top  V \nu}
\end{equation}
Thus we have:
\begin{equation}
    \frac{T(h-h_0)}{\|T(h-h_0)\|_2} = \sum_{i=1}^p \frac{\nu_i}{\sqrt{\nu^\top V\nu}} \E[x_i\mid z]
\end{equation}
Thus we can write $\frac{T(h-h_0)}{\|T(h-h_0)\|_2}$ as $\sum_{i=1}^p w_i f_i$, with $f_i\in \mcF_U$ and:
\begin{equation}
    \|w\|_1 = \frac{\|\nu\|_1}{\sqrt{\nu^\top V \nu}} \leq 2\sqrt{\frac{s}{\gamma}} + \frac{\delta_{n,\zeta}}{\lambda} \frac{1}{\|T(h-h_0)\|_2}.
\end{equation}
Thus: $\frac{T(h-h_0)}{\|T(h-h_0)\|_2} \in \spanF_{\kappa}(\mcF_U)$ for $\kappa=2\sqrt{\frac{s}{\gamma}} + \frac{\delta_{n,\zeta}}{\lambda} \frac{1}{\|T(h-h_0)\|_2}$.

Moreover, observe that by the triangle inequality:
\begin{equation}
    \|h_0\|_{\mcH} - \|\hat{h}\|_{\mcH} = \|\theta_0\|_1 - \|\hat{\theta}\|_1 \leq \|\theta_0-\hat{\theta}\|_1 = \|\nu\|_1 \leq 2 \sqrt{\frac{s}{\gamma} \nu^\top  V \nu}  + \delta_{n,\zeta}/\lambda
\end{equation}
Moreover, by standard results on the Rademacher complexity of linear function classes (see e.g. Lemma~26.11 of \citep{shalev2014understanding}), we have $\mcR(\mcH_B)\leq B\sqrt{\frac{2\log(2\, p)}{n}}\max_{x\in \mcX} \|x\|_{\infty}$ and $\mcR(\mcF_U)\leq U\sqrt{\frac{2\log(2\,p)}{n}}\max_{z\in \mcZ} \|z\|_{\infty}$ for $\mcF_U=\{z \to \ldot{\beta}{z}: \beta\in \R^p, \|\beta\|_1\leq U\}$. Thus invoking \Cref{thm:reg-main-error-2}:
\begin{align}
    \|T(\hat{h}-h_0)\|_2 \leq~& \left(2\sqrt{\frac{s}{\gamma}} + \frac{\delta_{n,\zeta}}{\lambda} \frac{1}{\|T(h-h_0)\|_2}\right)\cdot \left(2 (B+1) \sqrt{\frac{\log(2p)}{n}} + \delta_{n,\zeta} + \lambda \sqrt{\frac{s}{\gamma}} \|T(h-h_0)\|_2\right)
\end{align}
The right hand side is upper bounded by the sum of the following four terms:
\begin{align}
    Q_1 :=~& 2\sqrt{\frac{s}{\gamma}} \left(2(B+1) \sqrt{\frac{\log(2p)}{n}} + \delta_{n,\zeta}\right)\\
    Q_2 :=~& \left(\frac{\delta_{n,\zeta}}{\lambda} \frac{1}{\|T(h-h_0)\|_2}\right)\left(2 (B+1) \sqrt{\frac{\log(2p)}{n}} + \delta_{n,\zeta} \right)\\
    Q_3 :=~& 2 \lambda \frac{s}{\gamma} \|T(h-h_0)\|_2\\
    Q_4 :=~& \delta_{n,\zeta} \sqrt{\frac{s}{\gamma}}
\end{align}
If $\|T(h-h_0)\|_2 \geq \sqrt{\frac{s}{\gamma}} \delta_{n,\zeta}$ and setting $\lambda \leq \frac{\gamma}{8s}$, yields:
\begin{align}
    Q_2 \leq~& 8 \frac{1}{\lambda}\sqrt{\frac{\gamma}{s}}\left(2 (B+1) \sqrt{\frac{\log(2p)}{n}} + \delta_{n,\zeta} \right)\\
    Q_3 \leq~& \frac{1}{4} \|T(h-h_0)\|_2\\
\end{align}
Thus bringing $Q_3$ on the left-hand-side and dividing by $3/4$, we have:
\begin{equation}
    \|T(h-h_0)\|_2 \leq \frac{4}{3} (Q_1 + Q_2 + Q_4) = \frac{4}{3}\max\left\{\sqrt{\frac{s}{\gamma}}, \frac{1}{\lambda} \sqrt{\frac{\gamma}{s}}\right\} \left(20\, (B+1) \sqrt{\frac{\log(2p)}{n}} + 11 \delta_{n,\zeta}\right)
\end{equation}

The result for the case when $\sup_{z\in \mcZ} \|z\|_2\leq R$ and $\mcF_U=\{z\to \ldot{\beta}{z}: \|\beta\|_2\leq U\}$, follows along the exact same lines, but invoking the Lemma 26.10 of \citep{shalev2014understanding}, instead of Lemma 26.11, in order to get that $\mcR(\mcF_U) \leq \frac{U\, R}{\sqrt{n}}$.
\end{proof}

\subsection{Proof of Propositions \ref{prop:sparse-optimization-ell1} and \ref{prop:sparse-optimization-ell2}}

\begin{proposition}
Consider an online linear optimization algorithm over a convex strategy space $S$ and consider the OFTRL algorithm with a $1$-strongly convex regularizer with respect to some norm $\|\cdot\|$ on space $S$:
\begin{equation}
    f_t = \argmin_{f \in S} f^\top \left(\sum_{\tau\leq t} \ell_{\tau} + \ell_t\right) + \frac{1}{\eta} R(f)
\end{equation}
Let $\|\cdot\|_*$ denote the dual norm of $\|\cdot\|$ and $R=\sup_{f\in S} R(f) - \inf_{f\in S} R(f)$. Then for any $f^*\in S$:
\begin{equation}
    \sum_{t=1}^T (f_t-f^*)^\top \ell_t \leq \frac{R}{\eta} + \eta \sum_{t=1}^T \|\ell_t - \ell_{t-1}\|_* - \frac{1}{4\eta} \sum_{t=1}^T \|f_t - f_{t-1}\|^2
\end{equation}
\end{proposition}
\begin{proof}
The proof follows by observing that Proposition~7 in \cite{syrgkanis2015fast} holds verbatim for any convex strategy space $S$ and not necessarily the simplex.
\end{proof}

\begin{proposition}\label{prop:appendix-minimax}
Consider a minimax objective: $\min_{\theta\in \Theta} \max_{w\in W} \ell(\theta, w)$. Suppose that $\Theta, W$ are convex sets and that $\ell(\theta, w)$ is convex in $\theta$ for every $w$ and concave in $\theta$ for any $w$. Let $\|\cdot\|_\Theta$ and $\|\cdot\|_W$ be arbitrary norms in the corresponding spaces. Moreover, suppose that the following Lipschitzness properties are satisfied:
\begin{align}
    \forall \theta\in \Theta, w, w'\in W: \left\|\nabla_{\theta}\ell(\theta, w)  - \nabla_{\theta}\ell(\theta, w')\right\|_{\Theta, *} \leq L \|w-w'\|_W\\
    \forall w\in W, \theta, \theta'\in \Theta: \left\|\nabla_{w}\ell(\theta, w)  - \nabla_{w}\ell(\theta', w)\right\|_{W, *} \leq L \|\theta-\theta'\|_W
\end{align}
where $\|\cdot\|_{\Theta, *}$ and $\|\cdot\|_{W, *}$ correspond to the dual norms of $\|\cdot\|_{\Theta}, \|\cdot\|_W$. Consider the algorithm where at each iteration each player updates their strategy based on:
\begin{align}
    \theta_{t+1} =~& \argmin_{\theta\in \Theta} \theta^\top \left(\sum_{\tau\leq t} \nabla_{\theta}\ell(\theta_\tau, w_\tau) + \nabla_{\theta} \ell(\theta_t, w_t)\right) + \frac{1}{\eta} R_{\min}(\theta)\\
    w_{t+1} =~& \argmax_{w\in W} w^T \left(\sum_{\tau \leq t} \nabla_{w} \ell(\theta_\tau, w_\tau) + \nabla_w \ell(\theta_t, w_t)\right) - \frac{1}{\eta} R_{\max}(w)
\end{align}
such that $R_{\min}$ is $1$-strongly convex in the set $\Theta$ with respect to norm $\|\cdot\|_\Theta$ and $R_{\max}$ is $1$-strongly convex in the set $W$ with respect to norm $\|\cdot\|_W$ and with any step-size $\eta \leq \frac{1}{4L}$. Then the parameters $\bar{\theta} = \frac{1}{T} \sum_{t=1}^T \theta_t$ and $\bar{w}=\frac{1}{T}\sum_{t=1}^T w_t$ correspond to an $\frac{2 R_*}{\eta \cdot T}$-approximate equilibrium and hence $\bar{\theta}$ is a $\frac{4 R_*}{\eta T}$-approximate solution to the minimax objective, where $R$ is defined as:
\begin{equation}
    R_* := \max\left\{ \sup_{\theta\in \Theta} R_{\min}(\theta) - \inf_{\theta\in \Theta} R_{\min}(\theta), \sup_{w\in W} R_{\max}(w)-\inf_{w\in W} R_{\max}(w)\right\}
\end{equation}
\end{proposition}
\begin{proof}
The proposition is essentially a re-statement of Theorem~25 of \cite{syrgkanis2015fast} (which in turn is an adaptation of Lemma~4 of \cite{Rakhlin2013}), specialized to the case of the OFTRL algorithm and to the case of a two-player convex-concave zero-sum game, which implies that the if the sum of regrets of players is at most $\epsilon$, then the pair of average solutions corresponds to an $\epsilon$-equilibrium (see e.g. \cite{FREUND199979} and Lemma~4 of \cite{Rakhlin2013}).
\end{proof}

\paragraph{Proof of \Cref{prop:sparse-optimization-ell1}: $\ell_1$-ball adversary} Let $R_E(x)=\sum_{i=1}^{2p} x_i \log(x_i)$. For the space $\Theta:=\{\rho\in \R^{2p}: \rho \geq 0, \|\rho\|_1\leq B\}$, the entropic regularizer is $\frac{1}{B}$-strongly convex with respect to the $\ell_1$ norm and hence we can set $R_{\min}(\rho)=B\, R_{E}(\rho)$. Similarly, for the space $W:=\{w\in \R^{2p}: w\geq 0, \|w\|_1=1\}$, the entropic regularizer is $1$-strongly convex with respect to the $\ell_1$ norm and thus we can set $R_{\max}(w)=R_E(w)$. For this choice of regularizers, the update rules can be easily verified to have a closed form solution provided in \Cref{prop:sparse-optimization-ell1}, by writing the Lagrangian of each OFTRL optimization problem and invoking strong duality. Further, we can verify the lipschitzness conditions. Since the dual of the $\ell_1$ norm is the $\ell_{\infty}$ norm,  $\nabla_{\rho}\ell(\rho, w) = \E_n[vu^\top] w + \frac{\mu}{W}$ and thus:
\begin{align}
    \left\|\nabla_{\rho}\ell(\rho, w) - \nabla_{\rho}\ell(\rho, w')\right\|_{\infty} =\|\E_n[vu^\top] (w-w')\|_{\infty} \leq \|\E_n[vu^\top]\|_{\infty} \|w-w'\|_1\\
    \left\|\nabla_{w}\ell(\rho, w) - \nabla_{w}\ell(\rho', w)\right\|_{\infty} =\|\E_n[uv^\top] (\rho-\rho')\|_{\infty} \leq \|\E_n[vu^\top]\|_{\infty} \|\rho-\rho'\|_1\\
\end{align}
Thus we have $L=\|\E_n[uv^\top]\|_{\infty}$. Finally, observe that:
\begin{align}
    \sup_{\rho\in \Theta} B\, R_{E}(\rho) - \inf_{\rho\in \Theta} B\, R_E(\rho) =~& B^2 \log(B\vee 1) + B \log(2p)\\
    \sup_{w\in W} R_{E}(w) - \inf_{w\in W} R_E(w) =~& \log(2p)
\end{align}
Thus we can take $R_*=B^2 \log(B\vee 1) + (B+1) \log(2p)$. Thus if we set $\eta = \frac{1}{4\|\E_n[vu^\top]\|_{\infty}}$, then we have that after $T$ iterations, $\bar{\theta}=\bar{\rho}^+-\bar{\rho}^-$ is an $\epsilon(T)$-approximate solution to the minimax problem, with \begin{equation}
\epsilon(T)=16\|\E_n[vu^\top]\|_{\infty} \frac{4B^2 \log(B\vee 1) + (B+1) \log(2p)}{T}.
\end{equation}
Combining all the above with \Cref{prop:appendix-minimax} yields the proof of \Cref{prop:sparse-optimization-ell1}.

\paragraph{Proof of \Cref{prop:sparse-optimization-ell2}: $\ell_2$-ball adversary} For the case when $W:=\{\beta\in\R^{p}: \|\beta\|_2\leq U\}$, then we have that the squared norm regularizer $R_{\max}(\beta)=\frac{1}{2}\|\beta\|_2^2$ is $1$-strongly convex with respect to the $\ell_2$ norm and we can use $\|\cdot\|_W=\|\cdot\|_2$. The choice of $R_{\min}$ is the same as in the case of an $\ell_1$ adversary, as detailed in the previous paragraph. For this choice of regularizers, the update rules can be easily verified to have a closed form solution provided in \Cref{prop:sparse-optimization-ell2}, by writing the Lagrangian of each OFTRL optimization problem and invoking strong duality. Moreover, the Lipschitzness conditions become:
\begin{align}
    \left\|\nabla_{\rho}\ell(\rho, \beta) - \nabla_{\rho}\ell(\rho, \beta')\right\|_{\infty} =~&\|\E_n[vz^\top] (\beta-\beta')\|_{\infty} \leq \|\E_n[vz^\top]\|_{\infty, 2} \|\beta-\beta'\|_2\\
    \left\|\nabla_{\beta}\ell(\rho, \beta) - \nabla_{\beta}\ell(\rho', \beta)\right\|_{2} =~&\|\E_n[zv^\top] (\rho-\rho')\|_{2}

    \leq \|\E_n[zv^\top]\|_{2, \infty} \|\rho-\rho'\|_1
\end{align}
where $\|A\|_{\infty, 2} = \max_{i} \sqrt{\sum_{j} A_{ij}^2}$ and $\|A\|_{2,\infty}=\sqrt{\sum_{i} \max_{j} A_{ij}^2 }$. Thus we can take
\begin{align}
    L=~&\max\left\{\max_i \sqrt{\sum_j \E_n[v_i z_j]^2} + \sqrt{\sum_{i} \max_{j} \E_n[z_i v_j]^2}\right\}\\
    \leq~& \sqrt{\sum_{i} \max_j\E_n[z_i v_j]^2} = \|\E_n[zv^T]\|_{2,\infty}
\end{align}
Finally, we also have that:
\begin{equation}
    \sup_{\beta \in W} R_{\max}(\beta) - \inf_{\beta\in W}R_{\max}(\beta) \leq \frac{1}{2} U^2
\end{equation}
Thus we can take $R_*=B^2 \log(B\vee 1) + B\log(2p) + \frac{1}{2}U^2$. Thus if we set $\eta = \frac{1}{4\|\E_n[zv^\top]\|_{2,\infty}}$, then we have that after $T$ iterations, $\bar{\theta}=\bar{\rho}^+-\bar{\rho}^-$ is an $\epsilon(T)$-approximate solution to the minimax problem, with \begin{equation}
\epsilon(T)=16\|\E_n[zv^\top]\|_{2,\infty} \frac{4B^2 \log(B\vee 1) + B \log(2p) + U^2/2}{T}.
\end{equation}
Combining all the above with \Cref{prop:appendix-minimax} yields the proof of \Cref{prop:sparse-optimization-ell2}.

\section{Proofs from Section \ref{sec:ensemble} and \Cref{app:ensemble}}

\subsection{Proof of \Cref{thm:reduction-algorithm}}

Observe that we can view the minimax problem as the solution to a convex-concave zero-sum game, where the strategy of each player is a vector in an $n$-dimensional space, subject to complex constraints imposed by the corresponding hypothesis. In particular, let $A=\{(f(z_1), \ldots, f(z_n)): f\in \mcF\}$ and $B=\{(h(x_1), \ldots, h(z_n)): h\in \mcH\}$. Then the minimax problem can be phrased as:
\begin{equation}
    \min_{b\in B} \max_{a\in A} \frac{1}{n}\sum_i ((y_i - b_i)\, a_i - a_i^2) = \max_{b\in B} \min_{a\in A} \frac{1}{n} \sum_i (a_i^2 - (y_i - b_i)\, a_i)
\end{equation}
 Moreover, we will denote with $\ell(a, b):=\frac{1}{n}\sum_i (a_i^2 - (y_i - b_i)\, a_i)$, which is a loss that is concave (in fact linear) in $b$ and convex in $a$. Moreover, our assumption on $\mcF$ implies that $A$ is a convex set.

Then the algorithm described in the statement of the theorem corresponds to solving this zero-sum game via the following iterative algorithm: at every period $t=1, \ldots, T$, the adversary chooses a vector $a_t$ based on the the follow the leader (FTL) algorithm, i.e.:
\begin{equation}
    a_t = \argmin_{a\in A} \frac{1}{t-1}\sum_{\tau=1}^{t-1} \ell(a, b_{\tau})
\end{equation}
and the learner chooses $b_t$ by best-responding to the current test function, i.e.:
\begin{equation}
    b_t = \argmax_{b\in B} \ell(a_t, b)
\end{equation}
The equivalent stems from the following two observations: First, for the adversary we can re-write the FTL algorithm by completing the square as:
\begin{align}
   a_t =~& \argmin_{a\in A} \frac{1}{n} \sum_i \frac{1}{t-1}\sum_{\tau=1}^{t-1} (a_i^2 - (y_i - b_{it})\, a_i)\\
   =~& \argmin_{a\in A} \frac{1}{n} \sum_i \left(a_i^2 - \left(y_i - \frac{1}{t-1} \sum_{\tau=1}^{t-1} b_{it}\right) a_i\right)\\
   =~& \argmin_{a\in A} \frac{1}{n} \sum_i \left(a_i^2 - \frac{1}{2}\left(y_i - \frac{1}{t-1} \sum_{\tau=1}^{t-1} b_{it}\right)\right)^2
\end{align}
which then is equivalent to the oracle call described in the statement of the theorem. Second for the learner we have:
\begin{align}
    b_t =~& \argmax_{b\in B} \ell(a_t, b)\\
    =~& \argmax_{b\in B} \frac{1}{n} \sum_i b_i a_{it}\\
    =~& \argmax_{b\in B} \frac{1}{n} \sum_i b_i |a_{it}| \sign(a_{it})\\
    =~& \argmax_{b\in B} \frac{1}{n} \sum_i |a_{it}| \E_{z\sim \text{Bernoulli}(\frac{b_i + 1}{2})}[(2\,z_i-1)\, \sign(a_{it})] \\
    =~& \argmax_{b\in B} \frac{1}{n} \sum_i |a_{it}| \left(\Pr_{z\sim \text{Bernoulli}(\frac{b_i + 1}{2})}[(2\,z_i-1) = \sign(a_{it})] - \Pr_{z\sim \text{Bernoulli}(\frac{b_i + 1}{2})}[(2\,z_i-1) \neq \sign(a_{it})]\right) \\
    =~& \argmax_{b\in B} \frac{1}{n} \sum_i |a_{it}| \left(2\Pr_{z\sim \text{Bernoulli}(\frac{b_i + 1}{2})}[(2\,z_i-1) = \sign(a_{it})] - 1\right) \\
    =~& \argmax_{b\in B} \frac{1}{n} \sum_i |a_{it}| \Pr_{z\sim \text{Bernoulli}(\frac{b_i + 1}{2})}[(2\,z_i-1) = \sign(a_{it})]\\
    =~& \argmax_{b\in B} \frac{1}{n} \sum_i |a_{it}| \Pr_{z\sim \text{Bernoulli}(\frac{b_i + 1}{2})}\left[z_i = \frac{\sign(a_{it})+1}{2}\right] \\
    =~& \argmax_{b\in B} \frac{1}{n} \sum_i |a_{it}| \Pr_{z\sim \text{Bernoulli}(\frac{b_i + 1}{2})}\left[z_i = 1\{a_{it}>0\}\right]
\end{align}
which is exactly the oracle call described in the statement of the theorem.

Thus it remains to show that the vector $\bar{b} = \frac{1}{T} \sum_{t=1}^T b_t$ is a solution to the minimax problem, which would imply that the corresponding ensemble hypothesis $\bar{h} = \frac{1}{T}\sum_{t=1}^T h_t$ is also a solution to the empirical minimax problem.

To achieve this it suffices to show that the FTL algorithm is a no-regret algorithm for the adversary. Then we can invoke classic results on solving zero-sum games via no-regret dynamics \citep{FREUND199979}. Observe that the learner obviously has zero regret as it best-responds at each period. Thus if we show that the FTL algorithm has $\epsilon(T)$-regret after $T$ periods, then $\bar{b}$ is an $\epsilon(T)$-approximate solution to the minimax problem, invoking the results of \citep{FREUND199979}.

Hence, we now focus on the online learning problem that the adversary is facing and show that FTL is a no-regret algorithm with regret rate $\epsilon(T) = \frac{4log(T)}{T}$. We will begin by invoking Lemma 2.1 of \citep{shalev2007convex}, which states that the regret of the FTL algorithm is bounded by:
\begin{equation}
    \epsilon(T) \leq \frac{1}{T}\sum_{t=1}^T (\ell(a_t, b_{t}) - \ell(a_{t+1}, b_{t}))
\end{equation}
Thus it remains to bound the RHS.

Observe that the loss function $\ell(\cdot, b)$ is $\frac{2}{n}$-strongly convex with respect the $\|\cdot\|_2$ norm on the space $A$, since $a^\top \nabla_{aa}^2 \ell(a, b) a = \frac{2}{n}\|a\|_2^2$. Moreover, observe that the loss function $\ell(\cdot, b)$ is also $\frac{4}{\sqrt{n}}$-Lispchitz with respect to the $\|\cdot\|_2$ norm on the space $A$, since
\begin{align}
\nabla_{a_i} \ell(a, b) = \frac{1}{n} (2\, a_i - (y_i - b_i))
\end{align}
and therefore:
\begin{align}
    \|\nabla_a \ell(a, b)\|_2 = \sqrt{\frac{1}{n^2} \sum_i(y_i - b_i - 2\, a_i)^2} = \frac{1}{\sqrt{n}} \sqrt{\frac{1}{n} \sum_i(y_i - b_i - 2\, a_i)^2} \leq \frac{4}{\sqrt{n}}
\end{align}
In the last inequality we used the fact $|y_i|, |h(x_i)|, |f(z_i)|\leq 1$.

Since $\ell_t$ is $\frac{2}{n}$-strongly convex, we have that $L_t=\sum_{\tau=1}^t \ell(\cdot, b_\tau)$ is $\frac{2t}{n}$ strongly convex. Since $a_{t+1}$ is the minimizer of $L_t$ and the set $A$ is a convex set, we have by strong convexity and the first order condition that:
\begin{equation}
    L_t(a_t) \geq L_t(a_{t+1}) + \ldot{a_t - a_{t+1}}{\nabla_{a} L_t(a_{t+1})} + \frac{t}{n} \|a_t - a_{t+1}\|_2^2 \geq L_t(a_{t+1}) +\frac{t}{n} \|a_t - a_{t+1}\|_2^2
\end{equation}
Moreover, since $a_t$ is a minimizer of $L_{t-1}$ and invoking the first order condition, in a similar way as above, we have:
\begin{equation}
    L_{t-1}(a_{t+1}) \geq L_{t-1}(a_t) + \frac{t}{n} \|a_t - a_{t+1}\|_2^2
\end{equation}
Adding the two inequalities and re-arranging we get:
\begin{equation}
    \ell(a_t, b_t) - \ell(a_{t+1}, b_t) \geq \frac{2t}{n} \|a_t - a_{t+1}\|_2^2
\end{equation}
Invoking the lipschitzness of $\ell_t$:
\begin{equation}
    \frac{4}{\sqrt{n}} \|a_t - a_{t+1}\|_2 \geq \ell(a_t, b_t) - \ell(a_{t+1}, b_t) \geq \frac{2t}{n} \|a_t - a_{t+1}\|_2^2
\end{equation}
Thus we have:
\begin{equation}
    \|a_t - a_{t+1}\|_2 \leq \frac{2 \sqrt{n}}{t}
\end{equation}
Moreover, by lipschitzness of $\ell(\cdot, b)$, we have:
\begin{equation}
    \ell(a_t, b_t) - \ell(a_{t+1}, b_t)  \leq \frac{4}{\sqrt{n}} \|a_t - a_{t+1}\|_2 \leq \frac{8}{t}
\end{equation}
Thus we get:
\begin{equation}
    \epsilon(T) \leq \frac{8}{T} \sum_{t=1}^T \frac{1}{t} \leq \frac{8 (\log(T)+1)}{T}
\end{equation}