Estimation and inference for high-dimensional nonparametric additive
instrumental-variables regression

Ziang Niu<EMAIL_ADDRESS>
Applied Mathematics and Computational Science Graduate Group

of Pennsylvania

209 S 33rd St, Philadelphia

Pennsylvania 19104, USA
Yuwen Gu<EMAIL_ADDRESS>
Department of Statistics
University of Connecticut

215 Glenbrook Road, Storrs

Connecticut 06269, USA
Wei LiCorresponding author<EMAIL_ADDRESS>
Center for Applied Statistics and School of Statistics

Renmin University of China

59 Zhongguancun Street

Beijing 100872, China

The method of instrumental variables provides a fundamental and practical tool
for causal inference in many empirical studies where unmeasured confounding
between the treatments and the outcome is present. Modern data such as the
genetical genomics data from these studies are often high-dimensional. The
high-dimensional linear instrumental-variables regression has been considered
in the literature due to its simplicity albeit a true nonlinear relationship
may exist. We propose a more data-driven approach by considering the
nonparametric additive models between the instruments and the treatments while
keeping a linear model between the treatments and the outcome so that the
coefficients therein can directly bear causal interpretation. We provide a
two-stage framework for estimation and inference under this more general
setup. The group lasso regularization is first employed to select optimal
instruments from the high-dimensional additive models, and the outcome
variable is then regressed on the fitted values from the additive models to
identify and estimate important treatment effects. We provide non-asymptotic
analysis of the estimation error of the proposed estimator. A debiasing
procedure is further employed to yield valid inference. Extensive numerical
experiments show that our method can rival or outperform existing approaches
in the literature. We finally analyze the mouse obesity data and discuss new
findings from our method.

Causal inference; Group lasso; High-dimensional inference; Instrumental
variables; Nonparametric additive models.

§ INTRODUCTION

The method of instrumental variables has been extensively used in observational
studies to control for unmeasured confounding. If measurements of the treatments
and the outcome are confounded by unobserved variables, the casual effects of
the endogenous treatments can be identified if instrumental variables are
available. The instrumental variables need to be independent of the unmeasured
confounders and can only affect the outcome indirectly through the treatment

The method originates from early research on structural equations in
econometrics [Wright, 1928, Anderson and Rubin, 1949],

and has become increasingly popular in biostatistics and epidemiology with
Mendelian randomization as one of the most exciting recent applications

[Davey Smith and Ebrahim, 2003, Didelez and Sheehan, 2007, Lin et al., 2015, Burgess et al., 2017].

The general setting of this method involves modeling the instrument-treatment
and treatment-outcome relationships. The classical two-stage least squares approach
assumes linearity of both relationships and is widely used in the

low-dimensional setting. However, in many concurrent studies, data are often
high-dimensional. For example, gene expression data collected to identify
genetic variants associated with complex traits in genome-wide association
studies are usually high-dimensional. Many factors such as unmeasured
environmental conditions may induce spurious associations and distort the true
relationships between the gene expressions and the outcome of interest.
Nevertheless, the random assortment of the genes transferred from parents to
offspring resembles the use of randomization in experiments, and single
nucleotide polymorphisms (SNPs) therefore serve as natural instrumental
variables. The SNPs are also high-dimensional.

Recent developments of the instrumental-variables regression have introduced
regularization as a means to address the high dimensionality issue
[Gautier and Rose, 2011, Belloni et al., 2012, Fan and Liao, 2014, Cheng and Liao, 2015, Belloni et al., 2022].

For example, [Belloni et al., 2012] use the lasso to select optimal
instruments from a large pool when the number of treatments remains fixed or
low-dimensional. Various procedures using different types of regularization have
been proposed thereafter. See [Hansen and Kozbur, 2014] and
[Fan and Zhong, 2018], among others. Linear methods in which the
instruments and treatments are both high-dimensional have also been considered
[Lin et al., 2015, Zhu, 2018, Gold et al., 2020].
[Lin et al., 2015] demonstrate an application of the high-dimensional
linear instrumental-variables regression to genetic genomics. However, nonlinear
effects of the SNPs on the gene expressions are likely to exist as can be seen
from some recent articles that employ different kernel-based procedures to
capture possible nonlinear relationships
[Wang et al., 2015, Zhang and Ghosh, 2017, Zhan et al., 2017, Zhao et al., 2019]. While these
methods keep fully nonparametric forms in linking the gene expressions and SNPs,
they are not very effective when applied to the high-dimensional regime.
[Zhu, 2018] also considers the high-dimensional linear
instrumental-variables regression for peer effect estimation in econometrics.
Specifically, to analyze the effects of peers' output on a firm's production
output using panel data, the Research and Development expenditures of peer firms from a
previous period are treated as potential instrumental variables for the
endogenous treatments. Nevertheless, when the linear relationships are in
question, which likely are, the approach by [Zhu, 2018] may lead to
unignorable bias.

Specification of the outcome equation, either in a parametric or nonparametric
form, is often based on expert knowledge or domain theory. The treatment model,
however, can be more data-driven and should involve nonlinear relationships when
possible to reduce bias [Newey, 1990, Fan and Zhong, 2018]. To
better approximate the treatments using optimal instruments, a general
nonparametric model can be beneficial. A substantial body of the recent
literature on high-dimensional nonparametric estimation focuses on the additive
models <cit.>. To this end,

we consider the high-dimensional additive models to capture the nonlinear
effects of the large number of instrumental variables on the treatments. We keep
the linearity assumption for the outcome model so that its coefficients directly
bear causal interpretations. We allow the dimensions of both the instrumental
variables and the treatments to be larger than the sample size. Similar to the
regularized two-stage framework for the high-dimensional linear
instrumental-variables regression, our proposed procedure consists of a first
stage in which we use the group lasso to select important instruments to best
predict the treatments, and a second stage in which we employ lasso to regress
the outcome on the first-stage predictions to perform variable selection and
estimation. We provide rigorous non-asymptotic analysis of the estimator and
further employ a debiasing procedure to establish valid inference.

In contrast to existing methods in the literature, the present work has the
following favorable features and makes several contributions to the
high-dimensional instrumental-variables regression. Firstly, the proposed
procedure is more data-adaptive which allows possible nonlinear
instrument-treatment relationships under high dimensions. A few recent articles
from the machine learning literature adopt deep learning to better estimate the
instrumental-treatment relationships
[Hartford et al., 2017, Xu et al., 2020]. However, these methods typically
require the dimensions of the instruments and the treatments both be smaller
than the sample size, and are not directly applicable to the setting considered
in the present article. Secondly, for the high-dimensional additive models in
the first stage, we develop a probabilistic bound for the estimation error of
the group lasso estimator. Compared with existing work in this area
<cit.>, we explicitly derive the non-asymptotic
probabilistic bounds of the estimation errors, which may be of independent
interest. Similar probabilistic bounds for the estimation error of the second
stage are also provided. Lastly, we provide statistical inference for the causal
parameters of interest by leveraging the debiasing procedures under high
dimensionality. It is recognized that inference for high-dimensional models is
typically difficult even when endogeneity is not present
[Javanmard and Montanari, 2014, Zhang and Zhang, 2014, van de Geer et al., 2014].
[Gold et al., 2020] consider inference in the high-dimensional linear
instrumental-variables regression to deal with endogeneity. The present work
goes beyond that by establishing valid inference in the more flexible additive
models. Hence, our work enriches the literature on high-dimensional inference
that explicitly handles endogeneity.

§ THE SPARSE ADDITIVE INSTRUMENTAL-VARIABLES MODEL

Suppose we have $n$ independent and identically distributed observations from a
population of interest. Let $y_i$, $x_i$, and $z_i$ denote the $i$th
observations of the outcome, the $p\times 1$ vector of treatment variables, and
the $q\times 1$ vector of instrumental variables, respectively, where

Without loss of generality, assume the $x_i$'s and $y_i$'s are centered.
Consider the following joint modeling framework:

\begin{equation}\label{eq:joint-modeling}
(i=1,\ldots,n;\ \ell=1,\ldots,p),
\end{equation}

where $\eta_i\sim N(0,\sigma^{2}_{0})$ and $\varepsilon_{i\ell}\sim
N(0,\sigma^{2}_{\ell})$. Assume the treatment variables are endogeneous in the
sense that $E(\eta_i\mid x_{i\ell})\neq 0$, and the instrumental variables
satisfy $E(\eta_i\mid z_{ij})=E(\varepsilon_{i\ell}\mid z_{ij})=0$. The
$f_{j\ell}(\cdot)$'s are unknown smooth functions with compact support $[a,b]$,
where $a<b$. To ensure identifiability, assume $E\{f_{j\ell}(z_{ij})\}=0$ for
each $i$, $j$ and $\ell$. This is commonly assumed in the literature on additive
models. We also impose some smoothness conditions on the $f_{j\ell}(\cdot)$'s
and set the function class of consideration to a Hölder space $\mathcal{F}$.

For $j=1,\ldots,q$ and $\ell=1,\ldots,p$, the function $f_{j\ell}$ belongs to
$\mathcal{F}$, where

