# High-dimensional variable selection with
heterogeneous signals: A precise asymptotic perspective

Saptarshi Roy  Ambuj Tewari  Ziwei Zhu

###### Abstract

We study the problem of exact support recovery for high-dimensional sparse
linear regression when the signals are weak, rare and possibly heterogeneous.
Specifically, we fix the minimum signal magnitude at the information-theoretic
optimal rate and investigate the asymptotic selection accuracy of best subset
selection (BSS) and marginal screening (MS) procedures under independent
Gaussian design. Despite of the ideal setup, somewhat surprisingly, marginal
screening can fail to achieve exact recovery with probability converging to
one in the presence of heterogeneous signals, whereas BSS enjoys model
consistency whenever the minimum signal strength is above the information-
theoretic threshold. To mitigate the computational issue of BSS, we also
propose a surrogate two-stage algorithm called ETS (Estimate Then Screen)
based on iterative hard thresholding and gradient coordinate screening, and we
show that ETS shares exactly the same asymptotic optimality in terms of exact
recovery as BSS. Finally, we present a simulation study comparing ETS with
LASSO and marginal screening. The numerical results echo with our asymptotic
theory even for realistic values of the sample size, dimension and sparsity.

Keywords: Variable Selection; High-Dimensional Statistics; Model Consistency;
Heterogeneous Signals; Marginal Screening; Iterative Hard Thresholding.

##  1 Introduction

Consider $n$ independent observations $(x_{i},y_{i})_{i\in[n]}$ of a random
pair $(x,y)$ drawn from the following linear regression model:

$\begin{split}(x,\varepsilon)\sim\operatorname{\mathcal{P}}_{x}\times\operatorname{\mathcal{P}}_{\varepsilon},\\\
y=x^{\top}\beta+\varepsilon,\end{split}$ (1)

where $\operatorname{\mathcal{P}}_{x}$ is the $p$-dimensional isotropic
Gaussian distribution ${\sf N}_{p}(0,{\sf I}_{p})$, and
$\operatorname{\mathcal{P}}_{\varepsilon}$ is the standard Gaussian
distribution on ${\mathbb{R}}$. In matrix notation, the observations can be
represented as

$Y=X\beta+E,$

where
$Y=(y_{1},y_{2},\ldots,y_{n})^{\top},X=(x_{1},x_{2},\ldots,x_{n})^{\top}$ and
$E=(\varepsilon_{1},\varepsilon_{2},\ldots,\varepsilon_{n})^{\top}$. The
vector $\beta$ is unknown but sparse in the sense that
$\left\|\beta\right\|_{0}:=\sum_{j=1}^{p}\operatorname{\mathbbm{1}}(\beta_{j}\neq
0)=s$, which is much smaller than $p$. Denote by
$\operatorname{\mathcal{S}}(v)$ the set of non-zero coordinates of a vector
$v\in{\mathbb{R}}^{p}$. Lastly we denote by
$\operatorname{\mathbb{P}}_{\beta_{0}}(\cdot)$ and
$\mathbb{E}_{\beta_{0}}(\cdot)$ the probability measure and the expectation
with $\beta=\beta_{0}$ respectively. In this paper, we focus on the variable
selection problem, i.e., identifying $\operatorname{\mathcal{S}}(\beta)$. We
primarily use the 0-1 loss, i.e.,
$\operatorname{\mathbb{P}}_{\beta}(\operatorname{\mathcal{S}}(\hat{\beta})\neq\operatorname{\mathcal{S}}(\beta))$,
to assess the quality of the selected model
$\operatorname{\mathcal{S}}(\hat{\beta})$. The isotropic Gaussian design has
been widely used to conduct precise analysis of variable selection procedures
(Fletcher et al.,, 2009; Genovese et al.,, 2012; Ndaoud and Tsybakov,, 2020;
Su et al.,, 2017). Specifically, these works either derived the necessary and
sufficient condition for exact model recovery (Fletcher et al.,, 2009), or
established tight asymptotic bounds of model selection error (Genovese et
al.,, 2012; Ndaoud and Tsybakov,, 2020; Su et al.,, 2017). The isotropic
Gaussian design is also used in compressed sensing to generate a measurement
matrix (Candès et al.,, 2006; Candes and Tao,, 2006; Donoho,, 2006), so that
one can sense the sparse signals with few measurements of the high-dimensional
signal.

Recently there has been growing interest in the variable selection problem in
the presence of weak and rare signals (Genovese et al.,, 2012; Ji et al.,,
2012). This regime is ubiquitous in modern data analytics such as those in
Genome-Wide Association Study (GWAS). There the genes that exhibit detectable
association with the trait of interest can be extremely few with weak effects
(Consortium et al.,, 2007; Marttinen et al.,, 2013). Moreover, the number of
subjects $n$ typically ranges in thousands, while the number of features $p$
can range from tens of thousands to hundreds of thousands. Such a high
dimension further adds to the difficulty of identifying the weak signals. Weak
and rare signals also arise in multi-user detection problems (Arias-Castro et
al.,, 2011) where one typically uses linear model of the form (1). There the
$j$th column of $X$, denoted by $X_{j}$, is the channel impulse response for
user $j$. The signal received from user $j$ is $\beta_{j}X_{j}$. Thus
$\beta_{j}=0$ means that $j$th user is not sending any signal. It is a common
practice to model the mixing matrix $X$ as random with i.i.d. entries. Under
the presence of strong noise, one might be interested in knowing whether
information is being transmitted or not. Typically, in some applications it is
reasonable to assume that a very few numbers of users are sending signals.
Also due to strong noise environment the signals become quite weak, making
them harder to detect. Therefore, from an application point of view,
understanding variable selection in weak and rare signal regimes is crucial.
Despite its importance, the problem remains fairly underexplored. Typically
most of the popular methods including LASSO (Tibshirani,, 1996), SCAD (Fan and
Li,, 2001), MC+ (Zhang,, 2010), etc. have been extensively analyzed in terms
of 0-1 loss when the signals are uniformly strong (Zhao and Yu,, 2006; Guo et
al.,, 2015; Zhang,, 2010) in the sense that

$\min_{j\in\operatorname{\mathcal{S}}(\beta)}\left|\beta_{j}\right|\gg\left(\frac{\log
p}{n}\right)^{1/2}.$

However, very little is known about the variable selection problem in the
presence of weak and rare signals.

Besides the weakness and rarity of signals, heterogeneity in the signal
strength is another important feature of modern data applications that has not
yet received sufficient attention. One limitation of the existing literature
on variable selection in the weak and rare signals regime is that it typically
assumes that the true signals are homogeneous (Genovese et al.,, 2012; Ji et
al.,, 2012; Jin et al.,, 2014). Ji et al., (2012) refer to this setup as the
Asymptotically Rare and Weak (ARW) signal regime. Many popular approaches have
been shown to enjoy satisfactory variable selection properties under the ARW
regime. For instance, Genovese et al., (2012) showed that both LASSO and
marginal screening enjoy model consistency in terms of Hamming loss under
independent random design. Ji et al., (2012) and Jin et al., (2014)
investigated the same problem under sparsely correlated design. They proposed
two-stage screen and clean algorithms that also exhibit model consistency in
terms of Hamming loss. However, their theory heavily relies on homogeneous
signals and does not extend to the heterogeneous case that is of interest to
us. In reality, the ARW setup seldom occurs: the signals almost always have
different strengths (Li et al.,, 2019). It is not yet fully understood whether
one can achieve exact support recovery in the presence of weak, rare and
heterogeneous signals.

To underscore the contrasting effects of homogeneous and heterogeneous signals
in terms of exact model recovery, we study the variable selection property of
marginal screening (see Section 3). We show that under the presence of strong
heterogeneity in the signal, marginal screening fails to recover the exact
model with probability converging to 1, whereas under homogeneous signal it
can recover the exact model asymptotically (Genovese et al.,, 2012). It turns
out that due to heterogeneity, the spurious correlations become large and
create impediment to selecting the exact model. In correlated design, a
different problem known as unfaithfulness (Wasserman and Roeder,, 2009; Robins
et al.,, 2003) prevents marginal screening from achieving model consistency.
Specifically, due to “correlation cancellation”, marginal correlation between
$Y$ and $X_{j}$ becomes negligible even when $\beta_{j}$ is large and this
ultimately leads to false negatives. In this paper, we study independent
random design model in which correlation cancellation does not occur. Instead,
we identify a different source of problem under the presence of signal
heterogeneity that affects exact variable selection performance of marginal
screening.

On the computational side, modern methods like LASSO, SCAD, MC+ were initially
motivated as alternatives to Best Subset Selection (BSS). BSS is in general an
NP-hard optimization problem and was believed to be practically intractable
even for $p$ as small as $30$. Thanks to recent advancements in algorithms and
hardware, the optimal solution to the BSS problem can now be computed,
sometimes with approximations, for some practical settings. Jain et al.,
(2014) showed that a wide family of iterative hard thresholding (IHT)
algorithms can approximately solve the BSS problem, in the sense that they can
achieve similar goodness of fit with the best subset with slight violation of
the sparsity constraint. Liu and Foygel Barber, (2020) studied the optimal
thresholding operator for such iterative thresholding algorithms, which
manages to exploit fewer variables than IHT to achieve the same goodness fit
as BSS. Bertsimas et al., (2016) viewed the BSS problem through the lens of
mixed integer optimization (MIO) and showed that for $n$ in 100s and $p$ in
1000s, the MIO algorithm can obtain a near optimal solution reasonably fast.
More recently, Bertsimas and Van Parys, (2020) developed a new cutting plane
method that solves to provable optimality the Tikhonov-regularized (Tikhonov,,
1943) BSS problem with $n$ and $p$ in the 100s and 1000s respectively. A
recent work (Zhu et al.,, 2020) proposed an iterative splicing method called
Adaptive Best Subset Selection (ABESS) to solve the BSS problem. They also
showed that ABESS enjoys both statistical accuracy and polynomial
computational complexity when the design matrix satisfies sparse Reisz
condition and minimum signal strength is of order $\Omega\\{(s\log p\log\log
n/n)^{1/2}\\}$.

Given these recent advances in solving BSS, there has been growing
acknowledgement that BSS enjoys significant statistical superiority over the
aforementioned alternative methods. Bertsimas et al., (2016) and Bertsimas and
Van Parys, (2020) numerically demonstrated higher predictive power and lower
false discovery rate (FDR) respectively of the BSS solution compared to LASSO.
Fan et al., (2020) and Zhu and Wu, (2021) reported that the approximate BSS
solutions provided by IHT have much fewer false discoveries than LASSO, SCAD
and SIS, especially in the presence of highly correlated design. They also
theoretically showed that the model selection behavior of BSS does not
explicitly depend on the restricted eigenvalue condition for the design
(Bickel et al.,, 2009; Van De Geer and Bühlmann,, 2009), a condition which
appears unavoidable (assuming a standard computational complexity conjecture)
for any polynomial-time method (Zhang et al.,, 2014). This suggests that BSS
is robust against design collinearity in terms of model selection.

In this paper, we mainly focus on the _precise_ asymptotic bound, i.e., the
bound with the optimal constant, for the minimum signal strength that allows
BSS to achieve model consistency. Under a specific asymptotic setup, we show
that BSS achieves asymptotic exact recovery of the true model once the minimum
signal strength is above the information theoretic lower bound, meaning that
BSS is optimal in terms of the requirement on the signal strength. In
contrast, previous works such as Wainwright, (2009) and Fan et al., (2020)
analyze BSS from a non-asymptotic perspective: they show that BSS can achieve
model consistency under the optimal rate of the sample complexity. But their
analyses are tight only up to a multiplicative constants.

The rest of the paper is organized as follows. Section 2 introduces the
Asymptotically Ultra-Rare and Weak Minimum signal (AURWM) regime that
accommodates heterogeneous signal strengths. Section 3 shows that in the
presence of strong heterogeneity of the signal strength, marginal screening
procedures fail to achieve model consistency under the AURWM regime with
probability converging to 1. In Section 4, we derive the asymptotic minimax
0-1 loss under the AURWM regime and show that BSS is optimal in terms of the
requirement on the minimum signal strength. In Section 5, we propose a
computationally tractable two-stage algorithm that also enjoys model
consistency under essentially the same condition as BSS. Finally, in Section
6, we carry out simulation studies and numerically demonstrate the superiority
of our method over other competing methods.

Notation. Let ${\mathbb{R}}$ and ${\mathbb{R}}_{+}$ denote the set of real
numbers and the set of non-negative real numbers respectively. Denote by
${\mathbb{R}}^{p}$ the $p$-dimensional Euclidean space and by
${\mathbb{R}}^{p\times q}$ the space of real matrices of order $p\times q$.
For a positive integer $K$, denote by $[K]$ the set $\\{1,2,\ldots,K\\}$.

Regarding vectors and matrices, for a vector $v\in{\mathbb{R}}^{p}$, we denote
by $\left\|v\right\|_{2}$ the $\ell_{2}$-norm of $v$. We use ${\sf
I}_{p}\in{\mathbb{R}}^{p\times p}$ to denote the $p$-dimensional identity
matrix. For a matrix $A\in{\mathbb{R}}^{p\times p}$, we denote by $A_{j}$ and
$a_{j}$ the $j$th column and the transposed $j$th row of $A$ respectively. We
let ${\sf tr}(A)=\sum_{j=1}^{p}A_{jj}$ denote the trace of $A$.

Throughout the paper, let $O(\cdot)$ (respectively $\Omega(\cdot)$) denote the
standard big-O (respectively big-Omega) notation, i.e., we say
$a_{n}=O(b_{n})$ if there exists a universal constant $C>0$, such that
$a_{n}\leq Cb_{n}$ (respectively $a_{n}\geq Cb_{n}$) for all $n\in\mathbb{N}$.
Sometimes for notational convenience, we write $a_{n}\lesssim b_{n}$ in place
of $a_{n}=O(b_{n})$ and $a_{n}\gtrsim b_{n}$ in place of
$a_{n}=\Omega(b_{n})$. We write $a_{n}\sim b_{n}$ if $a_{n}=O(b_{n})$ and
$a_{n}=\Omega(b_{n})$. We denote by $\Omega_{\operatorname{\mathbb{P}}}$ the
big-Omega in probability: for a set of random variables $Z_{n}$ and a set of
constants $a_{n}$, $X_{n}=\Omega_{\operatorname{\mathbb{P}}}(a_{n})$ means
that for any $\varepsilon_{0}>0$, there exist $C_{\varepsilon_{0}}>0$ and
$n_{\varepsilon_{0}}\in\mathbb{N}$, both of which depend on $\varepsilon_{0}$,
such that

$\operatorname{\mathbb{P}}(\left|X_{n}/a_{n}\right|<C_{\varepsilon_{0}})\leq\varepsilon_{0},\quad\forall
n\geq n_{\varepsilon_{0}}.$

We use $\overset{\rm p}{\to}$ and $\overset{\rm d}{\to}$ to denote convergence
in probability and distribution respectively. Also we say $X\overset{\rm
d}{=}Y$ for two random variables $X,Y$ if their distributions are equal. We
denote by $\operatorname{\mathbbm{1}}(\cdot)$ the indicator function.

Finally, regarding probabilistic distributions, we use ${\sf N}(0,1)$ to
denote the standard Gaussian distribution. We use ${\sf N}_{p}(0,\Sigma)$ to
denote the $p$-dimensional Gaussian distribution with mean zero and variance-
covariance matrix $\Sigma\in{\mathbb{R}}^{p\times p}$. We denote by
${\sf{Ber}}(\pi)$ the Bernoulli distribution with success probability
$\pi\in[0,1]$.

##  2 Ultra rare and weak minimum signal regime

In this section, we focus on a specific asymptotic setup that allows
_heterogeneity_ among the sparse signals in high dimension. Throughout our
paper, we consider the following signal class:

${\mathcal{M}}_{s}^{a}:=\\{\beta\in{\mathbb{R}}^{p}:\|\beta\|_{0}=s,\min_{j\in\operatorname{\mathcal{S}}(\beta)}|\beta_{j}|\geq
a\\}.$

Here $a$ denotes the minimum signal strength of $\beta$. Note that the signal
class ${\mathcal{M}}_{s}^{a}$ only imposes a lower bound for the minimum
signal strength and thus allows arbitrarily large magnitudes across the true
signals. This implicitly accommodates heterogeneity in the signal, which is in
sharp contrast with the homogeneous signal setup considered by Genovese et
al., (2012).

Now we are in a position to introduce the Asymptotically Ultra Rare and Weak
Minimum signal regime (AURWM), in which we mainly consider the signal class
above with

$a=\left(\frac{2r\log p}{n}\right)^{1/2}\quad\text{and}\quad s=O(\log p).$ (2)

Besides, we set the sample size $n$ as