\begin{align*}
\mathcal{F}=\biggl\{f\colon|f^{(k_0)}(z')-f^{(k_0)}(z)|
\leq C|z'-z|^{\alpha_0},\ z,z'\in[a,b];\ \sup_{z\in[a,b]}|f(z)|\leq C_0\biggr\}
\end{align*}

with $d=k_0+\alpha_0>1.5$ and a universal constant $C_0>0$.

This assumption is common in nonparametric regression <cit.>. Other similar
assumptions such as existence of high-order continuous derivatives are also
widely adopted <cit.>.

In model (<ref>), we are mainly interested in estimating the
average treatment effects, $\beta$. The linear setting has been investigated by
[Lin et al., 2015] and [Zhu, 2018], where
$\sum_{j=1}^{q}f_{j\ell}(z_{ij})=z_{i}^{\T}\gamma_{\ell}$. Here, we relax the
linearity assumption and embrace the more general nonparametric additive form.
Define $J_{\ell}=\{j\colon{}f_{j\ell}\neq0\}$ for $\ell=1,\ldots,p$, and
$\L=\{\ell\colon\beta_{\ell}\neq0\}$. The sparsity assumption for the
high-dimensional additive model entails that $|J_{\ell}|\leq r$ for all $\ell$
and some positive integer $r$, where $|J_{\ell}|$ denotes the cardinality of the
set $J_{\ell}$. Similarly, we assume $s$-sparsity in the second stage with
$|\L|\leq s$, where

$s$ is a positive integer. To rewrite model (<ref>) in matrix
form, let $Y=(y_1,\ldots,y_n)^\T\in \R^n$,
$F_j=(F_{j1},\ldots,F_{jp})\in\R^{n\times p}$ with
$F_{j\ell}=\{f_{j\ell}(z_{1j}),\ldots,f_{j\ell}(z_{nj})\}^{\T}\in\R^{n}$, and
$\varepsilon=(\varepsilon_1,\ldots,\varepsilon_p)\in\R^{n\times p}$ with
\R^{n}$ for $j=1,\ldots,q$ and $\ell=1,\ldots,p$. Then
model (<ref>) can be rewritten as

\begin{align*}%\label{eqn:matrix-model}
Y=X\beta+\eta,\quad X=F+\varepsilon,
\end{align*}

where $F=\sum_{j=1}^q F_j$. To handle high dimensionality and endogeneity, we
allow $p,q\gg n$ and develop a two-stage penalized estimation framework.

In the first stage, we estimate each univariate function $f_{j\ell}$ via the
B-spline approximation. Let $a=\xi_{0}<\xi_{1}<\cdots<\xi_{K}<\xi_{K+1}=b$ be an
equal-distanced partition of $[a,b]$, where $K=[n^{\nu}]$ is a positive integer
for some $0<\nu<0.5$. Let $I_{Kt}=[\xi_{t},\xi_{t+1})$ for $t=0,\ldots,K-1$ and
$I_{KK}=[\xi_{K},\xi_{K+1}]$. Let $\{\phi_k(\cdot)\}_{k=1}^{m}$ be the
normalized B-splines such that each of them is (i) a polynomial function of
degree $L$ on $I_{Kt}$ for $t=0,\ldots,K$, and (ii) up to $(L-1)$ times
continuously differentiable on $[a,b]$, where $L>1$ is an integer and $m=K+L$.

A well-known property of such normalized basis functions is that
$0\leq\phi_k(z)\leq1$ and $\sum_{k=1}^{m}\phi_k(z)=1$ for any $z\in[a,b]$

Given the $z_{ij}$'s, let
$\psi_{kj}(\cdot)=\phi_{k}(\cdot)-n^{-1}\sum_{i=1}^{n}\phi_{k}(z_{ij})$ for
$i=1,\ldots,n$ and $j=1,\ldots,q$. We will denote
$\psi_{k}(\cdot)=\psi_{kj}(\cdot)$ when no confusion arises. Now approximate the
additive functions using

\begin{align}\label{eqn:approximation}
\gamma_{kj\ell}\psi_k(z_{ij})+\varepsilon_{i\ell}
\quad (i=1,\ldots,n;\ \ell=1,\ldots,p).
\end{align}

Let $U=( U_{1},\ldots, U_{q})\in\R^{n\times qm}$, where for each $i$ and $j$,
$U_{j}=( U_{1j},\ldots, U_{nj})^{\T}\in \R^{n\times m}$ and
$U_{ij}=\{\psi_{1}(z_{ij}),\ldots,\psi_{m}(z_{ij})\}^{\T}$. Further define the
parameter matrix $\Gamma=(\gamma_{1},\ldots, \gamma_{p})\in\R^{qm\times p}$,
where for each $j$ and $\ell$,
$\gamma_{\ell}=(\gamma_{1\ell}^\T,\ldots,\gamma_{q\ell}^\T)^{\T}\in \R^{qm}$ and
$\gamma_{j\ell}=(\gamma_{1j\ell},\ldots,\gamma_{mj\ell})^\T\in\R^{m}$. The
approximation in (<ref>) becomes

$X\approx U\Gamma+\varepsilon$.

For each $f_{j\ell}\in\mathcal{F}$, there exists
such that with probability at least $1-2(pqm)^{-2}$, the following holds

\begin{align*}
\sup_{z\in[a,b]}\left|f_{j\ell}(z)-\tilde{f}_{nj\ell}(z)\right|
\leq 2C_Lm^{-d}+2C_0\{\log(pqm)/n\}^{1/2},
\end{align*}

where $\tilde{f}_{nj\ell}(z)=\sum_{k=1}^{m}\bar{\gamma}_{kj\ell}\psi_k(z)$ and

$C_L$ is a universal constant depending only on $L$.

Lemma <ref> characterizes the approximation error of
the centered B-splines $\{\psi_k(\cdot)\}_{k=1}^{m}$ to each $f_{j\ell}(\cdot)$
with corresponding coefficients $\bar\gamma_{j\ell}$. Define
\R^{qm}$.

Lemma <ref> implies that an intermediate step of
recovering $f_{j\ell}$ is to estimate $\bar{\gamma}_{\ell}$ by considering the
following penalized problem:

\begin{align}\label{eq:solution_first}
\widehat{\gamma}_{\ell}=\mathop{\arg\min}\limits_{\gamma_{\ell}\in\mathbb{R}^{qm}}
\bigg\{\frac{1}{2n}\|X_{\ell}- U\gamma_{\ell}\|^{2}_{2}
\quad (\ell=1,\ldots,p),
\end{align}

where $X_{\ell}$ is the $\ell$th column of $X$ and $\lambda_{\ell}\ge0$ is a
tuning parameter. This is a group lasso problem [Yuan and Lin, 2006]

and is motivated by the fact that when $f_{j\ell}=0$, the vector
$\bar\gamma_{j\ell}=0$. After obtaining the predicted treatments
$\widehat{X}=(\widehat{X}_{1},\ldots,\widehat{X}_{p})$ with $\widehat{X}_{\ell}=
U\widehat{\gamma}_{\ell}$, we plug $\widehat{X}$ into the following lasso
problem to estimate $\beta$ in the second stage:

\begin{align}\label{eq:solution_second}
\widehat{\beta}=\mathop{\arg\min}\limits_{\beta\in\mathbb{R}^{p}}
\bigg\{\frac{1}{2n}\|Y-\widehat{X}\beta\|^{2}_{2}+\mu\|\beta\|_{1}\bigg\}
\end{align}

for some tuning parameter $\mu\ge0$. Estimation with high-dimensional predictors
has been a popular research topic in the past two decades. We note that the
above formulation is slightly different from the original lasso problem due to
the observed data being replaced by their estimations from the first stage. This
turns out to be more involved when showing the estimation consistency.

§ NON-ASYMPTOTIC ANALYSIS

We provide an estimation error bound for the first-stage group lasso problem.
Compared with existing results in this area
[Huang et al., 2010, Ravikumar et al., 2009], we make contributions by explicitly
deriving the non-asymptotic probability bound. Based on this bound, we establish
a similar error bound for the parameter of interest in the second-stage lasso

Define $\Sigma_U=E(U^{\T}U/n)$. We make the following assumptions.

Each instrumental variable $z_{ij}$ has a continuous density on $[a,b]$ and
the density is bounded away from zero and infinity.

There exists a universal constant $\rho\in(0,1)$ such that

\begin{align*}
\min\biggl\{\frac{\gamma^{\T}\Sigma_U\gamma}{\|\gamma_{J}\|_2^{2}}\colon
|J|\leq r,\gamma\in\mathbb{R}^{qm}\backslash\{0\},\sum_{j\in J^{c}}
\|\gamma_{j}\|_2\leq 3\sum_{j\in J}
\|\gamma_{j}\|_2\biggr\}\geq\frac{\rho}{m},
\end{align*}

where $J\subset\{1,\ldots,q\}$ is an index set, $J^c$ denotes its complement,
and $\gamma_J=(\gamma_{j\ell}\colon{}j\in J)^{\T}$.

Assumption <ref> is rather standard in the
high-dimensional additive models [Huang et al., 2010, Fan and Zhong, 2018].
Assumption <ref> is often called the group restricted eigenvalue
condition [Lounici et al., 2011, Lv et al., 2018, Lu et al., 2020]. This is a
natural extension of the restricted eigenvalue condition for the standard lasso
and Dantzig selector problems [Bickel et al., 2009]. When the
instrumental variables are independent, $\Sigma_U$ is a block diagonal matrix
with diagonals $\Sigma_j=E( U_j^{\T} U_j/n)$. It is well known that
$\lambda_{\min}(\Sigma_j)\geq c_{*}/m$ [Lian, 2012, Huang et al., 2010]
when each instrumental variable is uniformly distributed, where $c_{*}>0$ is a
constant depending on the smoothness of the B-splines. Denote $a_n=o(b_n)$ if
$\lim_{n\rightarrow\infty}a_n/b_n=0$, $a_n=O(b_n)$ if there exists a positive
constant $C_1$ such that $\limsup_{n\rightarrow\infty}a_n/b_n\leq C_1$, and
$a_n=\Theta(b_n)$ if there exists positive constants $C_2$ and $C_3$ such that
$C_2\leq \liminf_{n\rightarrow\infty}a_n/b_n\leq
\limsup_{n\rightarrow\infty}a_n/b_n\leq C_3$.

Suppose Assumptions <ref>–<ref> hold. There exist
positive constants $c_1$, $c_2$, and $c_3$ such that if

\begin{align*}
\lambda_{\max}=\max_{\ell}\lambda_{\ell}=\max
\biggl[c_1\sigma_{\max}\Bigl\{\frac{\log(pqm)}{n}\Bigr\}^{1/2},
\end{align*}

then for sufficiently large $n$, with probability at least $1-20(pqm)^{-1}$,
the regularized estimator $\widehat\gamma_{\ell}$ in (<ref>)

\begin{align*}
\max_{\ell}\bigg\Vert\sum_{j=1}^{q}F_{j\ell}-U\widehat{\gamma}_{\ell}
\bigg\Vert_2^2\leq \frac{50rmn\lambda_{\max}^2}{\rho},
\quad \max_{\ell}\sum_{j=1}^{q}\|\widehat{\gamma}_{j\ell}
\end{align*}

where $\sigma_{\max}=\max_{\ell}\sigma_{\ell}$,
$m=\Theta\{n^{1/(2d+1)}\}$, and $r^{2}=o[n /\{m^4\log(pqm)\}]$.

The performance of the group lasso estimator depends crucially on the eigen
behavior of the empirical covariance matrix $U^{\T}U/n$.

While it can be shown that the group restricted eigenvalue condition for the
empirical covariance matrix is satisfied under Assumption <ref>, this
does come with a price on the rate of the sparsity level, that is,
$r^{2}=o[n/\{m^4\log(pqm)\}]$. Similar requirements can be found in Corollary 1
of [Raskutti et al., 2010]. In view of the conditions of
Theorem <ref>, it is easy to verify that

Thus, to ensure the consistency of the average in-sample prediction, it is
required that $r^{3}=o\{n^{2d/(2d+1)}/\log(pqm)\}$, while for the estimation
consistency of the coefficients, $r^{4}=o\{n^{(2d-1)/(2d+1)}/\log(pqm)\}$ is
required. This is a more restrictive requirement than that in the standard
lasso, but it is expected due to the unspecified additive functional forms.

In contrast to Theorem 1 of [Huang et al., 2010] that only gives the convergence
rates, our result is completely non-asymptotic. Moreover, our
result allows the sparsity $r$ to grow with the sample size and dimension of the
data while this is not allowed in [Huang et al., 2010]. Guaranteed by Theorem 1 part
(i) of [Huang et al., 2010], we can directly compare the estimation consistency
result obtained here with part (ii) of their theorem, and when $r$ is a fixed
number, it is easy to show they are the same. Other aligned results include
[Ravikumar et al., 2009] and [Lu et al., 2020]. [Ravikumar et al., 2009]
obtain the out-of-sample risk consistency while both the explicit rate and the
in-sample error bound remain unclear. [Lu et al., 2020] consider a
kernel-sieve hybrid estimator and obtain a similar non-asymptotic bound.

To provide an estimation error bound for $\widehat{\beta}$ defined
by (<ref>), we make an extra assumption on the population
covariance matrix $\Sigma_{F}=E(F^\T F/n)$.

There exists a  constant $\kappa>0$ such that

\begin{align*}
\min\bigg\{\frac{\beta^{\T}\Sigma_F\beta}{\|\beta_{\mathcal{L}}\|_2^{2}}
\colon|\mathcal{L}|\leq s,\beta\in\mathbb{R}^{p}\backslash\{0\},
\sum_{\ell\in \mathcal{L}^{c}}|\beta_{\ell}|\leq
3\sum_{\ell\in \mathcal{L}}|\beta_{\ell}|\bigg\}\geq\kappa,
\end{align*}

where $\mathcal{L}\subset\{1,\ldots,p\}$ is an index set, $\mathcal{L}^{c}$
denotes its complement, and

Assumption <ref> is the restricted eigenvalue condition on
$\Sigma_F$ and is useful for deriving the error bounds in the second-stage lasso
problem. This assumption imposes some requirements on the covariance structures
of the treatment matrix $X$ and the noise variables $\varepsilon$. For example,
when $\mathrm{cov}(\varepsilon_{\ell},\varepsilon_{\ell'})=0$ for
$\ell\neq\ell'$ and the minimum eigenvalue of $\Sigma_{X}=E(X^{\T}X/n)$ is
larger than $\max_{\ell}\sigma_\ell$, the above condition immediately holds.

To provide an estimation error bound for $\widehat{\beta}$, we restrict the
parameter space of consideration to an $L_{1}$-ball $\|\beta\|_{1}\leq B$ for
some $B>0$. Similar technique has been frequently used in the literature
<cit.>. This restriction can be further
relaxed to the $L_{\infty}$-ball $\|\beta\|_{\infty}\leq{}B$, but it may lead
to a sacrifice of the convergence rate.

Suppose Assumptions <ref>–<ref> hold. Let the
regularization parameter $\lambda_{\max}$ be chosen as in
Theorem <ref>. Further assume $\lambda_{\max}$
satisfies $560C_0\lambda_{\max}(2rm/\rho)^{1/2}\leq\kappa^{2}/(4rs)$. If we
choose the second-stage regularization parameter as

\begin{align*}
\mu=2r\lambda_{\max}(7\sigma_0+8\sqrt{5}B\sigma_{\max}+30B)(2m/\rho)^{1/2},
\end{align*}

then with probability at least $1-234(pqm)^{-1}$, the estimator
$\widehat\beta$ in (<ref>) satisfies

\begin{align*}
\|\widehat{\beta}-\beta\|_{1}\leq\frac{64}{\kappa^{2}}s\mu,
\quad \|\widehat{X}(\widehat{\beta}-\beta)\|^{2}_{2}
\leq\frac{64}{\kappa^{2}}ns\mu^{2}.
\end{align*}

As far as we know, Theorem <ref> is the first to
present a non-asymptotic estimation error bound for the two-stage additive
model. Straightforward analysis shows that consistency is guaranteed if we take
$\mu^{2}=O\{r^{4}\log(pqm)/n^{2d/(2d+1)}\}$ and
$s^{2}r^{5}=o\{n^{2d/(2d+1)}/\log(pqm)\}$. When $r$ is fixed, we have
$s^{2}=o[n/\{m\log(pqm)\}]$. This almost recovers the sparsity in the
classical lasso setting when $d$ is large enough.

Since the two-stage linear model considered by [Lin et al., 2015] and
[Zhu, 2018] is a special case of our setting, the empirical results in
Section <ref> demonstrate similar performance between the two
models when the true relationship in both stages is linear.

§ INFERENCE

We develop a method to draw inference on each component of the outcome
regression parameters $\beta$. When endogeneity is absent, various methods via
debiasing the penalized estimator have been proposed to conduct valid inference
for high-dimensional models [Zhang and Zhang, 2014, van de Geer et al., 2014, Javanmard and Montanari, 2014]. A recent article by [Gold et al., 2020]
considers the endogeneity issue and adapts the parametric one-step update
procedure to construct confidence intervals for parameters in the
high-dimensional two-stage linear model. We extend the approach therein to draw
inference under the more general setting entailed by
model (<ref>).

A key step in deriving the debiased estimator is to utilize the conditional
moment restriction $E(\eta\mid Z)=0$. This equation entails the orthogonality
condition $E(\Gamma^{\T} U^\T\eta)=0$. Here $\Gamma$ is any fixed coefficient
matrix and will later be set to
$\overline\Gamma=(\bar{\gamma}_{1},\ldots,\bar{\gamma}_{p})^{\T}$. Let
$D=U\overline\Gamma$ and $d_i$ be the $i$th row of $ D$. By
Lemma <ref> and
Theorem <ref>, the estimate
$\widehat{D}=U\widehat\Gamma$ is a good estimate of $D$, where
$\widehat\Gamma=(\widehat\gamma_1,\ldots,\widehat\gamma_p)$. The orthogonality
condition then implies that the empirical counterpart based on the estimate
$\widehat{D}$ is approximately equal to zero, that is,

\begin{align*}
\end{align*}

where $E_n(w_i)=n^{-1}\sum_{i=1}^{n}w_i$ is the expectation with respect to the
empirical measure. The one-step update to the second-stage estimator
$\widehat{\beta}$ can thus be written as

\begin{align*}
\widetilde{\beta}=\widehat{\beta}-\widehat{\Omega}
E_n\{h(y_i,x_i,\widehat{ d}_i;\widehat\beta)\}
\end{align*}

where $\widehat{\Omega}$ is some estimate of $\Omega=\Sigma_F^{-1}$. Similar to
[Gold et al., 2020], we construct the estimator $\widehat{\Omega}$ from a
modification of the constrained $L_1$-minimization approach to sparse precision
matrix estimation proposed by [Cai et al., 2011]. The rows
$\widehat{\theta}_\ell$ of the estimator $\widehat{\Omega}$ are obtained as
solutions to the following program:

\begin{align}\label{eq:opt_precision}
\min_{\theta_\ell\in\mathbb{R}^p}\|\theta_\ell\|_1,\enskip
\text{subject to}\enskip\|\widehat{\Sigma}_F\theta_\ell-e_\ell\|_{\infty}
\leq\upsilon\quad (\ell=1,\ldots,p),
\end{align}

where $e_\ell$ is the $\ell$th canonical basis vector in $\R^{p}$ and
$\upsilon>0$ is the tolerance parameter. The following lemma characterizes a
decomposition of the one-step estimator $\widetilde{\beta}$.

The one-step estimator $\widetilde{\beta}=\widehat{\beta}+\widehat{\Omega}
\widehat{D}^{\T}(Y-X\widehat{\beta})/n$ satisfies
$\sqrt{n}(\widetilde{\beta}-\beta)=\Omega D^\T\eta/\sqrt{n}+\sum_{k=1}^4R_k$,

\begin{align*}
R_1 & =(\widehat{\Omega}-\Omega) D^\T\eta/\sqrt{n},
& R_2 & =\widehat{\Omega}(\widehat{ D}- D)^\T\eta/\sqrt{n},\\
R_3 & =\widehat{\Omega}\widehat{D}^\T(X-\widehat{D})
& R_4 & =\sqrt{n}(\widehat{\Omega}\widehat{\Sigma}_{F}-I)
\end{align*}

Lemma <ref> implies that to establish the asymptotic normality of each
component $\widetilde\beta_{\ell}$, it suffices to make sure each remainder term
$\|R_k\|_{\infty}=o_p(1)$, $k=1,2,3,4$. The $L_1$-bound on
$\widehat{\theta}_\ell-\theta_\ell$, which is needed for controlling the
remainder terms, becomes manageable when the following restriction on the
population precision matrix is imposed.

There exist some positive number $m_{\Omega}$, tolerance $b\in[0,1)$, and
generalized sparsity level $s_{\Omega}$ such that the population precision
matrix $\Omega\in\mathcal{U}(m_{\Omega},b,s_{\Omega})$, where

\begin{align*}
\mathcal{U}(m_{\Omega},b,s_{\Omega})
\colon\|\Omega\|_{1}\leq m_{\Omega};\max_{\ell\in\{1,\ldots,p\}}
\sum_{\ell'=1}^{p}|\theta_{\ell\ell'}|^{b}\leq s_{\Omega}\bigg\}
\end{align*}

and $\|\Omega\|_1=\sup_{\ell}\|\theta_\ell\|_1$.

We assume the event that the rows $\theta_\ell$ of $\Omega$ are feasible for the
minimization program (<ref>) has probability approaching one,
that is, $\pr(\|\Omega\widehat{\Sigma}_F-I\|_{\infty}\leq\upsilon)\rightarrow1$
as $n\rightarrow\infty$. The validity of such a requirement mainly depends on
the choice of the tolerance $\upsilon$. We give a theoretical choice of
$\upsilon$ in Theorem <ref> and provide some rate conditions
under which $\sqrt{n}(\widetilde{\beta}_\ell-\beta_\ell)/\omega_{\ell}$ is
asymptotically normal, where $\omega_\ell^{2}=\sigma_0^{2}\theta_{\ell\ell}$ and
$\theta_{\ell\ell}$ is the $\ell$th diagonal entry of $\Omega$.

Suppose Assumptions <ref>–<ref> and the
conditions of Theorems
<ref>–<ref> hold.
Assume each element $\theta_{\ell\ell}>\vartheta$ for some universal constant
$\vartheta>0$ and let
$\upsilon=36C_0m_{\Omega}\lambda_{\max}r(2rm/\rho)^{1/2}$. If the following
rate conditions hold

\begin{align*}
& \Big[m^{1/2}+\big\{\log(pqm)\big\}^{1/2}\Big]=o(1),\\
& \quad r^{7/2}s\big(m^{2}/n\big)^{1/2}\log(pqm)=o(1),
\end{align*}

then $\|R_k\|_{\infty}=o_{p}(1)$ as $n\rightarrow\infty$, where $k=1,2,3,4$.
Moreover, $\sqrt{n}(\widetilde{\beta}_\ell-\beta_\ell)/\omega_\ell$ converges in
distribution to the standard normal distribution.

The rates in Theorem <ref> are required to make the remainder
terms negligible. When the sparsity parameter $r$ in the first stage is fixed,
the requirement for the sparsity level of the second stage is
$s=o(n^{1/2}m/\log(pqm))$. This is almost the same as the requirement for the
debiased lasso: $s=o(n^{1/2}/\log(p))$ <cit.>. By
Theorem <ref>, we can construct a confidence interval for
$\beta_{\ell}$ if a consistent estimator $\widehat{\omega}_{\ell}$ is available.
Given $\ell\in\{1,\ldots,p\}$ and $\alpha\in(0,1)$, an asymptotic
$100(1-\alpha)\%$ confidence interval for $\beta_{\ell}$ is
\widetilde{\beta}_\ell+z_\alpha\widehat{\omega}_\ell/\sqrt{n}]$, where
$z_\alpha=\Phi^{-1}(1-\alpha/2)$ and $\Phi$ is the cumulative distribution
function of the standard normal distribution. The following theorem provides a
way to construct a consistent estimator $\widehat{\omega}_\ell$.

Suppose the conditions of Theorems
<ref>–<ref> hold. Define

\begin{align*}
\widehat{\omega}_\ell=\widehat{\sigma}_0\big(\widehat{\theta}_\ell^{\T}
\widehat{\Gamma}^{\T}U^{\T}U\widehat{\Gamma}
\widehat{\theta}_\ell/n\big)^{1/2},\enskip
\widehat{\sigma}_0=n^{-1/2}\|Y-X\widehat{\beta}\|_2.
\end{align*}

Then $\widehat\omega_{\ell}$ is a consistent estimator of $\omega_\ell$ for each

§ SIMULATION

We conduct simulation studies to evaluate the finite-sample performance of the
proposed methods. Our objective is to test both the estimation and inferential
procedures under a variety of experiments. For estimation purpose, we compare
our procedure with the classical one-stage  penalized least squares (PLS)
and the two-stage least-square method with lasso penalty (2SLS-L). For
inferential purpose, we compare our method with the updated two-stage lasso estimator (Up-2SLS-L)

proposed by
[Gold et al., 2020].

We first investigate the estimation performance under two design settings where
endogeneous treatments are generated from linear and nonlinear models,
respectively. In both settings, we take $p=q=600$ and vary $n$ from $100$ to
$2100$. Experiments for other values of $(p,q,n)$ are provided in the
Appendix. We generate the instrumental variables $z_i$ of the
$i$th observation from the multivariate normal distribution with zero mean and
covariance matrix $\Sigma_{Z}=\{(\Sigma_Z)_{jj'}\}$, where
$(\Sigma_Z)_{jj'}=0.2^{|j-j'|}$ for $j,j'=1,\ldots,q$. To generate the noise
vector $(\eta_i, \v_{i1},\ldots, \v_{ip})^{\T}$, we sample from another normal
distribution with zero mean and covariance matrix $\Sigma=(\Sigma_{\ell\ell'})$,
where $\Sigma_{\ell\ell'} =0.2^{|\ell-\ell'|}$ for $\ell,\ell'=2,\ldots,p+1$ and
$\Sigma_{11}=1$. We also set $\Sigma_{1\ell}$ for $\ell=2,\ldots,6$ and five
other random selected entries from the first column of $\Sigma$ to $0.3$. All
other entries are set to zero. We finally set $\Sigma_{1\ell}=\Sigma_{\ell1}$
for $\ell=2,\ldots,p+1$ to make $\Sigma$ symmetric. Note that the nonzero
$\Sigma_{1\ell}$'s induce endogeneity in the data. In the linear setting, we
generate the treatment variables $x_i$ according to $x_i=\Gamma^\T z_i+\v_i$,
where $\Gamma=(\gamma_{j\ell})\in\R^{q\times p}$ is a sparse coefficient matrix
obtained by sampling $r=5$ nonzero entries of each column from the uniform
distribution $U(0.75,1)$. In the nonlinear setting, we generate $x_i$ from the
following equation:

\begin{align*}
x_{i\ell} = \gamma_{1\ell}z_{i1}^{2} + \gamma_{2\ell}z_{i2}
+ \gamma_{3\ell}z_{i3}^{2}+\gamma_{4\ell}\sin(\pi z_{i4})
+\gamma_{5\ell}z_{i5}^{2}+\epsilon_{i\ell},\quad (\ell=1,\ldots,p).
\end{align*}

The sampling strategy for the coefficients $\gamma_{j\ell}$'s are the same as
that in the linear setting. We finally generate the outcome response according
to $y_i=x_i^\T\beta+\eta_i$, where the coefficient vector $\beta$ is generated
by sampling $s=5$ nonzero components from the uniform distribution over two
disjoint intervals $U\{(-1,-0.75)\cup(0.75,1)\}$.

In all simulations, we use the Bayesian information criterion to select the
first-stage tuning parameters $m$ and $\lambda_\ell$, and the five-fold cross
validation to select the second-stage regularization parameter $\mu$. We report
the $L_1$ error $\Vert\widehat\beta-\beta\Vert_1$ of each method based on one
hundred replications. The results are summarized in
Table <ref>. It can be seen that when the
instrument-treatment relationship is linear, our method is as good as 2SLS-L. As
the sample size increases, the $L_1$ errors of both 2SLS-L and our method decrease,
whereas that of PLS increases due to ignorance of endogeneity. When the
instrument-treatment relationship is nonlinear, the performance of our method is
similar to that in the linear setting. The $L_1$ error of our method is the
smallest in almost all settings and exhibits a decreasing trend as the sample
size increases. In contrast, neither PLS nor 2SLS-L shares such trend in their
performance under these settings. We can clearly see that the 2SLS-L method induces a large bias since it suffers losses from model misspecification. It is also implied that model misspecification
may have a heavier impact on 2SLS-L than on PLS. This may be related to the specific setting of
our simulations.
Overall, the results from
Table <ref> demonstrate the estimation consistency of our

$L_1$ estimation loss of each method averaged over 100 replications with standard deviation shown in parentheses for $p=600$.

2*Sample     3cLinear     3cNonlinear(lr)2-4(lr)5-7
Our method     2SLS-L     PLS     Our method     2SLS-L     PLS

100     1.26 (0.53)    2.52 (1.19)    0.86 (0.22)     2.96 (1.41)    1.19 (0.38)    1.25 (0.38)

300    0.50 (0.23)    0.59 (0.30)    0.51 (0.16)     0.74 (0.28)    1.27 (0.79)     0.79 (0.27)

600    0.34 (0.15)     0.27 (0.11)    0.46 (0.17)    0.43 (0.17)     1.73 (1.26)    0.89 (0.27)

900    0.25 (0.10)    0.23 (0.08)    0.52 (0.14)    0.32 (0.13)    2.90 (1.98)    1.10 (0.23)

1200   0.19 (0.09)     0.19 (0.08)    0.61 (0.18)    0.27 (0.13)    3.35 (2.36)    1.18 (0.18)

1500    0.17 (0.08)    0.18 (0.08)    0.70 (0.25)    0.24 (0.10)    4.17 (3.29)    1.27 (0.16)

1800    0.16 (0.07)    0.17 (0.08)    0.81 (0.29)    0.21 (0.09)    4.74 (4.03)    1.34 (0.17)

2100    0.14 (0.05)     0.16 (0.07)    1.09 (0.42)    0.21 (0.11)    5.61 (4.38)    1.43 (0.16)

Next, we evaluate the performance of the inferential procedure proposed in
Section <ref> and consider a more challenging nonlinear setting:

\begin{align*}
x_{i\ell}=-8\gamma_{1\ell}z_{i1}^{2}+\gamma_{2\ell}\sin(\pi z_{i2})
+\gamma_{5\ell}z_{i5}^{2}+\epsilon_{i\ell},~~ (\ell=1,\ldots,p).
\end{align*}

We calculate the 95% confidence interval of each element of $\beta$ based on
our method and the Up-2SLS-L based inferential procedure. The coverage probabilities
and interval lengths averaged over all elements are shown in
Table <ref>. Our confidence intervals have coverage
probabilities close to the nominal level of 0.95. In contrast, the
Up-2SLS-L based confidence intervals have coverage probabilities well
below the nominal level. Moreover, their intervals are much wider than ours
under all settings. These are expected as the Up-2SLS-L based inferential procedure
is proposed under the linear setting and may not
perform well under the nonlinear setting. The results from
Table <ref> validate our theory in
Section <ref>.

Coverage probabilities and lengths of the 95% confidence intervals
by our method and the Up-2SLS-L based inferential procedure. Numbers shown are
multiplied by one hundred.

2*Dimension     Sample
2cOur method
size     Coverage     Length     Coverage      Length

250   200     92.0    0.396     87.3    2.663

400    300     93.5    0.264    89.0    1.585

500   400    94.0    0.245    87.1    2.102

600   500    93.7    0.140    87.8    2.614

§ REAL DATA ANALYSIS

We further illustrate our proposed method by analyzing the mouse obesity data
described by [Wang et al., 2006]. This data set consists of genotype, gene
expression and clinical information about the F2 intercross mice. The genotypes
are characterized by the SNPs at an average density of 1.5 cM across the whole
genome and the gene expressions in the liver tissues of the mice are profiled by
microarrays. We are interested in the causal effect of the gene expressions on
the body weights of the mice. We
consider the data collected from $n=287$ (144 female and 143 male) mice with
$q=1250$ SNPs and $p=2816$ genes.

Since there are only three genotypes, a sparse high-dimensional linear model
between the SNPs and the gene expressions is often postulated. However, some
recent articles suspect that nonlinear effects may exist in the data
[Li et al., 2014, Zhang and Ghosh, 2017, Guha Majumdar et al., 2020]. This has motivated us to
consider the nonlinear setting. Before applying our method to the data, we
adjust for confounding induced by the sex of the mice. We first regress the body
weight on the sex and subtract the estimated effect from the body weight. We
then apply our proposed estimation method to the data with the adjusted body
weight as the outcome. We use the five-fold cross validation to select the
tuning parameters. The resultant model includes 28 genes. To increase the
stability and interpretability of our analysis, we apply the stability selection
approach proposed by [Meinshausen and Bühlmann, 2010] to compute the selection
probability of each gene over one hundred subsamples of size $\lfloor
n/2\rfloor$ for a sequence of values of the tuning parameter $\mu$. We set the threshold probability to 0.5 so that the number of genes selected is reasonable. The
results are shown in Table <ref>.

Stability selection: genes with selection probability exceeding 0.5.

Gene Name    Selection Probability     Gene Name     Selection Probability

Vwf     0.77     Krtap19-2     0.59

Akap12     0.63     Tmem184c     0.74

2010002N04Rik     0.84     Igfbp2     0.51

Slc43a1    0.76     Gstm2     0.91

Ccnl2     0.54     D14Abb1e    0.52

B4galnt4    0.71

We observe that
there are five genes: Igfbp2, Gstm2, Vwf, 2010002N04Rik, and Ccnl2 that are also selected in [Lin et al., 2015].
These overlapping genes are highly likely to be connected to obesity from a
biological point of view. This  indicates the effectiveness and stability of
our method in finding the risky genes.
Table <ref> also shows
our method identifies some other genes that were not previously found, possibly due to the nonparametric form we consider
in the first stage. In fact, some genes selected with high probability by our
method have been verified by many biological studies. In particular, Solute
Carrier Family 43 Member 1 (Slc43a1) is a Protein Coding gene which can encode
the amino acid transporters that are known to regulate the transmembrane
transport of phenylalanine. [Gill et al., 2010] found that the expression
of Slc43a1 in the fat mice group is quite different from that in the lean mice
group so Slc43a1 is a potential factor leading to obesity.

We also apply our inferential procedure to the data to quantify the uncertainty
associated with our estimation. We use the R package Flare to obtain the optimal
precision matrix from solving (<ref>). Then we construct the
confidence intervals based on the debiased estimator. Table <ref>
presents the causal effects of the genes on the body weight whose corresponding
confidence intervals do not contain zero. Note that the confidence intervals are
wide due to possibly low signal-to-noise ratio and small sample size of the
data. The result is generally consistent with that obtained by the stable
selection. The confidence intervals that are far away from zero include Vwf,
2010002N04Rik, Gstm2, Gp1d1, Slc43a1, Igfbp2, which are also shown to have high
selection probability in both Table <ref>. Moreover, several other genes from
Table <ref> that are shown to have significant causal effects on the
body weight are newly found and have been confirmed to have close biological
relation with obesity. For example, a recent study in [Wang et al., 2019] shows
that by silencing Anxa2, the obesity-induced insulin resistance is attenuated
and our result confirms such a positive relation. Cyp4f15 genes are known to
control the omega-hydroxylated fatty acids in the liver tissue and such acids can be
used for energy production [Hardwick et al., 2009]. Another
independent study shows that the downregulation of Cyp4f15 happens in the liver
tissue among the group of mice fed with high-fat diet [Gai et al., 2020].
The negative causal effect obtained in our result coincides with these findings.

$95\%$ confidence intervals for the causal effects of the genes on
the body weights of the mice. Shown are only the genes whose corresponding
intervals do not contain zero.

Gene Name     Confidence Interval      Gene Name     Confidence Interval

Anxa5             (0.010, 7.269)     Kif22       (0.615, 7.930)

Vwf               (0.500, 7.841)     Gstm2       (0.537, 8.231)

Aqp8              (0.066, 6.855)     Gpld1       ($-$7.448, $-$0.447)

Lamc1             (0.094, 5.877)     Slc43a1     ($-$6.641, $-$1.412)

Acot9             (0.056, 8.298)     Abca8a      ($-$7.152, $-$0.072)

Anxa2             (1.086, 9.331)     Cyp4f15     ($-$7.468, $-$0.250)

2010002N04Rik     (1.343, 8.240)     Igfbp2      ($-$6.451, $-$0.666)

Msr1              (0.004, 6.783)

§ DISCUSSION

Motivated by the data-driven modeling spirit, we develop a high-dimensional additive instrumental-variables regression method with a sound non-asymptotic analysis and valid inference procedure when both instruments and treatments are allowed to be high-dimensional.
There are a lot of directions that are worth further exploration. Firstly, while we estimate the nonparametric functions separately in the first-stage, it would be helpful to take into account the correlations among the
treatments and borrow information across those regressions.
Secondly, although it is convenient for causal interpretation when considering the linear outcome model, many applications in econometrics  and biostatistics have considered nonparametric method to model the relationships between treatments and outcome.
One may further relax our second-stage model setting by considering a
high-dimensional single-index or nonparametric additive outcome model [Radchenko, 2015, Huang et al., 2010].
Thirdly, it is of great interest to extend our results to handling other types of outcome data, for instance, binary or survival outcome. We plan to pursue
these and other related issues in future research.

We are grateful to Prof. Wei Lin for sharing the data. This research is
supported by the National Natural Science Foundation of China (Grant Nos.
12101607, 12071015), the Fundamental Research Funds for the Central Universities
and the Research Funds of Renmin University of China.

§ APPENDIX

We first define the notation that will be used throughout the Appendix.
Suppose $A\in\mathbb{R}^{m\times n}$. Let $\|A\|_2$ be the largest singular
value of $A$ and $\|A\|_{\infty}$ the largest absolute value of the entries of
$A$. Also, let $\|A\|_{1}=\max_{j\in [n]}\|A_{j}\|_1$, where $A_j$ is the $j$th
column of $A$ and $[n]=\{1,\ldots,n\}$. Denote by $\lambda_{\min}(A)$ and
$\lambda_{\max}(A)$ the minimum and maximum eigenvalues of $A$, respectively,
when $A$ is symmetric.

§ APPENDIX A: SOME USEFUL LEMMAS

In this appendix present some useful lemmas that will be used for proving the main results.

Suppose random variables $X_i,i=1,\ldots,n$ are independent, and $X_i$ has
mean $\mu_i$ and sub-Gaussian parameter $\sigma_i$. Then for all $t\geq0$, we

\begin{align*}
\pr\bigg\{\bigg|\frac{1}{n}\sum_{i=1}^{n}(X_i-\mu_i)\bigg|\geq t\bigg\}
\leq2\exp\bigg(-\frac{n^{2}t^{2}}{2\sum_{i=1}^{n}\sigma_i^{2}}\bigg).
\end{align*}

Specifically, if $X_i$ is bounded in $[a,b]$, we have

\begin{align*}
\pr\bigg\{\bigg|\frac{1}{n}\sum_{i=1}^{n}(X_i-\mu_i)\bigg|\geq
\end{align*}

Lemma <ref> follows from Lemma 5 of
[Stone, 1985], which is a well-known result on B-spline approximation
and has been frequently used in the additive models.

For each $f_{j\ell}\in\mathcal{F}$, there exists
$\Bar{f}_{nj\ell}=\sum_{k=1}^{m}\Bar{\gamma}_{kj\ell}\phi_k$ such that

$$\sup_{z\in[a,b]}|f_{j\ell}(z)-\Bar{f}_{nj\ell}(z)|\leq C_{L}m^{-d},$$

where $C_L>0$ is a constant depending only on the degree $L$ of the B-splines.

The next lemma shows the minimum and maximum eigenvalues of the spline matrix
have the same order. Denote the non-centered B-spline matrix by
$\tU=(\tU_1,\ldots,\tU_q)\in\mathbb{R}^{n\times{}qm},$ where
$\tU_j=(\tU_{1j},\ldots,\tU_{nj})^{\T}\in\mathbb{R}^{n\times m}$ and
$\tU_{ij}=\{\phi_{1}(Z_{ij}),\ldots,\phi_{m}(Z_{ij})\}^{\T} \in\mathbb{R}^{m}$.

For $j=1,\ldots,q$, we have

\begin{align*}
\frac{3c_*}{m}-2\|\mathbb{P}_j-\mathbb{P}^n_j\|_{\infty}
\leq\lambda_{\min}\Bigl(\frac{\tU_j^\T\tU_j}{n}\Bigr)
\leq\lambda_{\max}\Bigl(\frac{\tU_j^\T\tU_j}{n}\Bigr)
\leq \frac{c^*}{2m}+2\|\mathbb{P}_j-\mathbb{P}^n_j\|_{\infty},
\end{align*}

where $0<c_{*}<1<c^{*}$ and
z)-n^{-1}\sum_{i=1}^{n}\bm{1}\{Z_{ij}\leq z\}\bigr|$ with $\mathbb{P}_j$ and
$\mathbb{P}^{n}_j$ being the population and empirical distributions of
$Z_{j}$, respectively.

The following result gives a bound on $\|\mathbb{P}_j-\mathbb{P}^n_j\|_{\infty}$.

Let $\mathbb{P}$ be the distribution of a random variable $X$ and
$\mathbb{P}_n$ the empirical distribution based on $n$ i.i.d. copies,
$X_{1},\ldots,X_{n}$, of $X$. Then, we have

\begin{align*}
\left\|\mathbb{P}_n-\mathbb{P}\right\|_\infty\leq 8\{\log(n+1)/n\}^{1/2}+\delta
\end{align*}

with probability at least $1-\exp(-n\delta^{2}/2)$ for any $\delta\geq0$.

Combining Lemmas <ref> and <ref>, we
obtain the following result on the B-spline matrix.

If $8\{\log(n+1)/n\}^{1/2}+\{2/m^{3}\}^{1/2}\leq\min\{c_*/m,c^*/4m\}$, then
for each $j=1,\ldots,q$, with probability at least $1-\exp(-nm^{-3})$, we

\begin{align*}
\frac{c_*}{m}\leq\lambda_{\min}\Bigl(\frac{\tU_j^\T\tU_j}{n}\Bigr)
\leq\lambda_{\max}\Bigl(\frac{\tU_j^\T\tU_j}{n}\Bigr)\leq\frac{c^*}{m}.
\end{align*}

Apply the Glivenko–Cantelli Theorem (Lemma <ref>) with

Let $X_1,\ldots, X_n$ be i.i.d. random variables such that
$\left|X_i\right|\leq b$ almost surely for some constant $b>0$. Then for any

\begin{align*}
\pr\biggl(\bigg|\frac{1}{n}\sum_{i=1}^{n}\big\{X_i-\mathbb{E}(X_i)\big\}
\bigg|\leq\delta\biggr)\geq1-2\exp[-n\delta^{2}/\{2(b^{2}+\delta/3)\}].
\end{align*}

Suppose $X_1,\ldots, X_n$ are i.i.d. standard normal random variables. Then
for $t\geq0$,

\begin{align*}
\pr\bigg(\bigg|\frac{1}{n}\sum_{i=1}^{n}X_i^{2}-1\bigg|
\geq t\bigg)\leq 2\exp\big(-nt^{2}/8\big).
\end{align*}

See Example 2.11 of [Wainwright, 2019].

We present a key concentration result which will be used in the following

Let $\xi_{1},\ldots,\xi_{N}$ be i.i.d. standard normal random variables.
Moreover, let $v=(v_{1},\ldots,v_{N})\neq0,$
$\eta_{v}=\sum_{i=1}^{N}(\xi_{i}^{2}-1)v_{i}/(\sqrt{2}\|v\|_{2})$ and
$m(v)=\|v\|_{\infty}/\|v\|_{2}$. Then for all $x>0$, we have

\begin{align*}
\pr(|\eta_{v}|>x)\leq2\exp\bigg(-\frac{x^{2}}
\end{align*}

The following lemma validates the restricted eigenvalue condition given
Assumption <ref>.

Suppose Assumption <ref> holds. If
$160r\{m^{4}\log(pqm)/n\}^{1/2}\leq\rho/2$, with probability at least

\begin{align*}
\min\bigg\{\frac{\gama^{\T}\U^{\T}\U\gama}{n\|\gama_{J}\|_2^{2}}\colon
|J|\leq r,\gama\in\mathbb{R}^{qm}\backslash\{0\},
\sum_{j\in J^{c}}\|\gama_{j}\|_2\leq 3\sum_{j\in J}\|\gama_{j}\|_2
\bigg\}\geq\frac{\rho}{2m}.
\end{align*}

Note that

\begin{align*}
\Sigma-\frac{\U^{\T}\U}{n}
\frac{\big\{\tU-\mathbb{E}(\tU)\big\}^{\T}
\big\{\tU-\mathbb{E}(\tU)\big\}}{n}\Bigg]
\big\{\tU-\mathbb{E}(\tU)\big\}}{n}\\
\end{align*}

Let $\tilde{\Sigma}=\mathbb{E}[\{\tU-\mathbb{E}(\tU)\}^{\T}
\{\tU-\mathbb{E}(\tU)\}/n]$. Then, we have

\begin{align*}
\left\|\Sigma-\frac{\U^{\T}\U}{n}\right\|_{\infty}
\mathbb{E}(\tU)}{n}-\frac{\tU^{\T}
\mathbb{E}(\tU)}{n}\right\|_{\infty}\\
\eqqcolon{}T_1+T_2+T_3.
\end{align*}

For $T_2$, we can write

\begin{align*}
\mathbb{E}\{\phi_{k'}(Z_{ij'})\}-\sum_{i=1}^{n}
\frac{\phi_{k}(Z_{ij})}{n}\mathbb{E}\{\phi_{k'}(Z_{ij'})\}
\bigg|,\\
&\qquad\qquad \forall j,j'\in\{1,\ldots,q\},\ k,k'\in\{1,\ldots,m\},
\end{align*}

where $\tU_i^{\T}$ is the $i$th row of $\tU$. By the property of the B-spline
matrix, we have $0\leq\phi_k(z)\leq 1$ and
$0\leq\mathbb{E}\{\phi_{k}(Z_{ij})\}\leq1$. Therefore,

\begin{align*}
\leq\max_{j,k}\bigg|\mathbb{E}\{\phi_{k}(Z_{ij})\}
\end{align*}

Now apply Lemma <ref> with $t=\{4\log(pqm_n)/n\}^{1/2}$ and take
the union bound to obtain

\begin{align*}
\leq\bigg\{\frac{4\log(pqm_n)}{n}\bigg\}^{1/2},
\end{align*}

which holds with probability at least $1-2(qm)(pqm)^{-2}$. For $T_3$, the same
bound can be obtained with the same probability since

\begin{align*}
\phi_{k}(Z_{ij})\right\}\left\{\sum_{i=1}^{n}\phi_{k'}
\phi_{k}(Z_{ij})\right\}\mathbb{E}\{\phi_{k'}(Z_{ij'})\}}{n}\bigg|\\
\biggr|\cdot\biggl|\frac{\sum_{i=1}^{n}\phi_{k'}(Z_{ij'})}
&\leq \max_{j',k'}\bigg|\frac{\sum_{i=1}^{n}\phi_{k'}
\end{align*}

Now we bound $T_1$

\begin{align*}
\phi_{k'}(Z_{ij'})-\mathbb{E}\big\{\phi_{k}(Z_{ij})\big\}
\mathbb{E}\big\{\phi_{k'}(Z_{ij'})\big\}\big]\\
\end{align*}

where the last inequality holds because each entry of $\tilde{\Sigma}$ can be

\begin{align*}
\leq 2.
\end{align*}

Note $|\phi_{k}(Z_{ij})-\mathbb{E}\{\phi_{k}(Z_{ij})\}|\leq1,\forall
j\in\{1,\ldots,q\},k\in\{1,\ldots,m\}$, so we have

\begin{align*}
\big|\big[\phi_{k}(Z_{ij})-\mathbb{E}\{\phi_{k}(Z_{ij})\}\big]
\big[\phi_{k'}(Z_{ij'})-\mathbb{E}\{\phi_{k'}(Z_{ij'})\}\big]\big|
&\leq 1.
\end{align*}

Applying Lemma <ref> again with $t=\{8\log(pqm)/n\}^{1/2}$, we
obtain by the union bound argument over $j,j'\in\{1,\ldots,

\begin{align}\label{eq:bound-t123}
T_1\leq \{8\log(pqm)/n\}^{1/2}+2/n\leq 6\{\log(pqm)/n\}^{1/2},
\end{align}

which holds with probability at least $1-2(qm)^{2}(pqm)^{-4}$. Thus, we have

\begin{align*}
\bigg\|\Sigma-\frac{\U^{\T}\U}{n}\bigg\|_{\infty}
\leq T_1+T_2+T_3\leq 10\bigg\{\frac{\log(pqm)}{n}\bigg\}^{1/2},
\end{align*}

which holds with probability at least $1-6(qm)(pqm)^{-2}$. Consider $\gama$
such that $\sum_{j\in J^{c}}\|\gama_j\|_2\leq 3\sum_{j\in J}\|\gama_j\|_2$,
then by Assumption <ref>, we have