$n=\left\lfloor p^{k}\right\rfloor,\quad 0<k<1.$

The assumption that $s\lesssim\log p$ characterizes the ultra-rarity of the
signals, which is common in genetic studies such as GWAS (Yang et al.,, 2020).
Unless stated otherwise, from now on our statistical analysis follows the
scalings of $n,p,s,a$ in this AURWM regime. We say a support estimator
$\widehat{\operatorname{\mathcal{S}}}$ achieves asymptotic consistent recovery
in the AURWM regime if

$\lim_{p\to\infty}\sup_{\beta\in{\mathcal{M}}_{s}^{a}}\operatorname{\mathbb{P}}_{\beta}(\widehat{\operatorname{\mathcal{S}}}\neq\operatorname{\mathcal{S}}(\beta))=0.$
(3)

This paper mainly focuses on the criterion (3) to measure the quality of exact
recovery performance for an estimator $\widehat{\operatorname{\mathcal{S}}}$.

It is also worth mentioning that a relevant but different asymptotic setup is
studied by Genovese et al., (2012) and Ji et al., (2012). There the authors
assumed a Bayesian model such that all the signals are independent and
identicially distrbuted and that the sparsity $s\sim p^{1-\vartheta}$ for some
$\vartheta\in(0,1)$. Under such a setup they obtained asymptotically tight
phase transition boundaries with respect to Bayesian Hamming risk, which
partitions the $r\mathchar 45\relax\vartheta$ plane into three regions: (a)
Region of exact recovery, (b) Region of almost recovery, (c) Region of no
recovery. We skip the details of these results for brevity. The major
differences between their setup and ours are twofold: (1) They essentially
assume homogeneous signals; (2) They assume $s$ to grow in a polynomial
fashion with respect to $p$.

##  3 Marginal screening under heterogeneous signal

Marginal screening (MS) is one of the most widely used variable selection
methods in practice. It selects the variables with top absolute marginal
correlation with the response. Formally, for any $j\in[p]$, write
$\mu_{j}:=X_{j}^{\top}Y/n$. Given any possibly data-driven threshold
$\tau(X,Y)$, define the marginal screening estimator as follows:

$\widehat{\operatorname{\mathcal{S}}}_{\tau}:=\\{j\in[p]:\left|\mu_{j}\right|\geq\tau(X,Y)\\}.$
(4)

Note that $\mu_{j}$ is essentially equivalent to the marginal correlation
between $X_{j}$ and $Y$ because of isotropy of $X$. Marginal screening has
been applied in various fields for feature selection and dimension reduction,
including biomedicine (Huang et al.,, 2019; Lu,, 2005; Leisenring et al.,,
1997), survival data analysis (Hong et al.,, 2018; Li et al.,, 2016),
economics and econometrics (Wang et al.,, 2022; Huang et al.,, 2014).

Besides the broad applications, marginal screening has been shown to enjoy
some desirable statistical properties. Fan and Lv, (2008) established the sure
screening property of marginal screening under an ultra-high dimensional
setup, which serves as theoretical justification for MS to be used for
dimension reduction in many applications. Later, Genovese et al., (2012)
showed that MS enjoys the minimax optimal rate under Hamming loss with
homogeneous signals. Nevertheless, as mentioned in Section 1, precise
asymptotic characterization of the 0-1 loss of MS remains fairly underexplored
under high dimension, especially in the presence of heterogeneity in signal
strength.

###  3.1 Failure of MS in the AURWM regime

In this section, we study the 0-1 risk of the MS estimator. Define
$\operatorname{\mathcal{T}}:=\\{\widehat{\operatorname{\mathcal{S}}}_{\tau}\mid\tau:{\mathbb{R}}^{n\times
p}\times{\mathbb{R}}^{n}\to{\mathbb{R}}_{+}\\}$, which is the class of all
possible marginal screening estimators. Perhaps surprisingly, under the AURWM
regime, we show that MS fails to achieve exact model recovery in the minimax
sense.

###### Theorem 1

Under the AURWM regime, none of the MS estimators of form (4) can achieve
asymptotic exact recovery, i.e.,

$\lim_{p\to\infty}\inf_{\widehat{\operatorname{\mathcal{S}}}_{\tau}\in\operatorname{\mathcal{T}}}\sup_{\beta\in{\mathcal{M}}_{s}^{a}}\operatorname{\mathbb{P}}_{\beta}(\widehat{\operatorname{\mathcal{S}}}_{\tau}\neq\operatorname{\mathcal{S}}(\beta))=1.$

To understand the main message of this theorem, it is instructive to compare
it with the parallel result in Genovese et al., (2012) with homogeneous
signal. Specifically, Genovese et al., (2012) consider a Bayesian setup where
all the signal coefficients are independent and identically distributed
Bernoulli random variables (up to a universal constant). Under the AURWM
regime, $s=O(\log p)$, which implies that $\vartheta=1$ in Theorem 10 of
Genovese et al., (2012). Then Theorem 10 in Genovese et al., (2012) says that
when $r>1$, MS enjoys consistency in terms of Hamming risk and thus 0-1 risk
too. In contrast, when we broaden the signal class to ${\mathcal{M}}_{s}^{a}$
that embraces possibly heterogeneous signals, the same model consistency fails
to hold any more for MS as shown in Theorem 1. This comparison clearly reveals
the curse of signal heterogeneity on MS.

To see how signal heterogeneity hurts MS, for any $j\in[p]$, write $\mu_{j}$
as

$\mu_{j}=\frac{\beta_{j}\left\|X_{j}\right\|_{2}^{2}}{n}+\frac{X_{j}^{\top}(\sum_{\ell\neq
j}X_{\ell}\beta_{\ell}+E)}{n}=:\mu_{j}^{(1)}+\mu_{j}^{(2)}.$ (5)

Here $\mu^{(1)}_{j}=n^{-1}\beta_{j}\left\|X_{j}\right\|_{2}^{2}$ represents
the marginal contribution from $\beta_{j}$ to $\mu_{j}$, and $\mu_{j}^{(2)}$
represents the random error of $\mu_{j}$ due to the cross covariance between
$X_{j}$ and the other signals and noise. Suppose there are spiky signals among
$\\{\beta_{\ell}\\}_{\ell\neq j}$. Though
$\operatorname{\mathbb{E}}(\mu^{(2)}_{j})=0$ regardless of the magnitude of
$\beta$, the spiky signals may incur large variance of $\mu_{j}^{(2)}$ and
overwhelm the magnitude of $\mu^{(1)}_{j}$, which is the essential indicator
of the significance of $\beta_{j}$. Consequently, one cannot tell if
$\beta_{j}$ is a true variable based on only $\mu_{j}$ in the presence of
spiky signals.

##  4 Best subset selection

Now we shift our focus to BSS, one of the most classical variable selection
approaches. With the oracle knowledge of true sparsity $s$, BSS solves for

$\hat{\beta}_{\rm
best}\in{\arg\min}_{\beta\in{\mathbb{R}}^{p},\left\|\beta\right\|_{0}=s}\left\|Y-X\beta\right\|_{2}^{2}.$

Define
$P_{\operatorname{\mathcal{D}}}:=X_{\operatorname{\mathcal{D}}}(X_{\operatorname{\mathcal{D}}}^{\top}X_{\operatorname{\mathcal{D}}})^{-1}X_{\operatorname{\mathcal{D}}}^{\top}$,
which is the orthogonal projection operator onto the column space of
$X_{\operatorname{\mathcal{D}}}$. The BSS above can be alternatively viewed as
solving for

$\widehat{\operatorname{\mathcal{S}}}_{\rm
best}:={\arg\min}_{\operatorname{\mathcal{D}}\subseteq[p]:\left|\operatorname{\mathcal{D}}\right|=s}\frac{Y^{\top}({\sf
I}_{n}-P_{\operatorname{\mathcal{D}}})Y}{n}={\arg\max}_{\operatorname{\mathcal{D}}\subseteq[p]:\left|\operatorname{\mathcal{D}}\right|=s}\frac{Y^{\top}P_{\operatorname{\mathcal{D}}}Y}{n}.$
(6)

Using a union bound as in Wainwright, (2009) or Fan et al., (2020), one can
show that there exists a universal positive constant $\varphi$ (approximately
equal to 0.618) such that whenever $r>4/(1-\varphi)$, BSS achieves model
consistency, i.e.,

$\lim_{p\to\infty}\sup_{\beta\in{\mathcal{M}}_{s}^{a}}\operatorname{\mathbb{P}}_{\beta}(\widehat{\operatorname{\mathcal{S}}}_{\rm
best}\neq\operatorname{\mathcal{S}}(\beta))=0.$

We emphasize that the requirement on $r$ here is more stringent than needed:
we will show that BSS achieves model consistency whenever $r>1$, which turns
out to be necessary for any approach to obtain exact support recovery.

###  4.1 Exact support recovery of BSS

In the following theorem, we show that $r>1$ is sufficient for BSS to achieve
asymptotic exact recovery. Recall that $n=\left\lfloor p^{k}\right\rfloor$
with $0<k<1$.

###### Theorem 2

Let $r>1$ and write $\delta=r-1$. Then there exists a universal positive
constant $C_{0}$ such that whenever

$s<C_{0}\min\bigg{\\{}2k,\frac{4\delta}{1+\delta},\frac{\delta^{2}}{\\{(1+0.75\delta)^{1/2}+(1+0.5\delta)^{1/2}\\}^{2}}\bigg{\\}}\log
p,$

we have

$\lim_{p\to\infty}\sup_{\beta\in{\mathcal{M}}_{s}^{a}}\operatorname{\mathbb{P}}_{\beta}(\widehat{\operatorname{\mathcal{S}}}_{\rm
best}\neq\operatorname{\mathcal{S}}(\beta))=0.$

In order for BSS to achieve model consistency, we need to ensure that the
maximum spurious correlation, i.e., correlation between the spurious variables
and the response, is well controlled so that the best subset does not involve
any false discovery. One important ingredient of our analysis is the
asymptotic distribution of the maximum spurious correlation due to Fan et al.,
(2018), based on which we can derive the sharp constant in the minimum signal
strength for BSS to be model-consistent. It is worth emphasis that pursuing
the exact asymptotic distributions is crucial to obtain constant-sharp
results; typically, standard non-asymptotic analysis can only yield optimal
rates rather than optimal constants.

Note also that Theorem 2 requires $s$ to grow slowly. Given that we have at
least $\binom{p-s}{s}$ spurious models and that this number increases with
respect to $s$ when $s$ is small, a larger $s$ implies higher maximum spurious
correlation due to randomness and thus thinner chance for the best subset to
remain the true model.

In the next section, we show that $r>1$ is the weakest possible requirement on
the minimum signal strength for any approach to achieve asymptotic model
consistency.

###  4.2 Necessary condition for exact recovery

In this section, we show that under the AURWM regime, $r=1$ is the information
theoretic boundary of the minimum signal strength for exact model recovery. In
other words, if $r<1$, then no method can achieve model consistency. Towards
this end, we study the minimax 0-1 loss

$\inf_{\widehat{\operatorname{\mathcal{S}}}}\sup_{\beta\in{\mathcal{M}}_{s}^{a}}\operatorname{\mathbb{P}}_{\beta}(\widehat{\operatorname{\mathcal{S}}}\neq\operatorname{\mathcal{S}}(\beta)),$

where the infimum is taken over the class of all possible measurable functions
$\widehat{\operatorname{\mathcal{S}}}:(X,Y)\to\\{\operatorname{\mathcal{D}}\subseteq[p]:\left|\operatorname{\mathcal{D}}\right|=s\\}$.
The next theorem establishes a lower bound of the above minimax 0-1 loss.

######  Theorem 3 (Necessary Condition)

Under the AURWM regime, when $r<1$, there exists a universal positive constant
$c$ such that

$\lim_{p\to\infty}\inf_{\widehat{\operatorname{\mathcal{S}}}}\sup_{\beta\in{\mathcal{M}}_{a}^{s}}\operatorname{\mathbb{P}}_{\beta}(\widehat{\operatorname{\mathcal{S}}}\neq\operatorname{\mathcal{S}}(\beta))\geq
c.$

Theorem 3 suggests that $r\geq 1$ is a necessary condition for consistent
support recovery. Combining this with Theorem 2, we can see that BSS is
optimal in terms of requirement on the minimum signal strength to achieve
exact model recovery. The proof of Theorem 3 leverages Theorem 1 in Wang et
al., (2010) and detailed proof can be found in Section 9.1 of the
supplementary material.

##  5 Achieving statistical optimality with computational efficiency

In spite of the optimality of BSS in terms of model selection, its NP-hardness
seriously restricts its practical applicability. To address the computational
issue, we propose a two-stage algorithm called ETS (Estimate then Screen) that
combines the well-known IHT algorithm (Blumensath and Davies,, 2009; Jain et
al.,, 2014) with a follow-up coordinate screening step. We show that in the
AURWM regime, ETS enjoys the same selection optimality as BSS in terms of the
requirement on the minimum signal strength.

###  5.1 The ETS algorithm

In this section we introduce our ETS algorithm (Algorithm 2) in detail. Given
a partition parameter $0<\gamma<1$, ETS first splits the full sample
$(x_{i},Y_{i})_{i\in[n]}$ into two subsamples
$\operatorname{\mathcal{D}}_{1},\operatorname{\mathcal{D}}_{2}$ of respective
sizes $n_{1}=\left\lfloor\gamma n\right\rfloor$ and $n_{2}=n-n_{1}$. Then ETS
performs two main steps on these two sub-samples respectively. It applies IHT
to $\operatorname{\mathcal{D}}_{1}$ to obtain an estimator $\hat{\beta}^{\rm
iht}$ of the true signal vector $\beta$. Then ETS performs a coordinatewise
screening based on $\operatorname{\mathcal{D}}_{2}$ and $\hat{\beta}^{\rm
iht}$ to select the true variables. To be self-contained, in Algorithm 1 we
give the pseudo-code of the IHT algorithm. Let
$f:{\mathbb{R}}^{p}\to{\mathbb{R}},\theta\mapsto f(\theta)$ be a
differentiable objective function and $P_{\hat{s}}^{0}(\cdot)$ be the
projection operator that projects a vector onto the space of $\hat{s}$-sparse
vectors by selecting the $\hat{s}$ largest elements in magnitude.

Input: Objective function $f$, sparsity level $\hat{s}$, step size $h$ ;

$\beta^{(0)}=0$;

$t=0$ ;

while  _not converged_ do

$\beta^{(t+1)}=P_{\hat{s}}^{0}(\beta^{(t)}-h\nabla_{\theta}f(\beta^{t}))$;

$t\leftarrow t+1$

end while

Output: $\hat{\beta}^{\rm iht}=\beta^{(t)}$.

Algorithm 1 IHT

For $\ell\in\\{1,2\\}$, let $X^{(\ell)}\in{\mathbb{R}}^{n_{\ell}\times p}$ and
$Y^{(\ell)}\in{\mathbb{R}}^{n_{\ell}}$ denote the design matrix and the
response vector of the $\ell$th sub-sample respectively. ETS computes
$\hat{\beta}^{\rm iht}$ based on the first sub-sample $(X^{(1)},Y^{(1)})$ with
the square loss as the objective function:

$f_{n_{1}}(\theta;X^{(1)},Y^{(1)}):=\frac{1}{n_{1}}\|Y^{(1)}-X^{(1)}\theta\|_{2}^{2}.$
(7)

Next comes the screening step of ETS. For each $i\in[p]$, define

$\Delta_{i}:=\frac{X_{i}^{(2)\top}\left(Y^{(2)}-\sum_{j\neq
i}X_{j}^{(2)}\hat{\beta}^{\rm iht}_{j}\right)}{\|X_{i}^{(2)}\|_{2}}$ (8)

and

$\kappa_{\varsigma}(u):=\frac{a\left\|u\right\|_{2}}{2}+\frac{\varsigma^{2}\log
p}{a\left\|u\right\|_{2}},\quad\forall u\in{\mathbb{R}}^{n_{2}},$ (9)

where $\varsigma>0$ is specified later. ETS selects the $i$th variable if and
only if $\left|\Delta_{i}\right|>\kappa_{\varsigma}(X_{i}^{(2)})$. To see why
we can screen variables based on $\\{\Delta_{i}\\}_{i\in[p]}$, note that

$\Delta_{i}=\beta_{i}\|X_{i}^{(2)}\|_{2}+\frac{X_{i}^{(2)\top}\bigl{(}\sum_{j\neq
i}X_{j}^{(2)}(\beta_{j}-\hat{\beta}_{j}^{\rm
iht})+E\bigr{)}}{\|X_{i}^{(2)}\|_{2}}.$ (10)

Some probability argument shows that conditioned on
$\operatorname{\mathcal{D}}_{1}$ and $X_{i}^{(2)}$, $\Delta_{i}$ is
distributed as:

$\Delta_{i}\,\bigr{|}\,\big{(}\operatorname{\mathcal{D}}_{1},X_{i}^{(2)}\big{)}\overset{\rm
d}{=}\beta_{i}\big{\|}X_{i}^{(2)}\big{\|}_{2}+\bigg{\\{}1+\sum_{j\neq
i}(\beta_{j}-\hat{\beta}_{j}^{\rm iht})^{2}\bigg{\\}}^{1/2}g_{i},$

where $g_{i}\sim{\sf N}(0,1)$ and is independent of $X_{i}^{(2)}$. If IHT
performs well in the sense that $\|\hat{\beta}^{\rm iht}-\beta\|_{2}$ is
small, then for all $i\in\operatorname{\mathcal{S}}(\beta)$,
$\beta_{i}\|X_{i}^{(2)}\|_{2}$ becomes the dominant term in $\Delta_{i}$. In
contrast, for all $i\notin\operatorname{\mathcal{S}}(\beta)$,
$\beta_{i}\|X_{i}^{(2)}\|_{2}=0$ and we thus expect $\Delta_{i}$ to be small.
This suggests the existence of a threshold $t(\cdot)$ on
$(\Delta_{i})_{i\in[p]}$ that distinguishes the true support
$\operatorname{\mathcal{S}}(\beta)$ from the irrelevant variables. We follow
Ndaoud and Tsybakov, (2020) to choose the threshold function in (9), which is
shown to be a reasonable choice to identify the true variables.

For each $i\in[p]$, define
$\hat{\eta}_{i}(X,Y):=\operatorname{\mathbbm{1}}\\{\left|\Delta_{i}\right|>\kappa_{\varsigma}(X_{i}^{(2)})\\}$
and write

$\hat{\eta}(X,Y):=(\hat{\eta}_{1}(X,Y),\ldots,\hat{\eta}_{p}(X,Y))^{\top}.$

The selector $\hat{\eta}(X,Y)$ is the final estimate of the support
$\mathcal{S}(\beta)$ produced by the ETS algorithm. Algorithm 2 shows the
detailed steps of the ETS algorithm.

Input: Data $\mathcal{D}=\\{(x_{i},Y_{i})\\}_{i=1}^{n}$, sparsity level
$\hat{s}$, step size $h$, partition parameter $\gamma$ , threshold parameter
$\varsigma$;

1\. Randomly partition the whole dataset $\mathcal{D}$ into two disjoint
subsets $\mathcal{D}_{1}=(X^{(1)},Y^{(1)})$ and
$\mathcal{D}_{2}=(X^{(2)},Y^{(2)})$ ;

2\. Run IHT (Algorithm 1) with the objective function
$f_{n_{1}}(\theta;X^{(1)},Y^{(1)})$ in (7) and compute $\hat{\beta}_{\rm
iht}$;

3\. Construct the statistics $\\{\Delta_{i}\\}_{i=1}^{p}$ and thresholds
$\\{\kappa_{\varsigma}(X_{i}^{(2)})\\}_{i=1}^{p}$ using (8)-(9);

4\. Finally compute the selector $\hat{\eta}(X,Y)$;

Output: The selector $\hat{\eta}(X.Y)$.

Algorithm 2 ETS

Next, we define the binary decoder of the true support
$\operatorname{\mathcal{S}}(\beta)$ as
$\eta_{\beta}:=(\operatorname{\mathbbm{1}}\\{\beta_{1}\neq
0\\},\ldots,\operatorname{\mathbbm{1}}\\{\beta_{p}\neq 0\\})^{\top}$. The next
theorem shows that ETS can achieve exact recovery under suitable choices of
the tuning parameters.

###### Theorem 4

Assume $r>1$ and write $\delta=r-1$. Then there exist universal positive
constants $A_{1},A_{2}$ and suitable choices of tuning parameters $\hat{s}$,
$h$ such that with IHT iteration count
$t\geq\log\\{A_{1}(1+\|\beta\|_{2}^{2})/\delta\\}$,
$\gamma\in(0,\delta/(8+8\delta))$ and $\varsigma=(1+A_{2}\delta)^{1/2}$, we
have that
$\lim_{p\to\infty}\sup_{\beta\in{\mathcal{M}}_{a}^{s}}\operatorname{\mathbb{P}}_{\beta}(\hat{\eta}\neq\eta_{\beta})=0.$

Note that as the signal strength parameter $r$ approaches the information-
theoretic boundary, i.e., as $\delta$ approaches 0, ETS requires more
iterations to achieve model consistency. This is not surprising: intuitively,
weaker signals are harder to identify than strong ones.

Besides, ETS does not require the knowledge of the true sparsity $s$, but
requires the knowledge of $a$ in the second stage for accurate screening. If
the true sparsity $s$ is known, then we can enforce ETS to select exactly $s$
features as follows. Let $\left|\Delta\right|_{(m)}$ denote the $m$th largest
value of $\\{\left|\Delta_{i}\right|\\}_{i\in[p]}$. For each $i\in[p]$, define

$\hat{\eta}_{i}(X,Y;s)=\operatorname{\mathbbm{1}}\\{\left|\Delta_{i}\right|\geq\left|\Delta\right|_{(s)}\\}.$
(11)

Hence $\hat{\eta}(s):=(\hat{\eta}_{1}(s),\ldots,\hat{\eta}_{p}(s))^{\top}$
selects exactly $s$ features. The following corollary shows that under the
same conditions of Theorem 4, $\hat{\eta}(s)$ achieves model consistency.

###### Corollary 1

Assume $r>1$ and write $\delta=r-1$. If
$t\geq\log\\{A_{1}(1+\|\beta\|_{2}^{2})/\delta\\}$ with the same $A_{1}$ in
Theorem 4 and $\gamma\in(0,\delta/(8+8\delta))$, we have that

$\lim_{p\to\infty}\sup_{\beta\in{\mathcal{M}}_{s}^{a}}\operatorname{\mathbb{P}}_{\beta}(\hat{\eta}(s)\neq\eta_{\beta})=0.$
(12)

Detailed proofs of Theorem 4 and Corollary 12 can be found in Section 9.2 and
Section 9.3 of supplementary materials respectively. Besides, the sparsity
level $\hat{s}$ and the step-size $h$ are chosen based on certain geometric
properties of $f_{n_{1}}(\theta)$ and details can be found in Section 9.2 of
the supplementary materials.

##  6 Numerical experiments

In this section, we first numerically investigate the probability for MS to
achieve exact recovery of the true model with growing ambient dimension $p$
under both homogeneous and heterogeneous signal setup. Our results show that
while MS exhibits model consistency under the homogeneous signal regime, it
completely fails to do so under the heterogeneous signal regime, which is
consistent with Theorem 1. We then conduct simulation experiments to
demonstrate superiority of ETS over competing methods including LASSO and MS
as signal strength grows or signal heterogeneity grows. To this end, we
mention that we do not numerically compare exact BSS in this section mainly
due to computational issues. In most of our simulation setups we consider $p$
in thousands and BSS generally suffers from high computational costs in such
regimes, even with the help of modern optimization tools such as Gurobi
(Hastie et al.,, 2020). Instead, we focus on LASSO and ETS, both of which are
two different computational surrogates of the BSS problem.

###  6.1 Exact recovery performance of MS

In Figure 1, we demonstrate the asymptotics of MS under both homogeneous and
heterogeneous signal patterns. We consider $p\in\\{1000,2000,\ldots,8000\\}$
and signal strength parameter $r\in\\{2,3,4,5,6\\}$. We set $s=\left\lfloor
2\log p\right\rfloor$ and $n=\left\lfloor p^{0.9}\right\rfloor$. We let
$\tau(X,Y)$ in (4) be equal to the $s$th largest value of
$\\{\left|\mu_{1}\right|,\ldots,\left|\mu_{p}\right|\\}$, so that MS always
chooses a model of size $s$. For the homogeneous signal setup, we consider
$\beta$ with $\left\|\beta\right\|_{0}=s$ and $\beta_{j}=a$ for all
$j\in\operatorname{\mathcal{S}}(\beta)$, where $a$ is defined in (2). This
implies that the SNR varies between 0.19 and 2.15 across different choices of
$(r,p)$. For the heterogeneous signal setup, we consider $\beta$ with $(s-1)$
active coordinates equal to $a$ and one “spiky” coordinate equal to
$\\{10-(s-1)a^{2}\\}^{1/2}$. This ensures that the SNR is fixed at 10 for all
choices of $r,p$.

Figure 1(a) shows that under homogeneous signal MS is able to recover the
exact model with probability converging to 1 as $p$ grows. In contrast, Figure
1(b) shows that under heterogeneous signal MS never achieves exact model
recovery: plots for all values of $r$ are at level 0. Such a contrast
corroborates Theorem 1: signal spikes can give rise to substantial spurious
correlation and jeopardize the accuracy of MS.

(a) Proportion of exact recovery under weak homogeneous signal.
($0.19\leq\text{SNR}\leq 2.15$)

(b) Proportion of exact recovery under heterogeneous signal. (SNR = 10)

Figure 1: Asymptotics of MS with growing dimension $p$.

###  6.2 Effect of growing signal strength

Here we numerically compare the probability of exact support recovery of ETS
with those of LASSO and MS as signal strength parameter $r$ grows. We
investigate both homogeneous and heterogeneous signal patterns. We set
$p=2000$, $s\in\\{13,52\\}$ and $n=\left\lfloor p^{0.9}\right\rfloor=935$. We
set signal strength parameter $r\in\\{1.5,2,2.5,\ldots,9\\}$ in (2). The
support $\operatorname{\mathcal{S}}$ is chosen uniformly over all the size-$s$
subsets of $[p]$, and each support coordinate of $\beta$ is chosen as follows:

$\beta_{j}=(1-b_{j})(1+Z_{j}^{2}/n)^{1/2}a+b_{j}r^{1/2}\quad\forall
j\in\operatorname{\mathcal{S}},$

where $(Z_{j})_{j\in\operatorname{\mathcal{S}}}\overset{i.i.d.}{\sim}{\sf
N}(0,1)$, and where
$(b_{j})_{j\in\operatorname{\mathcal{S}}}\overset{i.i.d.}{\sim}{\sf{Ber}}(\pi)$
with $\pi\in\\{0,0.2\\}$. $\pi=0$ corresponds to the homogeneous signal
pattern, and $\pi=0.2$ corresponds to the heterogeneous signal pattern, where
spiky signals are present with probability $0.2$. Each entry $x_{ij}$ of the
design matrix $X$ is generated independently from ${\sf N}(0,1)$.

In this experiment, we grant all the approaches with the knowledge of $s$, so
that the comparison is fair. Using this oracle knowledge, we only look at the
solutions of the aforementioned three methods with sparsity exactly equal to
$s$. Specifically, for LASSO, we look at the solution path and select the
model of size exactly equal to $s$. For MS, we just select the top $s$
variables corresponding to the largest absolute values of $\mu$’s. For ETS, we
do not split data for estimation and screening separately; instead we use the
full data in both steps. Specifically, we replace $X^{(1)}$ and $Y^{(1)}$ with
$X$ and $Y$ respectively in (7) and replace $X^{(2)}$ and $Y^{(2)}$ with $X$
and $Y$ respectively in (8). We set gradient step size $h=0.5$ in IHT. We
choose projection size $\hat{s}$ by cross validation in terms of mean squared
prediction error. Lastly, for selecting exactly $s$ features we use (11) in
the screening stage of ETS. It is worthwhile to mention that from an
application point of view, incorporating data splitting in ETS is not
necessary as we are only interested in identifying the active signals, which
is akin to point estimation. Also, given the fact that $n\ll p$ in high
dimensional regime, using full sample in both the estimation and screening
step delivers greater sample efficiency and provides better inference.

Next, for each choice of $r$, we run LASSO, MS and ETS over 200 independent
Monte Carlo experiments to compute the empirical probability of exact
recovery. Figure 2 presents the results. We make the following important
observations:

  1. 1.

All the three methods enjoy higher chance of exact support recovery as the
signal strength grows;

  2. 2.

MS completely fails to achieve exact support recovery when $s$ becomes large
(compare panels (a) and (c)) or the signal becomes heterogeneous (compare
panels (a) and (b)).

  3. 3.

LASSO and ETS are insensitive to heterogeneity of the signal. However, LASSO
suffers from larger sparsity, while ETS is much more robust against it.

  4. 4.

Overall, ETS is the best among all the three methods in terms of exact support
recovery.

(a) Homogeneous signal pattern $(\pi=0,s=13)$

(b) Heterogeneous signal pattern $(\pi=0.2,s=13)$

(c) Homogeneous signal pattern $(\pi=0,s=52)$

(d) Heterogeneous signal pattern $(\pi=0.2,s=52)$

Figure 2: Plot of proportion of exact recovery for varying $r$.

###  6.3 Effect of growing heterogeneity

In this numerical experiment we study the effect of growing heterogeneity on
ETS, LASSO and MS. We set $p=2000,s=13$, $n=\left\lfloor
p^{0.9}\right\rfloor=935$ and $r\in\\{2,6\\}$ in (2). Next, we introduce
${\mathsf{n}_{\rm spike}}$, the number of “spiky” signals in
$\operatorname{\mathcal{S}}$. We vary ${\mathsf{n}_{\rm spike}}$ in
$\\{0\\}\cup[6]$. The case ${\mathsf{n}_{\rm spike}}=0$ corresponds to the
homogeneous signal setup where the true signals are set as $a$ uniformly. For
${\mathsf{n}_{\rm spike}}>0$, we randomly set $(s-{\mathsf{n}_{\rm spike}})$
signals in $\operatorname{\mathcal{S}}$ to be equal to $a$ and the remaining
signals to be equal to ${a_{\rm spike}}$, which is defined as

${a_{\rm spike}}:=\bigg{\\{}\frac{(2-sa^{2})}{{\mathsf{n}_{\rm
spike}}}+a^{2}\bigg{\\}}^{1/2}.$

Such a choice of ${a_{\rm spike}}$ ensures that the SNR always equals 2
whenever ${\mathsf{n}_{\rm spike}}>0$. We again perform ETS, LASSO and MS over
200 Monte Carlo simulations for each choice of $r$ and ${\mathsf{n}_{\rm
spike}}$ to obtain the empirical probability of exact support recovery.
Similarly to the previous sections, we assume that the true sparsity $s$ is
known and we apply the three methods in the same fashion as before.

Figure 3 shows again the detrimental effect of heterogeneity on MS in terms of
exact recovery. In both panels we see a significant drop in the proportion of
exact recovery for MS when ${\mathsf{n}_{\rm spike}}$ changes from 0 to 1.
This is consistent with the theory in Section 3. However, in Figure 3(b) we
see that the proportion of exact recovery is slowly increasing as
${\mathsf{n}_{\rm spike}}$ grows from 1 to 6. This is because as
${\mathsf{n}_{\rm spike}}$ increases, ${a_{\rm spike}}$ monotonically
decreases, so that the signals become more homogeneous. MS is then able to
recover the exact model more frequently. We do not see a similar phenomenon in
Figure 3(a) because ${a_{\rm spike}}$ is too large. Another important
observation is that while ETS and LASSO are both performing nearly perfectly
when $r=6$, ETS significantly outperforms both LASSO and MS when $r=2$.
Therefore, ETS is again the overall winner.

(a) $r=2$

(b) $r=6$

Figure 3: Plot of proportion of exact recovery with varying ${\mathsf{n}_{\rm
spike}}$.

##  7 Conclusion

In this paper, we study exact support recovery in high-dimensional sparse
linear regression with independent Gaussian design. We focus on the AURWM
regime that not only accommodates rare and weak signals as the ARW regime
does, but also allows heterogeneity in the signal strength. Our first
theoretical result (Theorem 1) shows that marginal screening fails to achieve
exact support recovery under the AURWM regime. The main reason is that the
presence of “spiky” signals increases the maximum spurious marginal
correlation, thereby blinding the marginal screening procedure to weak
signals. Therefore, one needs to be cautious with usage of marginal screening
for variable selection in practice.

In contrast, we show that BSS is robust to signal spikes and is able to
achieve model consistency under the AURWM regime with the optimal requirement
on signal strength (Theorem 2). The primary reason behind this is that unlike
MS, BSS takes into account multiple features simultaneously and thus selects
variables based on their capability of fitting the residualized responses
given the other variables rather than the responses themselves. Therefore,
spiky signals do not affect BSS: They are very likely to be in plausible
candidate models in the first place and their effect on the response has been
removed in the residualization procedure. Given the recent computational
advancements in solving BSS, our positive result on BSS makes it more
appealing from an application point of view.