\begin{align*}
& \geq\frac{\rho\|\gama_J\|_2^{2}}{m}-\|\gama\|_1^{2}
\bigg\|\Sigma-\frac{\U^{\T}\U}{n}\bigg\|_{\infty}
=\frac{\rho\|\gama_J\|_2^{2}}{m}-\Bigg(\sum_{j\in J}\|\gama_j\|_1
+\sum_{j\in J^{c}}\|\gama_{j}\|_1\Bigg)^{2}\bigg\|
\Sigma-\frac{\U^{\T}\U}{n}\bigg\|_{\infty}\\
\geq\frac{\rho\|\gama_J\|_2^{2}}{m}-m\Bigg(\sum_{j\in J}
\|\gama_j\|_2+\sum_{j\in J^{c}}\|\gama_j\|_2\Bigg)^{2}
\bigg\|\Sigma-\frac{\U^{\T}\U}{n}\bigg\|_{\infty}\\
\geq \frac{\rho\|\gama_J\|_2^{2}}{m}-m\Bigg(\sum_{j\in J}
\|\gama_j\|_2+3\sum_{j\in J}\|\gama_j\|_2\Bigg)^{2}
\bigg\|\Sigma-\frac{\U^{\T}\U}{n}\bigg\|_{\infty}\\
=\frac{\rho\|\gama_J\|_2^{2}}{m}-16m\Bigg(\sum_{j\in J}
\|\gama_j\|_2\Bigg)^{2}\bigg\|\Sigma
\geq \frac{\rho\|\gama_J\|_2^{2}}{m}-16rm
\sum_{j\in J}\|\gama_j\|_2^{2}\bigg\|\Sigma
= \frac{\rho\|\gama_J\|_2^{2}}{m}-16rm
\|\gama_{J}\|_2^{2}\bigg\|\Sigma
\end{align*}

It follows from (<ref>) that

\begin{align*}
\frac{\gama^{\T}\U^{\T}\U\gama}{n}\geq
\frac{\rho\|\gama_J\|_2^{2}}{m}-160rm_n\|\gama_{J}\|_2^{2}
\bigg\{\frac{\log(pqm)}{n}\bigg\}^{1/2}\geq
\frac{\|\gama_J\|_2^{2}\rho}{2m}
\end{align*}

holds with probability at least $1-6qm(pqm)^{-2}$, as long as
$160r\{m^{4}\log(pqm)/n\}^{1/2}\leq\rho/2$. This concludes the lemma.

The following lemma is crucial for proving
Theorem <ref>.

For every $\ell\in\{1,\ldots,p\}$, consider the random event
$\mathcal{A}_\ell=\cap_{j=1}^{q}\mathcal{A}_{j\ell}$, where

\begin{align*}
\mathcal{A}_{j\ell}=\bigg\{\frac{1}{n}\|\U^{\T}_{j}\varepsilon_{\ell}\|_{2}
\leq\frac{\lambda_{\ell}}{4}\bigg\}.
\end{align*}

If $\lambda_{\ell}\geq4c^{*}\sigma_{\ell}\{14\log(pqm)/n\}^{1/2}$, then

Note that

\begin{align*}
\pr(\mathcal{A}_{j\ell})=\pr\bigg(\frac{1}{n^{2}}\varepsilon_{\ell}^{\T}
\U_{j}\U_{j}^{\T}\varepsilon_{\ell}\leq\frac{\lambda^{2}_{\ell}}{16}\bigg)
\end{align*}

where $\xi_{1},\ldots,\xi_{n}$ are i.i.d standard normal random variables,
$\nu_{j,1},\ldots,\nu_{j,n}$ denote the eigenvalues of the matrix
$\U_{j}\U_{j}^{\T}/n$ among which the positive ones are the same as those of
$\Psi_{j}=\U_{j}^{\T}\U_{j}/n$, and

\begin{align*}
\end{align*}

We bound from the above the probability of $\mathcal{A}_{j}^{c}$ using Lemma
<ref>. Specifically, choose
$\nu=\nu_{j}=(\nu_{j,1},\ldots,\nu_{j,n}),x=x_{j}$ and
$m(\nu)=\|\Psi_{j}\|_{2}/\|\Psi_{j}\|_{F}=\|v_{j}\|_{\infty}/\|v_{j}\|_{2}$ to

\begin{align*}
\pr(\mathcal{A}^{c}_{j\ell})\leq 2\exp\bigg\{-\frac{x^{2}_{j\ell}}
\end{align*}

Now we find the appropriate $\lambda_{\ell}$ such that the above probability
approaches one as $p$ and $q$ increase with $n$. Let

\begin{align*}
\exp\bigg\{-\frac{x^{2}_{j\ell}}{2\big(1+\sqrt{2}x_{j\ell}
\|\Psi_{j}\|_2/\|\Psi_{j}\|_{\mathrm{F}}\big)}\bigg\}\leq (pqm)^{-2},
\end{align*}

which implies $x_{j\ell}^{2}\geq4\log(pqm)(1+\sqrt{2}x_{j\ell}
\|\Psi_{j}\|_2/\|\Psi\|_{F})$.

Equivalently, we need to have

\begin{align*}
\log(pqm)\}^{2}+4\log(pqm)\big]^{1/2},
\end{align*}

\begin{align*}
\lambda^{2}_{\ell}n/16\sigma_\ell^{2}\geq\mathrm{tr}(\Psi_{j})
\end{align*}

Therefore, a sufficient condition for the above is (by noting $(a+b)^{1/2}\leq
a^{1/2}+b^{1/2}$, $\forall a,b\geq0$)

\begin{align}\label{eq:lambda-bound}
\lambda^{2}_{\ell}n/16\sigma^{2}_\ell\geq \mathrm{tr}(\Psi_{j})
\end{align}

Therefore, when (<ref>) holds,
$\pr(\mathcal{A}_{j\ell}^{c})\leq2(pqm)^{-2}$. It suffices to find the
probability of (<ref>) when
$\lambda_{\ell}\geq4c^{*}\sigma_\ell \{14\log(pqm)/n\}^{1/2}$. By Lemma
<ref>, when $\{2\log(pqm)m^{3}/n\}^{1/2}\leq 1$, with
probability at least $1-\exp(-nm^{-3})\geq1-(pqm)^{-2}$, the following

\begin{align*}
\mathrm{tr}(\Psi_{j})=\mathrm{tr}\bigg(\frac{\U_j^{\T}\U_j}{n}\bigg)
\leq m\frac{c^{*}}{m}=c^{*}
\end{align*}

holds, where the first inequality follows from the positive semidefiniteness
of $\tU_j^{\T}11^{\T}\tU_j/n^{2}$. Similarly, $\|\Psi_{j}\|_{2}\leq c^{*}/m$
and $\|\Psi_{j}\|_{F}\leq\{\|\Psi_{j}\|_{2}\mathrm{tr}(\Psi_{j})\}^{1/2}\leq
c^{*}/\sqrt{m}$. Therefore, when $\lambda_{\ell}\geq4c^{*}\sigma_\ell
\{14\log(pqm)/n\}^{1/2}$,

\begin{align*}
\geq c^{*}+\frac{8c^{*}\log(pqm)}{m}
\end{align*}

holds with probability at least $1-(pqm)^{-2}$. As a result,
$\mathcal{A}_{j\ell}$ happens with probability at least $1-3(pqm)^{-2}$ when
$\lambda_{\ell}\geq4\sigma_\ell c^{*}\{14\log(pqm)/n\}^{1/2}$. The proof is
complete by taking the union bound argument.

Under the assumption of Lemma <ref>, the event

\begin{align*}
\mathcal{D}_\ell\coloneqq\Bigg\{\frac{1}{\sqrt{n}}\bigg\|
\U\bar{\gama}_{\ell}-\sum_{j=1}^{q}F_{j\ell}\bigg\|_{2}
\leq 4rC_Lm^{-d}+4r\big\{C_0^{2}\log(pqm)/n\big\}^{1/2}\Bigg\}
\end{align*}

happens with probability at least $1-2r(pqm)^{-2}$.

First, note that

\begin{align*}
\bar{\gamma}_{jk\ell}\psi_{k}(Z_{ij})\bigg\}^{2}\\
&=\frac{1}{n}\sum_{i=1}^{n}\bigg\{\sum_{j\in J_\ell}
f_{j\ell}(Z_{ij})-\sum_{j\in J_\ell}\sum_{k=1}^{m}
\bar{\gamma}_{jk\ell}\psi_{k}(Z_{ij})\bigg\}^{2}\\
&\leq \frac{1}{n}r\sum_{i=1}^{n}\sum_{j\in J_\ell}
\bigg\{f_{j\ell}(Z_{ij})-\sum_{k=1}^{m}
\bar{\gamma}_{jk\ell}\psi_{k}(Z_{ij})\bigg\}^{2}
\leq r\sum_{j\in J_\ell}\sup_{z\in[a,b]}
\big|f_{j\ell}(z)-\tilde{f}_{nj\ell}(z)\big|^{2}.
\end{align*}

Applying Lemma <ref> with a union bound argument,
with probability at least $1-2r(pqm)^{-2}$, we have

\begin{align*}
\frac{1}{\sqrt{n}}\bigg\|\U\bar{\gama}_{\ell}
-\sum_{j=1}^{q}F_{j\ell}\bigg\|_{2}\leq 4rC_Lm^{-d}
\end{align*}

This completes the proof.

The next lemma resembles Lemma <ref> but focuses on the classical
restricted eigenvalue condition.

Suppose Assumption <ref> holds. Then, with probability at
least $1-2(pqm)^{-2}$,

\begin{align*}
\min\bigg\{\frac{\bata^{\T}F^{\T}F\bata}{n\|\bata_{\mathcal{L}}\|_2^{2}}
\colon|\mathcal{L}|\leq s,\bata\in\mathbb{R}^{p}\backslash\{0\},
\sum_{\ell\in \mathcal{L}^{c}}|\beta_{\ell}|\leq
3\sum_{\ell\in \mathcal{L}}|\beta_{\ell}|\bigg\}\geq\frac{\kappa}{2}
\end{align*}

holds when $64sr^{2}C_0\{\log(pqm)/n\}^{1/2}\leq\kappa/2$ and

Consider $\bata$ satisfying $\sum_{\ell\in\mathcal{L}^{c}}|\beta_{\ell}|
\leq3\sum_{\ell\in\mathcal{L}}|\beta_{\ell}|$. Then,