However, it is worth mentioning that even with modern advances in
optimization, BSS suffers from high computational costs when the ambient
dimension is high. To address this issue, we propose a computationally
tractable two-stage method ETS that delivers essentially the same optimal
exact recovery performance as BSS (Theorem 4). Similar to BSS, ETS seeks for
the features that exhibit high explanation power for the residuals from the
model that excludes these features themselves (see (8)). Therefore, ETS is
robust to spiky signals. This fact together with the slowly growing sparsity
condition in (2) yields the optimal exact recovery accuracy of ETS.

Our work naturally raises several important questions for future research. One
question is whether similar optimality results hold for BSS and ETS when the
sparsity $s$ grows faster than $\log p$. The same question can also be asked
for correlated random design. Another direction of our interest is studying
the problem of exact recovery in a distributed setting where data are stored
at different places and communication between them is restricted.

##  8 Proofs

In this section we collect the proofs of some important results of this paper.
Specifically, we present the proofs of Theorem 1 and Theorem 2. The first one
shows that marginal screening fails to achieve exact recovery in the AURWM
regime, whereas BSS achieves model consistency under optimal requiremnt on the
signal strength.

###  8.1 Proof of Theorem 1

Consider a MS procedure
$\widehat{\operatorname{\mathcal{S}}}_{\tau}\in\operatorname{\mathcal{T}}$. To
start with note that

$\operatorname{\mathbb{P}}_{\beta}(\widehat{\operatorname{\mathcal{S}}}_{\tau}=\operatorname{\mathcal{S}}(\beta))=\operatorname{\mathbb{P}}_{\beta}\left(\max_{j\in
S^{c}}\left|\mu_{j}\right|<\tau(X,Y)\leq\min_{j\in
S}\left|\mu_{j}\right|\right)\leq\operatorname{\mathbb{P}}\left(\max_{j\in
S^{c}}\left|\mu_{j}\right|<\min_{j\in S}\left|\mu_{j}\right|\right).$

Recall that for for $j\in[p]$ we have

$\mu_{j}=\begin{cases}\beta_{j}\left\|X_{j}\right\|^{2}_{2}/n+\omega_{j}\left\|X_{j}\right\|_{2}g_{j}/n&\text{if
$j\in\operatorname{\mathcal{S}}(\beta)$}\\\
\omega_{j}\left\|X_{j}\right\|_{2}g_{j}/n&\text{if
$j\notin\operatorname{\mathcal{S}}(\beta)$},\end{cases}$

where $\omega_{j}^{2}=1+\sum_{k\neq j}\beta_{k}^{2}$ and
$g_{j}=X_{j}^{\top}(\sum_{k\neq
j}X_{k}\beta_{k}+E)/(\omega_{j}\left\|X_{j}\right\|_{2})\sim{\sf N}(0,1)$.
Thus we have

$\operatorname{\mathbb{P}}_{\beta}(\widehat{\operatorname{\mathcal{S}}}_{\tau}=S)\leq\operatorname{\mathbb{P}}_{\beta}\left(\max_{j\in
S^{c}}{\omega_{j}\left\|X_{j}\right\|_{2}\left|g_{j}\right|}/(2n\log
p)^{1/2}<\min_{j\in
S}|\beta_{j}\left\|X_{j}\right\|_{2}^{2}+\omega_{j}\left\|X_{j}\right\|_{2}g_{j}|/(2n\log
p)^{1/2}\right).$ (13)

Right hand side of Equation (13) does not depend on
$\widehat{\operatorname{\mathcal{S}}}_{\tau}$ hence the above inequality is
valid uniformly over the class $\operatorname{\mathcal{T}}$. Now choose a
sequence $\\{c_{p}\\}_{p=1}^{\infty}$ such that
$\lim_{p\to\infty}c_{p}^{2}/r\geq 1$ . Next construct a sequence of
$\beta^{(p)}$ in the following manner,

  * •

Consider the set $\operatorname{\mathcal{S}}_{0}=\\{1,\ldots,s\\}\subseteq[p]$
with $s=O(\log p)$.

  * •

Set $\beta_{1}^{(p)}=c_{p}$. For all other
$i\in\operatorname{\mathcal{S}}_{0}\setminus\\{1\\}$ set
$\beta_{i}^{(p)}=a=(2r(\log p)/n)^{1/2}$.

  * •

Set $\beta_{i}^{(p)}=0$ if $i\notin\operatorname{\mathcal{S}}_{0}$.

In this setup we have $\omega_{j}\sim(1+c_{p}^{2})^{1/2}$ for all $j\neq 1$.
Now fix $k_{0}\in\operatorname{\mathcal{S}}_{0}\setminus\\{1\\}$ (say
$k_{0}=2$). From Equation (13) it can be concluded that

$\displaystyle\sup_{\widehat{\operatorname{\mathcal{S}}}_{\tau}\in\operatorname{\mathcal{T}}}\operatorname{\mathbb{P}}_{\beta^{(p)}}(\widehat{\operatorname{\mathcal{S}}}_{\tau}=\operatorname{\mathcal{S}}_{0})$
$\displaystyle\leq\operatorname{\mathbb{P}}_{\beta^{(p)}}\left(\max_{j\in
S^{c}}{\omega_{j}\left\|X_{j}\right\|_{2}\left|g_{j}\right|}/(2n\log
p)^{1/2}<|\beta_{k_{0}}\left\|X_{k_{0}}\right\|_{2}^{2}+\omega_{k_{0}}\left\|X_{k_{0}}\right\|_{2}g_{k_{0}}|/(2n\log
p)^{1/2}\right).$

An elementary calculation shows that
$\\{g_{j}\\}_{j\in\operatorname{\mathcal{S}}_{0}^{c}}$ are i.i.d. standard
Gaussian. Also note that $\omega_{j}>(1+c_{p}^{2})^{1/2}$ for all
$j\in\operatorname{\mathcal{S}}_{0}^{c}$. Using these facts and lemma 3 from
Fletcher et al., (2009) we get

$\frac{1}{(1+c_{p}^{2})^{1/2}}\max_{j\in\operatorname{\mathcal{S}}_{0}^{c}}\omega_{j}\frac{\left\|X_{j}\right\|_{2}\left|g_{j}\right|}{(2n\log
p)^{1/2}}\geq\min_{j\in\operatorname{\mathcal{S}}_{0}^{c}}\frac{\left\|X_{j}\right\|_{2}}{n^{1/2}}\max_{j\in\operatorname{\mathcal{S}}_{0}^{c}}\frac{\left|g_{j}\right|}{(2\log
p)^{1/2}}\overset{\rm p}{\to}1.$

But recall that $\beta_{k_{0}}=\\{2r(\log p)/n\\}^{1/2}$ and thus we have the
following:

$\operatorname{\mathbb{P}}(|\beta_{k_{0}}\left\|X_{k_{0}}\right\|_{2}^{2}+\omega_{k_{0}}\left\|X_{k_{0}}\right\|_{2}g_{k_{0}}|/\\{(1+c_{p}^{2})(2n\log
p)\\}^{1/2}<1)\to 1.$

This tells that,

$\lim_{p\to\infty}\sup_{\widehat{\operatorname{\mathcal{S}}}_{\tau}\in\operatorname{\mathcal{T}}}\inf_{\beta\in{\mathcal{M}}_{s}^{a}}\operatorname{\mathbb{P}}_{\beta}(\widehat{\operatorname{\mathcal{S}}}_{\tau}=\operatorname{\mathcal{S}}(\beta))\leq\lim_{p\to\infty}\sup_{\widehat{\operatorname{\mathcal{S}}}_{\tau}\in\operatorname{\mathcal{T}}}\operatorname{\mathbb{P}}_{\beta^{(p)}}(\widehat{\operatorname{\mathcal{S}}}_{\tau}=\operatorname{\mathcal{S}}_{0})=0.$

This finishes the proof.

###  8.2 Proof of Theorem 2

In this proof we repararmetrize $\delta$ by $8\delta_{0}$ for algebraic
convenience. The main result can be salvaged by back substituting
$8\delta_{0}$ by $\delta$ in all the main equations in this section. Also, for
brevity of notation, in this proof we use
$\widehat{\operatorname{\mathcal{S}}}$ and $\operatorname{\mathcal{S}}$ to
denote that oracle BSS estimator and $\operatorname{\mathcal{S}}(\beta)$
respectively. Recall that BSS is defined as

$\widehat{\operatorname{\mathcal{S}}}={\arg\max}_{\operatorname{\mathcal{D}},\left|\operatorname{\mathcal{D}}\right|=s}\left\|P_{\operatorname{\mathcal{D}}}Y\right\|_{2}^{2}={\arg\min}_{\operatorname{\mathcal{D}}:\left|\operatorname{\mathcal{D}}\right|=s}\left\|({\sf
I}_{n}-P_{\operatorname{\mathcal{D}}})Y\right\|_{2}^{2}.$

Thus from the above definition we have the following equality:

$\displaystyle\operatorname{\mathbb{P}}(\widehat{\operatorname{\mathcal{S}}}\neq\operatorname{\mathcal{S}})$
$\displaystyle=\operatorname{\mathbb{P}}\left(\left\|P_{\operatorname{\mathcal{S}}}Y\right\|_{2}^{2}<\max_{\operatorname{\mathcal{D}}\neq\operatorname{\mathcal{S}}}\left\|P_{\operatorname{\mathcal{D}}}Y\right\|_{2}^{2}\right).$

Now we will try to understand how the quantity
$\left\|P_{\operatorname{\mathcal{S}}}Y\right\|_{2}^{2}$ behaves
asymptotically. First it is easy to see that
$P_{\operatorname{\mathcal{S}}}Y=\sum_{j\in\operatorname{\mathcal{S}}}X_{j}\beta_{j}+P_{\operatorname{\mathcal{S}}}E$.
Note that
$\sum_{j\in\operatorname{\mathcal{S}}}X_{j}\beta_{j}=\left\|\beta\right\|_{2}\tilde{\varepsilon}$
where $\tilde{\varepsilon}\sim{\sf N}_{n}(0,{\sf I}_{n})$ and independent of
the noise $z$. Hence we up with the following:

$\displaystyle\left\|P_{\operatorname{\mathcal{S}}}Y\right\|_{2}^{2}$
$\displaystyle=\Big{\|}\sum_{j\in
S}X_{j}\beta_{j}\Big{\|}_{2}^{2}+2E^{\top}P_{\operatorname{\mathcal{S}}}\left(\sum_{j\in
S}X_{j}\beta_{j}\right)+E^{\top}P_{\operatorname{\mathcal{S}}}E$
$\displaystyle=\Big{\|}\sum_{j\in
S}X_{j}\beta_{j}\Big{\|}_{2}^{2}+2E^{\top}\left(\sum_{j\in
S}X_{j}\beta_{j}\right)+E^{\top}P_{\operatorname{\mathcal{S}}}E$
$\displaystyle=\left\|\beta\right\|_{2}^{2}\left\|\tilde{\varepsilon}\right\|_{2}^{2}+2\left\|\beta\right\|_{2}E^{\top}\tilde{\varepsilon}+E^{\top}P_{\operatorname{\mathcal{S}}}E.$

Recall that $\left|\beta_{j}\right|\geq\\{2r(\log p)/n\\}^{1/2}$ for all
$j\in\operatorname{\mathcal{S}}$. This is presumably the hardest setup as
increasing signal strength can only decrease the error probability. Then
$\left\|\beta\right\|_{2}^{2}\geq(2rs\log p)/n$. also note that
$E^{\top}P_{\operatorname{\mathcal{S}}}E\sim\chi^{2}_{s}$. Hence we have ,

$\frac{\left\|P_{\operatorname{\mathcal{S}}}Y\right\|_{2}^{2}}{s\log p}\geq
2r\frac{\left\|\tilde{\varepsilon}\right\|_{2}^{2}}{n}+2\left(\frac{2r}{s\log
p}\right)^{1/2}\frac{\tilde{\varepsilon}^{\top}E}{n^{1/2}}+\frac{E^{\top}P_{\operatorname{\mathcal{S}}}E}{s\log
p}\overset{\rm p}{\longrightarrow}2r.$

Thus
$\lim_{p\to\infty}\operatorname{\mathbb{P}}(\widehat{\operatorname{\mathcal{S}}}\neq\operatorname{\mathcal{S}})\leq\operatorname{\mathbb{P}}\left(2r\leq\limsup_{p\to\infty}\max_{\operatorname{\mathcal{D}}\neq\operatorname{\mathcal{S}}}\left\|P_{\operatorname{\mathcal{D}}}Y\right\|_{2}^{2}/(s\log
p)\right)$. The limiting behaviour of the obtained maximal process turns out
to be very challenging to analyze and hence we do not directly study this
maximal process. Instead we focus on a related maximal process (will be
defined shortly) derived from the earlier one and we use the results from Fan
et al., (2018) to study its asymptotic behaviour. Now let us denote the set
$\operatorname{\mathcal{S}}\cap\operatorname{\mathcal{D}}$ by
$\operatorname{\mathcal{I}}_{0}$ and
$\operatorname{\mathcal{D}}\setminus\operatorname{\mathcal{S}}$ by
$\operatorname{\mathcal{I}}_{1}$, i.e.,
$\operatorname{\mathcal{D}}=\operatorname{\mathcal{I}}_{0}\cup\operatorname{\mathcal{I}}_{1}$.
Next define the class
$\operatorname{\mathcal{J}}_{\operatorname{\mathcal{I}}_{0}}=\\{\operatorname{\mathcal{I}}_{1}\subseteq[p]:\operatorname{\mathcal{I}}_{1}\cap
S=\emptyset,\left|\operatorname{\mathcal{I}}_{1}\cup\operatorname{\mathcal{I}}_{0}\right|=s\\}$
for each $\operatorname{\mathcal{I}}_{0}\subset S$. Note that
$0\leq\left|\operatorname{\mathcal{I}}_{0}\right|\leq s-1$ from the
construction (if $\left|\operatorname{\mathcal{I}}_{0}\right|=s$ then
$\operatorname{\mathcal{D}}=\operatorname{\mathcal{S}}$). The random variable
of interest can be rewritten as follows:

$\max_{\operatorname{\mathcal{D}}\neq\operatorname{\mathcal{S}}}\frac{\left\|P_{\operatorname{\mathcal{D}}}Y\right\|_{2}^{2}}{s\log
p}=\max_{\operatorname{\mathcal{I}}_{0}:\operatorname{\mathcal{I}}_{0}\subset\operatorname{\mathcal{S}}}\;\max_{\operatorname{\mathcal{I}}_{1}:\operatorname{\mathcal{I}}_{1}\in\operatorname{\mathcal{J}}_{\operatorname{\mathcal{I}}_{0}}}\frac{\left\|P_{\operatorname{\mathcal{I}}_{0}\cup\operatorname{\mathcal{I}}_{1}}Y\right\|_{2}^{2}}{s\log
p}.$

Using union bound we get,

$\operatorname{\mathbb{P}}(\widehat{\operatorname{\mathcal{S}}}\neq\operatorname{\mathcal{S}})\leq\sum_{\operatorname{\mathcal{I}}_{0}\subset\operatorname{\mathcal{S}}}\operatorname{\mathbb{P}}\left(\max_{\operatorname{\mathcal{I}}_{1}:\operatorname{\mathcal{I}}_{1}\in\operatorname{\mathcal{J}}_{\operatorname{\mathcal{I}}_{0}}}\frac{\left\|P_{\operatorname{\mathcal{I}}_{0}\cup\operatorname{\mathcal{I}}_{1}}Y\right\|_{2}^{2}}{s\log
p}>\frac{\left\|P_{\operatorname{\mathcal{S}}}Y\right\|_{2}^{2}}{s\log
p}\right).$ (14)

Now fix a subset $\operatorname{\mathcal{I}}_{0}$ of the true support
$\operatorname{\mathcal{S}}$. Similar to previous section define
$\tilde{Y}=Y-X_{\operatorname{\mathcal{I}}_{0}}\beta_{\operatorname{\mathcal{I}}_{0}}$
and this independent of the features in
$\operatorname{\mathcal{I}}_{0}\cup\operatorname{\mathcal{I}}_{1}$. Also we
have

$\displaystyle\left\|P_{\operatorname{\mathcal{I}}_{0}\cup\operatorname{\mathcal{I}}_{1}}Y\right\|_{2}^{2}=\left\|P_{\operatorname{\mathcal{I}}_{0}\cup\operatorname{\mathcal{I}}_{1}}\tilde{Y}\right\|_{2}^{2}+\left\|X_{\operatorname{\mathcal{I}}_{0}}\beta_{\operatorname{\mathcal{I}}_{0}}\right\|_{2}^{2}+2\beta_{\operatorname{\mathcal{I}}_{0}}^{\top}X_{\operatorname{\mathcal{I}}_{0}}^{\top}\tilde{Y},$
$\displaystyle\hskip
22.1931pt\left\|P_{\operatorname{\mathcal{S}}}Y\right\|_{2}^{2}=\left\|P_{\operatorname{\mathcal{S}}}\tilde{Y}\right\|_{2}^{2}+\left\|X_{\operatorname{\mathcal{I}}_{0}}\beta_{\operatorname{\mathcal{I}}_{0}}\right\|_{2}^{2}+2\beta_{\operatorname{\mathcal{I}}_{0}}^{\top}X_{\operatorname{\mathcal{I}}_{0}}^{\top}\tilde{Y}.$

Thus the summands in the right hand side of (14) can be written as the
probability of the event
$\\{\max_{\operatorname{\mathcal{I}}_{1}:\operatorname{\mathcal{I}}_{1}\in\operatorname{\mathcal{J}}_{\operatorname{\mathcal{I}}_{0}}}\left\|P_{\operatorname{\mathcal{I}}_{0}\cup\operatorname{\mathcal{I}}_{1}}g\right\|_{2}^{2}/(s\log
p)>\left\|P_{\operatorname{\mathcal{S}}}g\right\|_{2}^{2}/(s\log p)\\}$, where
$g:=(1+\left\|\beta_{\operatorname{\mathcal{S}}\setminus\operatorname{\mathcal{I}}_{0}}\right\|_{2}^{2})^{-1/2}\tilde{Y}$.
Note that $g\sim{\sf N}_{n}(0,{\sf I}_{n})$ and is independent of the features
in $\mathcal{D}$. Now fix a specific $\operatorname{\mathcal{I}}_{0}$. In the
analysis we encounter the maximal process

$\max_{\operatorname{\mathcal{I}}_{1}:\operatorname{\mathcal{I}}_{1}\in\operatorname{\mathcal{J}}_{\operatorname{\mathcal{I}}_{0}}}\frac{\left\|P_{\operatorname{\mathcal{I}}_{0}\cup\operatorname{\mathcal{I}}_{1}}g\right\|_{2}^{2}}{s\log
p},$

Now consider the set of indices
$F_{\operatorname{\mathcal{I}}_{0}}=(\\{1,\cdots,p\\}\setminus\operatorname{\mathcal{S}})\cup\operatorname{\mathcal{I}}_{0}$.
Hence it is easy to see that
$\tilde{p}:=\left|F_{\operatorname{\mathcal{I}}_{0}}\right|=p-s+\left|\operatorname{\mathcal{I}}_{0}\right|$.
Without loss of generality, let
$F_{\operatorname{\mathcal{I}}_{0}}=\\{1,\ldots,\tilde{p}\\}$. Also define
$\tilde{s}:=s-\left|\operatorname{\mathcal{I}}_{0}\right|$. Let the set
$\operatorname{\mathcal{V}}_{\operatorname{\mathcal{I}}_{0}}=\\{\alpha\in{\mathbb{R}}^{p}:\left\|\alpha\right\|_{0}=s,\left\|\alpha\right\|_{2}=1,\operatorname{\mathcal{I}}_{0}\subseteq\operatorname{\mathcal{S}}(\alpha),\alpha_{F_{\operatorname{\mathcal{I}}_{0}}^{c}}=0\\}$.
Here $\alpha_{J}$ denotes the sub-vector of $\alpha$ corresponding to the
indices in $J\subseteq[p]$. Next we will focus on the random variable,

$\widehat{L}_{n}:=\widehat{L}_{n}(\tilde{s},\tilde{p})=\sup_{\alpha\in\operatorname{\mathcal{V}}_{\operatorname{\mathcal{I}}_{0}}}\frac{1}{n^{1/2}}\sum_{i=1}^{n}\frac{\alpha^{\top}(g_{i}x_{i})}{(\alpha^{\top}\widehat{\boldsymbol{\Sigma}}_{n}\alpha)^{1/2}},$
(15)

here
$\widehat{\boldsymbol{\Sigma}}_{n}=\frac{1}{n}\sum_{i=1}^{n}x_{i}x_{i}^{\top}$.
Now recall that
$\operatorname{\mathcal{D}}=\operatorname{\mathcal{I}}_{0}\cup\operatorname{\mathcal{I}}_{1}$
for all $\operatorname{\mathcal{D}}\neq\operatorname{\mathcal{S}}$ with
$\left|\operatorname{\mathcal{D}}\right|=s$. To see the connection, first note
that the above optimization problem can be viewed as the following:

$\displaystyle\widehat{L}_{n}$
$\displaystyle=\max_{\operatorname{\mathcal{I}}_{1}\in\operatorname{\mathcal{J}}_{\operatorname{\mathcal{I}}_{0}}}\max_{\alpha\in\operatorname{\mathcal{V}}_{\operatorname{\mathcal{I}}_{0}\cup\operatorname{\mathcal{I}}_{1}}}\frac{\alpha_{\operatorname{\mathcal{D}}}^{\top}(\sum_{i=1}^{n}g_{i}x_{i,\operatorname{\mathcal{D}}}/n^{1/2})}{\\{\alpha_{\operatorname{\mathcal{D}}}^{\top}(\widehat{\boldsymbol{\Sigma}}_{n,\operatorname{\mathcal{D}}\operatorname{\mathcal{D}}})\alpha_{\operatorname{\mathcal{D}}}\\}^{1/2}}$
$\displaystyle=\max_{\operatorname{\mathcal{I}}_{1}\in\operatorname{\mathcal{J}}_{\operatorname{\mathcal{I}}_{0}}}\left\\{(\sum_{i=1}^{n}g_{i}x_{i,\operatorname{\mathcal{D}}}/n^{1/2})^{\top}\widehat{\boldsymbol{\Sigma}}^{-1}_{n,\operatorname{\mathcal{D}}\operatorname{\mathcal{D}}}(\sum_{i=1}^{n}g_{i}x_{i,\operatorname{\mathcal{D}}}/n^{1/2})\right\\}^{1/2}$
$\displaystyle=\max_{\operatorname{\mathcal{I}}_{1}\in\operatorname{\mathcal{J}}_{\operatorname{\mathcal{I}}_{0}}}\\{g^{\top}X_{\operatorname{\mathcal{D}}}(X_{\operatorname{\mathcal{D}}}^{\top}X_{\operatorname{\mathcal{D}}})^{-1}X_{\operatorname{\mathcal{D}}}^{\top}g\\}^{1/2}$
$\displaystyle=\max_{\operatorname{\mathcal{I}}_{1}\in\operatorname{\mathcal{J}}_{\operatorname{\mathcal{I}}_{0}}}\left\|P_{\operatorname{\mathcal{I}}_{0}\cup\operatorname{\mathcal{I}}_{1}}g\right\|_{2}.$

Thus it is essential to study the asymptotic property of $\widehat{L}_{n}$.
Now we define the standardized version of $\widehat{L}_{n}$ as follows

$L_{n}:=L_{n}(\tilde{s},\tilde{p})=\sup_{\alpha\in\operatorname{\mathcal{V}}_{\operatorname{\mathcal{I}}_{0}}}\frac{1}{n^{1/2}}\sum_{i=1}^{n}\alpha^{\top}(g_{i}x_{i}).$

Let $\boldsymbol{Z}=(Z_{1},\cdots,Z_{\tilde{p}})$ be $\tilde{p}-$variate
Gaussian random variable with covariance matrix $I_{\tilde{p}\times\tilde{p}}$
and define the random variable
$T^{*}:=T^{*}(\tilde{s},\tilde{p})=\sup_{\alpha\in\operatorname{\mathcal{V}}_{\operatorname{\mathcal{I}}_{0}}}\alpha_{F_{\operatorname{\mathcal{I}}_{0}}}^{\top}\boldsymbol{Z}$.

###### Lemma 1

There exists universal constants $K_{0},K_{1}$ such that for any
$\delta_{1}\in(0,K_{0}K_{1}]$,

$\left|L_{n}-T^{*}\right|\lesssim
n^{-1}c_{n}^{1/2}(\tilde{s},\tilde{p})+K_{0}K_{1}n^{-3/2}c_{n}^{2}(\tilde{s},\tilde{p})+\delta_{1}$
(16)

holds with probability at least $1-C\Delta_{n}(s,\tilde{p};\delta_{1})$ where
$c_{n}(s,\tilde{p})=s\log(e\tilde{p}/s)\vee\log n$ and

$\Delta_{n}(\tilde{s},\tilde{p};\delta_{1})=(K_{0}K_{1})^{3}\frac{\\{\tilde{s}b_{n}(\tilde{s},\tilde{p})\\}^{2}}{\delta_{1}^{3}n^{1/2}}+(K_{0}K_{1})^{4}\frac{\\{\tilde{s}b_{n}(\tilde{s},\tilde{p})\\}^{5}}{\delta_{1}^{4}n}$

with $b_{n}(\tilde{s},\tilde{p})=\log(\tilde{p}/\tilde{s})\vee\log n$.

###### Lemma 2

Assume that the sample size satisfies $n\geq C_{1}(K_{0}\vee
K_{1})^{4}c_{n}(\tilde{s},\tilde{p})$. then with probability at least
$1-C_{2}n^{-1/2}c_{n}^{1/2}(\tilde{s},\tilde{p})$,

$\left|\widehat{L}_{n}-L_{n}\right|\lesssim(K_{0}\vee
K_{1})^{2}K_{0}K_{1}n^{-1/2}c_{n}(\tilde{s},\tilde{p}),$ (17)

where $c_{n}(\tilde{s},\tilde{p})=\tilde{s}\log(e\tilde{p}/\tilde{s})\vee\log
n$.

Proof of the above two lemmas are omitted as it is in the same line of the
proofs of Fan et al., (2018). Now applying Lemma 1 and 2 with

$\delta_{1}=\delta_{1}(s,\tilde{p})=(K_{0}K_{1})^{3/4}\min[1,n^{-1/8}\\{\tilde{s}b_{n}(\tilde{s},\tilde{p})^{3/8}\\}]$

yields that with probability at least
$1-C(K_{0}K_{1})^{3/4}n^{-1/8}\\{sb_{n}(s,\tilde{p})\\}^{7/8}$,

$\left|\widehat{L}_{n}-T^{*}\right|\lesssim(K_{0}K_{1})^{3/4}n^{-1/8}\\{\tilde{s}b_{n}(\tilde{s},\tilde{p})\\}^{3/8}.$

Together with Lemma 2.3 from Chernozhukov et al., (2014) we can conclude that

$\sup_{t\in{\mathbb{R}}}\left|\operatorname{\mathbb{P}}(\widehat{L}_{n}\leq
t)-\operatorname{\mathbb{P}}(T^{*}\leq t)\right|\lesssim
C(K_{0}K_{1})^{3/4}n^{-1/8}\\{\tilde{s}b_{n}(\tilde{s},\tilde{p})\\}^{7/8}.$
(18)

Next by the definition of $T^{*}$ it follows that

$T^{*2}=\max_{\operatorname{\mathcal{I}}_{1}\in\operatorname{\mathcal{J}}_{\operatorname{\mathcal{I}}_{0}}}\left\|\boldsymbol{Z}_{\operatorname{\mathcal{I}}_{0}\cup\operatorname{\mathcal{I}}_{1}}\right\|_{2}=\sum_{j\in\operatorname{\mathcal{I}}_{0}}Z_{j}^{2}+\max_{\operatorname{\mathcal{I}}_{1}\in\operatorname{\mathcal{J}}_{\operatorname{\mathcal{I}}_{0}}}\sum_{k\in\operatorname{\mathcal{I}}_{1}}Z_{k}^{2}.$

Let $\boldsymbol{W}\sim{\sf N}_{(p-s)}(0,{\sf I}_{(p-s)})$ be a Gaussian
vector independent of $\boldsymbol{Z}$. Thus it follows that

$T^{*2}\overset{d}{=}\sum_{j\in\operatorname{\mathcal{I}}_{0}}Z_{j}^{2}+\sum_{k=p-2s+\left|\operatorname{\mathcal{I}}_{0}\right|+1}^{p-s}W^{2}_{(k:p-s)}\leq\sum_{j\in\operatorname{\mathcal{I}}_{0}}Z_{j}^{2}+(s-\left|\operatorname{\mathcal{I}}_{0}\right|)W^{2}_{(p-s:p-s)}.$

From Equation (18) it also follows that

$\sup_{t\geq 0}\left|\operatorname{\mathbb{P}}(\widehat{L}^{2}_{n}\leq
t)-\operatorname{\mathbb{P}}(T^{*2}\leq t)\right|\lesssim
C(K_{0}K_{1})^{3/4}n^{-1/8}\\{\tilde{s}b_{n}(\tilde{s},\tilde{p})\\}^{7/8}.$
(19)

Now from the assumption we have $r=1+8\delta_{0}$. Assume that

$s\leq
0.5\min\\{\delta_{0},\frac{2\delta_{0}^{2}}{\\{(1+6\delta_{0})^{1/2}+(1+4\delta_{0})^{1/2}\\}^{2}}\\}\log
p.$

Hence $\left|\operatorname{\mathcal{I}}_{0}\right|\leq s-1\leq\delta_{0}\log
p$. Thus we have,

$\displaystyle\quad\operatorname{\mathbb{P}}(T^{*2}>2(1+4\delta_{0})(s-\left|\operatorname{\mathcal{I}}_{0}\right|)\log
p)$
$\displaystyle=\operatorname{\mathbb{P}}(\sum_{j\in\operatorname{\mathcal{I}}_{0}}Z_{j}^{2}+(s-\left|\operatorname{\mathcal{I}}_{0}\right|)W^{2}_{(p-s:p-s)}>2(1+4\delta_{0})(s-\left|\operatorname{\mathcal{I}}_{0}\right|)\log
p)$
$\displaystyle\leq\operatorname{\mathbb{P}}\left(\sum_{j\in\operatorname{\mathcal{I}}_{0}}Z_{j}^{2}>4\delta_{0}(s-\left|\operatorname{\mathcal{I}}_{0}\right|)\log
p\right)+\operatorname{\mathbb{P}}\left((s-\left|\operatorname{\mathcal{I}}_{0}\right|)W^{2}_{(p-s:p-s)}>2(1+2\delta_{0})(s-\left|\operatorname{\mathcal{I}}_{0}\right|)\log
p\right)$
$\displaystyle\overset{(a)}{\leq}\operatorname{\mathbb{P}}\left(\frac{\sum_{j\in\operatorname{\mathcal{I}}_{0}}Z_{j}^{2}-\left|\operatorname{\mathcal{I}}_{0}\right|}{\left|\operatorname{\mathcal{I}}_{0}\right|}>(4\delta_{0}\log
p-\left|\operatorname{\mathcal{I}}_{0}\right|)/\left|\operatorname{\mathcal{I}}_{0}\right|\right)+\operatorname{\mathbb{P}}\left(W^{2}_{(p-s:p-s)}>2(1+2\delta_{0})\log
p\right)$
$\displaystyle\overset{(b)}{\leq}\operatorname{\mathbb{P}}\left(\frac{\sum_{j\in\operatorname{\mathcal{I}}_{0}}Z_{j}^{2}-\left|\operatorname{\mathcal{I}}_{0}\right|}{\left|\operatorname{\mathcal{I}}_{0}\right|}>(3\delta_{0}\log
p)/\left|\operatorname{\mathcal{I}}_{0}\right|\right)+\operatorname{\mathbb{P}}\left(W^{2}_{(p-s:p-s)}>2(1+2\delta_{0})\log
p\right)$ $\displaystyle\lesssim\exp(-0.75\delta_{0}\log
p)+(p-s)\operatorname{\mathbb{P}}\left(W^{2}_{1}>2(1+2\delta_{0})\log
p\right)$ $\displaystyle\lesssim
p^{-0.75\delta_{0}}+C\frac{p^{-2\delta_{0}}}{\sqrt{\log p}}.$

Inequality $(a)$ uses $s-\left|\operatorname{\mathcal{I}}_{0}\right|\geq 1$
and inequality $(b)$ uses
$\left|\operatorname{\mathcal{I}}_{0}\right|<s<\delta_{0}\log p$. Now define
the event
$\operatorname{\mathcal{E}}_{\operatorname{\mathcal{I}}_{0}}:=\\{\left\|P_{\operatorname{\mathcal{S}}}{g}\right\|_{2}^{2}/s\log
p>2(1+4\delta_{0})R_{\operatorname{\mathcal{I}}_{0}}\\}$ where
$R_{\operatorname{\mathcal{I}}_{0}}:=(s-\left|\operatorname{\mathcal{I}}_{0}\right|)/s$.
Recall that