\begin{align*}
&\geq \|\bata_{\mathcal{L}}\|_2^{2}\kappa-\bata^{\T}
\bigg(\Sigma_f-\frac{F^{\T}F}{n}\bigg)\bata
\geq \|\bata_{\mathcal{L}}\|_2^{2}\kappa-\|\bata\|_1^{2}
\bigg\|\Sigma_f-\frac{F^{\T}F}{n}\bigg\|_{\infty}\\
\Sigma_f-\frac{F^{\T}F}{n}\bigg\|_{\infty}\\
&\geq \|\bata_{\mathcal{L}}\|_2^{2}\kappa-\bigg(3\sum_{\ell\in\mathcal{L}}
\bigg\|\Sigma_f-\frac{F^{\T}F}{n}\bigg\|_{\infty}\\
& =\|\bata_{\mathcal{L}}\|_2^{2}\kappa-16\|\bata_{\mathcal{L}}\|_1^{2}
\bigg\|\Sigma_f-\frac{F^{\T}F}{n}\bigg\|_{\infty}
\geq \|\bata_{\mathcal{L}}\|_2^{2}\kappa-16s\|\bata_{\mathcal{L}}\|_2^{2}
\bigg\|\Sigma_f-\frac{F^{\T}F}{n}\bigg\|_{\infty}.
\end{align*}

We first bound $\|\Sigma_f-F^{\T}F/n\|_{\infty}$,

\begin{align*}
\bigg\|\Sigma_f-\frac{F^{\T}F}{n}\bigg\|_{\infty}
\big\{\sum_{j=1}^{q}f_{j\ell'}(Z_{ij})\big\}}{n}\Bigg|.
\end{align*}

As there are at most $r$ nonzero summands in $\sum_{j=1}^{q}f_{j\ell}(z)$ for
every $\ell\in\{1,\ldots,p\}$, we have

\begin{align*}
\bigg|\sum_{j=1}^{q}f_{j\ell}(z)\bigg|\leq rC_0,
\ \forall\ell\in\{1,\ldots,p\},
\end{align*}

which implies

\begin{align*}
\bigg|\bigg\{\sum_{j=1}^{q}f_{j\ell}(Z_j)\bigg\}
\bigg\{\sum_{j=1}^{q}f_{j\ell}(Z_j)\bigg\}\bigg|\leq r^{2}C_0^{2}.
\end{align*}

Now applying Lemma <ref> with $b=r^{2}C_0^{2}$ and
$\delta=4r^{2}C_0^{2}\{\log(pqm)/n\}^{1/2}$, with probability at least
$1-2\exp\big[-8r^{4}C_0^{4}\log(pqm)/ \big\{r^{4}C_0^{4}+\frac{4}{3}
r^{4}C_0^{4}(\log(pqm)/n)^{1/2}\big\}\big]$, we have

\begin{align*}
\bigg\}\bigg\{\sum_{j=1}^{q}f_{j\ell'}(Z_j)\bigg\}\bigg]
\big\}\big\{\sum_{j=1}^{q}f_{j\ell'}(Z_{ij})\big\}}{n}\Bigg|\\
\end{align*}

When $\{\log(pqm)/n\}^{1/2}\leq1/2$, we have

\begin{align*}
{\big\{r^{4}C_0^{4}+\frac{4}{3} r^{4}C_0^{4}
(\log(pqm)/n)^{1/2}\big\}}\bigg]\geq 1-2(pqm)^{-4}
\end{align*}

and by taking the union bound we get

\begin{align*}
\bigg\|\Sigma_f-\frac{F^{\T}F}{n}\bigg\|_{\infty}\leq
\end{align*}

which holds with probability at least $1-2p^{2}(pqm)^{-4}$. It follows that

\begin{align*}
\frac{\bata^{\T}F^{\T}F\bata}{n}\geq \|\bata_{\mathcal{L}}\|_2^{2}
\kappa-64s\|\bata_{\mathcal{L}}\|_2^{2}r^{2}C_0^{2}
\bigg\{\frac{\log(pqm)}{n}\bigg\}^{1/2}
\end{align*}

holds with probability at least $1-2p^{2}(pqm)^{-4}$. When
$64sr^{2}C_0^{2}\{\log(pqm)/n\}^{1/2}\leq\kappa/2$, we have
$\bata^{\T}F^{\T}F\bata/n\geq\kappa\|\bata_{\mathcal{L}}\|_2^{2}/2$ holds with
probability at least $1-2(pqm)^{-2}$. This concludes the lemma.

§ SECTION B: PROOF OF RESULTS IN SECTION <REF>

§.§ Proof of Lemma <ref> in the main text

By the assumption $\mathbb{E}\{f_{j\ell}(z_{ij})\}=0$ and Lemma
<ref>, we know

\begin{align*}
\bar{\gamma}_{jk\ell}\phi_{k}(Z_{ij})-\frac{1}{n}\sum_{i=1}^{n}
\qquad\qquad+\bigg|\frac{1}{n}\sum_{i=1}^{n}
&\leq C_Lm_n^{-d}+\frac{1}{n}\sum_{i=1}^{n}
\big|\Bar{f}_{nj\ell}(Z_{ij})-f_{j\ell}(Z_{ij})\big|
&\leq 2C_Lm_n^{-d}+\bigg|\frac{1}{n}\sum_{i=1}^{n}
\end{align*}

where the first inequality follows from the triangle inequality and the second
to the last inequalities follow from Lemma <ref>. Finally,
noting $|f(Z_{ij})|\leq C_0$ and applying Lemma <ref> with
$t=\{4C_0^{2}\log(pqm)/n\}^{1/2}$, we have

\begin{align*}
\bigg|\frac{1}{n}\sum_{i=1}^{n}f_{j\ell}(Z_{ij})
-\mathbb{E}(f_{j\ell})\bigg|\leq \{4C_0^{2}\log(pqm)/n\}^{1/2},
\end{align*}

which holds with probability at least $1-2(pqm)^{-2}$. This completes the

§ SECTION C: PROOF OF RESULTS IN SECTION <REF>

§.§ Proof of Theorem <ref> in the main text

Consider the $\ell$th optimization problem in (<ref>). By the
optimality of $\widehat{\gama}_\ell$ in (<ref>), we have

\begin{align*}
\frac{1}{n}\left\|\U\widehat{\gama}_\ell-\X_\ell\right\|_2^{2}
\leq\frac{1}{n}\left\|\U\gama_\ell-\X_\ell\right\|_2^{2}
\end{align*}

Recall that $\sum_{j=1}^{q}F_{j\ell}+\varepsilon_\ell=\X_\ell$. Therefore,

\begin{align*}
\|\widehat{\gama}_{j\ell}\|_2\\
\end{align*}

Thus, we obtain

\begin{align*}
\bigg(\varepsilon_\ell+\U\bar{\gama}_{\ell}-\sum_{j=1}^{q}F_{j\ell}\bigg)
\|\widehat{\gama}_{j\ell}\|_2\\
\bigg\|\sum_{j=1}^{q}F_{j\ell}-\U\bar{\gama}_\ell\bigg\|_2^{2}\\
\bigg(\varepsilon_\ell+\U\bar{\gama}_\ell-\sum_{j=1}^{q}F_{j\ell}\bigg)
\end{align*}

which implies that

\begin{align*}
\bigg(\varepsilon_\ell+\U\bar{\gama}_\ell
\sum_{j=1}^{q}\|\gama_{j\ell}\|_2\\
\bigg(\U\bar{\gama}_\ell-\sum_{j=1}^{q}F_{j\ell}\bigg)
\left\|\widehat{\gama}_{j\ell}-\gama_{j\ell}\right\|_2\\
\bigg(\U\bar{\gama}_{\ell}-\sum_{j=1}^{q}F_{j\ell}\bigg)\bigg\|_2
\|\widehat{\gama}_{j\ell}-\gama_{j\ell}\|_2
\end{align*}

where the last inequality holds due to the Cauchy–Schwarz inequality. By Lemma
<ref> and the choice of $\lambda_\ell$, we have

\begin{align*}
\|\widehat{\gama}_{j\ell}-\gama_{j\ell}\|_2\\
\frac{2\big\|\U\bar{\gama}_{\ell}-\sum_{j=1}^{q}F_{j\ell}\big\|_2}
\end{align*}

with probability at least $1-3q(pqm)^{-2}$. Then by Lemma
<ref>, we have

\begin{align*}
\|\widehat{\gama}_{j\ell}-\gama_{j\ell}\|_2\\
\frac{2\|\U\bar{\gama}_{\ell}-\sum_{j=1}^{q}F_{j\ell}\|_2}{n^{\frac12}}
\|\widehat{\gama}_{j\ell}-\gama_{j\ell}\|_2
\end{align*}

which holds with probability at least $1-4q(pqm)^{-2}$ (taking the union bound
here) when $1-\exp(-nm^{-3})\geq 1-(pqm)^{-2}$. Then by Lemma
<ref>, we have

\begin{align*}
\|\widehat{\gama}_{j\ell}-\gama_{j\ell}\|_2
\frac{c^{*}C_0^{2}\log(pqm)}{nm}\biggr\}^{1/2}\bigg]
\sum_{j=1}^{q}\|\widehat{\gama}_{j\ell}-\gama_{j\ell}\|_2\\
\|\widehat{\gama}_{j\ell}-\gama_{j\ell}\|_2
\end{align*}

which holds with probability at least $1-6q(pqm)^{-2}$ since $r\leq q$.

Setting $\gama_\ell=\bar{\gama}_\ell$, we have

\begin{align}
\frac{1}{n}\|\U\widehat{\gama}_\ell-\U\bar{\gama}_\ell\|_2^{2}
\leq 2\lambda_\ell\sum_{j=1}^{q}\|\widehat{\gama}_{j\ell}
\|\bar{\gama}_{j\ell}\|_2-2\lambda_\ell\sum_{j=1}^{q}
\|\widehat{\gama}_{j\ell}\|_2\\
\leq4\lambda_\ell\sum_{j\in J_\ell}
\|\bar{\gama}_{j\ell}-\widehat{\gama}_{j\ell}\|_2.
\end{align}

Therefore, we know $\lambda_{\ell}\sum_{j=1}^{q}
\|\widehat{\gama}_{j\ell}-\bar{\gama}_{j\ell}\|_2\leq
-\widehat{\gama}_{j\ell}\|_2$, which implies that

\begin{align*}
\sum_{j\in J^{c}_{\ell}}\|\widehat{\gama}_{j\ell}-\bar{\gama}_{j\ell}\|_2
\leq3\sum_{j\in J_\ell}\|\bar{\gama}_{j\ell}-\widehat{\gama}_{j\ell}\|_2
\end{align*}

holds with probability at least $1-6q(pqm)^{-2}$. Then, combine with Lemma
<ref> to obtain

\begin{align*}
\|\bar{\gama}_{J_{\ell}\ell}-\widehat{\gama}_{J_{\ell}\ell}\|_2
\leq\frac{\|\U\bar{\gama}_{\ell}-\U\widehat{\gama}_{\ell}\|_2}
\end{align*}

which holds with probability at least $1-12(qm)(pqm)^{-2}$. We also know from
inequality (<ref>) that

\begin{align*}
\frac{1}{n}\|\U\widehat{\gama}_\ell-\U\bar{\gama}_{\ell}\|_2^{2}
\leq4\lambda_\ell\sum_{j\in J_\ell}
\|\bar{\gama}_{j\ell}-\widehat{\gama}_{j\ell}\|_2
\leq4\lambda_\ell\sqrt{r}\|\bar{\gama}_{J_\ell}
\end{align*}

which holds with probability at least $1-6q(pqm)^{-2}$, where the second
inequality follows from the Cauchy–Schwarz inequality. Then, it follows that

\begin{align*}%\label{eq:estimation_bound}
\frac{\|\U\widehat{\gama}_\ell-\U\bar{\gama}_{\ell}\|_2}{n^{1/2}}
\leq4\lambda_{\ell}r^{1/2}\bigg(\frac{2m}{\rho}\bigg)^{1/2}
\end{align*}

holds with probability at least $1-18(qm)(pqm)^{-2}$. It holds with the same
probability that

\begin{align*}
\leq4\sum_{j\in J_\ell}\|\bar{\gama}_{j\ell}
\|\bar{\gama}_{J_\ell}-\widehat{\gama}_{J_\ell}\|_2\\
&\leq 4r^{1/2}\frac{\|\U\bar{\gama}_{\ell}
\bigg(\frac{2m}{\rho}\bigg)^{1/2}
\leq 32r\lambda_\ell\frac{m}{\rho}.
\end{align*}

Now we apply Lemma <ref> to obtain

\begin{align*}
\end{align*}

which holds with probability at least $1-20(qm)(pqm)^{-2}$, where the last
inequality follows from the definition of $\lambda_\ell$. The proof is complete
by taking the union bound over all $\ell$.

§.§ Proof of Theorem <ref>

Suppose Assumptions <ref>–<ref> hold. If the
regularization parameters satisfy

\begin{align}\label{eq:lambda_max}
\leq\frac{\kappa^{2}}{4rs},
\end{align}

then with probability at least $1-62(pqm)^{-1}$, the matrix
$\widehat{\X}=\U\widehat{\Gama}$ satisfies

\begin{align*}
\min\bigg\{\frac{\|\widehat{\X}\bata\|}{n^{1/2}\|\bata_{\L}\|}
\colon|\L|\leq s,\bata\in\mathbb{R}^{p}\backslash \{0\},
\|\bata_{\mathcal{L}^{c}}\|_{1}\leq 3\|\bata_{\mathcal{L}}\|_{1}\bigg\}
\geq\frac{\kappa}{2}
\end{align*}

when $n$ is sufficiently large.

For any subset $\L\subset\{1,\ldots,p\}$ with $|\L|\leq s$ and any
$\delta\in\mathbb{R}^{p}$ such that $\delta\neq0$ and
$\|\delta_{\L^{c}}\|_{1}\leq3\|\delta_{\L}\|_{1}$, we have

\begin{align*}
\frac{\delta^{\T}(F^{\T}F-\widehat{\X}^{\T}\widehat{\X})\delta}
\max_{1\leq\ell,\ell'\leq p}\bigg|\Bigl(\sum_{j=1}^{q}F_{j\ell'}\Bigr)^{\T}
\Bigl(\sum_{j=1}^{q}F_{j\ell}\Bigr)-\widehat{\X}^{\T}_{\ell'}
\widehat{\X}_{\ell}\bigg|}{n\|\delta_{\L}\|^{2}_{2}}.
\end{align*}

Since $\|\delta_{\L^{c}}\|_{1}\leq 3\|\delta_{\L}\|_{1}$, we have
\leq16\|\delta_{\L}\|^{2}_{1}\leq16s\|\delta_{\L}\|^{2}_{2}$, which implies that

\begin{align}\label{eq:entrywise_bound}
\frac{\delta^{\T}(F^{\T}F-\widehat{\X}^{\T}\widehat{\X})\delta}
{n\|\delta_{\L}\|^{2}_{2}}\leq\frac{16s\max_{1\leq \ell,\ell'\leq p}
\bigg|\Bigl(\sum_{j=1}^{q}F_{j\ell'}\Bigr)^{\T}
\Bigl(\sum_{j=1}^{q}F_{f\ell}\Bigr)
\end{align}

To bound the entrywise maximum, we write

\begin{align*}
\bigg(\sum_{j=1}^{q}F_{j\ell'}\bigg)
\bigg(\widehat{\X}_{\ell'}-\sum_{j=1}^{q}
\bigg(\widehat{\X}_{\ell}-\sum_{j=1}^{q}
\bigg(\widehat{\X}_{\ell'}
\end{align*}

For $T_1$, by the Cauchy–Schwarz inequality and Theorem
<ref>, we have

\begin{align*}
\end{align*}

which holds with probability at least $1-20(pqm)^{-1}$. For $T_2$, note that
$\forall\ell\in\{1,\ldots,p\}$ due to the sparsity assumption. Then, we have

\begin{align*}
\bigg\|\sum_{j=1}^{q}F_{j\ell'}\bigg\|_{2}
\leq5n\lambda_{\max}\bigg(\frac{2rm}{\rho}\bigg)^{1/2}rC_0
\end{align*}

which holds with probability at least $1-20(pqm)^{-1}$. Similarly, it holds with
the same probability that

\begin{align*}
\end{align*}

Therefore, we have

\begin{align*}
\max_{1\leq \ell,\ell'\leq p}\bigg|\bigg(\sum_{j=1}^{q}
\leq35n\lambda_{\max}rC_f\bigg(\frac{2rm}{\rho}\bigg)^{1/2},
\end{align*}

which holds with probability at least $1-60(pqm)^{-1}$ when
$\lambda_{\max}\{2m/\rho\}^{1/2}\leq r^{1/2} C_0$. We know from Lemma
<ref> that

\begin{align*}
\min\bigg\{\frac{\|\f\bata\|}{\sqrt{n}\|\bata_{\L}\|}
\colon|\L|\leq s,\bata\in\mathbb{R}^{p}\backslash \{0\},
\|\bata_{\L^{c}}\|_{1}\leq 3\|\bata_{\L}\|_{1}\bigg\}
\geq\frac{\sqrt{2}\kappa}{2}
\end{align*}

holds with probability at least $1-2(pqm)^{-2}$. It follows from
(<ref>) and (<ref>) that

\begin{align*}
\frac{\delta^{\T}(F^{\T}F-\widehat{\X}^{\T}\widehat{\X})\delta}
\bigg(\frac{2rm}{\rho}\bigg)^{1/2}\leq\frac{\kappa^{2}}{4}
\end{align*}

holds with probability at least $1-60(pqm)^{-1}$. Therefore, with probability at
least $1-62(pqm)^{-1}$, we get

\begin{align*}
\frac{\bata^{\T}(\widehat{\X}^{\T}\widehat{\X})\bata}
\end{align*}

This completes the proof.

Under Assumptions <ref>–<ref>, if the
regularization parameters $\lambda_{\ell}$'s are chosen as in Theorem
<ref> and

\begin{align*}
\mu=2\lambda_{\max}r\bigg(\frac{2m}{\rho}\bigg)^{1/2}
\big(7\sigma_0+8\sqrt{5}B\max_{\ell}\sigma_{\ell}+30B\big),
\end{align*}

then with probability at least $1-86(pqm)^{-1}$, the regularized estimator
$\widehat{\bata}$ of (<ref>) satisfies

\begin{align*}
\frac{1}{2n}\big\|\widehat{\X}(\widehat{\bata}-\bata)\big\|^{2}_{2}
\end{align*}

By the optimality of $\widehat{\bata}$, we have

\begin{align*}
\frac{1}{2n}\big\|\Y-\widehat{\X}\widehat{\bata}\big\|^{2}_{2}
\big\|\Y-\widehat{\X}\bata\big\|^{2}_{2}+\mu\|\bata\|_{1}.
\end{align*}

Substituting $\Y=\X\bata+\yita$, we have

\begin{align*}
\widehat{\X}(\widehat{\bata}-\bata)
\end{align*}

\begin{align*}
\|\Y-\widehat{\X}\bata\|^{2}_{2}
\end{align*}

It then follows that

\begin{align}
\frac{1}{2n}\|\widehat{\X}(\widehat{\bata}-\bata)\|^{2}_{2}
\widehat{\X}(\widehat{\bata}-\bata)\\
\leq\mu\|\bata\|_{1}-\mu\|\widehat{\bata}\|_{1}
\|\widehat{\bata}-\bata\|_{1}.
\end{align}

Next, we find a probability bound for the following event

\begin{align*}%\label{eq:inter_res}
\bigg\|\frac{1}{n}\widehat{\X}^{\T}\yita-\frac{1}{n}
\widehat{\X}^{\T}(\widehat{\X}-\X)\bata\bigg\|_{\infty}
\leq\frac{\mu}{2}.
\end{align*}

Substituting $\widehat{\X}=\U\widehat{\Gama}$ and $\X=F+\varepsilon$, we write

\begin{align*}
\big(\U\widehat{\Gama}-\X\big)\bata\\
\varepsilon\bata\\
\end{align*}

For $T_1$, it follows from Theorem <ref> that

\begin{align*}
\|T_1\|_\infty\leq\frac{1}{n} \max_{\ell}
\bigg\|\sum_{j=1}^{q}F_{j\ell}-\U\widehat{\gama}_\ell\bigg\|_2
\|\yita\|_2\leq5\lambda_{\max}r^{1/2}
\bigg(\frac{2m}{\rho}\bigg)^{1/2}\frac{\|\eta\|_2}{n^{1/2}}
\end{align*}

holds with probability at least $1-20(pqm)^{-1}$. For a Gaussian variable
$\eta_i$, $\eta_i^{2}$ is sub-exponential and follows the
$\sigma_{0}^{2}\chi^{2}$ distribution. Apply Lemma <ref> with
$t=\{16\log(pqm)/n\}^{1/2}$ to obtain

\begin{align*}
\bigg|\frac{1}{\sigma^{2}_0n}\|\eta\|^{2}_2-1\bigg|
\leq4\bigg\{\frac{\log(pqm)}{n}\bigg\}^{1/2},
\end{align*}

which holds with probability at least $1-2(pqm)^{-2}$. Therefore, we have

\begin{align*}
\leq5\lambda_{\max}r^{1/2}\bigg(\frac{2m}{\rho}\bigg)^{1/2}
\bigg(\frac{\|\eta\|_2^{2}}{n}\bigg)^{1/2}\\
\bigg[4\sigma_0^{2}\{\log(pqm)/n\}^{1/2}+\sigma_0^{2}\bigg]\\
\end{align*}

which holds with probability at least $1-22(pqm)^{-1}$ when
$\log(pqm)/n\leq1$. For $T_2$, we have

\begin{align*}
\|T_2\|_{\infty}=\frac{1}{n}\|\f^{\T}\yita\|_{\infty}
\leq\frac{1}{n}\max_{\ell}\bigg|\bigg(\sum_{j=1}^{q}
\end{align*}

By the tail bound of a standard normal random variable $X$

\begin{align*}
\pr\left(|X|\geq{}t\right)\leq2\bigg(
\frac{2}{\pi}\bigg)^{1/2}\frac{\exp(-t^{2}/2)}{t},
\end{align*}

we have

\begin{align*}
&\pr\left(\|T_2\|_{\infty}\geq t\right)
\leq\pr\bigg\{\frac{1}{n}\max_{\ell}\bigg|
\bigg(\sum_{j=1}^{q}F_{j\ell}\bigg)^{\T}
\yita\bigg|\geq t\bigg\}\\
\sum_{j=1}^{q}F_{j\ell}^{\T}\yita\bigg|\geq{}t\bigg)
=p\pr\bigg(\frac{1}{n}\bigg|\sum_{j\in J_\ell}
F_{j\ell}^{\T}\yita\bigg|\geq t\bigg)\\
&=p\pr\Bigg[\frac{\big|\sum_{j\in J_\ell}F_{j\ell}^{\T}
\yita\big|}{\big\{\sum_{i=1}^{n}\big(\sum_{j\in J_{\ell}}
\frac{nt}{\big\{\sum_{i=1}^{n}\big(\sum_{j\in J_{\ell}}
\exp\Bigg[\frac{-n^{2}t^{2}}{2\sigma_0^{2}
\sum_{i=1}^{n}\big\{\sum_{j\in J_{\ell} }
\frac{\sigma_0\bigg[\sum_{i=1}^{n}\big\{
\sum_{j\in J_{\ell}}f_{j\ell}(Z_{ij})\big\}^{2}\bigg]}{nt}.
\end{align*}

Note that

$\sum_{i=1}^{n}\bigl\{\sum_{j\in J_{\ell} } f_{j\ell}(Z_{ij})\bigr\}^{2}\leq

By setting $t=2C_0\{r^{2}\log(pqm)/n\}^{1/2}\sigma_0$, we obtain

\begin{align*}
\pr\bigg[\|T_2\|_{\infty}\geq 2C_0\{r^{2}\log(pqm)/n\}^{1/2}
\sigma_0\bigg]\leq p(pqm)^{-2}.
\end{align*}

Therefore, we can bound $T_2$ as

\begin{align*}
\|T_2\|_\infty\leq2C_0\sigma_0\bigg\{
\frac{r^{2}\log(pqm)}{n}\bigg\}^{1/2},
\end{align*}

which holds with probability at least $1-p(pqm)^{-2}$. Now for $T_3$, we have

\begin{align*}
\|T_3\|_{\infty}\leq \|\bata\|_1\frac{1}{n}
\big\|\big(\U\widehat{\Gama}-F\big)^{\T}
\varepsilon\big\|_{\infty}\leq{}B\frac{1}{n}
\max_{\ell,\ell'}\bigg\|\U\widehat{\gama}_\ell
\end{align*}

It follows from Theorem <ref> that

\begin{align*}
\max_{\ell}\bigg\|\U\widehat{\gama}_\ell-\sum_{j=1}^{q}
F_{j\ell}\bigg\|_2\leq 5\lambda_{\max}r^{1/2}
\bigg(\frac{2m}{\rho}\bigg)^{1/2}\sqrt{n},
\end{align*}

which implies that

\begin{align*}
\|T_3\|_{\infty}\leq5B\lambda_{\max}r^{1/2}
\bigg(\frac{2m}{\rho}\bigg)^{1/2}
\max_{\ell'}\frac{\|\varepsilon_{\ell'}\|_2}{\sqrt{n}}
\end{align*}

holds with probability at least $1-20(pqm)^{-1}$. Now applying Lemma
<ref> again with $t=\{16\log(pqm)/n\}^{1/2}$, we get

\begin{align*}
\left|\frac{\|\varepsilon_{\ell'}\|_2^{2}}{n\sigma_{\ell'}^{2}}
-1\right|\leq 4\bigg\{\frac{\log(pqm)}{n}\bigg\}^{1/2}
\end{align*}

with probability at least $1-2(pqm)^{-2}$, which implies

\begin{align*}
\frac{\|\varepsilon_{\ell'}\|_2}{n^{1/2}}\leq
\bigg[4\sigma_{\ell'}^{2}\{\log(pqm)/n\}^{1/2}
+\sigma_{\ell'}^{2}\bigg]\leq \sqrt{5}\sigma_{\ell'}
\end{align*}

with probability at least $1-2(pqm)^{-2}$ when $\log(pqm)/n\leq 1$. Therefore,
we have

\begin{align*}
\|T_3\|_{\infty}\leq 5\sqrt{5}\max_{\ell'}\sigma_{\ell'}B
\lambda_{\max}r^{1/2}\bigg(\frac{2m}{\rho}\bigg)^{1/2}
\bigg(\frac{10rm}{\rho}\bigg)^{1/2},
\end{align*}

which holds with probability at least $1-22(pqm)^{-1}$. Similar to $T_2$, we
can bound $T_4$ as

\begin{align*}
\|T_4\|_{\infty}= \frac{1}{n}\|F^{\T}\eps\bata\|_{\infty}
\leq\frac{\|\bata\|_1}{n}\left\|F^{\T}\eps\right\|_{\infty}
\leq\frac{B}{n}\big\|F^{\T}\eps\big\|_{\infty}\leq
\frac{B}{n}\max_{\ell,\ell'}\bigg|\bigg(\sum_{j=1}^{q}
\end{align*}

Therefore, we have

\begin{align*}
&\pr(\|T_4\|_{\infty}\geq t)
\leq p^{2}\pr\bigg\{\frac{B}{n}\bigg|\bigg(\sum_{j=1}^{q}
F_{j\ell}\bigg)^{\T}\varepsilon_{\ell'}\bigg|\geq t\bigg\}
=p^{2}\pr\bigg\{\frac{B}{n}\bigg|\bigg(\sum_{j\in J_\ell}
F_{j\ell}\bigg)^{\T}\varepsilon_{\ell'}\bigg|\geq t\bigg\}\\
& =p^{2}\pr\Bigg[\frac{\big|\big(\sum_{j\in J_\ell}F_{j\ell}\big)^{\T}
\varepsilon_{\ell'}\big|}{\big\{\sum_{i=1}^{n}
\big(\sum_{j\in J_{\ell} }f_{j\ell}(Z_{ij})\big)^{2}\big\}^{1/2}
\sigma_{\ell'}}\geq \frac{nt}{B\big\{\sum_{i=1}^{n}
\big(\sum_{j\in J_{\ell} }f_{j\ell}(Z_{ij})\big)^{2}
\big\}^{1/2}\sigma_{\ell'}}\Bigg]\\
\exp\Bigg[\frac{-n^{2}t^{2}}{2\sigma_{\ell'}^{2}B^{2}
\sum_{i}\big\{\sum_{j\in J_{\ell} }f_{j\ell}(Z_{ij})
\big\}^{2}}\Bigg]\frac{B\big[\sum_{i}
\big\{\sum_{j\in J_{\ell} }f_{j\ell}(Z_{ij})\big\}^{2}
\big]^{1/2}\sigma_{\ell'}}{nt}.
\end{align*}

Since $\sum_{i=1}^{n}\bigl\{\sum_{j\in{}J_{\ell}}f_{j\ell}(Z_{ij})\bigr\}^{2}
\leq{}r^{2}C_0^{2}n$, by setting
$t=\max_{\ell'}\sigma_{\ell'}BC_0\{8r^{2}\log(pqm)/n\}^{1/2}$, we obtain

\begin{align*}
\|T_4\|_{\infty}\leq \max_{\ell'}\sigma_{\ell'}BC_0
\{8r^{2}\log(pqm)/n\}^{1/2},
\end{align*}

which holds with probability at least $1-p^{2}(pqm)^{-4}$. Now we bound $T_5$

\begin{align*}
\leq\frac{\|\bata\|_1}{n}\big\|\big(\widehat{\Gama}^{\T}
\U^{\T}-F^{\T}\big)\big(\U\widehat{\Gama}-F\big)\big\|_{\infty}\\
\leq25B\lambda_{\max}^{2}r\frac{2m}{\rho},
\end{align*}

which holds with probability at least $1-20(pqm)^{-1}$. Finally, we bound $T_6$

\begin{align*}
\leq\frac{\|\bata\|_1}{n}\big\|F^{\T}
\big(\U\widehat{\Gama}-F\big)\big\|_{\infty}
\leq\frac{B}{n}\max_{\ell,\ell'}\bigg\|\sum_{j=1}^{q}
\bigg\|\U\widehat{\gama}_{\ell'}
\leq5Br\lambda_{\max}\bigg(\frac{2m}{\rho}\bigg)^{1/2},
\end{align*}

which holds with probability at least $1-20(pqm)^{-1}$. Therefore, we have

\begin{align*}
\bata\bigg\|_{\infty}\\
&\leq 5\sigma_0\lambda_{\max}\bigg(\frac{2rm}{\rho}\bigg)^{1/2}
\bigg(\frac{10rm}{n}\bigg)^{1/2}\\
\bigg\}^{1/2}+25B\lambda_{\max}^{2}r\frac{2m}{\rho}
\frac{10m}{n}\bigg)^{1/2}+30rB\lambda_{\max}
\bigg(\frac{2m}{\rho}\bigg)^{1/2}\\
\big(7\sigma_0+8\sqrt{5}B\max_{\ell}\sigma_{\ell}+30B\big),
\end{align*}

which holds with probability at least $1-86(pqm)^{-1}$ when
$\lambda_{\max}(2m/\rho)^{1/2}\leq1$, where we use the definition of
$\lambda_{\max}$ in the third inequality. This, together with
(<ref>), implies

\begin{align*}
\frac{1}{2n}\|\widehat{\X}(\widehat{\bata}-\bata)\|^{2}_{2}
\leq\mu\|\bata\|_{1}-\mu\|\widehat{\bata}\|_{1}
\end{align*}

Adding $\mu\|\widehat{\bata}-\bata\|^{2}_{2}/2$ to both sides yields

\begin{align*}
\leq\mu(\|\bata\|_{1}-\|\widehat{\bata}\|_{1}
\leq2\mu\|\widehat{\bata}_{\L}-\bata_{\L}\|_{1},
\end{align*}

where the last inequality follows from the triangle inequality. This completes
the proof.

Now we are ready to prove our main results. We note from Lemma
<ref> that with probability at least
$1-86(pqm)^{-1}$, we have

\begin{align}\label{eq:combination_estimation_pred}
\frac{1}{2n}\|\widehat{\X}(\widehat{\bata}-\bata)\|^{2}_{2}
\leq2\mu\|\widehat{\bata}_{\L}-\bata_{\L}\|_{1}\leq2\mu s^{1/2}
\|\widehat{\bata}_{\L}-\bata_{\L}\|_{2}
\end{align}

\begin{align*}
\frac{\mu}{2}\|\widehat{\bata}-\bata\|_{1}
\leq 2\mu\|\widehat{\bata}_{\L}-\bata_{\L}\|_{1}.
\end{align*}

By Lemma <ref>, with probability at least
$1-148(pqm)^{-1}$, we have

\begin{align}\label{eq:re_application}
\|\widehat{\bata}_{\L}-\bata_{\L}\|_{2}\leq
\frac{2\|\widehat{\X}(\widehat{\bata}-\bata)\|_{2}}
\end{align}

Combining (<ref>) and (<ref>),
we obtain

\begin{align*}
\|\widehat{\X}(\widehat{\bata}-\bata)\|^{2}_{2}
\leq\frac{64}{\kappa^{2}}ns\mu^{2}
\end{align*}

\begin{align*}
\|\widehat{\bata}-\bata\|_{1}\leq4\|\widehat{\bata}_{\L}-\bata_{\L}\|_{1}
\leq4s^{1/2}\|\widehat{\bata}_{\mathcal{L}}-\bata_{\mathcal{L}}\|_{2}
\leq\frac{64}{\kappa^{2}}s\mu,
\end{align*}

which hold with probability at least $1-234(pqm)^{-1}$.

§ SECTION D: PROOF OF RESULTS IN SECTION <REF>

§.§ Proof of Lemma <ref>

Note that

\begin{align*}
\widehat{\Sigma}_{\d}-I)(\bata-\widehat{\bata}).
\end{align*}

Now decompose the second term on the rightmost-hand side of the above display

\begin{align*}
\yita/n+\widehat{\Omega}(\widehat{\D}-\D)^\T\yita/n,
\end{align*}

which completes the proof.

§.§ Proof of Theorem <ref>

In this proof, we will temporarily assume $\|R_k\|_{\infty}=o_p(1),k=1,2,3,4$,
which will be elaborated in the next section. First, we show $\omega_\ell$ is
bounded away from zero for large $n$, that is, $\omega_{\ell}\geq c$ for some
constant $c>0$. This follows immediately since $\Omega_{\ell\ell}$ is lower
bounded by a constant. Define

\begin{align*}
\big(\theta_{\ell}^{\T}\Gama^{\T}
\tU_i^{\T}\big)\frac{\yita_i}{\omega_{\ell}},
\end{align*}

where $\tU_{i}$ is the $i$th row of $\tU$. Then, we have

\begin{align*}
\Gama^{\T}\tU_i^{\T}\frac{\yita_i}{\omega_{\ell}}
-\sum_{l=1}^{4}\frac{R_{\ell l}}{\omega_{\ell}}\bigg|\\
\theta_{\ell}^{\T}\Gama^{\T}\tU_i^{\T}
\frac{\yita_i}{\omega_{\ell}}-\frac{\theta_{\ell}^{\T}
\D^{\T}\yita}{\sqrt{n}\omega_{\ell}}\bigg|+\sum_{l=1}^{4}
\frac{\left\|R_{l}\right\|_{\infty}}{\min_{\ell}\omega_{\ell}}.
\end{align*}

As noted above, we will control the remainder terms $\|R_k\|_{\infty}$ in the
next section and for now we assume they are all $o_{p}(1)$. Now we have

\begin{align*}
\Gama^{\T}\tU_i^{\T}\frac{\yita_i}{\omega_{\ell}}
\omega_{\ell}}\bigg|
\theta_{\ell}^{\T}\Gama^{\T}\tU^{\T}11^{\T}\yita\bigg|
\leq\frac{1}{n^{1/2}\omega_{\ell}}\big|\theta_{\ell}^{\T}
\Gama^{\T}\tU^{\T}1\big|\frac{1}{n}|1^{\T}\yita|\\
\theta_{\ell}^{\T}\left(\Gama^{\T}\tU^{\T}-F^{\T}\right)
\bigg\}\frac{1}{n}\left|1^{\T}\yita\right|,
\end{align*}

where the first equality follows from the fact

\begin{align*}
\Gama^{\T}\tU_i^{\T}\frac{\yita_i}{\omega_{\ell}}
\theta_{\ell}^{\T}\Gama^{\T}\tU_i^{\T}\yita_i
\yita_i+\frac{1}{n}\theta_{\ell}^{\T}\Gama^{\T}
\tU^{\T}11^{\T}\yita\bigg|\\
\theta_{\ell}^{\T}\Gama^{\T}\tU^{\T}11^{\T}\yita\bigg|.
\end{align*}

Apply the Hölder's inequality to get

\begin{align*}
\Gama^{\T}\tU_i^{\T}\frac{\yita_i}{\omega_{\ell}}
\omega_{\ell}}\bigg|\\
\theta_\ell\right\|_{1}\max_{\ell}\bigg\|\sum_{j=1}^{q}
\Bigg)\frac{1}{n}\bigg|\sum_{i=1}^{n}\yita_i\bigg|\\
\frac{1}{n}\bigg|\sum_{i=1}^{n}\yita_i\bigg|\\
&\leq\Bigg\{m_{\Omega}\max_{\ell}\sum_{j\in J_{\ell}}
\sqrt{n}\sup_{z}\bigg|f_{j\ell}(z)-\sum_{k=1}^{m_n}
\phi_{k}(z)\bar{\gamma}_{kj\ell}\bigg|/\omega_{\ell}\\
&\qquad+m_{\Omega}\max_{\ell}\sum_{j\in J_{\ell}}\bigg|\sum_{i=1}^{n}
\frac{1}{n}\bigg|\sum_{i=1}^{n}\yita_i\bigg|\\
\bigg|/(n^{1/2}\omega_{\ell})\bigg\}\frac{1}{n}
\bigg|\sum_{i=1}^{n}\yita_i\bigg|.
\end{align*}

It then follows from the Gaussian tail probability that

\begin{align*}
\pr\bigg(\bigg|\frac{1}{n^{1/2}\sigma_0}\sum_{i=1}^{n}
\yita_i\bigg|\geq t \bigg)\leq 2\bigg(
\frac{2}{\pi}\bigg)^{1/2}\frac{\exp(-t^{2}/2)}{t}.
\end{align*}

Setting $t=2\{\log(n)\}^{1/2}$, we obtain with probability at least $1-n^{-1}$

\begin{align*}
\bigg|\frac{1}{n}\sum_{i=1}^{n}\yita_i\bigg|
\leq\sigma_0\bigg\{\frac{\log(n)}{n}\bigg\}^{1/2}.
\end{align*}

Similarly, by Lemma <ref>, we have

\begin{align*}
\bigg|\frac{1}{n}\sum_{i=1}^{n}f_{j\ell}(Z_{ij})
\bigg|\leq \big\{4C_0^{2}\log(pqm)/n\big\}^{1/2}
\end{align*}

with probability at least $1-2(pqm)^{-2}$. By the union bound argument, we

\begin{align*}
\left|\frac{1}{n^{1/2}}\sum_{i=1}^{n}\theta_{\ell}^{\T}
\Gama^{\T}\tU_i^{\T}\frac{\yita_i}{\omega_{\ell}}
\omega_{\ell}}\right|\leq \sigma_0\{\log(n)\}^{1/2}
\end{align*}

with probability at least $1-2(pqm)^{-1}-n^{-1}$. Since $\omega_{\ell}$ is
lower bounded by a constant when $n$ is large, we conclude

\begin{align*}
\bigg|\frac{1}{n^{1/2}}\sum_{i=1}^{n}\theta_{\ell}^{\T}
\Gama^{\T}\tU_i^{\T}\frac{\yita_i}{\omega_{\ell}}
\omega_{\ell}}\bigg|=o_{p}(1).
\end{align*}

It follows that

\begin{align*}
\bigg|T_{\ell,n}-n^{1/2}\frac{\tilde{\bata}_{\ell}
\end{align*}

Therefore, $T_{\ell,n}$ and
$\sqrt{n}(\tilde{\bata}_{\ell}-\bata_{\ell})/\omega_{\ell}$ share the same
weak limit. Now we show $T_{\ell,n}$ converges in distribution to the standard
normal. Note that

\begin{align*}
\Gama^{\T}\tU_i^{\T}\frac{\yita_i}{\omega_{\ell}}\\
\Gama^{\T}\tU_{i}^{\T}\yita_i}{\big\{\frac{1}{n}
\sum_{i=1}^{n}\theta_{\ell}^{\T}\Gama^{\T}\mathbb{E}
\big(\tU_i^{\T}\tU_i\big)\Gama\theta_{\ell}\big\}^{1/2}
\sigma_0}\frac{\big\{\frac{1}{n}\sum_{i=1}^{n}
\theta_{\ell}^{\T}\Gama^{\T}\mathbb{E}\big(\tU_i^{\T}
\tU_i\big)\Gama\theta_{\ell}\big\}^{1/2}\sigma_0}
\left(\Gama^{\T}\U^{\T}\U\Gama/n\right)
\theta_{\ell}\big\}^{1/2}}\\
&\qquad\times \frac{\sigma_0\big\{\theta_{\ell}^{\T}
\mathbb{E}\left(\Gama^{\T}\U^{\T}\U\Gama/n\right)
\theta_{\ell}\big\}^{1/2}}{\omega_\ell}.
\end{align*}

It follows that

\begin{align*}
\frac{\theta_{\ell}^{\T}\Gama^{\T}\tU_{i}^{\T}
\yita_i}{\big\{\frac{1}{n}\sum_{i=1}^{n}\theta_{\ell}^{\T}
\Gama^{\T}\mathbb{E}\big(\tU_i^{\T}\tU_i\big)
\Gama\theta_{\ell}\big\}^{1/2}\sigma_0}\leq t\bigg]\\
\sum_{i=1}^{n}\frac{\theta_{\ell}^{\T}\Gama^{\T}
\tU_{i}^{\T}\yita_i}{\big(\frac{1}{n}\sum_{i=1}^{n}
\theta_{\ell}^{\T}\Gama^{\T}\mathbb{E}\big(\tU_i^{\T}
\tU_i\big)\Gama\theta_{\ell}\big)^{1/2}\sigma_0}
\leq t\bigg|\Z \Bigg\}\Bigg]=\mathbb{E}\{\Phi(t)\}=\Phi(t)
\end{align*}

for all $t\in\mathbb{R}$, where $\Phi(t)$ denotes the CDF of the standard
normal distribution. Finally, we note that

\begin{align*}
\frac{\big\{\frac{1}{n}\sum_{i=1}^{n}\theta_{\ell}^{\T}
\Gama^{\T}\mathbb{E}\big(\tU_i^{\T}\tU_i\big)\Gama
\theta_{\ell}\big\}^{1/2}\sigma_0}{\sigma_0\big\{\theta_{\ell}^{\T}
\mathbb{E}\left(\Gama^{\T}\U^{\T}\U\Gama/n\right)
\theta_{\ell}\big\}^{1/2}}=\bigg(\frac{n}{n-1}\bigg)^{1/2},
\end{align*}

where we use the fact that
$\mathbb{E}(\U^{\T}\U/n)=\frac{n-1}{n}\mathbb{E}(\tU^{\T}\tU/n)$. Then, by the
triangular inequality, we have

\begin{align*}
\frac{\Gama^{\T}\tU^{\T}\tU\Gama}{n}\bigg)\bigg\|_{\infty}\\
\bigg\|_{\infty}+\bigg\|\frac{1}{n}\mathbb{E}
\bigg(\frac{F^{\T}F}{n}\bigg)\bigg\|_{\infty}\\
\frac{F^{\T}F}{n}-\frac{F^{\T}\tU\Gama}{n}
\tU^{\T}\tU\Gama}{n}\bigg)\bigg\|_{\infty}
\frac{F^{\T}F}{n}\bigg)\bigg\|_{\infty}\\
\frac{F^{\T}F}{n}-\frac{F^{\T}\tU\Gama}{n}
\bigg)\bigg\|_{\infty}+\frac{n-1}{n}
\bigg\|\mathbb{E}\bigg(\frac{F^{\T}\tU\Gama}{n}
\frac{F^{\T}F}{n}\bigg)\bigg\|_{\infty}.
\end{align*}

It then follows from the triangular and Jensen's inequalities that

\begin{align*}
\bigg\{\frac{\big(\sum_{j=1}^{q}F_{j\ell}^{\T}\big)
\big(\sum_{j=1}^{q}F_{j\ell'}-\tU\gama_{\ell'}
\big)}{n}\bigg\}\bigg|\\
\bigg\{\frac{\big(\tU\gama_{\ell}\big)^{\T}\big(
\sum_{j=1}^{q}F_{j\ell'}-\tU\gama_{\ell'}\big)}{n}\bigg\}\bigg|\\
\bigg\|\frac{\sum_{j=1}^{q}F_{j\ell}}{n^{1/2}}\bigg\|_{2}
\bigg\|\frac{\sum_{j=1}^{q}F_{j\ell'}}{n^{1/2}}\bigg\|_{2}\\
\bigg|\bigg\{\frac{\big(\sum_{j=1}^{q}F_{j\ell}^{\T}
\big)\big(\sum_{j=1}^{q}F_{j\ell'}-\tU\gama_{\ell'}
\big)}{n}\bigg\}\bigg|\\
\bigg|\frac{\big(\tU\gama_{\ell}\big)^{\T}
\big(\sum_{j=1}^{q}F_{j\ell'}-\tU\gama_{\ell'}\big)}{n}\bigg|\\
\bigg\|\frac{\sum_{j=1}^{q}F_{j\ell}}{n^{1/2}}
\bigg\|_{2}\bigg\|\frac{\sum_{j=1}^{q}
\end{align*}

Now we apply the Cauchy–Schwarz inequality to get

\begin{align*}
\frac{\big\|\sum_{j=1}^{q}F_{j\ell}\big\|_2}{n^{1/2}}
\bigg\|\frac{\sum_{j=1}^{q}F_{j\ell'}
\bigg\|\frac{\tU\gama_{\ell}-\sum_{j=1}^{q}F_{j\ell}
\frac{\sum_{j=1}^{q}F_{j\ell'}-\tU\gama_{\ell'}}
\frac{\sum_{j=1}^{q}F_{j\ell}}{n^{1/2}}\bigg\|_{2}
\bigg\|\frac{\sum_{j=1}^{q}F_{j\ell'}}{n^{1/2}}\bigg\|_{2}\\
\frac{\sum_{j=1}^{q}F_{j\ell'}-\tU\gama_{\ell'}}
\Bigg\{\bigg(\bigg\|\frac{\tU\gama_{\ell}
\bigg\|\bigg)\bigg\|\frac{\sum_{j=1}^{q}
\frac{\sum_{j=1}^{q}F_{j\ell}}{n^{1/2}}\bigg\|_{2}
\bigg\|\frac{\sum_{j=1}^{q}F_{j\ell'}}{n^{1/2}}\bigg\|_{2}\\
\sum_{j\in J_{\ell}}\bigg\|\frac{F_{j\ell}
\sum_{j\in J_{\ell}}\bigg\|\frac{\tU\gama_{\ell}
\end{align*}

Finally, we obtain

\begin{align*}
\phi_{k}(z)\bigg|+\frac{r^{2}C_0^{2}}{n}\\
\phi_{k}(z)\bigg|+rC_0\bigg\}\bigg\{r\sup_{z\in[a,b]}
\bigg|f_{j\ell}(z)-\sum_{k=1}^{m_n}
\bar{\gamma}_{jk\ell}\phi_{k}(z)\bigg|\bigg\},
\end{align*}

where we apply the Cauchy–Schwarz and triangular inequalities. Now by Lemma
<ref>, we have

\begin{align*}
\end{align*}

The above inequality holds when $C_Lm_n^{-d}\leq C_0$. It then follows that

\begin{align*}
\bigg)\theta_\ell-\theta_\ell^{\T}\mathbb{E}\bigg(
\frac{\Gama^{\T}\U^{\T}\U\Gama}{n}\bigg)\theta_\ell\bigg|
\leq\|\theta_\ell\|_1^{2}\bigg\|\mathbb{E}\bigg(
\frac{F^{\T}F}{n}-\frac{\Gama^{\T}\U^{\T}\U\Gama}{n}
\bigg)\bigg\|_{\infty}\\
\end{align*}

from which we know

\begin{align*}
\sigma_0\bigg\{\theta_\ell^{\T}\mathbb{E}\bigg(
\frac{\Gama^{\T}\U^{\T}\U\Gama}{n}\bigg)
\theta_\ell\bigg\}^{1/2}\rightarrow\omega_\ell.
\end{align*}

Next, we show, for a consistent estimator $\widehat{\omega}_{\ell}$,

\begin{align*}
\widehat{\omega}_\ell\rightsquigarrow\mathcal{N}(0,1).
\end{align*}

Since $\omega_{\ell}$ is lower bounded by a constant when $n$ is large, for a
consistent estimator $\widehat{\omega}_{\ell}$, we know
$\widehat{\omega}_\ell=\Omega_{p}(1)$. Therefore, we have

\begin{align*}
\left|\frac{\omega_\ell}{\widehat{\omega}_\ell}-1\right|=o_{p}(1).
\end{align*}

By Slutsky's theorem, we have

\begin{align*}
\frac{\omega_\ell}{\widehat{\omega}_\ell}\rightsquigarrow
\mathcal{N}(0,1).
\end{align*}

§.§ Proof of Lemma <ref>

Note that

\begin{align*}
\widehat{\omega}_\ell^{2}=\widehat{\sigma}_0^{2}
\widehat{\theta}_\ell^{\T}\frac{\widehat{\Gama}^{\T}
\U^{\T}\U\widehat{\Gama}}{n}\widehat{\theta}_\ell.
\end{align*}

We first prove $\widehat{\sigma}_0$ is a consistent estimator of $\sigma_0$ by
showing $\widehat{\sigma}_0^{2}\rightarrow\sigma_0^{2}$ in probability. Note

\begin{align*}
\big\|_2^{2}-\sigma_0^{2}\bigg|
\big\|\X\bata-\X\widehat{\bata}\big\|_2^{2}
\X\bata-\X\widehat{\bata}\big\|_2}{n^{1/2}}\\
\end{align*}

which holds with probability at least $1-234(pqm)^{-1}$ by applying Theorem
<ref>. From Theorem
<ref>, we have

\begin{align*}
\bigg|\frac{1}{n}\left\|\yita\right\|_2^{2}
-\sigma_0^{2}\bigg|\leq 4\sigma_0^{2}
\bigg\{\frac{\log(pqm)}{n}\bigg\}^{1/2}
\end{align*}

with probability at least $1-2(pqm)^{-2}$, which implies

\begin{align*}
\left|\widehat{\sigma}^{2}_0-\sigma_0^{2}\right|
\leq4\sigma_0^{2}\bigg\{\frac{\log(pqm)}{n}\bigg\}^{1/2}
\sigma_0\frac{8}{\kappa}s^{1/2}\mu
\end{align*}

holds with probability at least $1-236(pqm)^{-1}$ when $\log(pqm_n)/n\leq 1$.
Then we consider

\begin{align*}
\U^{\T}\U\widehat{\Gama}}{n}\widehat{\theta}_\ell
\U^{\T}\U\widehat{\Gama}}{n}\widehat{\theta}_\ell
\widehat{\theta}_\ell+\widehat{\theta}_{\ell}^{\T}
\frac{F^{\T}F}{n}\widehat{\theta}_\ell-\theta_\ell^{\T}
\frac{\Gama^{\T}\mathbb{E}(\U^{\T}\U)\Gama}{n}\theta_\ell\bigg|\\
\U^{\T}\U\widehat{\Gama}}{n}\widehat{\theta}_\ell
\widehat{\theta}_\ell\bigg|+\bigg|\theta_{\ell}^{\T}\frac{F^{\T}F}{n}
\theta_\ell-\theta_\ell^{\T}\frac{F^{\T}F}{n}\theta_\ell\\
\widehat{\theta}_\ell-\theta_\ell^{\T}\frac{\Gama^{\T}
\mathbb{E}(\U^{\T}\U)\Gama}{n}\theta_\ell\bigg|\\
\frac{\widehat{\Gama}^{\T}\U^{\T}\U\widehat{\Gama}}{n}
\frac{F^{\T}F}{n}\theta_{\ell}-\widehat{\theta}_\ell^{\T}
\frac{F^{\T}F}{n}\widehat{\theta}_\ell\bigg|\\
\Gama}{n}\theta_\ell\bigg|\\
\end{align*}

For $T_1$, we know that

\begin{align*}
\U^{\T}\U\widehat{\Gama}}{n}-\frac{F^{\T}F}{n}
\bigg\|_{\infty}\leq 35m_{\Omega}^{2}\lambda_{\max}rC_0
\bigg\{\frac{2rm}{\rho}\bigg\}^{1/2}
\end{align*}

with probability at least $1-60(pqm)^{-1}$ when
$\lambda_{\max}\{2m/\rho\}^{1/2}\leq\sqrt{r}C_0$, where we use the definition
of $\widehat{\theta}$ in the first inequality and Lemma
<ref> in the second inequality. For $T_2$, we have

\begin{align*}
T_2\leq \bigg|\theta_{\ell}^{\T}\frac{F^{\T}F}{n}\theta_{\ell}
\leq\left\|\theta_\ell\right\|_1\big\|\theta_\ell-\widehat{\theta}_{\ell}
\big\|_1\bigg\|\frac{F^{\T}F}{n}\bigg\|_{\infty}\leq
\end{align*}

where we use the fact that

\begin{align*}
\bigg|\bigg(\sum_{j=1}^{q}f_{j\ell}(Z_j)\bigg)
\bigg(\sum_{j=1}^{q}f_{j\ell}(Z_j)\bigg)\bigg|\leq r^{2}C_0^{2}.
\end{align*}

Applying Lemma <ref>, we have

\begin{align*}
T_2\leq 2c_bm_{\Omega}s_{\Omega}r^{2}C_0^{2}(2m_{\Omega}\upsilon)^{1-b}
\end{align*}

with probability at least $1-62(pqm)^{-1}$, where
$\upsilon=36m_{\Omega}\lambda_{\max}rC_0\{2rm/\rho\}^{1/2}, 0\leq b<1$. For
$T_3$, we have

\begin{align*}
\Gama}{n}\bigg)\theta_{\ell}\bigg|\\
\Gama}{n}\bigg\|_{\infty}\\
\end{align*}

which follows from Theorem <ref>. Therefore, we have

\begin{align*}
\U^{\T}\U\widehat{\Gama}}{n}\widehat{\theta}_\ell
\Gama}{n}\theta_\ell\bigg|\\
\end{align*}

Now we have

\begin{align*}
\frac{\widehat{\Gama}^{\T}\U^{\T}\U\widehat{\Gama}}{n}
\widehat{\theta}_\ell-\sigma_0^{2}\theta_{\ell}^{\T}
\mathbb{E}\left(\Gama^{\T}\U^{\T}\U\Gama/n
\right)\theta_{\ell}\bigg|\\
\frac{\widehat{\Gama}^{\T}\U^{\T}\U\widehat{\Gama}}{n}
\widehat{\theta}_\ell-\sigma_0^{2}\widehat{\theta}_\ell^{\T}
\frac{\widehat{\Gama}^{\T}\U^{\T}\U\widehat{\Gama}}{n}
\widehat{\theta}_\ell+\sigma_0^{2}\widehat{\theta}_\ell^{\T}
\frac{\widehat{\Gama}^{\T}\U^{\T}\U\widehat{\Gama}}{n}
\widehat{\theta}_\ell-\sigma_0^{2}\theta_{\ell}^{\T}
\mathbb{E}\left(\Gama^{\T}\U^{\T}\U\Gama/n
\right)\theta_{\ell}\bigg|\\
\widehat{\theta}_\ell^{\T}\frac{\widehat{\Gama}^{\T}\U^{\T}
\U\widehat{\Gama}}{n}\widehat{\theta}_\ell\bigg|
\U^{\T}\U\widehat{\Gama}}{n}\widehat{\theta}_\ell
\U\Gama/n\right)\theta_{\ell}\bigg|\sigma_0^{2}\\
\bigg\}^{1/2}+\frac{64}{\kappa^{2}}s\mu^{2}
\bigg|\widehat{\theta}_\ell^{\T}\frac{\widehat{\Gama}^{\T}
\U^{\T}\U\widehat{\Gama}}{n}\widehat{\theta}_{\ell}\bigg|\\
\bigg\{\frac{2rm}{\rho}\bigg\}^{1/2}
\end{align*}

We also have

\begin{align*}
\U^{\T}\U\widehat{\Gama}}{n}\widehat{\theta}_{\ell}\right|
\leq\left\|\widehat{\theta}_\ell\right\|_{1}^{2}
\left\|\frac{\widehat{\Gama}^{\T}\U^{\T}\U\widehat{\Gama}}{n}
\U^{\T}\U\widehat{\Gama}}{n}-\frac{F^{\T}F}{n}
\right\|_{\infty}+\left\|\frac{F^{\T}F}{n}
\right\|_{\infty}\right)\\
\frac{2rm}{\rho}\bigg\}^{1/2}+m_{\Omega}r^{2}C_0^{2}.
\end{align*}

Therefore, we get

\begin{align*}
\big|\widehat{\omega}^{2}_{\ell}-\omega^{2}_{\ell}\big|
\frac{8}{\kappa}s^{1/2}\mu\bigg]\\
\bigg\{\frac{2rm}{\rho}\bigg\}^{1/2}
\bigg\{\frac{2rm}{\rho}\bigg\}^{1/2}
\end{align*}

§ SECTION E: CONTROL OF REMAINDER TERMS

We derive specific conditions to make sure the remainder terms satisfy
$\|R_k\|_{\infty}=o_p(1)$, $k=1,2,3,4$. We first find the probability of the
event $\|\Omega\widehat{\Sigma}_d-I\|_{\infty}=o_p(1)$.

Suppose Assumption <ref> holds, then

\begin{align*}
\big\|\Omega\widehat{\Sigma}_d-I\big\|_{\infty}\leq
\end{align*}

holds with probability at least $1-62(pqm)^{-1}$, when
$\lambda_{\max}\{2m/\rho\}^{1/2}\leq r^{1/2}C_0$.

Note that

\begin{align*}
\U\widehat{\Gama}}{n}-\mathbb{E}\bigg(\frac{F^{\T}F}
\U\widehat{\Gama}}{n}-\mathbb{E}\bigg(
\frac{F^{\T}F}{n}\bigg)\bigg\|_{\infty}\\
&\leq m_{\Omega}\bigg\|\Sigma_{f}-\frac{F^{\T}F}{n}
\widehat{\X}}{n}\bigg\|_{\infty}\\
\bigg\|_{\infty}+\bigg\|\frac{F^{\T}F}{n}
\bigg\|_{\infty}\bigg).
\end{align*}

From Lemma <ref>, we have

\begin{align*}
\bigg\|\Sigma_{f}-\frac{F^{\T}F}{n}\bigg\|_{\infty}\leq
\end{align*}

with probability at least $1-2p^{2}(pqm)^{-4}$. Also, from Lemma
<ref>, we have

\begin{align*}
\bigg\|\frac{F^{\T}F}{n}-\frac{\widehat{\X}^{\T}\widehat{\X}}{n}
\bigg\|_{\infty}\leq 35\lambda_{\max}rC_0\bigg(\frac{2rm}{\rho}\bigg)
\end{align*}

with probability at least $1-60(pqm)^{-1}$, when
$\lambda_{\max}(2m/\rho)^{1/2}\leq r^{1/2}C_0$. Therefore, we have

\begin{align*}
\big\|\Omega\widehat{\Sigma}_d-I\big\|_{\infty}\leq
\end{align*}

with probability at least $1-62(pqm)^{-1}$ when
$\lambda_{\max}(2m/\rho)^{1/2}\leq r^{1/2}C_0$.

By Lemma <ref>, $\upsilon$ can be chosen as
$36m_{\Omega}\lambda_{\max}rC_0\sqrt{2rm/\rho}$. The next result provides a
probabilistic bound for $\|\U^{\T}\yita/n\|_{\infty}$.

\begin{align*}
\pr\bigg[\frac{\|\U^{\T}\yita\|_{\infty}}{n}\leq
\end{align*}

Write $\|\U^{\T}\yita/n\|_{\infty}$ as
\left|\U_{jk}^{\T}\yita\right|$, where $\U_{jk}$ is the $k$th column of spline
matrix $\U_{j}$. Recall the Gaussian tail bound

\begin{align*}
\pr\left(|X|\geq t \right)\leq2\bigg(
\frac{2}{\pi}\bigg)^{1/2}\frac{\exp(-t^{2}/2)}{t}.
\end{align*}

Now we have

\begin{align*}
&\pr\bigg(\frac{\|\U^{\T}\yita\|_{\infty}}{n}\geq t\bigg)
\leq{}qm_n\pr\bigg(\frac{\big|\U_{jk}^{\T}\yita\big|}{n}\geq{}t\bigg)\\
\big\{\sum_{i=1}^{n}\psi_{k}^{2}(Z_{ij})\big\}^{1/2}}\geq
\frac{nt}{\sigma_0\big\{\sum_{i=1}^{n}\psi_{k}^{2}(Z_{ij})\big\}^{1/2}}\Bigg]\\
\big\{\sum_{i=1}^{n}\psi_{k}^{2}(Z_{ij})\big\}^{1/2}}{nt}.
\end{align*}

Since $|\psi_{k}(Z_{ij})|\leq 2$, setting $t=4\sigma_0\{\log(pqm/n)\}^{1/2}$,
we obtain

\begin{align*}
\pr\bigg[\frac{\|\U^{\T}\yita\|_{\infty}}{n}\geq
\end{align*}

Suppose (i) $\|\Omega\|_{1}$ is bounded above by a constant
$m_{\Omega}<\infty$, and (ii) $\widehat{\Omega}$ is an estimate of $\Omega$
with rows $\widehat{\theta}_\ell$ obtained as solutions to
(<ref>). Then

\begin{align*}
\|\widehat{\Omega}-\Omega\|_{\infty}\leq 2m_{\Omega}\upsilon
\end{align*}

holds with probability at least $1-62(pqm)^{-1}$.

Under the conditions of this Lemma, we know that

\begin{align*}
\big\|\Omega\widehat{\Sigma}_d-I\big\|_{\infty}\leq
\upsilon,\ \big\|\widehat{\Omega}\widehat{\Sigma}_{d}
\end{align*}

Then, we have

\begin{align*}
\bigg(\frac{F^{\T}F}{n}\bigg)\bigg\}\Omega
\end{align*}

By Hölder's inequality, we have

\begin{align*}
\widehat{\Omega}\big(I-\widehat{\Sigma}_d\Omega\big)
\leq\big\|I-\widehat{\Sigma}_d\Omega\big\|_{\infty}
\big\|\widehat{\Omega}^{\T}\big\|_{1}\leq
\big\|I-\widehat{\Sigma}_d\Omega\big\|_{\infty}
\big\|\Omega^{\T}\big\|_{1}\leq m_{\Omega}\upsilon
\end{align*}

\begin{align*}
\big\|\big(I-\widehat{\Omega}\widehat{\Sigma}_d\big)
\Omega\big\|_{\infty}\leq\|\Omega\|_1
\big\|I-\widehat{\Omega}\widehat{\Sigma}_d\big\|_{\infty}
\leq{}m_\Omega\upsilon.
\end{align*}

Therefore, we have

\begin{align*}
\big\|\Omega-\widehat{\Omega}\big\|_{\infty}\leq 2m_{\Omega}\upsilon.
\end{align*}

We have shown in Lemma <ref> that

\begin{align*}
\big\|\Omega\widehat{\Sigma}_d-I\big\|_{\infty}\leq
\end{align*}

with probability at least $1-62(pqm)^{-1}$. Therefore, the proof is complete
by choosing $\upsilon=36m_{\Omega}\lambda_{\max}rC_0(2rm/\rho)^{1/2}$.

Suppose, in addition to the conditions of Lemma <ref>,
$\Omega\in\mathcal{U}(m_\Omega,b,s_{\Omega})$. Then,

\begin{align*}
\|\widehat{\theta}_\ell-\theta_\ell\|_1\leq
\end{align*}

for each $\ell\in\{1,\ldots,p\}$, where $c_b=1+2^{1-b}+3^{1-b}$.

See the proof of Theorem 6 in [Cai et al., 2011].

Now we examine each remainder term $R_k$ in the subsequent subsections.

§.§ Control of $R_1$

Suppose the conditions of Theorem <ref> and Lemma
<ref> hold. If

\begin{align*}
\sigma_0rC_Lm_n^{-d}+14C_0\sigma_0\bigg\{\frac{r^{2}\log(pqm)}
\end{align*}

then $\|R_1\|_{\infty}=o(1)$ with probability at least $1-5(pqm)^{-1}$.

Note that

\begin{align*}
\|(\Gama^{\T}\U^{\T}-F^{\T})\yita/n^{1/2}\|_{\infty}
\big(\|(\Gama^{\T}\U^{\T}-F^{\T})\yita/n\|_{\infty}
\end{align*}

As is shown in the proof of Theorem <ref>,
$\|F^{\T}\yita/n\|_{\infty}\leq2C_0\sigma_0\{r^{2}\log(pqm)/n\}^{1/2}$ holds
with probability at least $1-p(pqm)^{-2}$. Applying Lemma
<ref> and the $\chi^{2}$ concentration bound, we

\begin{align*}
\leq\frac{1}{n^{1/2}}\max_{\ell}\bigg\|\U
\bar{\gama}_{\ell}-\sum_{j=1}^{q}F_{j\ell}\bigg\|_{2}
\max_{\ell'}\frac{\|\yita\|_2}{n^{1/2}}\\
\end{align*}

with probability at least $1-4r(pqm)^{-2}$. It follows that

\begin{align*}
\bigg[4\sqrt{5}\sigma_0rC_Lm_n^{-d}+14C_0\sigma_0
\bigg\{\frac{r^{2}\log(pqm)}{n}\bigg\}^{1/2}\bigg]
\end{align*}

holds with probability at least $1-5(pqm)^{-1}$.

§.§ Control of $R_2$

Under the same conditions of Lemma <ref>, if

\begin{align*}
\frac{m_n^{3/2}}{\rho}=o(1),
\end{align*}

then $\|R_2\|_{\infty}=o(1)$ with probability at least $1-19(pqm)^{-1}$.

Note that

\begin{align*}
\U^{\T}-\Gama^{\T}\U^{\T}\big)\yita/n^{1/2}\\
\|\widehat{\Gama}-\Gama\|_1\|\U^{\T}\yita/n\|_{\infty}\\
\|\widehat{\gama}_\ell-\bar{\gama}_\ell\|_{1}
\|\U^{\T}\yita/n\|_{\infty}.
\end{align*}

Recall from the proof of Theorem <ref>

\begin{align*}
\max_{\ell=1,\ldots,p}\|\widehat{\gama}_\ell-\bar{\gama}_\ell\|_{1}
\leq\max_{\ell}m^{1/2}\sum_{j=1}^{q}\|\widehat{\gama}_{j\ell}
-\bar{\gama}_{j\ell}\|_{2}\leq{}m^{1/2} 32r\lambda_{\max}\frac{m}{\rho}
\end{align*}

holds with probability at least $1-18(pqm)^{-1}$. Then, by Lemma
<ref>, we obtain

\begin{align*}
\|\U^{\T}\yita/n\|_{\infty}\leq4\bigg\{\frac{\log(pqm)}{n}\bigg\}^{1/2}
\end{align*}

with probability at least $1-(pqm)^{-1}$. It follows that

\begin{align*}
\|R_2\|_{\infty}\leq{}m_{\Omega}128r\lambda_{\max}
\{\log(pqm)\}^{1/2}\frac{m_n^{3/2}}{\rho}
\end{align*}

holds with with probability at least $1-19(pqm)^{-1}$.

§.§ Control of $R_3$

Under the same conditions of Lemma <ref>, if

\begin{align*}
\big(30rC_0+16\max_{\ell}\sigma_\ell\big)\frac{64}{\kappa^{2}}s\mu=o(1),
\end{align*}

then $\|R_3\|_{\infty}=o(1)$ with probability at least $1-277(pqm)^{-1}$.

We first apply Hölder's inequality to get

\begin{align*}
\|R_3\|_{\infty}\leq n^{1/2}\|\widehat{\Omega}^{\T}\|_1
\|\widehat{\D}^{\T}(\X-\widehat{\D})/n\|_{\infty}
\|\widehat{\bata}-\bata\|_1.
\end{align*}

Note that

\begin{align*}
\end{align*}

For the first term on the right-hand side above, we have

\begin{align*}
\leq\|F^{\T}(F-\widehat{\D})/n\|_{\infty}+\|(\widehat{\Gama}^{\T}
\U^{\T}-F^{\T})(F-\widehat{\D})/n\|_{\infty}\\
\big\|_{2}}{n^{1/2}}\bigg\|\frac{\sum_{j=1}^{q}F_{j\ell}
\leq30\lambda_{\max}r^{3/2}C_0\bigg(\frac{2m}{\rho}\bigg)^{1/2},
\end{align*}

which holds with probability at least $1-20(pqm)^{-1}$ when
$\lambda_{\max}r^{1/2}(2m/\rho)^{1/2}\leq C_0$. For the second term, we have

\begin{align*}
\frac{\varepsilon_{\ell'}}{n^{1/2}}\bigg\|_2+\|F^{\T}\E/n\|_{\infty}\\
& \leq 5\sqrt{5}\lambda_{\max}\max_{\ell'}\sigma_{\ell'}r^{1/2}
\bigg(\frac{2m}{\rho}\bigg)^{1/2}+\|F^{\T}\E/n\|_{\infty}\\
\bigg(\frac{2m}{\rho}\bigg)^{1/2}+\max_{\ell'}
\sigma_{\ell'}C_0\{8r^{2}\log(pqm)/n\}^{1/2}\\
\bigg(\frac{2m}{\rho}\bigg)^{1/2}
\end{align*}

with probability at least $1-23(pqm)^{-1}$, where we use the arguments in the
proof of Theorem <ref> when bounding
$\|T_3\|_{\infty}$ and $\|T_4\|_{\infty}$. It follows that

\begin{align*}
\big\|\widehat{\D}^{\T}(\X-\widehat{\D})/n\big\|_{\infty}\leq
\lambda_{\max}\bigg(\frac{2rm}{\rho}\bigg)^{1/2}
\big(30rC_0+16\max_{\ell}\sigma_\ell\big)
\end{align*}

holds with probability at leasst $1-43(pqm)^{-1}$. Finally, we have

\begin{align*}
&\|R_3\|_{\infty}\leq n^{1/2}\|\Omega^{\T}\|_1\|\widehat{\D}^{\T}
\big(30rC_0+16\max_{\ell}\sigma_\ell\big)\frac{64}{\kappa^{2}}s\mu
\end{align*}

with probability at least $1-277(pqm)^{-1}$.

§.§ Control of $R_4$

Under the same conditions of Lemma <ref>, if

\begin{align*}
\end{align*}

then $\|R_4\|_{\infty}=o(1)$ with probability at least $1-296(pqm)^{-1}$.

By Hölder's inequality, we have

\begin{align*}
\|R_4\|_{\infty}\leq n^{1/2}\|\widehat{\Omega}
\widehat{\Sigma}_d-I\|_\infty\|\widehat{\bata}-\bata\|_1.
\end{align*}

Applying Lemma <ref>, we obtain

\begin{align*}
\|R_4\|_{\infty}\leq n^{1/2}36m_{\Omega}\lambda_{\max}rC_{0}
\bigg(\frac{2rm}{\rho}\bigg)^{1/2}\|\widehat{\bata}-\bata\|_1,
\end{align*}

which holds with probability at least $1-62(pqm)^{-1}$ when
$\lambda_{\max}(2m/\rho)^{1/2}\leq r^{1/2}C_0$. Apply Theorem
<ref> to get

\begin{align*}
\|R_4\|_{\infty}\leq n^{1/2}36m_{\Omega}\lambda_{\max}rC_{0}
\bigg(\frac{2rm}{\rho}\bigg)^{1/2}\frac{64}{\kappa^{2}}s\mu,
\end{align*}

which holds with probability at least $1-296(pqm)^{-1}$.

§ SECTION F: ADDITIONAL EXPERIMENTS

In this section, we present more experiments to demonstrate the finite sample
performance of our proposed estimator when the dimensions of the treatment and
instrumental variables are large.

$L_1$ error of each method averaged over one hundred replications with standard
deviation shown in parentheses for $p=100$.

Sample     3cLinear     3cNonlinear(lr)2-4(lr)5-7
size     Our method     2SLS-L     PLS     Our method     2SLS-L     PLS

100     1.15 (0.42)    1.48 (0.62)     0.75 (0.27)     1.94 (0.98)    1.55 (0.74)     0.94 (0.31)

200     0.62 (0.33)    0.62 (0.32)    0.69 (0.28)     1.06 (0.53)    1.87 (0.81)     0.84 (0.32)

400    0.47 (0.26)     0.41 (0.21)    0.82 (0.34)    0.65 (0.29)     2.55 (0.95)    1.06 (0.28)

600    0.38 (0.18)    0.35 (0.16)    1.20 (0.51)     0.45 (0.26)    2.95 (1.11)    1.27 (0.26)

800    0.29 (0.14)    0.28 (0.13)    1.36 (0.42)     0.38 (0.17)    3.22 (1.52)    1.38 (0.29)

1000    0.23 (0.10)     0.23 (0.10)    1.44 (0.38)    0.30 (0.13)    3.49 (1.52)    1.43 (0.26)

1200    0.23 (0.10)    0.24 (0.10)    1.68 (0.37)    0.27 (0.15)    3.44 (1.28)    1.46 (0.21)

1400    0.19 (0.08)    0.20 (0.09)    1.75 (0.38)    0.24 (0.13)    3.83 (1.75)    1.55 (0.21)

$L_1$ error of each method averaged over one hundred replications with
standard deviation shown in parentheses for $p=200$.

Sample     3cLinear     3cNonlinear(lr)2-4(lr)5-7
size     Our method     2SLS-L     PLS     Our method     2SLS-lL     PLS

100     1.14 (0.52)    0.52 (0.76)     0.80 (0.23)     2.53 (1.47)    1.21 (0.50)     1.02 (0.31)

200     0.65 (0.28)    0.74 (0.38)    0.74 (0.38)     0.99 (0.41)    1.35 (0.80)     0.84 (0.29)

400    0.42 (0.20)     0.42 (0.20)    0.64 (0.23)     0.62 (0.24)     2.03 (1.02)    0.93 (0.29)

600    0.30 (0.16)    0.28 (0.14)    0.73 (0.31)     0.43 (0.19)    2.49 (1.16)    1.09 (0.27)

800    0.26 (0.11)    0.25 (0.10)    0.92 (0.35)     0.35 (0.15)    3.42 (1.75)    1.25 (0.24)

1000    0.21 (0.09)     0.20 (0.08)    1.14 (0.43)    0.30 (0.15)    3.42 (1.65)    1.35 (0.22)

1200    0.21 (0.10)    0.21 (0.10)    1.36 (0.47)    0.29 (0.15)    3.92 (2.32)    1.39 (0.24)

1400    0.18 (0.07)    0.19 (0.08)    1.58 (0.47)    0.24 (0.12)    4.12 (2.43)    1.47 (0.21)

$L_1$ error of each method averaged over one hundred replications with
standard deviation shown in the parentheses when $p=400$.

Sample     3cLinear     3cNonlinear(lr)2-4(lr)5-7
size     Our method     2SLS-L     PLS     Our method     2SLS-L     PLS

100     1.23 (0.68)    2.05 (0.86)     0.91 (0.34)     2.69 (1.24)    1.08 (0.45)     1.14 (0.36)

200     0.74 (0.60)    1.00 (0.50)    0.57 (0.21)     1.11 (0.64)    1.11 (0.63)     0.83 (0.29)

400    0.39 (0.16)     0.36 (0.15)    0.51 (0.18)     0.57 (0.26)     1.62 (0.98)    0.86 (0.28)

600    0.35 (0.15)    0.29 (0.11)    0.55 (0.21)     0.43 (0.17)    1.75 (1.04)    1.00 (0.22)

800    0.27 (0.13)    0.23 (0.11)    0.56 (0.19)     0.37 (0.18)    2.81 (1.50)    1.14 (0.24)

1000    0.24 (0.11)     0.23 (0.09)    0.56 (0.19)    0.31 (0.12)    3.20 (2.11)    1.19 (0.22)

1200    0.20 (0.08)    0.20 (0.08)    0.84 (0.37)    0.23 (0.09)    3.47 (1.90)    1.29 (0.22)

1400    0.19 (0.07)    0.18 (0.08)    0.91 (0.41)    0.22 (0.10)    4.12 (2.98)    1.36 (0.20)

[Anderson and Rubin, 1949]
Theodore W Anderson and Herman Rubin.
Estimation of the parameters of a single equation in a complete
system of stochastic equations.
The Annals of Mathematical Statistics, 200
(1):0 46–63, 1949.

[Belloni et al., 2012]
Alexandre Belloni, Daniel Chen, Victor Chernozhukov, and Christian Hansen.
Sparse models and methods for optimal instruments with an application
to eminent domain.
Econometrica, 800 (6):0 2369–2429, 2012.

[Belloni et al., 2022]
Alexandre Belloni, Christian Hansen, and Whitney Newey.
High-dimensional linear models with many endogenous variables.
Journal of Econometrics, 2022.

[Bickel et al., 2009]
Peter J Bickel, Ya’acov Ritov, and Alexandre B Tsybakov.
Simultaneous analysis of Lasso and Dantzig selector.
The Annals of Statistics, 370 (4):0
1705–1732, 2009.

[Burgess et al., 2017]
Stephen Burgess, Dylan S Small, and Simon G Thompson.
A review of instrumental variable estimators for Mendelian
Statistical Methods in Medical Research, 260
(5):0 2333–2355, 2017.

[Cai et al., 2011]
Tony Cai, Weidong Liu, and Xi Luo.
A constrained $\ell_1$ minimization approach to sparse precision
matrix estimation.
Journal of the American Statistical Association, 1060
(494):0 594–607, 2011.

[Cheng and Liao, 2015]
Xu Cheng and Zhipeng Liao.
Select the valid and relevant moments: An information-based LASSO
for GMM with many moments.
Journal of Econometrics, 1860 (2):0 443–464,

[Davey Smith and Ebrahim, 2003]
George Davey Smith and Shah Ebrahim.
“Mendelian randomization": can genetic epidemiology contribute to
understanding environmental determinants of disease?
International Journal of Epidemiology, 320
(1):0 1–22, 2003.

[Didelez and Sheehan, 2007]
Vanessa Didelez and Nuala Sheehan.
Mendelian randomization as an instrumental variable approach to
causal inference.
Statistical Methods in Medical Research, 160
(4):0 309–330, 2007.

[Fan and Liao, 2014]
Jianqing Fan and Yuan Liao.
Endogeneity in high dimensions.
The Annals of Statistics, 420 (3):0 872–917,

[Fan and Zhong, 2018]
Qingliang Fan and Wei Zhong.
Nonparametric additive instrumental variable estimator: A group
shrinkage estimation perspective.
Journal of Business & Economic Statistics, 360
(3):0 388–399, 2018.

[Fan et al., 2015]
Yingying Fan, Gareth M James, Peter Radchenko, et al.
Functional additive regression.
The Annals of Statistics, 430 (5):0
2296–2325, 2015.

[Gai et al., 2020]
Zhibo Gai, Ting Gui, Irina Alecu, Museer A Lone, Thorsten Hornemann, Qingfa
Chen, Michele Visentin, Christian Hiller, Stephanie Hausler, and Gerd A
Farnesoid X receptor activation induces the degradation of
hepatotoxic 1-deoxysphingolipids in non-alcoholic fatty liver disease.
Liver International, 400 (4):0 844–859, 2020.

[Gautier and Rose, 2011]
Eric Gautier and Christiern Rose.
High-dimensional instrumental variables regression and confidence
arXiv preprint arXiv:1105.2454, 2011.

[Gill et al., 2010]
Ryan Gill, Somnath Datta, and Susmita Datta.
A statistical framework for differential network analysis from
microarray data.
BMC Bioinformatics, 110 (1):0 1–10, 2010.

[Gold et al., 2020]
David Gold, Johannes Lederer, and Jing Tao.
Inference for high-dimensional instrumental variables regression.
Journal of Econometrics, 2170 (1):0 79–111,

[Guha Majumdar et al., 2020]
Sayanti Guha Majumdar, Anil Rai, and Dwijesh C Mishra.
Integrated framework for selection of additive and nonadditive
genetic markers for genomic selection.
Journal of Computational Biology, 270 (6):0
845–855, 2020.

[Hansen and Kozbur, 2014]
Christian Hansen and Damian Kozbur.
Instrumental variables estimation with many weak instruments using
regularized JIVE.
Journal of Econometrics, 1820 (2):0 290–308,

[Hardwick et al., 2009]
James P Hardwick, Douglas Osei-Hyiaman, Homer Wiland, Mohamed A Abdelmegeed,
and Byoung-Joon Song.
PPAR/RXR regulation of fatty acid metabolism and fatty acid
$\omega$-hydroxylase (CYP4) isozymes: implications for prevention of
lipotoxicity in fatty liver disease.
PPAR Research, 2009, 2009.

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

[Horowitz and Mammen, 2004]
Joel L Horowitz and Enno Mammen.
Nonparametric estimation of an additive model with a link function.
The Annals of Statistics, 320 (6):0
2412–2443, 2004.

[Huang et al., 2010]
Jian Huang, Joel L. Horowitz, and Fengrong Wei.
Variable selection in nonparametric additive models.
The Annals of Statistics, 380 (4):0
2282–2313, 08 2010.
URL <https://doi.org/10.1214/09-AOS781>.

[Javanmard and Montanari, 2014]
Adel Javanmard and Andrea Montanari.
Confidence intervals and hypothesis testing for high-dimensional
The Journal of Machine Learning Research, 150
(1):0 2869–2909, 2014.

[Li et al., 2014]
Jiahan Li, Jun Dan, Chunlei Li, and Rongling Wu.
A model-free approach for detecting interactions in genetic
association studies.
Briefings in Bioinformatics, 150 (6):0
1057–1068, 2014.

[Lian, 2012]
Heng Lian.
Semiparametric estimation of additive quantile regression models by
two-fold penalty.
Journal of Business & Economic Statistics, 300
(3):0 337–350, 2012.

[Lin et al., 2015]
Wei Lin, Rui Feng, and Hongzhe Li.
Regularization methods for high-dimensional instrumental variables
regression with an application to genetical genomics.
Journal of the American Statistical Association, 1100
(509):0 270–288, 2015.

[Lounici et al., 2011]
Karim Lounici, Massimiliano Pontil, Sara van de Geer, Alexandre B Tsybakov,
et al.
Oracle inequalities and optimal inference under group sparsity.
The Annals of Statistics, 390 (4):0
2164–2204, 2011.

[Lu et al., 2020]
Junwei Lu, Mladen Kolar, and Han Liu.
Kernel meets sieve: Post-regularization confidence bands for sparse
additive model.
Journal of the American Statistical Association, 1150
(532):0 2084–2099, 2020.

[Lv et al., 2018]
Shaogao Lv, Jiakun Jiang, Fanyin Zhou, Jian Huang, and Huazhen Lin.
Estimating high-dimensional additive Cox model with time-dependent
covariate processes.
Scandinavian Journal of Statistics, 450 (4):0
900–922, 2018.

[Meinshausen and Bühlmann, 2010]
Nicolai Meinshausen and Peter Bühlmann.
Stability selection.
Journal of the Royal Statistical Society: Series B (Statistical
Methodology), 720 (4):0 417–473, 2010.

[Newey, 1990]
Whitney K Newey.
Efficient instrumental variables estimation of nonlinear models.
Econometrica, pages 809–837, 1990.

[Radchenko, 2015]
Peter Radchenko.
High dimensional single index models.
Journal of Multivariate Analysis, 139:0 266–282,

[Raskutti et al., 2010]
Garvesh Raskutti, Martin J Wainwright, and Bin Yu.
Restricted eigenvalue properties for correlated Gaussian designs.
The Journal of Machine Learning Research, 11:0
2241–2259, 2010.

[Ravikumar et al., 2009]
Pradeep Ravikumar, John Lafferty, Han Liu, and Larry Wasserman.
Sparse additive models.
Journal of the Royal Statistical Society: Series B (Statistical
Methodology), 710 (5):0 1009–1030, 2009.

[Schumaker, 2007]
Larry Schumaker.
Spline Functions: Basic Theory.
Cambridge University Press, 2007.

[Shen et al., 1998]
Xiaotong Shen, DA Wolfe, and S Zhou.
Local asymptotics for regression splines and confidence regions.
The Annals of Statistics, 260 (5):0
1760–1782, 1998.

[Stone, 1985]
Charles J Stone.
Additive regression and other nonparametric models.
The Annals of Statistics, pages 689–705, 1985.

[van de Geer et al., 2014]
Sara van de Geer, Peter Bühlmann, Ya’acov Ritov, Ruben Dezeure, et al.
On asymptotically optimal confidence regions and tests for
high-dimensional models.
The Annals of Statistics, 420 (3):0
1166–1202, 2014.

[Wainwright, 2019]
Martin J. Wainwright.
High-Dimensional Statistics: A Non-Asymptotic Viewpoint.
Cambridge Series in Statistical and Probabilistic Mathematics.
Cambridge University Press, 2019.

[Wang et al., 2006]
Susanna Wang, Nadir Yehya, Eric E Schadt, Hui Wang, Thomas A Drake, and
Aldons J Lusis.
Genetic and genomic analysis of a fat mass trait with complex
inheritance reveals marked sex specificity.
PLoS Genetics, 20 (2):0 e15, 2006.

[Wang et al., 2015]
Xuefeng Wang, Eric P Xing, and Daniel J Schaid.
Kernel methods for large-scale genomic data analysis.
Briefings in Bioinformatics, 160 (2):0
183–192, 2015.

[Wang et al., 2019]
Yong Wang, Yun-Sheng Cheng, Xiao-Qiang Yin, Gang Yu, and Ben-Li Jia.
Anxa2 gene silencing attenuates obesity-induced insulin resistance by
suppressing the nf-$\kappa$b signaling pathway.
American Journal of Physiology-Cell Physiology, 3160
(2):0 C223–C234, 2019.

[Wright, 1928]
Philip G Wright.
Tariff on Animal and Vegetable Oils.
Macmillan, New York, 1928.

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

[Yuan and Lin, 2006]
Ming Yuan and Yi Lin.
Model selection and estimation in regression with grouped variables.
Journal of the Royal Statistical Society: Series B (Statistical
Methodology), 680 (1):0 49–67, 2006.

[Zhan et al., 2017]
Xiang Zhan, Anna Plantinga, Ni Zhao, and Michael C Wu.
A fast small-sample kernel independence test for microbiome
community-level association analysis.
Biometrics, 730 (4):0 1453–1463, 2017.

[Zhang and Zhang, 2014]
Cun-Hui Zhang and Stephanie S Zhang.
Confidence intervals for low dimensional parameters in high
dimensional linear models.
Journal of the Royal Statistical Society: Series B: Statistical
Methodology, pages 217–242, 2014.

[Zhang and Ghosh, 2017]
Weiming Zhang and Debashis Ghosh.
On the use of kernel machines for Mendelian randomization.
Quantitative Biology, 50 (4):0 368–379, 2017.

[Zhao et al., 2019]
Ni Zhao, Haoyu Zhang, Jennifer J Clark, Arnab Maity, and Michael C Wu.
Composite kernel machine regression based on likelihood ratio test
for joint testing of genetic and gene-environment interaction effect.
Biometrics, 750 (2):0 625–637, 2019.

[Zhu, 2018]
Ying Zhu.
Sparse linear models and $l_1$-regularized 2SLS with
high-dimensional endogenous regressors and instruments.
Journal of Econometrics, 2020 (2):0 196–213,