$\begin{split}\frac{\left\|P_{\operatorname{\mathcal{S}}}{g}\right\|_{2}^{2}}{s\log
p}&\geq\frac{\left\|\frac{\sum_{j\in\operatorname{\mathcal{S}}\setminus\operatorname{\mathcal{I}}_{0}}X_{j}\beta_{j}+P_{\operatorname{\mathcal{S}}}E}{(1+\|\beta_{\operatorname{\mathcal{S}}\setminus\operatorname{\mathcal{I}}_{0}}\|_{2}^{2})^{1/2}}\right\|_{2}^{2}}{s\log
p}\\\
&\geq\frac{\left\\{\frac{\left\|\sum_{j\in\operatorname{\mathcal{S}}\setminus\operatorname{\mathcal{I}}_{0}}X_{j}\beta_{j}\right\|_{2}}{(1+\|\beta_{\operatorname{\mathcal{S}}\setminus\operatorname{\mathcal{I}}_{0}}\|_{2}^{2})^{1/2}}-\frac{\left\|P_{\operatorname{\mathcal{S}}}E\right\|_{2}}{(1+\|\beta_{\operatorname{\mathcal{S}}\setminus\operatorname{\mathcal{I}}_{0}}\|_{2}^{2})^{1/2}}\right\\}^{2}}{s\log
p}=(T_{1}^{1/2}-T_{2}^{1/2})^{2}.\end{split}$ (20)

where

$T_{1}:=\frac{\left\|\sum_{j\in\operatorname{\mathcal{S}}\setminus\operatorname{\mathcal{I}}_{0}}X_{j}\beta_{j}\right\|_{2}^{2}}{(1+\left\|\beta_{\operatorname{\mathcal{S}}\setminus\operatorname{\mathcal{I}}_{0}}\right\|_{2}^{2})s\log
p}\geq\frac{2rR_{\operatorname{\mathcal{I}}_{0}}}{1+2r(s-\left|\operatorname{\mathcal{I}}_{0}\right|)(\log
p)/n}\frac{V_{n}}{n},$

and

$V_{n}:=\frac{\|\sum_{j\in\operatorname{\mathcal{S}}\setminus\operatorname{\mathcal{I}}_{0}}X_{j}\beta_{j}\|_{2}^{2}}{\|\beta_{\operatorname{\mathcal{S}}\setminus\operatorname{\mathcal{I}}_{0}}\|_{2}^{2}}$

is an $\chi^{2}_{n}$ random variable. Also we have

$T_{2}:=\frac{\left\|P_{\operatorname{\mathcal{S}}}E\right\|^{2}_{2}}{(1+\left\|\beta_{\operatorname{\mathcal{S}}\setminus\operatorname{\mathcal{I}}_{0}}\right\|_{2}^{2})s\log
p}\leq V_{s}/s\log p$

where $V_{s}:=\left\|P_{\operatorname{\mathcal{S}}}E\right\|_{2}^{2}$ is an
$\chi^{2}_{s}$ random variable independent of $X_{S}$. Next we state the
following simple algebraic relationship:

$(1+6\delta_{0})^{1/2}-\frac{\delta_{0}}{(1+6\delta_{0})^{1/2}+(1+4\delta_{0})^{1/2}}\geq(1+4\delta_{0})^{1/2}.$

In light of Equation (20) and using the above algebraic inequality we have the
following:

$\operatorname{\mathcal{E}}_{\operatorname{\mathcal{I}}_{0}}^{c}\subseteq\\{T_{1}\leq
2(1+6\delta_{0})R_{\operatorname{\mathcal{I}}_{0}}\\}\bigcup\\{T_{2}\geq
2\delta_{0}^{2}\\{(1+6\delta_{0})^{1/2}+(1+4\delta_{0})^{1/2}\\}^{-2}R_{\operatorname{\mathcal{I}}_{0}}\\}$

Next we have

$\operatorname{\mathbb{P}}(T_{1}\leq
2(1+6\delta_{0})R_{\operatorname{\mathcal{I}}_{0}})\leq\operatorname{\mathbb{P}}\left(\frac{V_{n}}{n}\leq\frac{1+6\delta_{0}}{1+8\delta_{0}}(1+2rs\log
p/n)\right).$

Now choose large $n$ such that $(1+6\delta_{0})(1+2rs\log
p/n))<(1+7\delta_{0})$. Then for large $n$ we have,

$\operatorname{\mathbb{P}}(T_{1}\leq
2(1+6\delta_{0})R_{\operatorname{\mathcal{I}}_{0}})\leq\operatorname{\mathbb{P}}\left(\left|V_{n}/n-1\right|\geq\frac{\delta_{0}}{1+8\delta_{0}}\right)\lesssim\exp\left(-C^{*}\frac{\delta_{0}}{1+8\delta_{0}}n\right),$

where $C^{*}$ is an universal constant. Now we analyze the quantity $T_{2}$.
We have

$\displaystyle\operatorname{\mathbb{P}}(T_{2}\geq\frac{2\delta_{0}^{2}}{\\{(1+6\delta_{0})^{1/2}+(1+4\delta_{0})^{1/2}\\}^{2}}R_{\operatorname{\mathcal{I}}_{0}})$
$\displaystyle\leq\operatorname{\mathbb{P}}(V_{s}/s\geq\frac{2\delta_{0}^{2}}{\\{(1+6\delta_{0})^{1/2}+(1+4\delta_{0})^{1/2}\\}^{2}}R_{\operatorname{\mathcal{I}}_{0}}\log
p)$
$\displaystyle\leq\operatorname{\mathbb{P}}(V_{s}\geq\frac{2\delta_{0}^{2}}{\\{(1+6\delta_{0})^{1/2}+(1+4\delta_{0})^{1/2}\\}^{2}}\log
p)$ $\displaystyle\leq\operatorname{\mathbb{P}}(\left|V_{s}/s-1\right|\geq
0.5\frac{2\delta_{0}^{2}}{\\{(1+6\delta_{0})^{1/2}+(1+4\delta_{0})^{1/2}\\}^{2}}(\log
p)/s)$
$\displaystyle\leq\exp(-C^{\prime}\frac{2\delta_{0}^{2}}{\\{(1+6\delta_{0})^{1/2}+(1+4\delta_{0})^{1/2}\\}^{2}}\log
p)$
$\displaystyle=p^{-C^{\prime}\frac{2\delta_{0}^{2}}{\\{(1+6\delta_{0})^{1/2}+(1+4\delta_{0})^{1/2}\\}^{2}}}.$

Ultimately it shows that
$\operatorname{\mathbb{P}}(\operatorname{\mathcal{E}}_{\operatorname{\mathcal{I}}_{0}}^{c})\lesssim\exp\left(-C^{*}\frac{\delta_{0}}{1+8\delta_{0}}n\right)+p^{-C^{\prime}\frac{2\delta_{0}^{2}}{\\{(1+6\delta_{0})^{1/2}+(1+4\delta_{0})^{1/2}\\}^{2}}}$.
Now we are ready to show that the error probability goes to $0$.

$\displaystyle\operatorname{\mathbb{P}}_{\beta}(\widehat{\operatorname{\mathcal{S}}}\neq\operatorname{\mathcal{S}})\leq\sum_{\operatorname{\mathcal{I}}_{0}\subset
S}\operatorname{\mathbb{P}}_{\beta}\left(\max_{\operatorname{\mathcal{I}}_{1}:\operatorname{\mathcal{I}}_{1}\in\operatorname{\mathcal{J}}_{\operatorname{\mathcal{I}}_{0}}}\frac{\left\|P_{\operatorname{\mathcal{I}}_{0}\cup\operatorname{\mathcal{I}}_{1}}Y\right\|_{2}^{2}}{s\log
p}>\frac{\left\|P_{\operatorname{\mathcal{S}}}Y\right\|_{2}^{2}}{s\log
p}\right)$
$\displaystyle\leq\sum_{k=0}^{s-1}\sum_{\operatorname{\mathcal{I}}_{0}:\left|\operatorname{\mathcal{I}}_{0}\right|=k}\operatorname{\mathbb{P}}_{\beta}\left(\max_{\operatorname{\mathcal{I}}_{1}:\operatorname{\mathcal{I}}_{1}\in\operatorname{\mathcal{J}}_{\operatorname{\mathcal{I}}_{0}}}\frac{\left\|P_{\operatorname{\mathcal{I}}_{0}\cup\operatorname{\mathcal{I}}_{1}}Y\right\|_{2}^{2}}{s\log
p}>\frac{\left\|P_{\operatorname{\mathcal{S}}}Y\right\|_{2}^{2}}{s\log
p}\right)$
$\displaystyle\leq\sum_{k=0}^{s-1}\sum_{\operatorname{\mathcal{I}}_{0}:\left|\operatorname{\mathcal{I}}_{0}\right|=k}\operatorname{\mathbb{P}}_{\beta}\left(\max_{\operatorname{\mathcal{I}}_{1}:\operatorname{\mathcal{I}}_{1}\in\operatorname{\mathcal{J}}_{\operatorname{\mathcal{I}}_{0}}}\frac{\left\|P_{\operatorname{\mathcal{I}}_{0}\cup\operatorname{\mathcal{I}}_{1}}Y\right\|_{2}^{2}}{s\log
p}>\frac{\left\|P_{\operatorname{\mathcal{S}}}Y\right\|_{2}^{2}}{s\log
p},\operatorname{\mathcal{E}}_{\operatorname{\mathcal{I}}_{0}}\right)+\operatorname{\mathbb{P}}(\operatorname{\mathcal{E}}_{\operatorname{\mathcal{I}}_{0}}^{c})$
$\displaystyle\leq\sum_{k=0}^{s-1}\sum_{\operatorname{\mathcal{I}}_{0}:\left|\operatorname{\mathcal{I}}_{0}\right|=k}\operatorname{\mathbb{P}}_{\beta}\left(\max_{\operatorname{\mathcal{I}}_{1}:\operatorname{\mathcal{I}}_{1}\in\operatorname{\mathcal{J}}_{\operatorname{\mathcal{I}}_{0}}}\frac{\left\|P_{\operatorname{\mathcal{I}}_{0}\cup\operatorname{\mathcal{I}}_{1}}g\right\|_{2}^{2}}{s\log
p}>\frac{\left\|P_{\operatorname{\mathcal{S}}}g\right\|_{2}^{2}}{s\log
p},\operatorname{\mathcal{E}}_{\operatorname{\mathcal{I}}_{0}}\right)+\operatorname{\mathbb{P}}(\operatorname{\mathcal{E}}_{\operatorname{\mathcal{I}}_{0}}^{c})$
$\displaystyle\leq\sum_{k=0}^{s-1}\sum_{\operatorname{\mathcal{I}}_{0}:\left|\operatorname{\mathcal{I}}_{0}\right|=k}\operatorname{\mathbb{P}}_{\beta}\left(\max_{\operatorname{\mathcal{I}}_{1}:\operatorname{\mathcal{I}}_{1}\in\operatorname{\mathcal{J}}_{\operatorname{\mathcal{I}}_{0}}}\frac{\left\|P_{\operatorname{\mathcal{I}}_{0}\cup\operatorname{\mathcal{I}}_{1}}g\right\|_{2}^{2}}{s\log
p}>2(1+4\delta_{0})R_{\operatorname{\mathcal{I}}_{0}}\right)+\operatorname{\mathbb{P}}(\operatorname{\mathcal{E}}_{\operatorname{\mathcal{I}}_{0}}^{c})$
$\displaystyle\overset{(a)}{\lesssim}\sum_{k=0}^{s-1}\sum_{\operatorname{\mathcal{I}}_{0}:\left|\operatorname{\mathcal{I}}_{0}\right|=k}\left[\operatorname{\mathbb{P}}\left(T^{*2}>2(1+4\delta_{0})sR_{\operatorname{\mathcal{I}}_{0}}\log
p\right)+\operatorname{\mathbb{P}}(\operatorname{\mathcal{E}}_{\operatorname{\mathcal{I}}_{0}}^{c})+C(K_{0}K_{1})^{3/4}n^{-1/8}\\{sb_{n}(s,p)\\}^{7/8}\right]$
$\displaystyle\lesssim\sum_{k=0}^{s-1}\sum_{\operatorname{\mathcal{I}}_{0}:\left|\operatorname{\mathcal{I}}_{0}\right|=k}p^{-0.75\delta_{0}}+C\frac{p^{-2\delta_{0}}}{\sqrt{\log
p}}+\exp\left(-C^{*}\frac{\delta_{0}}{1+8\delta_{0}}n\right)+p^{-C^{\prime}\frac{2\delta_{0}^{2}}{\\{(1+6\delta_{0})^{1/2}+(1+4\delta_{0})^{1/2}\\}^{2}}}$
$\displaystyle+n^{-1/8}\\{sb_{n}(s,p)\\}^{7/8}$
$\displaystyle\lesssim\sum_{k=0}^{s-1}\binom{s}{k}\left[C_{\delta_{0}}\frac{p^{-2\delta_{0}}}{\sqrt{\log
p}}+\exp\left(-C^{*}\frac{\delta_{0}}{1+8\delta_{0}}n\right)+p^{-C^{\prime}\frac{2\delta_{0}^{2}}{\\{(1+6\delta_{0})^{1/2}+(1+4\delta_{0})^{1/2}\\}^{2}}}+n^{-1/8}\\{sb_{n}(s,p)\\}^{7/8}\right]$
$\displaystyle\lesssim
2^{s}\left[C_{\delta_{0}}\frac{p^{-2\delta_{0}}}{\sqrt{\log
p}}+\exp\left(-C^{*}\frac{\delta_{0}}{1+8\delta_{0}}n\right)+p^{-C^{\prime}\frac{2\delta_{0}^{2}}{\\{(1+6\delta_{0})^{1/2}+(1+4\delta_{0})^{1/2}\\}^{2}}}+n^{-1/8}\\{sb_{n}(s,p)\\}^{7/8}\right].$

Inequality $(a)$ uses $\tilde{s}b_{n}(\tilde{s},\tilde{p})\leq sb_{n}(s,p)$
for large $p$. Thus if

$s\lesssim\left(\frac{\delta_{0}}{1+8\delta_{0}}\wedge\frac{2\delta_{0}^{2}}{\\{(1+6\delta_{0})^{1/2}+(1+4\delta_{0})^{1/2}\\}^{2}}\wedge\frac{k}{16}\right)\log
p$

, then error probability goes to 0 uniformly over
$\beta\in{\mathcal{M}}_{s}^{a}$.

## References

  * Agarwal et al., (2012) Agarwal, A., Negahban, S., and Wainwright, M. J. (2012).  Fast global convergence of gradient methods for high-dimensional statistical recovery.  The Annals of Statistics, 40(5):2452–2482.
  * Arias-Castro et al., (2011) Arias-Castro, E., Candès, E. J., and Plan, Y. (2011).  Global testing under sparse alternatives: ANOVA, multiple comparisons and the higher criticism.  The Annals of Statistics, 39:2533–2556.
  * Bertsimas et al., (2016) Bertsimas, D., King, A., and Mazumder, R. (2016).  Best subset selection via a modern optimization lens.  The Annals of Statistics, 44(2):813–852.
  * Bertsimas and Van Parys, (2020) Bertsimas, D. and Van Parys, B. (2020).  Sparse high-dimensional regression: Exact scalable algorithms and phase transitions.  The Annals of Statistics, 48(1):300–323.
  * Bickel et al., (2009) Bickel, P. J., Ritov, Y., and Tsybakov, A. B. (2009).  Simultaneous analysis of lasso and dantzig selector.  The Annals of Statistics, 37(4):1705–1732.
  * Blumensath and Davies, (2009) Blumensath, T. and Davies, M. E. (2009).  Iterative hard thresholding for compressed sensing.  Applied and Computational Harmonic Analysis, 27(3):265–274.
  * Candès et al., (2006) Candès, E. J., Romberg, J., and Tao, T. (2006).  Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information.  IEEE Transactions on information theory, 52(2):489–509.
  * Candes and Tao, (2006) Candes, E. J. and Tao, T. (2006).  Near-optimal signal recovery from random projections: Universal encoding strategies?  IEEE transactions on information theory, 52(12):5406–5425.
  * Chernozhukov et al., (2014) Chernozhukov, V., Chetverikov, D., and Kato, K. (2014).  Gaussian approximation of suprema of empirical processes.  The Annals of Statistics, 42(4):1564–1597.
  * Consortium et al., (2007) Consortium, W. T. C. C. et al. (2007).  Genome-wide association study of 14,000 cases of seven common diseases and 3,000 shared controls.  Nature, 447(7145):661.
  * Donoho, (2006) Donoho, D. L. (2006).  Compressed sensing.  IEEE Transactions on information theory, 52(4):1289–1306.
  * Fan et al., (2020) Fan, J., Guo, Y., and Zhu, Z. (2020).  When is best subset selection the” best”?  arXiv preprint arXiv:2007.01478.
  * Fan and Li, (2001) Fan, J. and Li, R. (2001).  Variable selection via nonconcave penalized likelihood and its oracle properties.  Journal of the American Statistical Association, 96(456):1348–1360.
  * Fan and Lv, (2008) Fan, J. and Lv, J. (2008).  Sure independence screening for ultrahigh dimensional feature space.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 70(5):849–911.
  * Fan et al., (2018) Fan, J., Shao, Q.-M., and Zhou, W.-X. (2018).  Are discoveries spurious? Distributions of maximum spurious correlations and their applications.  The Annals of Statistics, 46(3):989.
  * Fletcher et al., (2009) Fletcher, A. K., Rangan, S., and Goyal, V. K. (2009).  Necessary and sufficient conditions for sparsity pattern recovery.  IEEE Transactions on Information Theory, 55(12):5758–5772.
  * Genovese et al., (2012) Genovese, C. R., Jin, J., Wasserman, L., and Yao, Z. (2012).  A comparison of the lasso and marginal regression.  Journal of Machine Learning Research, 13:2107–2143.
  * Guo et al., (2015) Guo, X., Zhang, H., Wang, Y., and Wu, J.-L. (2015).  Model selection and estimation in high dimensional regression models with group SCAD.  Statistics & Probability Letters, 103:86–92.
  * Hastie et al., (2020) Hastie, T., Tibshirani, R., and Tibshirani, R. (2020).  Best subset, forward stepwise or lasso? analysis and recommendations based on extensive comparisons.  Statistical Science, 35(4):579–592.
  * Hong et al., (2018) Hong, H. G., Kang, J., and Li, Y. (2018).  Conditional screening for ultra-high dimensional covariates with survival outcomes.  Lifetime Data Analysis, 24(1):45–71.
  * Huang et al., (2014) Huang, D., Li, R., and Wang, H. (2014).  Feature screening for ultrahigh dimensional categorical data with applications.  Journal of Business & Economic Statistics, 32(2):237–244.
  * Huang et al., (2019) Huang, T.-J., McKeague, I. W., and Qian, M. (2019).  Marginal screening for high-dimensional predictors of survival outcomes.  Statistica Sinica, 29(4):2105.
  * Jain et al., (2014) Jain, P., Tewari, A., and Kar, P. (2014).  On iterative hard thresholding methods for high-dimensional M-estimation.  In Advances in Neural Information Processing Systems, volume 27.
  * Ji et al., (2012) Ji, P., Jin, J., et al. (2012).  Ups delivers optimal phase diagram in high-dimensional variable selection.  The Annals of Statistics, 40(1):73–103.
  * Jin et al., (2014) Jin, J., Zhang, C.-H., and Zhang, Q. (2014).  Optimality of graphlet screening in high dimensional variable selection.  Journal of Machine Learning Research, 15(1):2723–2772.
  * Leisenring et al., (1997) Leisenring, W., Pepe, M., and Longton, G. (1997).  A marginal regression modelling framework for evaluating medical diagnostic tests.  Statistics in Medicine, 16(11):1263—1281.
  * Li et al., (2016) Li, J., Zheng, Q., Peng, L., and Huang, Z. (2016).  Survival impact index and ultrahigh-dimensional model-free screening with survival outcomes.  Biometrics, 72(4):1145–1154.
  * Li et al., (2019) Li, Y., Hong, H. G., Ahmed, S. E., and Li, Y. (2019).  Weak signals in high-dimensional regression: Detection, estimation and prediction.  Applied Stochastic Models in Business and Industry, 35(2):283–298.
  * Liu and Foygel Barber, (2020) Liu, H. and Foygel Barber, R. (2020).  Between hard and soft thresholding: optimal iterative thresholding algorithms.  Information and Inference: A Journal of the IMA, 9(4):899–933.
  * Lu, (2005) Lu, W. (2005).  Marginal regression of multivariate event times based on linear transformation models.  Lifetime Data Analysis, 11(3):389–404.
  * Marttinen et al., (2013) Marttinen, P., Gillberg, J., Havulinna, A., Corander, J., and Kaski, S. (2013).  Genome-wide association studies with high-dimensional phenotypes.  Statistical Applications in Gnetics and Molecular Biology, 12(4):413–431.
  * Ndaoud and Tsybakov, (2020) Ndaoud, M. and Tsybakov, A. B. (2020).  Optimal variable selection and adaptive noisy compressed sensing.  IEEE Transactions on Information Theory, 66(4):2517–2532.
  * Robins et al., (2003) Robins, J. M., Scheines, R., Spirtes, P., and Wasserman, L. (2003).  Uniform consistency in causal inference.  Biometrika, 90(3):491–515.
  * Su et al., (2017) Su, W., Bogdan, M., and Candes, E. (2017).  False discoveries occur early on the lasso path.  The Annals of statistics, pages 2133–2150.
  * Tibshirani, (1996) Tibshirani, R. (1996).  Regression shrinkage and selection via the lasso.  Journal of the Royal Statistical Society: Series B (Methodological), 58(1):267–288.
  * Tikhonov, (1943) Tikhonov, A. (1943).  On the stability of inverse problems.  Proceedings of the USSR Academy of Sciences, 39:195–198.
  * Van De Geer and Bühlmann, (2009) Van De Geer, S. A. and Bühlmann, P. (2009).  On the conditions used to prove oracle results for the lasso.  Electronic Journal of Statistics, 3:1360–1392.
  * Wainwright, (2009) Wainwright, M. J. (2009).  Information-theoretic limits on sparsity recovery in the high-dimensional and noisy setting.  IEEE Transactions on Information Theory, 55(12):5728–5741.
  * Wang et al., (2022) Wang, C. D., Chen, Z., Lian, Y., and Chen, M. (2022).  Asset selection based on high frequency sharpe ratio.  Journal of Econometrics, to appear.
  * Wang et al., (2010) Wang, W., Wainwright, M. J., and Ramchandran, K. (2010).  Information-theoretic limits on sparse signal recovery: Dense versus sparse measurement matrices.  IEEE Transactions on Information Theory, 56(6):2967–2979.
  * Wasserman and Roeder, (2009) Wasserman, L. and Roeder, K. (2009).  High dimensional variable selection.  The Annals of Statistics, 37(5A):2178.
  * Yang et al., (2020) Yang, S., Wen, J., Eckert, S. T., Wang, Y., Liu, D. J., Wu, R., Li, R., and Zhan, X. (2020).  Prioritizing genetic variants in gwas with lasso using permutation-assisted tuning.  Bioinformatics, 36(12):3811–3817.
  * Zhang, (2010) Zhang, C.-H. (2010).  Nearly unbiased variable selection under minimax concave penalty.  The Annals of Statistics, 38(2):894–942.
  * Zhang et al., (2014) Zhang, Y., Wainwright, M. J., and Jordan, M. I. (2014).  Lower bounds on the performance of polynomial-time algorithms for sparse linear regression.  In Conference on Learning Theory, pages 921–948. PMLR.
  * Zhao and Yu, (2006) Zhao, P. and Yu, B. (2006).  On model selection consistency of lasso.  The Journal of Machine Learning Research, 7:2541–2563.
  * Zhu et al., (2020) Zhu, J., Wen, C., Zhu, J., Zhang, H., and Wang, X. (2020).  A polynomial algorithm for best-subset selection problem.  Proceedings of the National Academy of Sciences, 117(52):33117–33123.
  * Zhu and Wu, (2021) Zhu, Z. and Wu, S. (2021).  On the early solution path of best subset selection.  arXiv preprint arXiv:2107.06939.

##  9 Appendix

###  9.1 Proof of Theorem 3

We first present a result form Wang et al., (2010) which is gives us necessary
condition for asymptotic exact recovery.

######  Theorem 5 (Wang et al., (2010))

Consider the model (1) with the design matrix $X\in{\mathbb{R}}^{n\times p}$
be drawn with i.i.d elements from any distribution with zero mean and unit
variance. Let
$a:=\min_{j\in\operatorname{\mathcal{S}}(\beta)}\left|\beta_{j}\right|$, i.e.,
it denote the minimum signal strength of $\beta$. Define the function

$f_{m}(p,s,a):=\frac{\log\binom{p-s+m}{m}-1}{\frac{1}{2}\log\left(1+ma^{2}(1-\frac{m}{p-s+m})\right)},\quad
1\leq m\leq s.$

Then $n\geq\max\\{f_{1}(p,s,a),\ldots,f_{s}(p,s,a),s\\}$ is necessary for
asymptotic exact recovery.

In the light of Theorem 5 the proof of Theorem 3 is follows immediately. To
see this note that if $r<1$ then there exists $\alpha\in(0,1)$ such that
$r=1-\alpha$. Also recall that $a=\\{2r(\log p)/n\\}^{1/2},s=O(\log p)$ and
$n=\left\lfloor p^{k}\right\rfloor$. Note that
$f_{1}(p,s,a)/n\sim\frac{1}{1-\alpha}$. This shows that asymptotically the
necessary condition in above theorem is violated and hence $r\geq 1$ is
necessary.

###  9.2 Proof of Theorem 4

We first introduce some standard assumptions analyzing ETS.

######  Definition 1 (RSC property)

A differentiable function $F:{\mathbb{R}}^{p}\to{\mathbb{R}}$ is said to
satisfy restricted strong convexity (RSC) at sparsity level $s=s_{1}+s_{2}$
with strong convexity constraint $\ell_{s}$ if the following holds for all
$\theta_{1},\theta_{2}$ s.t. $\left\|\theta_{1}\right\|_{0}\leq s_{1}$ and
$\left\|\theta_{2}\right\|_{0}\leq s_{2}$:

$F(\theta_{1})-F(\theta_{2})\geq\left\langle\theta_{1}-\theta_{2},\nabla_{\theta}F(\theta_{2})\right\rangle+\frac{\ell_{s}}{2}\left\|\theta_{1}-\theta_{2}\right\|_{2}^{2}.$

######  Definition 2 (RSS property)

A differentiable function $F:{\mathbb{R}}^{p}\to{\mathbb{R}}$ is said to
satisfy restricted strong smoothness (RSS) at sparsity level $s=s_{1}+s_{2}$
with strong smoothness constraint $L_{s}$ if the following holds for all
$\theta_{1},\theta_{2}$ s.t. $\left\|\theta_{1}\right\|_{0}\leq s_{1}$ and
$\left\|\theta_{2}\right\|_{0}\leq s_{2}$:

$F(\theta_{1})-F(\theta_{2})\leq\left\langle\theta_{1}-\theta_{2},\nabla_{\theta}F(\theta_{2})\right\rangle+\frac{L_{s}}{2}\left\|\theta_{1}-\theta_{2}\right\|_{2}^{2}.$

Now we quote an important theorem from Jain et al., (2014) which quantifies
the sub-optimality gap of Algorithm 1.

######  Theorem 6 (Jain et al., (2014))

Let $F$ has RSC and RSS parameters given by $\ell_{2\hat{s}+s}(F)=\alpha$ and
$L_{2\hat{s}+\hat{\pi}}(F)=L$ respectively. Call Algorithm 1 with $\hat{s}\geq
32L^{2}\ell^{-2}s$ and $h=2/(3L)$. Also let
$\hat{\beta}={\arg\min}_{\theta,\left\|\theta\right\|_{0}\leq s}F(\theta)$.
Then $t{\rm th}$ iterate of Algorithm 1 for
$t=O(L\ell^{-1}\log(F(\beta^{(0)})/\epsilon))$ satisfies:

$F(\beta^{(t)})-F(\hat{\beta})\leq\epsilon.$

In our setup, the observations $\\{x_{i}\\}_{i=1}^{n}$ are coming from i.i.d.
mean zero isotropic Gaussian distribution. Thus, lemma 6 from Agarwal et al.,
(2012) immediately tells that RSC and RSS at any sparsity level $m$ hold for
$f_{n_{1}}(\cdot)$ with probability at least $1-\exp(-c_{0}n_{1})$ with
$\ell_{m}=\frac{1}{2}-c_{1}(m\log p)/n_{1}$ and $L_{m}=2+c_{1}(m\log
p)/n_{1}$, where $c_{0},c_{1}$ are universal constants. Now set $m=2\hat{s}+s$
and recall that $n_{1}\sim\gamma p^{k}$. If $n_{1}>4c_{1}(2\hat{s}+s)\log p$
then we have $\ell_{m}\geq 1/4$ and $L_{m}\leq 9/4$, which means that
$L_{m}/(9\ell_{m})\leq 1$. Thus to apply Theorem 6 it is enough to choose
$\hat{s}=2592s$. Also by the assumption on $n$ for large $p$ we have
$n_{1}>4c_{1}(2\hat{s}+s)\log p$. Let
$f_{n_{1}}(\theta):=n_{1}^{-1}\|Y^{(1)}-X^{(1)}\theta\|_{2}^{2}$ for
$\theta\in{\mathbb{R}}^{p}$. Note that
$f_{n_{1}}(0)=n_{1}^{-1}\|Y^{(1)}\|_{2}^{2}\overset{\rm
d}{=}(1+\left\|\beta\right\|_{2}^{2})V_{n_{1}}/n_{1}$, where $V_{n_{1}}$ is
chi-square random variable with $n_{1}$ degrees of freedom. Also by
Bernstein’s type inequality it follows that
$\left|(V_{n_{1}}/n_{1})-1\right|\leq 1/2$ with probability at least
$1-\exp(-c_{4}n_{1})$, where $c_{4}$ is a universal positive constant. Thus if
$t=O(L_{m}\ell_{m}^{-1}\log((1+\left\|\beta\right\|_{2}^{2})/\epsilon))$, then
we have $f_{n_{1}}(\beta^{(t)})-f_{n_{1}}(\hat{\beta})\leq\epsilon$. Thus by
Theorem 3 of Jain et al., (2014) it follows that with probability at lest
$1-\exp(-c_{0}n_{1})-\exp(-c_{4}n_{1})-c_{2}p^{-c_{3}}$ ($c_{2},c_{3}$ are
universal constants) we have

$\|\beta^{(t)}-\beta\|_{2}\leq C\left\\{\left(\frac{s\log
p}{n_{1}}\right)^{1/2}+\epsilon^{1/2}\right\\},$

for some universal positive constant $C$. In practice one can use the full
data in both steps of ETS but for theoretical convenience we use data
splitting in the algorithm. The key idea here is to use one subsample for
obtaining $\hat{\beta}^{\rm iht}$ and the remaining subsample in the screening
step. Due to this $\hat{\beta}^{\rm iht}$ becomes independent of the second
subsample, thus making the analysis of ETS easy.

For notational brevity we write $\eta$ instead of $\eta_{\beta}$. Define the
event $\operatorname{\mathcal{H}}=\left\\{\|\hat{\beta}^{\rm
iht}-\beta\|_{2}\leq\xi_{p}\right\\}$, where $\xi_{p}=C\left\\{(\frac{s\log
p}{n_{1}})^{1/2}+\epsilon^{1/2}\right\\}\leq(A\epsilon)^{1/2}$ (for some
universal constant $A>0$) and $\hat{\beta}^{\rm iht}$ is based on the
subsample $\operatorname{\mathcal{D}}_{1}$. From Theorem 3 of Jain et al.,
(2014), it follows that
$\operatorname{\mathbb{P}}(\operatorname{\mathcal{H}}^{c})\leq\exp(-c_{0}n_{1})+\exp(-c_{4}n_{1})+c_{2}p^{-c_{3}}$.

Next, for algebraic convenience we again reparametrize $\delta$ as
$8\delta_{0}$ and set $\epsilon=6\delta_{0}/A,\varsigma=(1+A\epsilon)^{1/2}$.
Now note that for any $\beta\in{\mathcal{M}}_{a}^{s}$, we have

$\displaystyle\operatorname{\mathbb{P}}_{\beta}(\hat{\eta}\neq\eta|\operatorname{\mathcal{D}}_{1})$
$\displaystyle\leq\sum_{j:\beta_{j}=0}\operatorname{\mathbb{P}}(\hat{\eta}_{j}=1,\operatorname{\mathcal{H}}|\operatorname{\mathcal{D}}_{1})+\sum_{j:\beta_{j}\neq
0}\operatorname{\mathbb{P}}_{\beta}(\hat{\eta}_{j}\neq
1,\operatorname{\mathcal{H}}|\operatorname{\mathcal{D}}_{1})+\operatorname{\mathbb{P}}(\operatorname{\mathcal{H}}^{c}|\operatorname{\mathcal{D}}_{1})$
$\displaystyle=\sum_{j:\beta_{j}=0}\operatorname{\mathbb{P}}_{\beta}(\left|\Delta_{j}\right|>\kappa_{\varsigma}(X_{j}^{(2)}),\operatorname{\mathcal{H}}|\operatorname{\mathcal{D}}_{1})+\sum_{j:\beta_{j}\neq
0}\operatorname{\mathbb{P}}_{\beta}(\left|\Delta_{j}\right|\leq\kappa_{\varsigma}(X_{j}^{(2)}),\operatorname{\mathcal{H}}|\operatorname{\mathcal{D}}_{1})+\operatorname{\mathbb{P}}(\operatorname{\mathcal{H}}^{c}|\operatorname{\mathcal{D}}_{1}),$

Using the fact that conditionally on $\hat{\beta}$ and $X_{j}$, the random
variable $\Delta_{j}$ has the same distribution as (10), we conclude that for
all $j\notin\operatorname{\mathcal{S}}(\beta)$ ,

$\displaystyle\operatorname{\mathbb{P}}(\eta_{j}=1,\operatorname{\mathcal{H}}|\operatorname{\mathcal{D}}_{1})$
$\displaystyle\leq\operatorname{\mathbb{P}}\left((1+\xi_{p}^{2})^{1/2}\left|g_{j}\right|>\frac{a\left\|X_{j}^{(2)}\right\|_{2}}{2}+\frac{(1+A\epsilon)\log
p}{a\left\|X_{j}^{(2)}\right\|_{2}},\operatorname{\mathcal{H}}\big{|}\operatorname{\mathcal{D}}_{1}\right)$
$\displaystyle\leq\operatorname{\mathbb{P}}\left((1+A\epsilon)^{1/2}\left|g_{j}\right|>\frac{a\left\|X_{j}^{(2)}\right\|_{2}}{2}+\frac{(1+A\epsilon)\log
p}{a\left\|X_{j}^{(2)}\right\|_{2}},\operatorname{\mathcal{H}}\big{|}\operatorname{\mathcal{D}}_{1}\right)$
$\displaystyle=2\operatorname{\mathbb{E}}\left\\{\overline{\Phi}\left(\frac{a\left\|X_{j}^{(2)}\right\|_{2}}{2(1+A\epsilon)^{1/2}}+\frac{(1+A\epsilon)^{1/2}\log
p}{a\left\|X_{j}^{(2)}\right\|_{2}}\right)\right\\}.$

Here $\overline{\Phi}(\cdot)$ denotes the survival function of standard
Gaussian random variable. Now note that for each $j$ we have
$\left\|X_{j}\right\|_{2}^{2}\overset{\rm d}{=}V_{n_{2}}$, where $V_{n_{2}}$
is a chi-squared random variable with $n_{2}$ degrees of freedom. Thus we have

$\operatorname{\mathbb{P}}(\eta_{j}=1,\operatorname{\mathcal{H}}|\operatorname{\mathcal{D}}_{1})\leq
2\operatorname{\mathbb{E}}\left\\{\overline{\Phi}\left(\frac{aV_{n_{2}}^{1/2}}{2(1+A\epsilon)^{1/2}}+\frac{(1+A\epsilon)^{1/2}\log
p}{aV_{n_{2}}^{1/2}}\right)\right\\}.$

Analogous argument and the fact that $\left|\beta_{j}\right|\geq a$ for all
$\beta_{j}\neq 0$, leads to the fact that for all
$j\in\operatorname{\mathcal{S}}(\beta)$,

$\operatorname{\mathbb{P}}(\eta_{j}\neq
1,\operatorname{\mathcal{H}}|\operatorname{\mathcal{D}}_{1})\leq
2\mathbb{E}\left\\{\overline{\Phi}\left(\max\left\\{\frac{aV_{n_{2}}^{1/2}}{2(1+A\epsilon)^{1/2}}-\frac{(1+A\epsilon)^{1/2}\log
p}{aV_{n_{2}}^{1/2}},0\right\\}\right)\right\\}.$

Now recall that $\epsilon=\frac{6\delta_{0}}{A}$ and
$\gamma\in(0,\frac{\delta_{0}}{1+8\delta_{0}})$. With this choice of tuning
parameters it is easy to see that
$r(1-\gamma)/(1+A\epsilon)\geq\frac{1+7\delta_{0}}{1+6\delta_{0}}>1$ and hence
as $p\to\infty$ we have

$\displaystyle W_{n_{2}}$ $\displaystyle:=\frac{1}{(\log
p)^{1/2}}\left(\frac{aV_{n_{2}}^{1/2}}{2(1+A\epsilon)^{1/2}}-\frac{(1+A\epsilon)^{1/2}\log
p}{aV_{n_{2}}^{1/2}}\right)$ $\displaystyle\overset{\rm
p}{\longrightarrow}\frac{1}{(2r)^{1/2}}\left\\{r\left(\frac{1-\gamma}{1+A\epsilon}\right)^{1/2}-\left(\frac{1+A\epsilon}{1-\gamma}\right)^{1/2}\right\\}>0.$

The above display uses the fact that $n_{2}/n\to 1-\gamma$ and
$V_{n_{2}}/n_{2}\overset{\rm p}{\to}1$ as $p\to\infty$. Next let is define the
following quantity $q$:

$q:=q(A,\epsilon,\delta_{0},\gamma)=\frac{1}{\\{2(1+8\delta_{0})\\}^{1/2}}\left\\{(1+8\delta_{0})\left(\frac{1-\gamma}{1+A\epsilon}\right)^{1/2}-\left(\frac{1+A\epsilon}{1-\gamma}\right)^{1/2}\right\\}.$

Due to choice of $\epsilon$ and $\gamma$ it is easy to show $q>0$. Now define
the event $G_{n_{2}}:=\\{W_{n_{2}}>q/2\\}$. Before we proceed it is useful to
note the following:

$W_{n_{2}}=\frac{1}{(2r)^{1/2}}\left(\frac{r\\{V_{n_{2}}/(n(1-\gamma))\\}^{1/2}(1-\gamma)^{1/2}}{(1+A\epsilon)^{1/2}}-\frac{(1+A\epsilon)^{1/2}}{\\{V_{n_{2}}/(n(1-\gamma))\\}^{1/2}(1-\gamma)^{1/2}}\right).$

Next define the function

$H(u):=\frac{1}{(2+16\delta_{0})^{1/2}}\left\\{u(1+8\delta_{0})\left(\frac{1-\gamma}{1+A\epsilon}\right)^{1/2}-\frac{1}{u}\left(\frac{1+A\epsilon}{1-\gamma}\right)^{1/2}\right\\},\quad
u>0.$

As $r=1+8\delta_{0}$ we have
$W_{n_{2}}=H(\\{V_{n_{2}}/(n(1-\gamma))\\}^{1/2})$. It is also easy to see
that $H(\cdot)$ is strictly increasing function on $(0,\infty)$ and $H(1)=q$.
Hence $\lambda_{\delta_{0}}:=H^{-1}(q/2)\in(0,1)$. Now
$G_{n_{2}}^{c}=\\{W_{n_{2}}\leq
q/2\\}\subseteq\\{H(\\{V_{n_{2}}/(n(1-\gamma))\\}^{1/2})\leq q/2\\}$. Thus a
straight forward calculation shows that that,

$\operatorname{\mathbb{P}}(G_{n_{2}}^{c})\leq\operatorname{\mathbb{P}}\left(\frac{V_{n_{2}}}{n_{2}}\leq\frac{n(1-\gamma)}{n_{2}}\lambda_{\delta_{0}}^{2}\right).$

Choose $p$ large enough such that $n_{2}/n>\lambda_{\delta_{0}}(1-\gamma)$ and
hence we have,

$\operatorname{\mathbb{P}}(G_{n_{2}}^{c})\leq\operatorname{\mathbb{P}}\left(\frac{V_{n_{2}}}{n_{2}}\leq\lambda_{\delta_{0}}\right)\lesssim\exp(-K_{\delta_{0}}n_{2}),$

where $K_{\delta_{0}}$ is constant depending on $\delta_{0}$. Note that
$\overline{\Phi}(t)\leq e^{-t^{2}/2}$ for all $t>0$. Using this fact we have
the following:

$\operatorname{\mathbb{P}}(\eta_{j}=1,\operatorname{\mathcal{H}}|\operatorname{\mathcal{D}}_{1})\leq\operatorname{\mathbb{E}}\left[\exp\left\\{-\left(1+\frac{W_{n_{2}}^{2}}{2}\right)\log
p\right\\}\operatorname{\mathbbm{1}}_{G_{n_{2}}}\right]+\operatorname{\mathbb{P}}(G_{n_{2}}^{c})\lesssim
p^{-(1+q^{2}/8)}+\exp(-K_{\delta_{0}}n_{2}),$

for all $j\notin\operatorname{\mathcal{S}}(\beta)$. Similarly,

$\operatorname{\mathbb{P}}(\eta_{j}\neq
1,\operatorname{\mathcal{H}}|\operatorname{\mathcal{D}}_{1})\lesssim
p^{-q^{2}/8}+\exp(-K_{\delta_{0}}n_{2}),\quad\forall
j\in\operatorname{\mathcal{S}}(\beta).$

Thus marginalizing out $\operatorname{\mathcal{D}}_{1}$ and summing over all
$j$ we get,

$\sup_{\beta\in{\mathcal{M}}_{a}^{s}}\operatorname{\mathbb{P}}_{\beta}(\hat{\eta}\neq\eta)\lesssim
p^{-q^{2}/8}\log
p+p\exp(-K_{\delta_{0}}n_{2})+\exp(-c_{0}n_{1})+\exp(-c_{4}n_{1})+c_{2}p^{-c_{3}}.$

Now taking $p\to\infty$ we have the result.

###  9.3 Proof of Corollary 12

Similar to previous proofs, we reparametrize $\delta$ by $8\delta_{0}$ and set
$\epsilon=6\delta_{0}/A,\varsigma=(1+A\epsilon)^{1/2}$. Now note that, it is
enough to prove the following:

$\lim_{p\to\infty}\inf_{\beta\in{\mathcal{M}}_{s}^{a}}\operatorname{\mathbb{P}}_{\beta}\left(\max_{j\notin\operatorname{\mathcal{S}}(\beta)}\left|\Delta_{j}\right|<\min_{j\in\operatorname{\mathcal{S}}(\beta)}\left|\Delta_{j}\right|\right)\to
1$

as $p\to\infty$. To this end first define the following quantity:

$t_{p}:=\frac{\left(\frac{2rn_{2}\log
p}{n}\right)^{1/2}}{2}+\frac{\varsigma^{2}\log p}{\left(\frac{2rn_{2}\log
p}{n}\right)^{1/2}}.$

We will show that
$\lim_{p\to\infty}\inf_{\beta\in{\mathcal{M}}_{s}^{a}}\operatorname{\mathbb{P}}_{\beta}\left(\min_{j\in\operatorname{\mathcal{S}}(\beta)}\left|\Delta_{j}\right|>t_{p},\max_{j\notin\operatorname{\mathcal{S}}(\beta)}\left|\Delta_{j}\right|\leq
t_{p}\right)\to 1$ as $p\to\infty$. For convenience let us define the events
$G_{\min}:=\\{\min_{j\in\operatorname{\mathcal{S}}(\beta)}\left|\Delta_{j}\right|>t_{p}\\}$
and
$G_{\max}:=\\{\max_{j\notin\operatorname{\mathcal{S}}(\beta)}\left|\Delta_{j}\right|\leq
t_{p}\\}$. Let $\operatorname{\mathcal{H}}$ be the event as defined in Section
9.2. First we will analyze $\operatorname{\mathbb{P}}_{\beta}(G_{\min}^{c})$.
Note that
$\operatorname{\mathbb{P}}_{\beta}(G_{\min}^{c})\leq\operatorname{\mathbb{P}}_{\beta}(G_{\min}^{c}\cap\operatorname{\mathcal{H}})+\operatorname{\mathbb{P}}_{\beta}(\operatorname{\mathcal{H}}^{c})$.
Now the second term goes to $0$ uniformly over
$\beta\in{\mathcal{M}}_{s}^{a}$. Also using Equation (10) under the event
$\operatorname{\mathcal{H}}$ we get

$\displaystyle\sup_{\beta\in{\mathcal{M}}_{s}^{a}}\operatorname{\mathbb{P}}_{\beta}(G_{\min}^{c}\cap\operatorname{\mathcal{H}})$
$\displaystyle\leq\sup_{\beta\in{\mathcal{M}}_{s}^{a}}\operatorname{\mathbb{P}}_{\beta}\left(\min_{j\in\operatorname{\mathcal{S}}(\beta)}\left|\beta_{j}\|X_{j}^{(2)}\|_{2}+(1+A\epsilon)^{1/2}g_{j}\right|\leq
t_{p}\right)$
$\displaystyle\leq\sup_{\beta\in{\mathcal{M}}_{s}^{a}}\operatorname{\mathbb{P}}_{\beta}\left(\max_{j\in\operatorname{\mathcal{S}}(\beta)}\frac{\left|g_{j}\right|}{(\log
p)^{1/2}}\geq\frac{1}{(1+A\epsilon)^{1/2}(\log
p)^{1/2}}\left\\{a\min_{j\in\operatorname{\mathcal{S}}(\beta)}\|X_{j}^{(2)}\|_{2}-t_{p}\right\\}\right)$
$\displaystyle\leq\sup_{\beta\in{\mathcal{M}}_{s}^{a}}\operatorname{\mathbb{P}}_{\beta}\left(\max_{j\in\operatorname{\mathcal{S}}(\beta)}\frac{\left|g_{j}\right|}{(\log
p)^{1/2}}\geq\frac{1}{(1+A\epsilon)^{1/2}(\log
p)^{1/2}}\left\\{a\min_{j\in[p]}\|X_{j}^{(2)}\|_{2}-t_{p}\right\\}\right)$

where $\\{g_{j}\\}_{j\in\operatorname{\mathcal{S}}(\beta)}$ are non i.i.d.
standard Gaussian. Note that
$\left|\operatorname{\mathcal{S}}(\beta)\right|=O(\log p)$. Hence
$\max_{j\in\operatorname{\mathcal{S}}(\beta)}\left|g_{j}\right|=O_{\operatorname{\mathbb{P}}}(\log\log
p)$, which tells that

$\max_{j\in\operatorname{\mathcal{S}}(\beta)}\frac{\left|g_{j}\right|}{(\log
p)^{1/2}}\overset{\rm p}{\to}0.$

Also using lemma 3 from Fletcher et al., (2009) we have

$\frac{1}{(1+A\epsilon)^{1/2}(\log
p)^{1/2}}\left(a\min_{j\in[p]}\|X_{j}^{(2)}\|_{2}-t_{p}\right)\overset{\rm
p}{\to}\frac{1}{(2r)^{1/2}}\left\\{r\left(\frac{1-\gamma}{1+A\epsilon}\right)^{1/2}-\left(\frac{1+A\epsilon}{1-\gamma}\right)^{1/2}\right\\}.$

The right hand side of the above display is at least
$q(A,\epsilon,\delta_{0},\gamma)$ (defined in Section 9.2) which is strictly
positive. Again for compactness we use $q$ instead of
$q(A,\epsilon,\delta_{0},\gamma)$. The above display motivates us to define
the following event:

$\operatorname{\mathcal{E}}_{p}=\left\\{\frac{1}{(1+A\epsilon)^{1/2}(\log
p)^{1/2}}\left(a\min_{j\in[p]}\|X_{j}^{(2)}\|_{2}-t_{p}\right)\geq
q/2\right\\},$

and it follows that
$\operatorname{\mathbb{P}}(\operatorname{\mathcal{E}}_{p}^{c})\to 0$ as
$p\to\infty$. This leads to the following inequality:

$\displaystyle\sup_{\beta\in{\mathcal{M}}_{s}^{a}}\operatorname{\mathbb{P}}_{\beta}(G_{\min}^{c}\cap\operatorname{\mathcal{H}})\leq$
$\displaystyle\sup_{\beta\in{\mathcal{M}}_{s}^{a}}\operatorname{\mathbb{P}}_{\beta}\left(\max_{j\in\operatorname{\mathcal{S}}(\beta)}\frac{\left|g_{j}\right|}{(\log
p)^{1/2}}\geq
q/2\right)+\operatorname{\mathbb{P}}(\operatorname{\mathcal{E}}_{p}^{c})$
$\displaystyle\lesssim p^{-q^{2}/8}\log
p+\operatorname{\mathbb{P}}(\operatorname{\mathcal{E}}_{p}^{c})\to 0.$

Thus we have
$\sup_{\beta\in{\mathcal{M}}_{a}^{a}}\operatorname{\mathbb{P}}_{\beta}(G_{\min}^{c})\to
0$. Similarly it can be shown that
$\sup_{\beta\in{\mathcal{M}}_{s}^{a}}\operatorname{\mathbb{P}}_{\beta}(G_{\max}^{c})\to
0$ as $p\to\infty$. These two claims together completes the proof.