# Understanding Best Subset Selection: A Tale of Two C(omplex)ities

Saptarshi Roy ⋆ Ambuj Tewari ⋆ Ziwei Zhu †
⋆ University of Michigan  Ann Arbor  USA
† Radix Trading  Chicago  USA

###### Abstract

For decades, best subset selection (BSS) has eluded statisticians mainly due
to its computational bottleneck. However, until recently, modern computational
breakthroughs have rekindled theoretical interest in BSS and have led to new
findings. Recently, Guo et al., (2020) showed that the model selection
performance of BSS is governed by a margin quantity that is robust to the
design dependence, unlike modern methods such as LASSO, SCAD, MCP, etc.
Motivated by their theoretical results, in this paper, we also study the
variable selection properties of best subset selection for high-dimensional
sparse linear regression setup. We show that apart from the identifiability
margin, the following two complexity measures play a fundamental role in
characterizing the margin condition for model consistency: (a) complexity of
_residualized features_ , (b) complexity of _spurious projections_. In
particular, we establish a simple margin condition that depends only on the
identifiability margin and the dominating one of the two complexity measures.
Furthermore, we show that a margin condition depending on similar margin
quantity and complexity measures is also necessary for model consistency of
BSS. For a broader understanding, we also consider some simple illustrative
examples to demonstrate the variation in the complexity measures that refines
our theoretical understanding of the model selection performance of BSS under
different correlation structures.

Keywords: Variable Selection; High-Dimensional Statistics; Model Consistency.

##  1 Introduction

Variable selection in high-dimensional sparse regression has been one of the
central topics in statistical research over the past few decades. Consider $n$
independent and identically distributed (i.i.d.) observations
$\\{(\mathbf{x}_{i},y_{i})\\}_{i=1}^{n}$ following the linear model:

$y_{i}=\mathbf{x}_{i}^{\top}{\boldsymbol{\beta}}+\varepsilon_{i},\quad
i\in\\{1,\ldots,n\\},$ (1.1)

where $\\{\mathbf{x}_{i}\\}_{i\in[n]}$ are fixed $p$-dimensional feature
vectors, $\\{\varepsilon_{i}\\}_{i\in[n]}$ are i.i.d. mean-zero $\sigma$-sub-
Gaussian noise, i.e.,
$\mathbb{E}\exp(\lambda\varepsilon_{i})\leq\exp(\lambda^{2}\sigma^{2}/2)$ for
all $\lambda\in\mathbb{R}$ and $i\in[n]$, and the signal vector
${\boldsymbol{\beta}}\in\mathbb{R}^{p}$ is unknown but is assumed to have a
sparse support. In matrix notation, the observations can be represented as

$\mathbf{y}=\mathbf{X}{\boldsymbol{\beta}}+{\boldsymbol{\varepsilon}},$

where $\mathbf{y}=(y_{1},\ldots,y_{n})^{\top}$,
$\mathbf{X}=(\mathbf{x}_{1},\ldots,\mathbf{x}_{n})^{\top}$, and
${\boldsymbol{\varepsilon}}=(\varepsilon_{1},\ldots,\varepsilon_{n})^{\top}$.
We consider the standard _high-dimensional sparse_ setup where $n<p$, and
possibly $n\ll p$, and the vector ${\boldsymbol{\beta}}$ is sparse in the
sense that
$\left\|{\boldsymbol{\beta}}\right\|_{0}:=\sum_{j=1}^{p}{\mathbbm{1}}(\beta_{j}\neq
0)=s$, which is much smaller than $p$. In this paper, we focus on the variable
selection problem, i.e., identifying the active set
$\mathcal{S}:=\\{j:\beta_{j}\neq 0\\}$. We primarily use the 0-1 loss, i.e.,
${\mathbb{P}}(\widehat{\mathcal{S}}\neq\mathcal{S})$, to assess the quality of
the selected model $\widehat{\mathcal{S}}$.

One of the well-studied methods for variable selection in high-dimensional
sparse regression is to penalize the empirical risk by model complexity,
thereby encouraging sparse solutions. Specifically, consider

$\widehat{{\boldsymbol{\beta}}}^{\rm
pen}:={\arg\min}_{{\boldsymbol{\beta}}\in\mathbb{R}^{p}}\,\mathcal{L}({\boldsymbol{\beta}})+{\rm
pen}_{\lambda}({\boldsymbol{\beta}}),$

were $\mathcal{L}({\boldsymbol{\beta}})$ is a loss function and ${\rm
pen}_{\lambda}({\boldsymbol{\beta}})$ is the penalization term that controls
the model complexity. Classical methods such as AIC (Akaike,, 1974, 1998), BIC
(Schwarz,, 1978), Mallow’s $C_{p}$ (Mallows,, 2000) use model complexity as
penalty term, i.e., $\ell_{0}$-norm of the regression coefficient, to penalize
the negative log-likelihood. Although these methods enjoy nice sampling
properties (Barron et al.,, 1999; Zhang and Zhang,, 2012), such $\ell_{0}$
regularized methods are known to suffer from huge computational bottleneck
(Foster et al.,, 2015). This motivated a whole generation of statisticians to
develop alternative penalization methods such as LASSO (Tibshirani,, 1996),
SCAD (Fan and Li,, 2001), MC+ (Zhang,, 2010), and many others that have both
strong statistical guarantee and computational expediency.

However, after recent computational advancements in solving BSS (Bertsimas et
al.,, 2016; Bertsimas and Parys,, 2020; Zhu et al.,, 2020), there has been
growing acknowledgment that BSS enjoys significant statistical superiority
over its computational surrogates and has inevitably motivated statisticians
to investigate the variable selection properties of BSS. For example, through
extensive simulations, Hastie et al., (2020) show that BSS performs better
than LASSO in high signal-to-noise ratio regime in terms of the prediction
risk. Guo et al., (2020) 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). Specifically, they show that BSS is robust to design
collinearity. Under a certain asymptotic regime and independent design, Roy et
al., (2022) further established information-theoretic optimality of BSS in
terms of precise constants for the signal strength parameter under weak and
heterogeneous signal regime.

In this paper, we also study the variable selection property of BSS and
identify novel quantities that are fundamental to understanding the model
consistency of BSS. Specifically, we take the geometric alignment of the
feature vectors $\\{\mathbf{X}_{j}\\}_{j\in[p]}$ into consideration to produce
a more refined analysis of BSS, and show that on top of a certain
identifiability margin (Guo et al.,, 2020), the following two geometric
quantities also control the model selection performance of BSS: (a) Geometric
complexity of the space of _residualized signals_ , and (b) Geometric
complexity of _spurious projections_. We show the explicit dependence of these
two complexity measures in our main results and demonstrate the interplay
between the margin condition and the underlying geometric structure of the
features through some illustrative examples. In the process, we also point out
the existence of a design that is more favorable to BSS than the orthogonal
design, which is commonly believed to be the easiest case for model selection.
To the best of our knowledge, this is the first work that identifies the
underlying geometric complexity of the feature space as a governing force
behind the performance of BSS.

The rest of the paper is organized as follows. In Section 2 we discuss the
preliminaries of BSS. Section 3 is devoted to the discussion of the key
quantities, i.e., identifiability margin and the two complexities. In
particular, Section 3.1-3.3 carefully introduce the notion of identifiability
margin discussed in Guo et al., (2020) and the two novel complexity measures.
In Section 3.4, we build intuition for understanding the effect of these two
complexities with varying correlation. In Section 4, we present both
sufficient (Section 4.1) and necessary (Section 4.3) conditions for model
consistency of BSS. We also partially extend our result to GLM models and
present a similar sufficiency result for model consistency in Section S2 of
the supplementary material.

##### Notation.

Let ${\mathbb{R}}$ denote the set of real numbers. 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
$\mathbb{I}_{p}\in{\mathbb{R}}^{p\times p}$ to denote the $p$-dimensional
identity matrix, and $\mathbf{1}_{p}\in\mathbb{R}^{p}$ denotes the
$p$-dimensional vector with all entries equal to 1. For a matrix
$\mathbf{A}\in{\mathbb{R}}^{p\times p}$, we denote by $\mathbf{A}_{j}$ and
$\mathbf{a}_{j}$ the $j$th column and the transposed $j$th row of $\mathbf{A}$
respectively. We use $\text{col}(\mathbf{A})$ and
$\text{col}(\mathbf{A})^{\perp}$ to denote the columnspace of $\mathbf{A}$ and
its orthogonal complement respectively.

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}\asymp b_{n}$ if $a_{n}=O(b_{n})$ and
$a_{n}=\Omega(b_{n})$. Finally, we write $a_{n}\sim b_{n}$ if
$\lim_{n\to\infty}a_{n}/b_{n}=1$, and $a_{n}=o(b_{n})$ if
$\lim_{n\to\infty}a_{n}/b_{n}=0$.

##  2 Best subset selection

We briefly review the preliminaries of best subset selection (BSS), one of the
most classical variable selection approaches. For a given sparsity level
$\widehat{s}$, BSS solves for

$\hat{{\boldsymbol{\beta}}}_{\rm
best}(\widehat{s}):={\arg\min}_{{\boldsymbol{\beta}}\in\mathbb{R}^{p},\left\|{\boldsymbol{\beta}}\right\|_{0}\leq\widehat{s}}\left\|\mathbf{y}-\mathbf{X}{\boldsymbol{\beta}}\right\|_{2}^{2}.$

For model selection purposes, we can choose the best fitting model to be
$\widehat{\mathcal{S}}_{\rm
best}(\widehat{s}):=\\{j:[\hat{{\boldsymbol{\beta}}}_{\rm
best}(\widehat{s})]_{j}\neq 0\\}$. For a subset $\mathcal{D}\subseteq[p]$,
define the matrix
$\mathbf{X}_{\mathcal{D}}:=(\mathbf{X}_{j};j\in\mathcal{D})$. Let
$\mathbf{P}_{\mathcal{D}}:=\mathbf{X}_{\mathcal{D}}(\mathbf{X}_{\mathcal{D}}^{\top}\mathbf{X}_{\mathcal{D}})^{-1}\mathbf{X}_{\mathcal{D}}^{\top}$
be orthogonal projection operator onto the column space of
$\mathbf{X}_{\mathcal{D}}$. Also, define the corresponding residual sum of
squares (RSS) for model $\mathcal{D}$ as

$R_{\mathcal{D}}:=\mathbf{y}^{\top}(\mathbb{I}_{n}-\mathbf{P}_{\mathcal{D}})\mathbf{y}.$

With this notation, the $\widehat{\mathcal{S}}_{\rm best}(\widehat{s})$ can be
alternatively written as

$\widehat{\mathcal{S}}_{\rm
best}(\widehat{s}):={\arg\min}_{\mathcal{D}\subseteq[p]:\left|\mathcal{D}\right|\leq\widehat{s}}\,R_{\mathcal{D}}.$
(2.1)

Given any candidate model $\mathcal{D}\subset[p]$, we can rewrite the model
(1.1) as

$\mathbf{y}=\mathbf{X}_{\mathcal{S}}{\boldsymbol{\beta}}_{\mathcal{S}}+{\boldsymbol{\varepsilon}}=\mathbf{P}_{\mathcal{D}}\mathbf{X}_{\mathcal{S}}{\boldsymbol{\beta}}_{\mathcal{S}}+(\mathbb{I}_{n}-\mathbf{P}_{\mathcal{D}})\mathbf{X}_{\mathcal{S}\setminus\mathcal{D}}{\boldsymbol{\beta}}_{\mathcal{S}\setminus\mathcal{D}}+{\boldsymbol{\varepsilon}}.$

The term
$(\mathbb{I}_{n}-\mathbf{P}_{\mathcal{D}})\mathbf{X}_{\mathcal{S}\setminus\mathcal{D}}{\boldsymbol{\beta}}_{\mathcal{S}\setminus\mathcal{D}}$
is the residual part of the signal that can not be linearly explained by
$\mathbf{X}_{\mathcal{D}}$. We refer to this part as the residualized signals.
We can thus measure the discrimination between the true model $\mathcal{S}$
and a different candidate model $\mathcal{D}$ through the quantity
$n^{-1}\left\|(\mathbb{I}_{n}-\mathbf{P}_{\mathcal{D}})\mathbf{X}_{\mathcal{S}\setminus\mathcal{D}}{\boldsymbol{\beta}}_{\mathcal{S}\setminus\mathcal{D}}\right\|_{2}^{2}$.

Let $\widehat{{\boldsymbol{\Sigma}}}:=n^{-1}\mathbf{X}^{\top}\mathbf{X}$ be
the sample covariance matrix and for any two sets
$\mathcal{D}_{1},\mathcal{D}_{2}\subset[p]$,
$\widehat{{\boldsymbol{\Sigma}}}_{\mathcal{D}_{1},\mathcal{D}_{2}}$ denotes
the submatrix of ${\boldsymbol{\Sigma}}$ with row indices in $\mathcal{D}_{1}$
and column indices in $\mathcal{D}_{2}$. Next, we define the collection
$\mathscr{A}_{\widehat{s}}:=\\{\mathcal{D}\subset[p]:\mathcal{D}\neq\mathcal{S},\left|\mathcal{D}\right|=\widehat{s}\\}$,
and for $\mathcal{D}\in\mathscr{A}_{\widehat{s}}$ write

$\Gamma(\mathcal{D})=\widehat{{\boldsymbol{\Sigma}}}_{\mathcal{S}\setminus\mathcal{D},\mathcal{S}\setminus\mathcal{D}}-\widehat{{\boldsymbol{\Sigma}}}_{\mathcal{S}\setminus\mathcal{D},\mathcal{D}}\widehat{{\boldsymbol{\Sigma}}}_{\mathcal{D},\mathcal{D}}^{-1}\widehat{{\boldsymbol{\Sigma}}}_{\mathcal{D},\mathcal{S}\setminus\mathcal{D}}.$

Then, it follows that
$n^{-1}\left\|(\mathbb{I}_{n}-\mathbf{P}_{\mathcal{D}})\mathbf{X}_{\mathcal{S}\setminus\mathcal{D}}{\boldsymbol{\beta}}_{\mathcal{S}\setminus\mathcal{D}}\right\|_{2}^{2}={\boldsymbol{\beta}}_{\mathcal{S}\setminus\mathcal{D}}^{\top}\Gamma(\mathcal{D}){\boldsymbol{\beta}}_{\mathcal{S}\setminus\mathcal{D}}$.
Intuitively, if
${\boldsymbol{\beta}}_{\mathcal{S}\setminus\mathcal{D}}^{\top}\Gamma(\mathcal{D}){\boldsymbol{\beta}}_{\mathcal{S}\setminus\mathcal{D}}$
is very close to zero, then there exists
$\mathbf{b}\in\mathbb{R}^{\left|\mathcal{D}\right|}$ such that
$\mathbf{X}_{\mathcal{S}}{\boldsymbol{\beta}}_{\mathcal{S}}\approx\mathbf{X}_{\mathcal{D}}\mathbf{b}$.
Hence, $\mathcal{S}$ and $\mathcal{D}$ have similar linear explanatory power,
and the true model $\mathcal{S}$ becomes practically indistinguishable from
$\mathcal{D}$. In fact, the following lemma shows that
${\boldsymbol{\beta}}_{\mathcal{S}\setminus\mathcal{D}}^{\top}\Gamma(\mathcal{D}){\boldsymbol{\beta}}_{\mathcal{S}\setminus\mathcal{D}}$
needs to be at least bounded away from 0 for all
$\mathcal{D}\in\mathscr{A}_{\widehat{s}}$ in order to make $\mathcal{S}$
identifiable.

######  Lemma 1.

If there exists a $\mathcal{D}\in\mathscr{A}_{\widehat{s}}$ such that
${\boldsymbol{\beta}}_{\mathcal{S}\setminus\mathcal{D}}^{\top}\Gamma(\mathcal{D}){\boldsymbol{\beta}}_{\mathcal{S}\setminus\mathcal{D}}=0$,
then there exists $\mathbf{b}\in\mathbb{R}^{\widehat{s}}$ such that
$\mathbf{X}_{\mathcal{S}}{\boldsymbol{\beta}}_{\mathcal{S}}=\mathbf{X}_{\mathcal{D}}\mathbf{b}$.
Hence, both $\mathbf{X}_{\mathcal{S}}{\boldsymbol{\beta}}_{\mathcal{S}}$ and
$\mathbf{X}_{\mathcal{D}}\mathbf{b}$ generates the same probability
distribution for $\mathbf{y}$, and $\mathcal{S}$ becomes non-identifiable.

Now we are ready to introduce the identifiability margin that characterizes
the _model discriminative power_ of BSS and the two complexity measures.

##  3 Identifiability margin and two complexities

###  3.1 Identifiability margin

The discussion in Section 2 motivates us to define the following
_identifiablity margin_ :

$\tau_{*}(\widehat{s}):=\min_{\mathcal{D}\in\mathscr{A}_{\widehat{s}}}\frac{{\boldsymbol{\beta}}_{\mathcal{S}\setminus\mathcal{D}}^{\top}\Gamma(\mathcal{D}){\boldsymbol{\beta}}_{\mathcal{S}\setminus\mathcal{D}}}{\left|\mathcal{D}\setminus\mathcal{S}\right|}.$
(3.1)

If we define
$\mathscr{A}_{\widehat{s},k}:=\\{\mathcal{D}\in\mathscr{A}_{\widehat{s}}:\left|\mathcal{D}\setminus\mathcal{S}\right|=k\\}$,
then the above can be rewritten as

$\tau_{*}(\widehat{s})=\min_{k\in[\widehat{s}]}\min_{\mathcal{D}\in\mathscr{A}_{\widehat{s},k}}\frac{{\boldsymbol{\beta}}_{\mathcal{S}\setminus\mathcal{D}}^{\top}\Gamma(\mathcal{D}){\boldsymbol{\beta}}_{\mathcal{D}\setminus\mathcal{S}}}{k}.$

As mentioned earlier, the quantity $\tau_{*}(\widehat{s})$ captures the model
discriminative power of BSS. To add more perspective, note that if the
features are highly correlated among themselves then it is expected that
$\tau_{*}(\widehat{s})$ is very close to $0$. Hence, any candidate model
$\mathcal{D}$ is practically indistinguishable from the true model
$\mathcal{S}$ which in turn makes the problem of exact model recovery harder.
On the contrary, if the features are uncorrelated then $\tau_{*}(\widehat{s})$
becomes bounded away from 0 making the true model $\mathcal{S}$ easily
recoverable. For example, Guo et al., (2020) showed that under the condition

$\tau_{*}(s)\gtrsim\sigma^{2}\frac{\log p}{n},$ (3.2)

BSS is able to achieve model consistency. From now on, unless otherwise
mentioned, we will assume that $\tau_{*}(\widehat{s})>0$ to avoid the non-
identifiability issue as pointed out in Lemma 1. .

Next, we will shift focus on the underlying geometric structure of two spaces
that governs the difficulty of the BSS problem (2.1). We essentially identify
the complexities of two types of sets that control the hardness of BSS: (i)
the set of residualized signals, (ii) the set of spurious projections. We
discuss these two sets and the associated complexities in detail below.

###  3.2  Complexity of residualized signals

We start with the definition of the residualized signal. For a candidate model
$\mathcal{D}\in\mathscr{A}_{\widehat{s}}$, define

${\boldsymbol{\gamma}}_{\mathcal{D}}:=n^{-1/2}(\mathbb{I}_{n}-\mathbf{P}_{\mathcal{D}})\mathbf{X}_{\mathcal{S}\setminus\mathcal{D}}{\boldsymbol{\beta}}_{\mathcal{S}\setminus\mathcal{D}},$

and the corresponding unit vector
$\widehat{{\boldsymbol{\gamma}}}_{\mathcal{D}}:={\boldsymbol{\gamma}}_{\mathcal{D}}/\left\|{\boldsymbol{\gamma}}_{\mathcal{D}}\right\|_{2}$.
Note that $\widehat{{\boldsymbol{\gamma}}}_{\mathcal{D}}$ is well-defined as
$\left\|{\boldsymbol{\gamma}}_{\mathcal{D}}\right\|_{2}^{2}\geq\tau_{*}(\widehat{s})>0$.
As mentioned before, ${\boldsymbol{\gamma}}_{\mathcal{D}}$ represents the part
of the signal that can not be linearly explained by the features in model
$\mathcal{D}$. Note that the margin condition (3.2) basically tells that the
vectors ${\boldsymbol{\gamma}}_{\mathcal{D}}$ are well bounded away from the
origin. However, this property does not quite capture the degree of their
radial spread in $\mathbb{R}^{n}$. It may happen that despite being well
bounded away from the origin, the vectors are clustered along one common unit
direction. To capture this notion of separation within the vectors
$\\{{\boldsymbol{\gamma}}_{\mathcal{D}}\\}_{\mathcal{D}\in\mathscr{A}_{\widehat{s}}}$,
we also need to capture the spatial alignment of their corresponding unit
vectors
$\\{\widehat{{\boldsymbol{\gamma}}}_{\mathcal{D}}\\}_{\mathcal{D}\in\mathscr{A}_{\widehat{s}}}$.
This motivates us to consider the geometric complexities of this set of unit
vectors. Specifically, for a set $\mathcal{I}\subset\mathcal{S}$, we define

$\mathcal{T}_{\mathcal{I}}^{(\widehat{s})}:=\\{\widehat{{\boldsymbol{\gamma}}}_{\mathcal{D}}:\mathcal{D}\in\mathscr{A}_{\widehat{s}},\mathcal{S}\cap\mathcal{D}=\mathcal{I}\\}\subseteq\mathbb{R}^{n},$

which is the set of all the normalized forms of the residualized signals
corresponding to the models $\mathcal{D}\in\mathscr{A}_{\widehat{s}}$ with
$\mathcal{I}$ as the common part with true model $\mathcal{S}$. To capture the
complexity of these spaces, we look at the scaled entropy integral

$\mathscr{E}_{\mathcal{T}_{\mathcal{I}}^{(\widehat{s})}}:=\frac{\int_{0}^{\infty}\sqrt{\log\mathcal{N}(\mathcal{T}_{\mathcal{I}}^{(\widehat{s})},\left\|\cdot\right\|_{2},\varepsilon)}\;{\rm
d}\varepsilon}{\sqrt{\log|\mathcal{T}_{\mathcal{I}}^{(\widehat{s})}|}}.$

The numerator in the above display is commonly known as entropy integral which
captures the topological complexity of
$\mathcal{T}_{\mathcal{I}}^{(\widehat{s})}$. In literature, this quantity has
a connection to the well-known Talagrand’s complexity (Talagrand,, 2005),
which often comes up in controlling the expectation of the supremum of
Gaussian processes (Krahmer et al.,, 2014; Lifshits,, 1995; Adler et al.,,
2007). In this paper, we look at the above scaled version of the entropy
integral which allows us to compare the quantity with the diameter and minimum
pairwise distance between the elements of the set
$\mathcal{T}_{\mathcal{I}}^{(\widehat{s})}$. To elaborate on this point,
define the diameter and minimum pairwise distance of
$\mathcal{T}_{\mathcal{I}}^{(\widehat{s})}$ as follows:

${\sf
D}_{\mathcal{T}_{\mathcal{I}}^{(\widehat{s})}}:=\max_{\mathbf{u},\mathbf{v}\in\mathcal{T}_{\mathcal{I}}^{(\widehat{s})}}\left\|\mathbf{u}-\mathbf{v}\right\|_{2},\quad\text{and}\quad{\sf
d}_{\mathcal{T}_{\mathcal{I}}^{(\widehat{s})}}:=\min_{\mathbf{u},\mathbf{v}\in\mathcal{T}_{\mathcal{I}}^{(\widehat{s})}}\left\|\mathbf{u}-\mathbf{v}\right\|_{2}.$

Now notice the following two simple facts:

$\log\mathcal{N}(\mathcal{T}_{\mathcal{I}}^{(\widehat{s})},\left\|\cdot\right\|_{2},{\sf
D}_{\mathcal{T}_{\mathcal{I}}})=0,\quad\log\mathcal{N}(\mathcal{T}_{\mathcal{I}}^{(\widehat{s})},\left\|\cdot\right\|_{2},{\sf
d}_{\mathcal{T}_{\mathcal{I}}})=\log\left|\mathcal{T}_{\mathcal{I}}^{(\widehat{s})}\right|.$

Noting that
$\log\mathcal{N}(\mathcal{T}_{\mathcal{I}}^{(\widehat{s})},\left\|\cdot\right\|_{2},\delta)$
is a decreasing function over $\delta$, we finally get

${\sf
d}_{\mathcal{T}_{\mathcal{I}}^{(\widehat{s})}}\leq\mathscr{E}_{\mathcal{T}_{\mathcal{I}}^{(\widehat{s})}}\leq{\sf
D}_{\mathcal{T}_{\mathcal{I}}^{(\widehat{s})}}.$

This shows that the quantity
$\mathscr{E}_{\mathcal{T}_{\mathcal{I}}^{(\widehat{s})}}$ roughly captures the
average separation of the set $\mathcal{T}_{\mathcal{I}}^{(\widehat{s})}$. In
our subsequent discussion, we will show that
$\mathscr{E}_{\mathcal{T}_{\mathcal{I}}^{(\widehat{s})}}$ heavily influences
the margin condition for exact recovery. This is indeed an important
observation, as the complexity of the set of residualized signals depends
heavily on the association between the features. For example, they can differ
vastly for highly correlated designs compared to almost uncorrelated designs.
Hence, the effect of $\mathscr{E}_{\mathcal{T}_{\mathcal{I}}^{(\widehat{s})}}$
on the exact model recovery also varies significantly across different classes
of distributions and leads to sharper margin conditions for exact model
recovery.

###  3.3 Complexity of spurious projections

In this section, we will introduce the space of projection operators that also
controls the level of difficulty of the true model recovery. Similar to the
previous section, for a fixed set $\mathcal{I}\subset\mathcal{S}$, we consider
the set

$\mathcal{G}_{\mathcal{I}}^{(\widehat{s})}:=\\{\mathbf{P}_{\mathcal{D}}-\mathbf{P}_{\mathcal{I}}:\mathcal{D}\in\mathscr{A}_{\widehat{s}},\mathcal{S}\cap\mathcal{D}=\mathcal{I}\\}\subseteq\mathbb{R}^{n\times
n}.$

It is a well-known fact that every projection operator of the form
$\mathbf{P}_{\mathcal{D}}-\mathbf{P}_{\mathcal{I}}\in\mathcal{G}_{\mathcal{I}}^{(\widehat{s})}$
has a one-to-one correspondence with the subspace
$\text{col}(\mathbf{X}_{\mathcal{D}})\cap\text{col}(\mathbf{X}_{\mathcal{I}})^{\perp}$.
Thus, $\mathcal{G}_{\mathcal{I}}^{(\widehat{s})}$ can be thought of as the
collection of all linear subspaces of the form
$\text{col}(\mathbf{X}_{\mathcal{D}})\cap\text{col}(\mathbf{X}_{\mathcal{I}})^{\perp}$,
which is essentially the set of spurious features that can not be linearly
explained by the set of features in model $\mathcal{I}\subset\mathcal{S}$. To
capture the proper measure of complexity of the set
$\mathcal{G}_{\mathcal{I}}^{(\widehat{s})}$, it is crucial to induce the space
of projection operators with a proper metric. It turns out that Grassmannian
distance is the correct distance to consider in this context. Specifically,
for two linear subspaces $L_{1},L_{2}$ we look at their _maximum sin-theta
distance_ :

${\rm
d}(L_{1},L_{2}):=\left\|{\boldsymbol{\Pi}}_{L_{1}}-{\boldsymbol{\Pi}}_{L_{2}}\right\|_{\rm
op},$

where ${\boldsymbol{\Pi}}_{L_{1}},{\boldsymbol{\Pi}}_{L_{2}}$ are the
orthogonal projection operators of $L_{1},L_{2}$ respectively. It turns out
that ${\rm d}(L_{1},L_{2})$ evaluates the trigonometric sine function at the
maximum principal angle between the subspaces $L_{1}$ and $L_{2}$. We point
the readers to Ye and Lim, (2016) for a more detailed discussion on this
topic.

Under this distance, we define the scaled entropy integral as

$\mathscr{E}_{\mathcal{G}_{\mathcal{I}}^{(\widehat{s})}}:=\frac{\int_{0}^{\infty}\sqrt{\log\mathcal{N}(\mathcal{G}_{\mathcal{I}}^{(\widehat{s})},\left\|\cdot\right\|_{{\rm
op}},\varepsilon)}\;{\rm
d}\varepsilon}{\sqrt{\log|\mathcal{G}_{\mathcal{I}}^{(\widehat{s})}|}}.$

Note that, unlike $\mathscr{E}_{\mathcal{T}_{\mathcal{I}}^{(\widehat{s})}}$,
the complexity measure
$\mathscr{E}_{\mathcal{G}_{\mathcal{I}}^{(\widehat{s})}}$ has no dependence on
${\boldsymbol{\beta}}$ or the residualized signal. Thus,
$\mathscr{E}_{\mathcal{G}_{\mathcal{I}}^{(\widehat{s})}}$ roughly captures the
geometric complexity of only the spurious features. In fact, via a similar
argument as in Section 3.2, it can be shown that ${\sf
d}_{\mathcal{G}_{\mathcal{I}}^{(\widehat{s})}}\leq\mathscr{E}_{\mathcal{G}_{\mathcal{I}}^{(\widehat{s})}}\leq{\sf
D}_{\mathcal{G}_{\mathcal{I}}^{(\widehat{s})}},$ where

${\sf
d}_{\mathcal{G}_{\mathcal{I}}^{(\widehat{s})}}:=\min_{\mathbf{U},\mathbf{V}\in\mathcal{G}_{\mathcal{I}}^{(\widehat{s})}}\left\|\mathbf{U}-\mathbf{V}\right\|_{\rm
op},\quad{\sf
D}_{\mathcal{G}_{\mathcal{I}}^{(\widehat{s})}}:=\max_{\mathbf{U},\mathbf{V}\in\mathcal{G}_{\mathcal{I}}^{(\widehat{s})}}\left\|\mathbf{U}-\mathbf{V}\right\|_{\rm
op}.$

Thus, $\mathscr{E}_{\mathcal{G}_{\mathcal{I}}^{(s)}}$ only captures the
separability in the set of subspaces generated by the spurious features. The
main motivation behind considering such quantity is to capture the influence
of the effective size of the set
$\\{\mathcal{G}_{\mathcal{I}}^{(\widehat{s})}\\}_{\mathcal{I}\subset\mathcal{S}}$
in the analysis of BSS. A naive union bound only uses
$\left|\mathcal{G}_{\mathcal{I}}^{(\widehat{s})}\right|=\binom{p-\widehat{s}}{\widehat{s}-\left|\mathcal{I}\right|}$
as a measure of complexity of the set
$\mathcal{G}_{\mathcal{I}}^{(\widehat{s})}$. This is rather loose, as the
effective complexity of the set is much smaller if
$\mathscr{E}_{\mathcal{G}_{\mathcal{I}}^{(\widehat{s})}}$ is small. Thus,
taking $\mathscr{E}_{\mathcal{G}_{\mathcal{I}}^{(\widehat{s})}}$ into account
unravels a broader picture of the effect incurred by the underlying geometry
of the feature space.

###  3.4 Correlation and complexities

From the discussion on the two complexities, it is quite evident that both of
the complexity measures heavily rely on the alignment of the feature vectors
$\\{\mathbf{X}_{j}:j\in[p]\\}$, which directly depends on the correlation
structure among the features in the model. Below, we discuss how these two
types of complexities may vary with correlation among the features.

#### Correlation and spurious projection operators

We first focus on the set $\mathcal{G}_{\mathcal{I}}^{(\widehat{s})}$, as it
is relatively easy to understand its behavior across different correlation
structures. Recall that for a fixed choice of $\mathcal{I}$, the set
$\mathcal{G}_{\mathcal{I}}^{(\widehat{s})}$ is the collection of all the
projection operators of the form
$\mathbf{P}_{\mathcal{D}}-\mathbf{P}_{\mathcal{I}}$ for all
$\mathcal{D}\in\mathscr{A}_{\widehat{s}}$, which can be thought of as the
collection of different subspaces generated by the spurious features. If the
spurious features are highly correlated then it is expected that these
subspaces are essentially indistinguishable from each other, i.e., the mutual
distance between the projection operators
$\\{\mathbf{P}_{\mathcal{D}}-\mathbf{P}_{\mathcal{I}}\\}_{\mathcal{D}\in\mathscr{A}_{\widehat{s}}}$
is significantly smaller compared to the case when they are weakly correlated.
As an example, let us consider the equi-correlated Gaussian design, i.e., the
row vectors $\\{x_{i}\\}_{i\in[n]}$ of $\mathbf{X}$ in (1.1) follows i.i.d.
mean-zero Gaussian distribution with covariance matrix

${\boldsymbol{\Sigma}}=(1-r)\mathbb{I}_{p}+r\mathbf{1}_{p}\mathbf{1}_{p}^{\top}.$

For the sake of simplicity, we also assume that the true model is a singleton
set. In particular, we consider $\mathcal{S}=\\{1\\}$ and set $\widehat{s}=1$.
Also, note that in this case $\mathscr{A}_{\widehat{s}}=\\{j\in[p]:j\neq 1\\}$
and $\mathcal{I}=\varnothing$. Under this setup, we have
$\mathcal{G}_{\varnothing}^{(1)}=\\{\mathbf{X}_{j}\mathbf{X}_{j}^{\top}/\left\|\mathbf{X}_{j}\right\|_{2}^{2}:j\notin\mathcal{S}\\}$
and $n^{-1}\left\|\mathbf{X}_{j}-\mathbf{X}_{k}\right\|_{2}^{2}\approx
2(1-r),$ for all $j,k\neq 1$. If $r$ is very close to 1 in the above display,
then it follows that the vectors $\\{\mathbf{X}_{j}/\sqrt{n}\\}_{j\neq 1}$ are
extremely clustered towards each other, and as a result, the spurious
projection operators are also very close to each other in operator norm. Due
to this, the complexity measure
$\mathscr{E}_{\mathcal{G}_{\varnothing}^{(1)}}$ becomes extremely small and
the subspaces become almost indistinguishable. In contrast, when the features
are approximately uncorrelated, i.e., $r\approx 0$, the scaled features
$\\{\mathbf{X}_{j}/\sqrt{n}\\}_{j\neq 1}$ are roughly orthogonal. In that case
the

$n^{-1}\left\|\mathbf{X}_{j}-\mathbf{X}_{k}\right\|_{2}^{2}\approx
2,\quad\text{for all $j,k\neq 1$}.$

This suggests that the linear spans generated by each of the set of features
$\\{n^{-1/2}\mathbf{X}_{j}\\}_{j\neq 1}$ are well separated and
$\mathscr{E}_{\mathcal{G}_{\varnothing}^{(1)}}$ is well bounded away from
zero. Thus, it follows that the features are well spread out in
$\mathbb{R}^{n}$. This phenomenon indicates that a higher correlation may aid
the model recovery chance for BSS by reducing the search space over the
features. As we will see in our subsequent discussion in Section 4.2, the
correlation between noise variables can significantly help BSS to identify the
correct model. Specifically, we construct an example where the true variables
are uncorrelated with the noise variables and show that a high correlation
among noise variables helps BSS to identify the correct model. The intuition
is that under the presence of correlation, the diversity of the elements in
$\mathcal{G}_{\mathcal{I}}^{(\widehat{s})}$ gets reduced as
$\mathscr{E}_{\mathcal{G}_{\mathcal{I}}^{(\widehat{s})}}$ becomes small. Thus,
BSS needs to search on a comparatively smaller feature space rather than
searching over all possible
$\binom{p-\widehat{s}}{\widehat{s}-\left|\mathcal{I}\right|}$ models, which in
turn aids the probability of finding the correct model out of the other
candidate ones. Thus, the smaller complexity of
$\mathcal{G}_{\mathcal{I}}^{(\widehat{s})}$ counteracts the adverse effect of
correlation to some degree, and it may improve the model recovery performance
of BSS.

#### Correlation and residualized signals

Now we shift our focus to understanding the behavior of the set of normalized
residualized signals denoted by $\mathcal{T}_{\mathcal{I}}^{(\widehat{s})}$.
Recall that for a fixed $\mathcal{I}$, the set
$\mathcal{T}_{\mathcal{I}}^{(\widehat{s})}$ denotes the collection of all the
unit vectors $\widehat{{\boldsymbol{\gamma}}}_{\mathcal{D}}$ (defined in
Section 3.2) such that $\mathcal{D}\cap\mathcal{S}=\mathcal{I}$. Similar to
$\mathcal{G}_{\mathcal{I}}^{(\widehat{s})}$, the complexity of the set
$\mathcal{T}_{\mathcal{I}}^{(\widehat{s})}$ also depends on the correlation
structure among the features. To elaborate more on this, we revisit the
example of equi-correlated Gaussian design with correlation parameter $r$ and
$\mathcal{S}=\\{1\\}$. We denote by $\mathbf{P}_{j}$ the orthogonal projection
operator onto the span of $\mathbf{X}_{j}$, i.e.,
$\mathbf{P}_{j}=\mathbf{X}_{j}\mathbf{X}_{j}^{\top}/\left\|\mathbf{X}_{j}\right\|_{2}^{2}.$
Similar to the previous section, in this case also the set
$\mathcal{T}_{\varnothing}^{(1)}$ consists of the scaled residualized signals
that take the following form for large $n$ with probability:

$\widehat{{\boldsymbol{\gamma}}}_{j}=\frac{(\mathbb{I}_{n}-\mathbf{P}_{j})\mathbf{X}_{1}}{\left\|(\mathbb{I}_{n}-\mathbf{P}_{j})\mathbf{X}_{1}\right\|_{2}}\approx\frac{\mathbf{X}_{1}-r\mathbf{X}_{j}}{\left\|\mathbf{X}_{1}-r\mathbf{X}_{j}\right\|_{2}},\quad\text{for
all $j\neq 1$.}$

Also, note that

$\widehat{\gamma}_{j}^{\top}\widehat{\gamma}_{k}\approx\frac{1-2r^{2}+r^{3}}{1-r^{2}}=:f(r).$

Since, $f(r)$ is a strictly decreasing function on $[0,1)$, and
$\left\|\widehat{{\boldsymbol{\gamma}}}_{j}-\widehat{{\boldsymbol{\gamma}}}_{k}\right\|_{2}^{2}=2(1-\widehat{{\boldsymbol{\gamma}}}_{j}^{\top}\widehat{{\boldsymbol{\gamma}}}_{k})$,
it follows that ${\sf d}_{\mathcal{T}_{\varnothing}}\geq 1/2$ when $r$ is very
close to 1. On the contrary, when $r\approx 0$, the above display suggests
that ${\sf D}_{\mathcal{T}_{\varnothing}^{(1)}}\approx 0$, i.e., for
uncorrelated design, the complexity
$\mathscr{E}_{\mathcal{T}_{\varnothing}^{(1)}}$ of the set
$\mathcal{T}_{\varnothing}^{(1)}$ is smaller compared to the highly correlated
case which is in sharp contrast with the behavior of
$\mathscr{E}_{\mathcal{G}_{\varnothing}^{(1)}}$.

However, it is worth pointing out that the above property of
$\mathscr{E}_{\mathcal{T}_{\varnothing}^{(1)}}$ is very specific to the above
considered model. There may exist a correlated structure where higher
correlation among noise variables does not increase
$\mathscr{E}_{\mathcal{T}_{\varnothing}^{(1)}}$ (see Section 4.2), and
improves the chance of identifying the correct model via BSS. However,
understanding such a phenomenon for a more general design could be
significantly more challenging.

##  4 Theoretical properties of BSS

###  4.1 Model selection consistency of BSS under known sparsity

This section illustrates the interaction between the identifiability margin
(3.1) and the two complexities that characterize the sufficient condition for
the exact model recovery. From here on, we assume that the true sparsity is
known, i.e., we set $\widehat{s}=s$ in (2.1), and BSS searches the best model
out of all possible models of size $s$. We now introduce a technical
assumption that essentially prevents the noisy features from becoming highly
correlated with the true features:

######  Assumption 1.

The design matrix $\mathbf{X}$ enjoys the following property:

$\min_{\mathcal{I}\subset\mathcal{S}}\mathscr{E}_{\mathcal{G}_{\mathcal{I}}^{(s)}}>\\{\log(ep)\\}^{-1/2}.$

The above assumption ensures that the noisy features are distinguishable
enough from the active features in order for BSS to identify the active
features. To see this, consider the case when the noise variables are highly
correlated with the true features $\\{\mathbf{X}_{j}\\}_{j\in\mathcal{S}}$. In
this case, the projection operator
$\mathbf{P}_{\mathcal{D}}-\mathbf{P}_{\mathcal{I}}$ can be written as
$(\mathbb{I}_{n}-\mathbf{P}_{\mathcal{I}})\mathbf{P}_{\mathcal{D}}$ for all
$\mathcal{D}\in\mathcal{G}_{\mathcal{I}}^{(s)}$, whenever
$\mathcal{I}\neq\varnothing$. As the features in
$\\{\mathbf{X}_{j}:j\in\mathcal{D}\setminus\mathcal{S}\\}$ are highly
correlated with $\mathbf{X}_{\mathcal{I}}$, it follows that
$\left\|\mathbf{P}_{\mathcal{D}}-\mathbf{P}_{\mathcal{I}}\right\|_{{\rm
op}}\approx 0$ and by triangle inequality it follows that
$\left\|\mathbf{P}_{\mathcal{D}}-\mathbf{P}_{\mathcal{D}^{\prime}}\right\|_{{\rm
op}}\approx 0$ for any two candidate models $\mathcal{D}$ and
$\mathcal{D}^{\prime}$ such that
$\mathcal{D}\cap\mathcal{S}=\mathcal{D}^{\prime}\cap\mathcal{S}=\mathcal{I}$.
Thus, Assumption 1 gets rid of such cases by indirectly controlling the
correlation between the active features and noisy features. Secondly, the
assumption also enforces diversity among the noise variables in the following
sense: If the features $\\{\mathbf{X}_{j}:j\notin\mathcal{S}^{c}\\}$ are too
similar to each other, then also
$\mathscr{E}_{\mathcal{G}_{\mathcal{I}}^{(s)}}$ shrinks towards 0. Thus,
Assumption 1 prevents the noise variables from becoming extremely correlated
with each other.

Assumptions with similar spirits are fairly common in the literature on high-
dimensional statistics. For example, the well-known Sparse Riesz Condition
(SRC) (Zhang and Huang,, 2008) assumes that there exist positive numbers
$\kappa_{-},\kappa_{+}$ and $\Psi\geq 1$ such that

$\kappa_{-}\leq\frac{\left\|\mathbf{X}\mathbf{v}\right\|_{2}^{2}}{n}\leq\kappa_{+},\quad\text{for
all
$\mathbf{v}\in\\{\mathbf{u}\in\mathbb{R}^{p}:\left\|\mathbf{u}\right\|_{2}=1,\left\|\mathbf{u}\right\|_{0}\leq\Psi
s\\}$.}$ (4.1)

The above SRC condition controls the maximum and minimum eigenvalues of all
the models of size $s$, which essentially prevents the features to become
extremely correlated with each other. In comparison, Assumption 1 is much
weaker than SRC condition in two aspects. First, unlike the SRC, Assumption 1
imposes conditions only over $(2^{s}-2)$ models, whereas SRC imposes
conditions on $\Omega((p/s)^{\left\lfloor\Psi s\right\rfloor})$ many models.
Second, the lower bound requirement in Assumption 1 is rather weak as the
bound decays with increasing ambient dimension and allows a higher degree of
correlation among the features. In other words, SRC condition (4.1) implies
the condition in Assumption 1, and we formalize this claim in the following
proposition.

######  Proposition 1.

Let the columns of $\mathbf{X}$ be normalized, i.e.,
$\left\|\mathbf{X}_{j}\right\|_{2}=\sqrt{n}$. Also, assume that there exist
positive constants $\kappa_{-},\kappa_{+}$ such that the SRC condition (4.1)
holds with $\Psi=2$. Then the condition in Assumption 1 also holds for large
enough $p$, i.e.,
$\min_{\mathcal{I}\subset\mathcal{S}}\mathscr{E}_{\mathcal{G}_{\mathcal{I}}^{(s)}}\geq\kappa_{-}/\kappa_{+}\gg\\{\log(ep)\\}^{-1/2}$.
Furthermore, the implication in the other direction is not true in general.

Now we are ready to state our main sufficiency result.

######  Theorem 4.1 (Sufficiency).

Under Assumption 1, there exists a positive universal constant $C_{0}$ such
that for any $0\leq\eta<1$, whenever the identifiability margin $\tau_{*}(s)$
satisfies

$\displaystyle\frac{\tau_{*}(s)}{\sigma^{2}}\geq$ (4.2)
$\displaystyle\frac{C_{0}}{(1-\eta)^{2}}\left[\max\left\\{\max_{\mathcal{I}\subset\mathcal{S}}\mathscr{E}_{\mathcal{T}_{\mathcal{I}}^{(s)}}^{2},\max_{\mathcal{I}\subset\mathcal{S}}\mathscr{E}_{\mathcal{G}_{\mathcal{I}}^{(s)}}^{2}\right\\}+\sqrt{\frac{\log(es)\vee\log\log(ep)}{\log(ep)}}\right]\frac{\log(ep)}{n},$

we have

$\left\\{\widehat{\mathcal{S}}:|\widehat{\mathcal{S}}|=s,R_{\widehat{\mathcal{S}}}\leq\min_{\mathcal{D}\in\mathscr{A}_{s}}R_{\mathcal{D}}+n\eta\tau_{*}(s)\right\\}=\\{\mathcal{S}\\},$

with probability at least $1-O(\\{s\vee\log p\\}^{-1})$. In particular,
setting $\eta=0$, we have
$\mathcal{S}={\arg\min}_{\mathcal{D}\in\mathscr{A}_{s}}R_{\mathcal{D}}$ with
high probability.

The proof of the above theorem is present Section S1.3 of the supplementary
material. The above theorem gives a sufficient condition for BSS to achieve
model consistency.

Moreover, note that the margin condition (4.2) involves the identifiability
margin $\tau_{*}(s)$ and the two complexities associated with the sets of
residualized signals and spurious projection operators. This condition reveals
an interesting interplay between the identifiability margin and the two
complexities. To highlight this phenomenon, it is instructive to consider the
case when the true model $\mathcal{S}=\\{1\\}$ and $\mathbf{X}_{1}$ is
orthogonal to the spurious features $\\{\mathbf{X}_{j}\\}_{j\neq 1}$ but the
spurious features may be extremely correlated to each other. As mentioned in
the independent block design example in Section 4.2, in this case, both of the
two complexities are small for higher correlation among the spurious features,
whereas $\tau_{*}(s)$ remains roughly unaffected by the strength of
correlation. Thus, the margin condition (4.2) becomes less stringent with
increasing strength of correlation, and the performance of BSS should improve.
To illustrate this phenomenon, we consider a simulation setup with
$p=2000,n=500$, and $s=1$. We generate $\mathbf{X}$ from independent Gaussian
block design mentioned in Section 4.2 with the cross-correlation $c=0$, and
$r\in[0,1)$ being the correlation within the noise variables. Thus, $r=0$
corresponds to the independent Gaussian design. We set
${\boldsymbol{\beta}}=(0.1,0,\ldots,0)^{\top}\in\mathbb{R}^{p}$, and the
errors $\\{\varepsilon_{i}\\}_{i\in[n]}$ are generated in i.i.d. fashion from
${\sf N}(0,1)$. Finally, the response $\mathbf{y}$ is generated according to
model (1.1). Assuming $s$ is known, we use ABESS (Zhu et al.,, 2020) as a fast
computational surrogate for BSS. The left panel of Figure 1 shows that the
mean model recovery rate of ABESS (across 20 independent runs) increases as
the correlation between the noise variables increases to 1, which validates
the findings in Theorem 4.1. The right panel of Figure 1 also shows that a
similar phenomenon is true even for $s>1$.

Figure 1: Model recovery rate of ABESS under independent block design.

On the other hand, as mentioned in Remark 1, under equicorrelated model with
$\mathcal{S}=\\{1\\}$ and high correlation,
$\mathscr{E}_{\mathcal{T}_{\varnothing}^{(1)}}$ remains strictly bounded away
from 0 and it dominates $\mathscr{E}_{\mathcal{G}_{\varnothing}^{(1)}}$.
However, the identifiability margin $\tau_{*}(s)$ becomes very small due to
the high correlation between the true and noise variables. Hence, the margin
condition (4.2) becomes harder to satisfy with increasing correlation. In the
case of independent design, it turns out
$\mathscr{E}_{\mathcal{G}_{\varnothing}^{(1)}}$ is the dominating complexity
measure. This is not surprising as under independent design, the features are
more spread out in the feature space compared to correlated design, whereas
the residualized signals are more concentrated towards a single unit
direction, making $\mathscr{E}_{\mathcal{T}_{\varnothing}^{(1)}}$ smaller
compared to $\mathscr{E}_{\mathcal{G}_{\varnothing}^{(1)}}$.

The above discussion shows that apart from the quantity $\tau_{*}(s)$, the
complexity of residualized signals and the complexity of spurious projection
operators also play a decisive role in the margin condition of the best subset
selection problem. Specifically, the set with higher complexity characterizes
the margin condition in Theorem 4.1.

We can further represent the condition (4.2) in terms of the diameter of the
sets $\mathcal{T}_{\mathcal{I}}^{(s)}$ and $\mathcal{G}_{\mathcal{I}}^{(s)}$.
To see this, recall that
$\mathscr{E}_{\mathcal{T}_{\mathcal{I}}^{(s)}}\leq{\sf
D}_{\mathcal{T}_{\mathcal{I}}^{(s)}}$ and
$\mathscr{E}_{\mathcal{G}_{\mathcal{I}}^{(s)}}\leq{\sf
D}_{\mathcal{G}_{\mathcal{I}}^{(s)}}$ for all $\mathcal{I}\subset\mathcal{S}$.
Under the light of this fact, we have the following corollary:

######  Corollary 1.

Let the condition in Assumption 1 hold. Then there exists a positive universal
constant $C_{0}$ such that for any $0\leq\eta<1$, whenever the identifiability
margin $\tau_{*}(s)$ satisfies

$\displaystyle\frac{\tau_{*}(s)}{\sigma^{2}}\geq$ (4.3)
$\displaystyle\frac{C_{0}}{(1-\eta)^{2}}\left[\max\left\\{\max_{\mathcal{I}\subset\mathcal{S}}{\sf
D}_{\mathcal{T}_{\mathcal{I}}^{(s)}}^{2},\max_{\mathcal{I}\subset\mathcal{S}}{\sf
D}_{\mathcal{G}_{\mathcal{I}}^{(s)}}^{2}\right\\}+\sqrt{\frac{\log(es)\vee\log\log(ep)}{\log(ep)}}\right]\frac{\log(ep)}{n},$

we have

$\left\\{\widehat{\mathcal{S}}:|\widehat{\mathcal{S}}|=s,R_{\widehat{\mathcal{S}}}\leq\min_{\mathcal{D}\in\mathscr{A}_{s}}R_{\mathcal{D}}+n\eta\tau_{*}(s)\right\\}=\\{\mathcal{S}\\},$

with probability at least $1-O(\\{s\vee\log p\\}^{-1})$. In particular,
setting $\eta=0$, we have
$\mathcal{S}={\arg\min}_{\mathcal{D}\in\mathscr{A}_{s}}R_{\mathcal{D}}$ with
high probability.

Corollary 1 essentially conveys the same message as Theorem 4.1, only under a
slightly stronger margin condition (4.3). However, in some cases, it could be
comparatively easier to give theoretical guarantees on the diameters ${\sf
D}_{\mathcal{T}_{\mathcal{I}}^{(s)}},{\sf
D}_{\mathcal{G}_{\mathcal{I}}^{(s)}}$ rather than their corresponding
complexity measures
$\mathscr{E}_{\mathcal{T}_{\mathcal{I}}^{(s)}},\mathscr{E}_{\mathcal{G}_{\mathcal{I}}^{(s)}}$
respectively. Now in the next section, we will discuss a few illustrative
examples to further elaborate on the effects of two complexities.

###  4.2 Illustrative examples

In this section well will discuss a few illustrative examples to highlight the
effect complexities of the two spaces described in Section 3.2 and Section
3.3.

#### Block design with a single active feature

Consider the model (1.1) where the rows of $\mathbf{X}$ are independently
generated from $p$-dimensional multivariate Gaussian distribution with mean-
zero and variance-covariance matrix

${\boldsymbol{\Sigma}}=\begin{pmatrix}1&c\boldsymbol{1}_{p-1}^{\top}\\\
c\boldsymbol{1}_{p-1}&(1-r)\mathbb{I}_{p-1}+r\boldsymbol{1}_{p-1}\boldsymbol{1}_{p-1}^{\top}\end{pmatrix},$

where $c\in[0,0.997],r\in[0,1)$. We need to further impose the restriction

$c^{2}<r+\frac{1-r}{p-1}$

to ensure positive definiteness of ${\boldsymbol{\Sigma}}$. In this case, we
also set the true model $\mathcal{S}=\\{1\\}$ and the noise variance
$\sigma=1$. Recall that in this case the sets of residualized signals and
spurious projection operators are denoted by $\mathcal{T}_{\varnothing}^{(1)}$
and $\mathcal{G}_{\varnothing}^{(1)}$ respectively. Under this setup, we have
the following lemma:

######  Lemma 2.

Assume that $\log p=o(n)$. Then under the above setup, there exist universal
positive constants $C,L,M$ such that the followings are true with
$\varepsilon_{n,p}=C\\{(\log p)/n\\}^{1/2}$:

  1. (a)

For large enough $n,p$ we have

$\displaystyle{\mathbb{P}}\left[\left\\{1+\varepsilon_{n,p}-\frac{(c-\varepsilon_{n,p})^{2}}{1+\varepsilon_{n,p}}\right\\}\geq\frac{\tau_{*}(1)}{\beta_{1}^{2}}\geq\left\\{1-\varepsilon_{n,p}-\frac{(c+\varepsilon_{n,p})^{2}}{1-\varepsilon_{n,p}}\right\\}\right]$
$\displaystyle=1+o(1/p).$

  2. (b)

For large enough $n,p$ we have

$\displaystyle{\mathbb{P}}\left[\max\left\\{\frac{2c^{2}(1-r)}{1-c^{2}}-L\varepsilon_{n,p},0\right\\}\leq{\sf
d}_{\mathcal{T}_{\varnothing}^{(1)}}^{2}\leq{\sf
D}_{\mathcal{T}_{\varnothing}^{(1)}}^{2}\leq\frac{2c^{2}(1-r)}{1-c^{2}}+L\varepsilon_{n,p}\right]$
$\displaystyle=1+o(1/p).$

  3. (c)

For large enough $n,p$ we have

$\displaystyle{\mathbb{P}}\left[\max\left\\{(1-r^{2})-M\varepsilon_{n,p},0\right\\}\leq{\sf
d}_{\mathcal{G}_{\varnothing}^{(1)}}^{2}\leq{\sf
D}_{\mathcal{G}_{\varnothing}^{(}1)}^{2}\leq(1-r^{2})+M\varepsilon_{n,p}\right]$
$\displaystyle=1+o(1/p).$

From part (b) and (c) of the above lemma, it follows that the complexity
$\mathscr{E}_{\mathcal{T}_{\varnothing}^{(1)}}\approx 0$ when $c=0$. For any
fixed $c>0$ and $r\in[0,1)$, we have

$\mathscr{E}_{\mathcal{T}_{\varnothing}^{(1)}}^{2}\sim\frac{2c^{2}(1-r)}{1-c^{2}},\quad\text{and}\quad\mathscr{E}_{\mathcal{G}_{\varnothing}^{(1)}}^{2}\sim(1-r^{2})\quad\text{for
large $n,p$}.$ (4.4)

A detailed derivation of the result is present in Section S1.5 of the
supplementary material. Left panel of Figure 2 shows the partition of
$c\mathchar 45\relax r$ plane based on the dominating complexity. It is
worthwhile to note that a high value $r$, i.e., a high correlation among the
noise variables results into a smaller value of the complexity terms in (4.2).
However, Lemma 2(a) suggest that $\tau_{*}(1)/\beta_{1}^{2}\sim(1-c^{2})$,
i.e., higher correlation between true and noise variables shrinks the margin
quantity $\tau_{*}(1)$ towards 0. This suggests that a smaller value of $c$
and a higher value $r$ is more favorable to BSS than other possible choices of
$(r,c)$. We now discuss these phenomena through some selected examples.

Figure 2: (left) Partition of $c\mathchar 45\relax r$ plane showing dominating
regions for the two complexities. The color gradient indicates the value of
$\mathscr{E}_{\mathcal{T}_{\varnothing}^{(1)}}-\mathscr{E}_{\mathcal{G}_{\varnothing}^{(1)}}$.
(right) The plot of two complexities for varying $r$ under equicorrelated
design.

##### Independent design:

In this case $c=r=0$. In this case (4.4) suggest that
$\mathscr{E}_{\mathcal{G}_{\varnothing}^{(1)}}^{2}\approx 1$. On the other
hand, we see that $\mathscr{E}^{2}_{\mathcal{T}_{\varnothing}^{(1)}}\approx
0$. Thus, the complexity of spurious projections is dominant in this case.
Also, in this case, $\tau_{*}(1)\approx\beta_{1}^{2}$ which suggests that
higher signal strength results in a better performance in terms of model
selection.

##### Independent block design:

In this case, we set $c=0$ and we vary $r$ in $(0,1)$. Note that (4.4) tells
that $\mathscr{E}_{\mathcal{T}_{\varnothing}^{(1)}}^{2}\approx 0$, and
$\mathscr{E}_{\mathcal{G}_{\varnothing}^{(1)}}^{2}$ has a decreasing trend
with $r\in(0,1)$. This suggests that the independent block design with high
value of $r$ is more favorable for BSS to identify the true model compared to
the independent random design model in the previous example. Finally, noting
the fact that $\tau_{*}(1)\approx\beta_{1}^{2}$, we can conclude that for high
values of $r$, the sufficient condition in Theorem 4.1 becomes less stringent.

##### Equicorrelated design:

Here we set $c=r$ and vary $r$ in the interval $[0,1)$. Let $r_{0}$ be the
unique positive solution to the following equation:

$\frac{2r^{2}}{1+r}-(1-r^{2})=0.$

Calculation shows that $r_{0}\approx 0.675$. Using (4.4), it follows that for
$r\in[0,r_{0})$ the complexity of spurious projection operators is dominating,
i.e.,
$\mathscr{E}_{\mathcal{G}_{\varnothing}^{(1)}}^{2}>\mathscr{E}_{\mathcal{T}_{\varnothing}^{(1)}}^{2}.$
In contrast, for $r\in(r_{0},1)$, we have the complexity of the residualized
signals to be dominating, i.e.,
$\mathscr{E}_{\mathcal{T}_{\varnothing}^{(1)}}^{2}>\mathscr{E}_{\mathcal{G}_{\varnothing}^{(1)}}^{2}.$
Right panel of Figure 2 indicates the phase transition between the two
complexities. Since the identifiability margin $\tau_{*}(1)$ roughly behaves
like $\beta_{1}^{2}(1-r^{2})$, the margin quantity becomes very small for high
value of $r$. Hence, for model consistency, we need a high value for
$\beta_{1}^{2}$.

######  Remark 1.

In the example of equi-correlated design with $\mathcal{S}=\\{1\\}$, the
effect of correlation parameter $r$ on the complexities
$\mathscr{E}_{\mathcal{T}_{\varnothing}^{(1)}}$ and
$\mathscr{E}_{\mathcal{G}_{\varnothing}^{(1)}}$ are complementary to each
other. In the case of the set of residualized signals, increasing correlation
among the features increases the overall complexity of the set
$\mathcal{T}_{\varnothing}^{(1)}$ and vice versa. In contrast, higher
correlation decreases the complexity
$\mathscr{E}_{\mathcal{G}_{\varnothing}^{(1)}}$, thus shrinking the effective
size of $\mathcal{G}_{\varnothing}^{(1)}$. Thus, in this case, the two
complexities act as two opposing forces in the margin condition (4.2).

###  4.3 Necessary condition

One question that arises from the preceding discussion is whether the margin
condition in Theorem 4.1 is necessary for model consistency or not.
Specifically, it is natural to ask whether the complexities of residualized
signals and spurious projections also characterize the necessary margin
condition. In this section, we show that a condition very similar to (4.2) is
necessary for model consistency of BSS, which is also governed by a similar
margin quantity and complexity measures.

For $j_{0}\in\mathcal{S}$, we define the set
$\mathscr{C}_{j_{0}}:=\\{\mathcal{D}:\mathcal{S}\setminus\mathcal{D}=\\{j_{0}\\},\left|\mathcal{D}\right|=s\\}\subset\mathscr{A}_{s,1}$.
We consider the maximum _leave-one-out_ identifiability margin for
$j_{0}\in\mathcal{S}$ as

$\widehat{\tau}(s):=\max_{j_{0}\in\mathcal{S}}\max_{\mathcal{D}\in\mathscr{C}_{j_{0}}}\frac{{\boldsymbol{\beta}}_{\mathcal{S}\setminus\mathcal{D}}^{\top}\Gamma(\mathcal{D}){\boldsymbol{\beta}}_{\mathcal{S}\setminus\mathcal{D}}}{\left|\mathcal{D}\setminus\mathcal{S}\right|}=\max_{j_{0}\in\mathcal{S}}\max_{\mathcal{D}\in\mathscr{C}_{j_{0}}}\Gamma(\mathcal{D})\beta_{j_{0}}^{2}.$
(4.5)

Consider the set $\mathcal{I}_{0}:=\mathcal{S}\setminus\\{j_{0}\\}$ for a
fixed index $j_{0}\in\mathcal{S}$. We capture the complexity of
$\mathcal{T}_{\mathcal{I}_{0}}^{(s)}$ through the following quantity:

$\mathscr{E}^{*}_{\mathcal{T}_{\mathcal{I}_{0}}^{(s)}}:=\frac{\sup_{\delta>0}\frac{\delta}{2}\sqrt{\log\mathcal{M}(\delta,\\{\widehat{{\boldsymbol{\gamma}}}_{\mathcal{I}_{0}\cup\\{j\\}}\\}_{j\in\mathcal{S}^{c}},\left\|\cdot\right\|_{2})}}{\sqrt{\log|\mathcal{T}_{\mathcal{I}_{0}}^{(s)}|}}.$
(4.6)

The above display immediately shows that
$\mathscr{E}^{*}_{\mathcal{T}_{\mathcal{I}_{0}}}\geq{\sf
d}_{\mathcal{T}_{\mathcal{I}_{0}}^{(s)}}/2.$ Also, from the property of
packing and covering number, it follows that

$\mathcal{M}(\delta,\\{\widehat{{\boldsymbol{\gamma}}}_{\mathcal{I}_{0}\cup\\{j\\}}\\}_{j\in\mathcal{S}^{c}},\left\|\cdot\right\|_{2})\leq\mathcal{N}(\delta/2,\\{\widehat{{\boldsymbol{\gamma}}}_{\mathcal{I}_{0}\cup\\{j\\}}\\}_{j\in\mathcal{S}^{c}},\left\|\cdot\right\|_{2}).$

As
$\mathcal{N}(\delta/2,\\{\widehat{{\boldsymbol{\gamma}}}_{\mathcal{I}_{0}\cup\\{j\\}}\\}_{j\in\mathcal{S}^{c}},\left\|\cdot\right\|_{2})$
is a decreasing function over $\delta\in(0,\infty)$, we have the follwing
inequality:

$\sup_{\delta>0}\frac{\delta}{2}\sqrt{\log\mathcal{N}(\delta/2,\\{\widehat{{\boldsymbol{\gamma}}}_{\mathcal{I}_{0}\cup\\{j\\}}\\}_{j\in\mathcal{S}^{c}},\left\|\cdot\right\|_{2})}\leq\int_{0}^{\infty}\sqrt{\log\mathcal{N}(\varepsilon,\\{\widehat{{\boldsymbol{\gamma}}}_{\mathcal{I}_{0}\cup\\{j\\}}\\}_{j\in\mathcal{S}^{c}},\left\|\cdot\right\|_{2})}\;d\varepsilon.$

The above inequality further shows that
$\mathscr{E}^{*}_{\mathcal{T}_{\mathcal{I}_{0}}^{(s)}}\leq\mathscr{E}_{\mathcal{T}_{\mathcal{I}_{0}}^{(s)}}\leq{\sf
D}_{\mathcal{T}_{\mathcal{I}_{0}}^{(s)}}$. Hence, similar to
$\mathscr{E}_{\mathcal{T}_{\mathcal{I}_{0}}^{(s)}}$, the alternative
complexity measure $\mathscr{E}^{*}_{\mathcal{T}_{\mathcal{I}_{0}}^{(s)}}$
also captures the average separation among the elements in
$\mathcal{T}_{\mathcal{I}_{0}}^{(s)}$.

Next, we focus on the set $\mathcal{G}_{\mathcal{I}_{0}}^{(s)}$ which is the
collection of all the spurious projection operators of the form
$\mathbf{P}_{\mathcal{D}}-\mathbf{P}_{\mathcal{I}_{0}}$ for all
$\mathcal{D}\in\mathscr{C}_{j_{0}}$. If
$\mathcal{D}=\mathcal{I}_{0}\cup\\{j\\}$ for some $j\in\mathcal{S}^{c}$, then
the corresponding spurious projection operator takes the form

$\mathbf{P}_{\mathcal{D}}-\mathbf{P}_{\mathcal{I}_{0}}=\widehat{\mathbf{u}}_{j}\widehat{\mathbf{u}}_{j}^{\top},$
(4.7)

where $\widehat{\mathbf{u}}_{j}$ denotes the unit vector along the
residualized feature vector
$\mathbf{u}_{j}:=(\mathbb{I}_{n}-\mathbf{P}_{\mathcal{I}_{0}})\mathbf{X}_{j}.$
Thus, the above display basically shows that the
$\mathbf{P}_{\mathcal{D}}-\mathbf{P}_{\mathcal{I}_{0}}$ is the orthogonal
projection operator onto the linear span generated by the residualized feature
$\mathbf{u}_{j}$. Similar to (4.6), we define the complexity measure of
$\mathcal{G}_{\mathcal{I}_{0}}^{(s)}$ as

$\mathscr{E}^{*}_{\mathcal{G}_{\mathcal{I}_{0}}^{(s)}}:=\frac{\sup_{\delta>0}\frac{\delta}{2}\sqrt{\log\mathcal{M}(\delta,\mathcal{G}_{\mathcal{I}_{0}}^{(s)},\left\|\cdot\right\|_{\rm
op})}}{\sqrt{\log|\mathcal{G}_{\mathcal{I}_{0}}^{(s)}|}}.$ (4.8)

By a similar argument, it also follows that ${\sf
d}_{\mathcal{G}_{\mathcal{I}_{0}}^{(s)}}/2\leq\mathscr{E}^{*}_{\mathcal{G}_{\mathcal{I}_{0}}^{(s)}}\leq{\sf
D}_{\mathcal{G}_{\mathcal{I}_{0}}^{(s)}}.$ Hence, combining the above
observation with (4.7), it also follows that
$\mathscr{E}^{*}_{\mathcal{G}_{\mathcal{I}_{0}}^{(s)}}$ captures the angular
separation among the elusive features
$\\{\mathbf{u}_{j}\\}_{j\in\mathcal{S}^{c}}$.

Next, we introduce some technical assumptions that are crucial for our
theoretical analysis of the necessity result.

######  Assumption 2.

The complexities of the $\mathcal{G}_{\mathcal{I}_{0}}^{(s)}$ and
$\mathcal{T}_{\mathcal{I}_{0}}^{(s)}$ are not too small, i.e.,

${\mathscr{E}^{*}}_{\mathcal{G}_{\mathcal{I}_{0}}^{(s)}}^{2}>16\\{\log(ep)\\}^{-1},\quad\text{and}\quad{\mathscr{E}^{*}}_{\mathcal{T}_{\mathcal{I}_{0}}^{(s)}}^{2}>16\\{\log(ep)\\}^{-1}$

for all $\mathcal{I}_{0}\subset\mathcal{S}$ and
$\left|\mathcal{I}_{0}\right|=s-1$.

Assumption 2 combined with the observation (4.7) essentially tells that the
set of elusive features $\\{\widehat{\mathbf{u}}_{j}\\}_{j\in\mathcal{S}^{c}}$
and the scaled spurious signals
$\\{\widehat{{\boldsymbol{\gamma}}}_{\mathcal{D}}\\}_{\mathcal{D}\in\mathscr{C}_{j_{0}}}$
are not too identical with each other, as
$\mathscr{E}^{*}_{\mathcal{G}_{\mathcal{I}_{0}}^{(s)}}$ and
$\mathscr{E}^{*}_{\mathcal{T}_{\mathcal{I}_{0}}^{(s)}}$ would be typically
small otherwise. Thus, Assumption 2 induces diversity in
$\mathcal{T}_{\mathcal{I}_{0}}^{(s)}$ and
$\mathcal{G}_{\mathcal{I}_{0}}^{(s)}$.

######  Condition 1.

There exists a constant $\alpha\in(0,1)$ such that
${\mathscr{E}^{*}}_{\mathcal{T}_{\mathcal{I}_{0}}^{(s)}}/\mathscr{E}_{\mathcal{T}_{\mathcal{I}_{0}}^{(s)}}\in(\alpha,1)$.

The condition basically tells that the set
$\mathcal{T}_{\mathcal{I}_{0}}^{(s)}$ has a somewhat regular geometric shape
in the sense that both the lower and upper complexity are of the same order.
This essentially implies that minimal separation and maximal separation of the
set $\mathcal{T}_{\mathcal{I}_{0}}^{(s)}$ are of the same order. .

Now we present our theorem on the necessary condition for model consistency of
BSS.

######  Theorem 4.2 (Necessity).

Assume ${\boldsymbol{\varepsilon}}\sim{\sf N}(0,\sigma^{2}\mathbb{I}_{n})$,
$p>16e^{3}$ and $s<p/2$. Also, let the Assumption 2 hold and write
$\mathcal{J}=\\{\mathcal{I}\subset\mathcal{S}:\left|\mathcal{I}\right|=s-1\\}$.
Then the followings are true:

  1. (a)

If
${\mathscr{E}^{*}}_{\mathcal{G}_{\mathcal{I}_{0}}^{(s)}}\notin({\mathscr{E}^{*}}_{\mathcal{T}_{\mathcal{I}_{0}}^{(s)}},\mathscr{E}_{\mathcal{T}_{\mathcal{I}_{0}}^{(s)}})$
for all $\mathcal{I}_{0}\in\mathcal{J}$, then there exists a universal
constant $C_{1}>0$ such that

$\widehat{\tau}(s)\leq
C_{1}\max\left\\{\max_{\mathcal{I}_{0}\in\mathcal{J}}{\mathscr{E}^{*}}_{\mathcal{T}_{\mathcal{I}_{0}}^{(s)}}^{2},\max_{\mathcal{I}_{0}\in\mathcal{J}}{\mathscr{E}^{*}}_{\mathcal{G}_{\mathcal{I}_{0}}^{(s)}}^{2}\right\\}\frac{\sigma^{2}\log(ep)}{n}$

implies that

${\mathbb{P}}(\widehat{\mathcal{S}}_{\rm
best}(s)\neq\mathcal{S})\geq\frac{1}{10}.$

  2. (b)

If there exists $\mathcal{I}_{\\#}\in\mathcal{J}$ such that
${\mathscr{E}^{*}}_{\mathcal{G}_{\mathcal{I}_{\\#}}^{(s)}}\in({\mathscr{E}^{*}}_{\mathcal{T}_{\mathcal{I}_{\\#}}^{(s)}},\mathscr{E}_{\mathcal{T}_{\mathcal{I}_{\\#}}^{(s)}})$,
then under Condition 1, there exists a constant $C_{\alpha}$ depending on
$\alpha$, such that

$\widehat{\tau}(s)\leq
C_{\alpha}\max\left\\{\max_{\mathcal{I}_{0}\in\mathcal{J}}{\mathscr{E}^{*}}_{\mathcal{T}_{\mathcal{I}_{0}}^{(s)}}^{2},\max_{\mathcal{I}_{0}\in\mathcal{J}}{\mathscr{E}^{*}}_{\mathcal{G}_{\mathcal{I}_{0}}^{(s)}}^{2}\right\\}\frac{\sigma^{2}\log(ep)}{n}$

implies that

${\mathbb{P}}(\widehat{\mathcal{S}}_{\rm
best}(s)\neq\mathcal{S})\geq\frac{1}{10}.$

The detailed proof can be found in Section S1.4 of the supplementary material.
The above theorem essentially says that if the maximum leave-one-out margin
$\widehat{\tau}(s)\lesssim\sigma^{2}(\log p)/n$ then the BSS fails to achieve
model consistency with positive probability. However, the interesting part of
the above theorem is to understand the effect of the term involving complexity
measures. Similar to Theorem 4.1, here also we see that the dominating
complexity characterizes the necessary condition for model consistency.
However, we reiterate a few major differences between the above theorem and
Theorem 4.1. First, Theorem 4.1 needs $\tau_{*}(s)$ to be lower bounded, which
is much stronger than the required condition on $\widehat{\tau}(s)$ in Theorem
4.2. Second, Theorem 4.2 involves the alternative complexity measures
${\mathscr{E}^{*}}_{\mathcal{T}_{\mathcal{I}_{0}}^{(s)}}$ and
${\mathscr{E}^{*}}_{\mathcal{G}_{\mathcal{I}_{0}}^{(s)}}$, which are typically
smaller than the complexity measures used in Theorem 4.1. Third, the resulting
complexity in Theorem 4.1 involves the maximum over all possible subsets of
$\mathcal{S}$, whereas Theorem 4.2 involves the maximum only over the subsets
of $\mathcal{S}$ of size $s-1$. These three facts are the main reasons that
the requirement in Theorem 4.2 is weaker compared to the margin condition
(4.2). Nonetheless, Theorem 4.2 is still interesting as it shows that the two
types of complexities are indeed important quantities to understand the model
selection performance of BSS.

Theorem 4.2 can also be stated in terms of the diameter and minimum
separability of the sets $\mathcal{T}_{\mathcal{I}_{0}}^{(s)}$ and
$\mathcal{G}_{\mathcal{I}_{0}}^{(s)}$. Recall that
${\mathscr{E}^{*}}_{\mathcal{T}_{\mathcal{I}_{0}}^{(s)}}\gtrsim{\sf
d}_{\mathcal{T}_{\mathcal{I}_{0}}^{(s)}}$ and
${\mathscr{E}^{*}}_{\mathcal{G}_{\mathcal{I}_{0}}^{(s)}}\gtrsim{\sf
d}_{\mathcal{G}_{\mathcal{I}_{0}}^{(s)}}$. Hence, it follows that under the
same conditions in Theorem 4.2, the margin condition

$\widehat{\tau}(s)\gtrsim\max\left\\{\max_{\mathcal{I}_{0}\in\mathcal{J}}{\sf
d}_{\mathcal{T}_{\mathcal{I}_{0}}^{(s)}}^{2},\max_{\mathcal{I}_{0}\in\mathcal{J}}{\sf
d}_{\mathcal{G}_{\mathcal{I}_{0}}^{(s)}}^{2}\right\\}\frac{\sigma^{2}\log(ep)}{n}$

is necessary for model consistency of BSS.

##  5 Conclusion

In this paper, we establish the sufficient and (nearly) necessary conditions
for BSS to achieve model consistency in a high-dimensional linear regression
setup. Apart from the identifiability margin, we show that the geometric
complexity of the residualized signals and spurious projections based on the
entropy number and packing numbers also play a crucial role in characterizing
the margin condition for model consistency of BSS. In particular, we establish
that the dominating complexity among the two plays a decisive role in the
margin condition. We also highlight the variation in these complexity measures
under different correlation strengths between the features through some simple
illustrative examples. Moreover, in the supplementary material, we extend the
results in Theorem 4.1 to the high-dimensional sparse generalized linear
models (Section S2). However, it is an open problem to find the analogs of the
two complexities in more general settings, e.g., the low-rank matrix
regression problem or multi-tasking regression problem.

### Supplementary Materials

The supplementary material contains the extension of Theorem 1 to the
generalized linear model and the proofs of the main results.

## References

  * Adamczak, (2015) Adamczak, R. (2015).  A note on the hanson-wright inequality for random vectors with dependencies.  Electron. Commun. Probab., 20:1–13.
  * Adler et al., (2007) Adler, R. J., Taylor, J. E., et al. (2007).  Random fields and geometry, volume 80.  Springer.
  * Akaike, (1974) Akaike, H. (1974).  A new look at the statistical model identification.  IEEE transactions on automatic control, 19(6):716–723.
  * Akaike, (1998) Akaike, H. (1998).  Information theory and an extension of the maximum likelihood principle.  In Selected papers of hirotugu akaike, pages 199–213. Springer.
  * Banerjee et al., (2019) Banerjee, A., Gu, Q., Sivakumar, V., and Wu, S. Z. (2019).  Random quadratic forms with dependence: Applications to restricted isometry and beyond.  Advances in Neural Information Processing Systems, 32.
  * Barron et al., (1999) Barron, A., Birgé, L., and Massart, P. (1999).  Risk bounds for model selection via penalization.  Probability theory and related fields, 113(3):301–413.
  * Bertsimas et al., (2016) Bertsimas, D., King, A., and Mazumder, R. (2016).  Best subset selection via a modern optimization lens.  Ann. Statist., 44(2):813–852.
  * Bertsimas and Parys, (2020) Bertsimas, D. and Parys, B. V. (2020).  Sparse high-dimensional regression: Exact scalable algorithms and phase transitions.  Ann. Statist., 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.  Ann. Statist., 37(4):1705–1732.
  * Boucheron et al., (2013) Boucheron, S., Lugosi, G., and Massart, P. (2013).  Concentration inequalities: A nonasymptotic theory of independence.  Oxford university press.
  * Fan and Li, (2001) Fan, J. and Li, R. (2001).  Variable selection via nonconcave penalized likelihood and its oracle properties.  J. Am. Stat. Assoc., 96(456):1348–1360.
  * Foster et al., (2015) Foster, D., Karloff, H., and Thaler, J. (2015).  Variable selection is hard.  In Conference on Learning Theory, pages 696–709. PMLR.
  * Gordon, (1941) Gordon, R. D. (1941).  Values of mills’ ratio of area to bounding ordinate and of the normal probability integral for large values of the argument.  Ann. Math. Stat., 12(3):364–366.
  * Guo et al., (2020) Guo, Y., Zhu, Z., and Fan, J. (2020).  Best subset selection is robust against design dependence.  arXiv preprint arXiv:2007.01478.
  * 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.  Stat. Sci., 35(4):579–592.
  * Krahmer et al., (2014) Krahmer, F., Mendelson, S., and Rauhut, H. (2014).  Suprema of chaos processes and the restricted isometry property.  Commun. Pure. Appl. Math., 67(11):1877–1904.
  * Laurent and Massart, (2000) Laurent, B. and Massart, P. (2000).  Adaptive estimation of a quadratic functional by model selection.  Ann. Stat., pages 1302–1338.
  * Lifshits, (1995) Lifshits, M. A. (1995).  Gaussian random functions, volume 322.  Springer Science & Business Media.
  * Mallows, (2000) Mallows, C. L. (2000).  Some comments on cp.  Technometrics, 42(1):87–94.
  * Pijyan et al., (2020) Pijyan, A., Zheng, Q., Hong, H. G., and Li, Y. (2020).  Consistent estimation of generalized linear models with high dimensional predictors via stepwise regression.  Entropy, 22(9).
  * Roy et al., (2022) Roy, S., Tewari, A., and Zhu, Z. (2022).  High-dimensional variable selection with heterogeneous signals: A precise asymptotic perspective.  arXiv preprint arXiv:2201.01508.
  * Rudelson and Vershynin, (2013) Rudelson, M. and Vershynin, R. (2013).  Hanson-Wright inequality and sub-gaussian concentration.  Electron. Commun. Probab., 18:1 – 9.
  * Schwarz, (1978) Schwarz, G. (1978).  Estimating the dimension of a model.  Ann. Statist., pages 461–464.
  * Talagrand, (2005) Talagrand, M. (2005).  The generic chaining: upper and lower bounds of stochastic processes.  Springer Science & Business Media.
  * Tibshirani, (1996) Tibshirani, R. (1996).  Regression shrinkage and selection via the lasso.  J. R. Stat. Soc. Ser. B, 58(1):267–288.
  * 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.  Electron. J. Stat., 3:1360–1392.
  * Wainwright, (2019) Wainwright, M. J. (2019).  High-Dimensional Statistics: A Non-Asymptotic Viewpoint.  Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press.
  * Ye and Lim, (2016) Ye, K. and Lim, L.-H. (2016).  Schubert varieties and distances between subspaces of different dimensions.  SIAM Journal on Matrix Analysis and Applications, 37(3):1176–1197.
  * Zhang, (2010) Zhang, C.-H. (2010).  Nearly unbiased variable selection under minimax concave penalty.  Ann. Statist., 38(2):894–942.
  * Zhang and Huang, (2008) Zhang, C.-H. and Huang, J. (2008).  The sparsity and bias of the lasso selection in high-dimensional linear regression.  Ann. Statist., 36(4):1567–1594.
  * Zhang and Zhang, (2012) Zhang, C.-H. and Zhang, T. (2012).  A general theory of concave regularization for high-dimensional sparse estimation problems.  Stat. Sci., 27(4):576–593.
  * 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.
  * Zheng et al., (2020) Zheng, Q., Hong, H. G., and Li, Y. (2020).  Building generalized linear models with ultrahigh dimensional features: A sequentially conditional approach.  Biometrics, 76(1):47–60.
  * Zhu et al., (2020) Zhu, J., Wen, C., Zhu, J., Zhang, H., and Wang, X. (2020).  A polynomial algorithm for best-subset selection problem.  Proc. Natl. Acad. Sci. U.S.A., 117(52):33117–33123.

##  S1 Proof of main results under linear model

###  S1.1 Proof of Lemma 1

First note that
${\boldsymbol{\beta}}^{\top}_{\mathcal{S}\setminus\mathcal{D}}\Gamma(\mathcal{D}){\boldsymbol{\beta}}_{\mathcal{S}\setminus\mathcal{D}}=0\Leftrightarrow(\mathbb{I}_{n}-\mathbf{P}_{\mathcal{D}})\mathbf{X}_{\mathcal{S}\setminus\mathcal{D}}{\boldsymbol{\beta}}_{\mathcal{S}\setminus\mathcal{D}}=0$.
This shows that
$\mathbf{X}_{\mathcal{S}\setminus\mathcal{D}}{\boldsymbol{\beta}}_{\mathcal{S}\setminus\mathcal{D}}\in\text{col}(\mathbf{X}_{\mathcal{D}})$.
Thus, we have
$\mathbf{X}_{\mathcal{S}}{\boldsymbol{\beta}}_{\mathcal{S}}=\mathbf{X}_{\mathcal{S}\setminus\mathcal{D}}{\boldsymbol{\beta}}_{\mathcal{S}\setminus\mathcal{D}}+\mathbf{X}_{\mathcal{S}\cap\mathcal{D}}{\boldsymbol{\beta}}_{\mathcal{S}\cap\mathcal{D}}\in\text{col}(\mathbf{X}_{\mathcal{D}})$.
This finishes the proof.

###  S1.2 Proof of Proposition 1

In this section, we will show that the SRC condition (4.7) is strictly
stronger than the condition in Assumption 1. Recall that the features are
normalized, i.e., $\left\|\mathbf{X}_{j}\right\|_{2}=\sqrt{n}$ for all
$j\in[p]$. Now, we will prove the proposition.

###### Proof.

#### SRC implies Assumption 1:

For a set $\mathcal{I}\subset\mathcal{S}$, define
$\mathscr{A}_{\mathcal{I}}:=\\{\mathcal{D}\in\mathscr{A}_{s}:\mathcal{S}\cap\mathcal{D}=\mathcal{I}\\}$.
Now recall that $\mathscr{E}_{\mathcal{G}_{\mathcal{I}}^{(s)}}\gtrsim{\sf
d}_{\mathcal{G}_{\mathcal{I}}^{(s)}}$ for large $p$. Thus, it suffices to show
that ${\sf d}_{\mathcal{G}_{\mathcal{I}}^{(s)}}$ is large for all choices of
$\mathcal{I}\subset\mathcal{S}$. Let
$\mathcal{D}_{1},\mathcal{D}_{2}\in\mathscr{A}_{\mathcal{I}}$ and write
$\mathcal{M}=\mathcal{D}_{1}\cap\mathcal{D}_{2}$. Let
$m=\left|\mathcal{M}\right|$ and consider the two subspaces
$L_{1}=\text{col}(\mathbf{X}_{\mathcal{D}_{1}})\cap\text{col}(\mathbf{X}_{\mathcal{M}})^{\perp}$
and
$L_{2}=\text{col}(\mathbf{X}_{\mathcal{D}_{2}})\cap\text{col}(\mathbf{X}_{\mathcal{M}})^{\perp}$.
Let $\\{{\boldsymbol{\xi}}_{j}\\}_{j=1}^{m}$ be an orthonormal basis of
$\mathcal{M}$. Let $\\{{\boldsymbol{\alpha}}_{j}\\}_{j=1}^{s-m}$ be an
orthonormal basis of $L_{1}$ and $\\{{\boldsymbol{\delta}}_{j}\\}_{j=1}^{s-m}$
be the orthonormal basis of $L_{2}$ such that

$\theta_{j}:=\angle({\boldsymbol{\alpha}}_{j},{\boldsymbol{\delta}}_{j}),\quad
j\in[k],$

are the principal angles between $L_{1}$ and $L_{2}$ in decreasing order. Now,
we construct the matrix $\mathbf{Z}$ in the following way:

$\mathbf{Z}=\left[\mathbf{X}_{\mathcal{D}_{1}\setminus\mathcal{D}_{2}}\mid\mathbf{X}_{\mathcal{M}}\mid\mathbf{X}_{\mathcal{D}_{2}\setminus\mathcal{D}_{1}}\right].$

There exists matrix $\mathbf{u},\mathbf{v}\in\mathbb{R}^{s-m}$ and
$\mathbf{w}_{1},\mathbf{w}_{2}\in\mathbb{R}^{m}$ such that

${\boldsymbol{\alpha}}_{1}=\mathbf{X}_{\mathcal{D}_{1}\setminus\mathcal{D}_{2}}\mathbf{u}+\mathbf{X}_{\mathcal{M}}\mathbf{w}_{1}\quad\text{and}\quad{\boldsymbol{\delta}}_{1}=\mathbf{X}_{\mathcal{D}_{2}\setminus\mathcal{D}_{1}}\mathbf{v}+\mathbf{X}_{\mathcal{M}}\mathbf{w}_{2}.$

As ${\boldsymbol{\alpha}}_{1}\perp\text{col}(\mathbf{X}_{\mathcal{M}})$, we
have

$1={\boldsymbol{\alpha}}_{1}^{\top}{\boldsymbol{\alpha}}_{1}={\boldsymbol{\alpha}}_{1}^{\top}\mathbf{X}_{\mathcal{D}_{1}\setminus\mathcal{D}_{2}}\mathbf{u}\leq\sqrt{n\kappa_{+}}\left\|\mathbf{u}\right\|_{2}\Rightarrow\left\|\mathbf{u}\right\|_{2}^{2}\geq
1/(n\kappa_{+}).$

By a similar argument, we have $\left\|\mathbf{v}\right\|_{2}^{2}\geq
1/(n\kappa_{+})$. Define the vectors
${\boldsymbol{\eta}}:=(\mathbf{u}^{\top},(\mathbf{w}_{1}-\mathbf{w}_{2})^{\top},\mathbf{v}^{\top})^{\top}$.
Hence,
$\left\|{\boldsymbol{\eta}}\right\|_{2}^{2}\geq\left\|\mathbf{u}\right\|_{2}^{2}+\left\|\mathbf{v}\right\|_{2}^{2}\geq
2/(n\kappa_{+})$. Due to SRC condition (4.7), we have

$\left\|\mathbf{Z}{\boldsymbol{\eta}}\right\|_{2}^{2}\geq
n\left\|{\boldsymbol{\eta}}\right\|_{2}^{2}\kappa_{-}\geq
2(\kappa_{-}/\kappa_{+}).$ (S1.1)

$\displaystyle\left\|\mathbf{Z}{\boldsymbol{\eta}}\right\|_{2}^{2}$
$\displaystyle=\left\|{\boldsymbol{\alpha}}_{1}-{\boldsymbol{\delta}}_{1}\right\|_{2}^{2}$
(S1.2) $\displaystyle=2(1-\sqrt{1-\sin^{2}\theta_{1}})$ $\displaystyle\leq
2\sin\theta_{1},$

where the last inequality follows from the fact that $1-x\leq\sqrt{1-x^{2}}$
for all $x\in[0,1]$. Combining (S1.1) and (S1.2), we have

$\left\|\mathbf{P}_{\mathcal{D}_{1}}-\mathbf{P}_{\mathcal{D}_{2}}\right\|_{\rm
op}\geq\frac{\kappa_{-}}{\kappa_{+}}.$

The above display shows that ${\sf
d}_{\mathcal{G}_{\mathcal{I}}^{(s)}}\gtrsim(\kappa_{-}/\kappa_{+})\gg\\{\log(ep)\\}^{-1/2}$
for all $\mathcal{I}\subset\mathcal{S}$. Hence, the claim follows.

#### Assumption 1 does not imply SRC:

In this case, assume $\mathcal{S}=\\{1\\}$ and $\mathbf{e}_{j}$ be the $j$th
canonical basis in $\mathbb{R}^{p}$. Under this setup, Assumption 1 becomes

$\mathscr{E}_{\mathcal{G}_{\varnothing}^{(1)}}>\\{\log(ep)\\}^{-1/2}.$ (S1.3)

Now assume that

$\frac{2}{\log(ep)}\leq\frac{\left\|\mathbf{X}_{j}-\mathbf{X}_{j^{\prime}}\right\|_{2}^{2}}{n}\leq\frac{3}{\log(ep)},\quad\text{for
all $j,j^{\prime}\in[p]$}.$

Then, for large $p$, the condition in (S1.3) holds but SRC fails with the
choice of $\mathbf{v}=(\mathbf{e}_{j}-\mathbf{e}_{j^{\prime}})/\sqrt{2}$. ∎

###  S1.3 Proof of Theorem 1

Recall that
${\boldsymbol{\mu}}=\mathbf{X}_{\mathcal{S}}{\boldsymbol{\beta}}_{\mathcal{S}}$
,
${\boldsymbol{\gamma}}_{\mathcal{D}}=n^{-1/2}(\mathbb{I}_{n}-\mathbf{P}_{\mathcal{D}}){\boldsymbol{\mu}}$
and

$\Gamma(\mathcal{D})=\widehat{{\boldsymbol{\Sigma}}}_{\mathcal{S}\setminus\mathcal{D},\mathcal{S}\setminus\mathcal{D}}-\widehat{{\boldsymbol{\Sigma}}}_{\mathcal{S}\setminus\mathcal{D},\mathcal{D}}\widehat{{\boldsymbol{\Sigma}}}_{\mathcal{D},\mathcal{D}}^{-1}\widehat{{\boldsymbol{\Sigma}}}_{\mathcal{D},\mathcal{S}\setminus\mathcal{D}}.$
(S1.4)

Note that for $\mathcal{D}\in\mathscr{A}_{s,k}$ and $0\leq\eta<1$ we have the
following:

$\displaystyle
n^{-1}(R_{\mathcal{D}}-R_{\mathcal{S}})=n^{-1}\\{\mathbf{y}^{\top}(\mathbb{I}_{n}-\mathbf{P}_{\mathcal{D}})\mathbf{y}-\mathbf{y}^{\top}(\mathbb{I}_{n}-\mathbf{P}_{\mathcal{S}})\mathbf{y}\\}$
(S1.5)
$\displaystyle=n^{-1}\left\\{(\mathbf{X}_{\mathcal{S}\setminus\mathcal{D}}{\boldsymbol{\beta}}_{\mathcal{S}\setminus\mathcal{D}}+{\boldsymbol{\varepsilon}})^{\top}(\mathbb{I}_{n}-\mathbf{P}_{\mathcal{D}})(\mathbf{X}_{\mathcal{S}\setminus\mathcal{D}}{\boldsymbol{\beta}}_{\mathcal{S}\setminus\mathcal{D}}+{\boldsymbol{\varepsilon}})-{\boldsymbol{\varepsilon}}^{\top}(\mathbb{I}_{n}-\mathbf{P}_{\mathcal{S}}){\boldsymbol{\varepsilon}}\right\\}$
$\displaystyle=\eta{\boldsymbol{\beta}}_{\mathcal{S}\setminus\mathcal{D}}^{\top}\Gamma(\mathcal{D}){\boldsymbol{\beta}}_{\mathcal{S}\setminus\mathcal{D}}++2^{-1}(1-\eta){\boldsymbol{\beta}}_{\mathcal{S}\setminus\mathcal{D}}^{\top}\Gamma(\mathcal{D}){\boldsymbol{\beta}}_{\mathcal{S}\setminus\mathcal{D}}-2\left\\{n^{-1}(\mathbb{I}_{n}-\mathbf{P}_{\mathcal{D}})\mathbf{X}_{\mathcal{S}\setminus\mathcal{D}}{\boldsymbol{\beta}}_{\mathcal{S}\setminus\mathcal{D}}\right\\}^{\top}(-{\boldsymbol{\varepsilon}})$
$\displaystyle\quad+2^{-1}(1-\eta){\boldsymbol{\beta}}_{\mathcal{S}\setminus\mathcal{D}}^{\top}\Gamma(\mathcal{D}){\boldsymbol{\beta}}_{\mathcal{S}\setminus\mathcal{D}}-n^{-1}{\boldsymbol{\varepsilon}}^{\top}(\mathbf{P}_{\mathcal{D}}-\mathbf{P}_{\mathcal{S}}){\boldsymbol{\varepsilon}}.$

Also, let
$\widetilde{{\boldsymbol{\varepsilon}}}:=(-{\boldsymbol{\varepsilon}})$. Now,
$\mathcal{E}$ be an event under which the following happens:

$\left\\{\widehat{\mathcal{S}}:|\widehat{\mathcal{S}}|=s,\min_{\mathcal{S}\in\mathscr{A}_{s}}R_{\widehat{\mathcal{S}}}\leq
R_{\mathcal{S}}+n\eta\tau_{*}(s)\right\\}=\\{\mathcal{S}\\}.$

Define the set
$\mathscr{A}_{\mathcal{I}}:=\\{\mathcal{D}\in\mathscr{A}_{s}:\mathcal{S}\cap\mathcal{D}=\mathcal{I}\\}$.
We also set $\left|\mathcal{I}\right|=s-k$ for $k\in[s]$. Then we have
$\mathscr{A}_{\mathcal{I}}\subset\mathscr{A}_{s,k}$. By union bound we have
the following:

${\mathbb{P}}(\mathcal{E}^{c})\leq\sum_{k=1}^{s}\sum_{\mathcal{I}\subset\mathcal{S}:\left|\mathcal{I}\right|=s-k}{\mathbb{P}}\left\\{\min_{\mathcal{D}\in\mathscr{A}_{\mathcal{I}}}n^{-1}(R_{\mathcal{D}}-R_{\mathcal{S}})<\eta\tau_{*}(s)\right\\}.$
(S1.6)

Thus, under the light of equation (S1.5) it is sufficient to show the
following with high probability:

$\max_{\mathcal{D}\in\mathscr{A}_{\mathcal{I}}}\widehat{{\boldsymbol{\gamma}}}_{\mathcal{D}}^{\top}\widetilde{{\boldsymbol{\varepsilon}}}<\frac{n^{1/2}(1-\eta)}{4}\min_{\mathcal{D}\in\mathscr{A}_{\mathcal{I}}}\left\|{\boldsymbol{\gamma}}_{\mathcal{D}}\right\|_{2},$
(S1.7)

$\max_{\mathcal{D}\in\mathscr{A}_{\mathcal{I}}}n^{-1}\left\\{{\boldsymbol{\varepsilon}}^{\top}(\mathbf{P}_{\mathcal{D}}-\mathbf{P}_{\mathcal{S}}){\boldsymbol{\varepsilon}}\right\\}<\frac{1-\eta}{2}\min_{\mathcal{D}\in\mathscr{A}_{\mathcal{I}}}\left\|{\boldsymbol{\gamma}}_{\mathcal{D}}\right\|_{2}^{2},$
(S1.8)

for every $k\in[s]$. We will analyze the above events separately. We recall
the two important sets below:

$\mathcal{T}_{\mathcal{I}}^{(s)}:=\\{\widehat{{\boldsymbol{\gamma}}}_{\mathcal{D}}:\mathcal{D}\in\mathscr{A}_{s},\mathcal{D}\cap\mathcal{S}=\mathcal{I}\\},\;\text{and}\;\;\mathcal{G}_{\mathcal{I}}^{(s)}:=\\{\mathbf{P}_{\mathcal{D}}-\mathbf{P}_{\mathcal{I}}:\mathcal{D}\in\mathscr{A}_{s},\mathcal{D}\cap\mathcal{S}=\mathcal{I}\\}.$

To reduce notational cluttering, we will drop the $s$ in the superscript, and
use $\mathcal{T}_{\mathcal{I}}$ and $\mathcal{G}_{\mathcal{I}}$ to denote the
above sets.

##### Linear term:

Let
$f_{\mathcal{D}}:=\widehat{{\boldsymbol{\gamma}}}_{\mathcal{D}}^{\top}\widetilde{{\boldsymbol{\varepsilon}}}$
and
$\left\|f\right\|:=\max_{\mathcal{D}\in\mathscr{A}_{\mathcal{I}}}f_{\mathcal{D}}$.
Since ${\sf D}_{\mathcal{T}_{\mathcal{I}}}\leq\sqrt{2}$, Theorem 5.36 of
Wainwright, (2019)tells that there exists a constant $A_{1}>0$ such that

${\mathbb{P}}\left\\{\left\|f\right\|\geq
A_{1}\sigma(\mathscr{E}_{\mathcal{T}_{\mathcal{I}}}\sqrt{k\log(ep)}+u)\right\\}\leq
3\exp\left(-\frac{u^{2}}{2}\right),$ (S1.9)

for all $u>0$. Setting $u=2c_{\mathcal{T}}\sqrt{k\log(ep)}$ in Equation (S1.9)
we get

${\mathbb{P}}(\left\|f\right\|\geq
A_{1}\sigma\mathscr{E}_{\mathcal{T}_{\mathcal{I}}}\sqrt{k\log(ep)}+2A_{1}c_{\mathcal{T}}\sigma\sqrt{k\log(ep)})\leq
3(ep)^{-2c_{\mathcal{T}}^{2}k}.$ (S1.10)

Writing $A_{1}$ as $c_{1}$, we get

${\mathbb{P}}\left\\{\max_{\mathcal{D}\in\mathscr{A}_{\mathcal{I}}}\widehat{{\boldsymbol{\gamma}}}_{\mathcal{D}}^{\top}\widetilde{{\boldsymbol{\varepsilon}}}\geq
c_{1}(\mathscr{E}_{\mathcal{T}_{\mathcal{I}}}+2c_{\mathcal{T}})\sigma\sqrt{k\log(ep)}\right\\}\leq
3(ep)^{-2c_{\mathcal{T}}^{2}k}.$ (S1.11)

##### Quadratic term:

Here we study the quadratic supremum process in Equation (S1.8). First, define
the two projection operators
$\mathbf{P}_{\mathcal{D}\mid\mathcal{I}}=\mathbf{P}_{\mathcal{D}}-\mathbf{P}_{\mathcal{I}}$
and
$\mathbf{P}_{\mathcal{S}\mid\mathcal{I}}:=\mathbf{P}_{\mathcal{S}}-\mathbf{P}_{\mathcal{I}}$.
For any number $c_{\mathcal{G}}\in(\\{\log(ep)\\}^{-1},1)$, by union bound we
have,

$\displaystyle{\mathbb{P}}\left\\{n^{-1}\max_{\mathcal{D}\in\mathscr{A}_{\mathcal{I}}}{\boldsymbol{\varepsilon}}^{\top}(\mathbf{P}_{\mathcal{D}}-\mathbf{P}_{\mathcal{S}}){\boldsymbol{\varepsilon}}>\sigma^{2}u+\sigma^{2}c_{\mathcal{G}}u_{0}\right\\}$
(S1.12)
$\displaystyle{\mathbb{P}}\left\\{n^{-1}\max_{\mathcal{D}\in\mathscr{A}_{\mathcal{I}}}{\boldsymbol{\varepsilon}}^{\top}(\mathbf{P}_{\mathcal{D}\mid\mathcal{I}}-\mathbf{P}_{\mathcal{S}\mid\mathcal{I}}){\boldsymbol{\varepsilon}}>\sigma^{2}u+\sigma^{2}c_{\mathcal{G}}u_{0}\right\\}$
$\displaystyle\leq{\mathbb{P}}\left\\{n^{-1}(k\sigma^{2}-{\boldsymbol{\varepsilon}}^{\top}\mathbf{P}_{\mathcal{S}\mid\mathcal{I}}{\boldsymbol{\varepsilon}})>\sigma^{2}u_{0}c_{\mathcal{G}}\right\\}+{\mathbb{P}}\left\\{n^{-1}\max_{\mathcal{D}\in\mathscr{A}_{\mathcal{I}}}({\boldsymbol{\varepsilon}}^{\top}\mathbf{P}_{\mathcal{D}\mid\mathcal{I}}{\boldsymbol{\varepsilon}}-k\sigma^{2})>\sigma^{2}u\right\\}.$

Also, note that
$\mathbb{E}({\boldsymbol{\varepsilon}}^{\top}\mathbf{P}_{\mathcal{D}\mid\mathcal{I}}{\boldsymbol{\varepsilon}})\leq
k\sigma^{2}$ and recall that
$\mathscr{E}_{\mathcal{G}_{\mathcal{I}}}>\\{\log(ep)\\}^{-1/2}$ for all
$\mathcal{I}\subset\mathcal{S}$. This shows that
$\sqrt{k}\leq\mathscr{E}_{\mathcal{G}_{\mathcal{I}}}\sqrt{k\log(ep)}$.
Furthermore, by the properties of projection matrices, we have $d_{\rm
op}(\mathcal{G}_{\mathcal{I}})=1$ and
$d_{F}(\mathcal{G}_{\mathcal{I}})=\sqrt{k}$ (defined in Section S4). Also, it
follows that the quantities $M,V$ and $U$ (defined in Theorem S4.1) have the
following properties:

$M\leq 2\mathscr{E}_{\mathcal{G}_{\mathcal{I}}}^{2}k\log(ep),\quad V\leq
2\sqrt{k\log(ep)},\quad\text{and}\quad U=1.$

Using these facts and Theorem S4.1, we get that there exists a universal
positive constants $A_{2},A_{3}$, such that for
$t=A_{3}c_{\mathcal{G}}k\log(ep)$, we get

${\mathbb{P}}\left\\{\max_{\mathcal{D}\in\mathscr{A}_{\mathcal{I}}}{\boldsymbol{\varepsilon}}^{\top}\mathbf{P}_{\mathcal{D}\mid\mathcal{I}}{\boldsymbol{\varepsilon}}\geq
A_{2}\sigma^{2}(\mathscr{E}_{\mathcal{G}_{\mathcal{I}}}^{2}+c_{\mathcal{G}})k\log(ep)\right\\}\leq(ep)^{-2c_{\mathcal{G}}^{2}k}.$
(S1.13)

Due to Theorem 1.1 of Rudelson and Vershynin, (2013), setting
$u_{0}=k\log(ep)/(2n)$ we can show that there exists an absolute constant
$A_{4}>0$ such that

$\displaystyle{\mathbb{P}}\left\\{n^{-1}\left|{\boldsymbol{\varepsilon}}^{\top}\mathbf{P}_{\mathcal{S}\mid\mathcal{I}}{\boldsymbol{\varepsilon}}-k\sigma^{2}\right|>\frac{c_{\mathcal{G}}\sigma^{2}k\log(ep)}{2n}\right\\}$
$\displaystyle\leq 2\exp\left\\{-A_{4}c_{\mathcal{G}}k\log(ep)\right\\}$
(S1.14) $\displaystyle=2(ep)^{-A_{4}c_{\mathcal{G}}k},$

Combining Equation (S1.12), (S1.13) and Equation (S1.14) yields

${\mathbb{P}}\left\\{n^{-1}\max_{\mathcal{D}\in\mathscr{A}_{\mathcal{I}}}{\boldsymbol{\varepsilon}}^{\top}(\mathbf{P}_{\mathcal{D}}-\mathbf{P}_{\mathcal{S}}){\boldsymbol{\varepsilon}}>c_{2}(\mathscr{E}_{\mathcal{G}_{\mathcal{I}}}^{2}+c_{\mathcal{G}})\sigma^{2}\frac{k\log(ep)}{n}\right\\}\leq(ep)^{-2c_{\mathcal{G}}^{2}k}+2(ep)^{-A_{4}c_{\mathcal{G}}k},$
(S1.15)

where $c_{2}$ is a universal constant. Now, if we have

$\displaystyle\tau_{*}(s)$
$\displaystyle\triangleq\min_{\mathcal{D}\neq\mathcal{S}}\frac{{\boldsymbol{\beta}}_{\mathcal{S}\setminus\mathcal{D}}^{\top}\Gamma(\mathcal{D}){\boldsymbol{\beta}}_{\mathcal{S}\setminus\mathcal{D}}}{\left|\mathcal{D}\setminus\mathcal{S}\right|}$
(S1.16)
$\displaystyle\geq\frac{64}{(1-\eta)^{2}}\max\left\\{c_{1}\max_{\mathcal{I}\subset\mathcal{S}}(\mathscr{E}_{\mathcal{T}_{\mathcal{I}}}+2c_{\mathcal{T}})^{2},c_{2}\max_{\mathcal{I}\subset\mathcal{S}}(\mathscr{E}_{\mathcal{G}_{\mathcal{I}}}^{2}+c_{\mathcal{G}})\right\\}\frac{\sigma^{2}\log(ep)}{n},$

it will ensure that (S1.7) and (S1.8) hold with high probability. Finally,
using (S1.11) and (S1.15), we have

$\displaystyle{\mathbb{P}}(\mathcal{E}^{c})$
$\displaystyle\leq\sum_{k=1}^{s}\sum_{\mathcal{I}\subset\mathcal{S}:\left|\mathcal{I}\right|=s-k}{\mathbb{P}}\left\\{\min_{\mathcal{D}\in\mathscr{A}_{\mathcal{I}}}n^{-1}(R_{\mathcal{D}}-R_{\mathcal{S}})<\eta\tau_{*}(s)\right\\}$
$\displaystyle\lesssim\sum_{k=1}^{s}\sum_{\mathcal{I}\subset\mathcal{S}:\left|\mathcal{I}\right|=s-k}\left\\{(ep)^{-2c_{\mathcal{T}}^{2}k}+(ep)^{-2c_{\mathcal{G}}^{2}k}+(ep)^{-A_{4}c_{\mathcal{G}}k}\right\\}$
$\displaystyle\lesssim\sum_{k=1}^{s}\binom{s}{k}\left\\{(ep)^{-2c_{\mathcal{T}}^{2}k}+(ep)^{-2c_{\mathcal{G}}^{2}k}+(ep)^{-A_{4}c_{\mathcal{G}}k}\right\\}$
$\displaystyle\lesssim\sum_{k=1}^{s}(es)^{k}\left\\{(ep)^{-2c_{\mathcal{T}}^{2}k}+(ep)^{-2c_{\mathcal{G}}^{2}k}+(ep)^{-A_{4}c_{\mathcal{G}}k}\right\\}$
$\displaystyle\lesssim\sum_{k=1}^{s}\exp\left[-k\\{2c_{\mathcal{T}}^{2}\log(ep)-\log(es)\\}\right]+\exp[-k\\{2c_{\mathcal{G}}^{2}\log(ep)-\log(es)\\}]$
$\displaystyle\quad\quad+\exp\left[-k\\{A_{4}c_{\mathcal{G}}\log(ep)-\log(es)\\}\right].$

Now, setting $c_{\mathcal{T}}=\sqrt{\\{\log(es)\vee\log\log(ep)\\}/\log(ep)}$
and $c_{\mathcal{G}}=(2\vee A_{4}^{-1})c_{\mathcal{T}}$ in the above display,
and using the identity $(a+b)^{2}\leq 2(a^{2}+b^{2})$ we can conclude that the
following is sufficient to hold (S1.16):

$\tau_{*}(s)\geq\frac{64}{(1-\eta)^{2}}\max\\{8c_{1},c_{2}(2\vee
A_{4}^{-1})\\}\left[\max\left\\{\max_{\mathcal{I}\subset\mathcal{S}}\mathscr{E}_{\mathcal{T}_{\mathcal{I}}}^{2},\max_{\mathcal{I}\subset\mathcal{S}}\mathscr{E}_{\mathcal{G}_{\mathcal{I}}}^{2}\right\\}+c_{\mathcal{T}}\right]\frac{\log(ep)}{n}.$

and renaming the absolute constant $64\max\\{8c_{1},c_{2}(2\vee
A_{4}^{-1})\\}$ as $C_{0}$.

###  S1.4 Proof of Theorem 2

###### Proof.

First, we will show that
$\left\|{\boldsymbol{\gamma}}_{\mathcal{D}}\right\|_{2}$ has to be well
bounded away from 0 for every
$\mathcal{D}\in\cup_{j_{0}\in\mathcal{S}}\mathscr{C}_{j_{0}}$. Again, to
reduce notational cluttering, we drop the $s$ in the superscript and use
$\mathcal{T}_{\mathcal{I}_{0}}$ and $\mathcal{G}_{\mathcal{I}_{0}}$ to denote
the sets of scaled residualized signals and spurious projections respectively.

####  Ruling out the case
$\min_{\mathcal{D}\in\cup_{j_{0}\in\mathcal{S}}\mathscr{C}_{j_{0}}}\left\|{\boldsymbol{\gamma}}_{\mathcal{D}}\right\|_{2}\leq\sigma/\sqrt{n}$
:

Let
$\min_{\mathcal{D}\in\cup_{j_{0}\in\mathcal{S}}\mathscr{C}_{j_{0}}}\left\|{\boldsymbol{\gamma}}_{\mathcal{D}}\right\|_{2}\leq\sigma/\sqrt{n}$,
i.e, there exists $\mathcal{D}\in\mathscr{C}_{j_{0}}$ for some
$j_{0}\in\mathcal{S}$ such that

$\left\|{\boldsymbol{\gamma}}_{\mathcal{D}}\right\|_{2}\leq\sigma/\sqrt{n}.$

Recall that $\mathbf{y}\sim{\sf
N}(\mathbf{X}_{\mathcal{S}}{\boldsymbol{\beta}}_{\mathcal{S}}^{*},\sigma^{2}\mathbb{I}_{n})$
and define
$\mathbf{w}:=\mathbf{y}-\mathbf{X}_{\mathcal{S}}{\boldsymbol{\beta}}_{\mathcal{S}}^{*}$.
Hence, we have

$\displaystyle{\mathbb{P}}\left(R_{\mathcal{D}}-R_{\mathcal{S}}<0\right)$
$\displaystyle={\mathbb{P}}\left(\left\|\mathbf{y}-\mathbf{P}_{\mathcal{D}}\mathbf{X}_{\mathcal{S}}{\boldsymbol{\beta}}_{\mathcal{S}}^{*}\right\|_{2}^{2}<\left\|\mathbf{y}-\mathbf{X}_{\mathcal{S}}{\boldsymbol{\beta}}_{\mathcal{S}}^{*}\right\|_{2}^{2}\right)$
$\displaystyle={\mathbb{P}}\left(\left\|\mathbf{w}+n^{1/2}{\boldsymbol{\gamma}}_{\mathcal{D}}\right\|_{2}^{2}<\left\|\mathbf{w}\right\|_{2}^{2}\right)$
$\displaystyle\geq\frac{1}{2}\frac{e^{-0.5}}{\sqrt{2\pi}}>0.1\quad(\text{By
Lemma \ref{lemma: similar Gaussians are hard to distinguish}}),$

i.e.,
${\mathbb{P}}\left(\mathcal{S}\notin{\arg\min}_{\mathcal{D}:\left|\mathcal{D}\right|=s}R_{\mathcal{D}}\right)$
is strictly bounded away from 0. Hence, BSS can not recover the true model.
Hence, we rule out this case.

####  Under the case
$\min_{\mathcal{D}\in\cup_{j_{0}\in\mathcal{S}}\mathscr{C}_{j_{0}}}\left\|{\boldsymbol{\gamma}}_{\mathcal{D}}\right\|_{2}>\sigma/\sqrt{n}$
:

We fix a $j_{0}\in\mathcal{S}$.

#### First decomposition:

$\displaystyle\min_{\mathcal{D}\in\mathscr{C}_{j_{0}}}n^{-1}(R_{\mathcal{D}}-R_{\mathcal{S}})$
$\displaystyle\leq\min_{\mathcal{D}\in\mathscr{C}_{j_{0}}}\\{{\boldsymbol{\beta}}_{\mathcal{S}\setminus\mathcal{D}}^{\top}\Gamma(\mathcal{D}){\boldsymbol{\beta}}_{\mathcal{S}\setminus\mathcal{D}}-2n^{-1/2}{\boldsymbol{\gamma}}_{\mathcal{D}}^{\top}\widetilde{{\boldsymbol{\varepsilon}}}-n^{-1}{\boldsymbol{\varepsilon}}^{\top}(\mathbf{P}_{\mathcal{D}}-\mathbf{P}_{\mathcal{S}}){\boldsymbol{\varepsilon}}\\}$
(S1.17)
$\displaystyle\leq\min_{\mathcal{D}\in\mathscr{C}_{j_{0}}}\left[\left\|{\boldsymbol{\gamma}}_{\mathcal{D}}\right\|_{2}\left\\{\left\|{\boldsymbol{\gamma}}_{\mathcal{D}}\right\|_{2}-2n^{-1/2}\widehat{{\boldsymbol{\gamma}}}_{\mathcal{D}}^{\top}\widetilde{{\boldsymbol{\varepsilon}}}\right\\}-n^{-1}{\boldsymbol{\varepsilon}}^{\top}(\mathbf{P}_{\mathcal{D}}-\mathbf{P}_{\mathcal{I}_{0}}){\boldsymbol{\varepsilon}}\right]$
$\displaystyle\quad+n^{-1}{\boldsymbol{\varepsilon}}^{\top}(\mathbf{P}_{\mathcal{S}}-\mathbf{P}_{\mathcal{I}_{0}}){\boldsymbol{\varepsilon}}$
$\displaystyle\leq\min_{\mathcal{D}\in\mathscr{C}_{j_{0}}}\left[\left\|{\boldsymbol{\gamma}}_{\mathcal{D}}\right\|_{2}\left\\{\left\|{\boldsymbol{\gamma}}_{\mathcal{D}}\right\|_{2}-2n^{-1/2}\widehat{{\boldsymbol{\gamma}}}_{\mathcal{D}}^{\top}\widetilde{{\boldsymbol{\varepsilon}}}\right\\}\right]\quad+n^{-1}{\boldsymbol{\varepsilon}}^{\top}(\mathbf{P}_{\mathcal{S}}-\mathbf{P}_{\mathcal{I}_{0}}){\boldsymbol{\varepsilon}}$

We start with the quadratic term in the right hand side of the above display.
First note that
${\boldsymbol{\varepsilon}}^{\top}(\mathbf{P}_{\mathcal{S}}-\mathbf{P}_{\mathcal{I}_{0}}){\boldsymbol{\varepsilon}}/\sigma^{2}=(\widehat{\mathbf{u}}_{j_{0}}^{\top}{\boldsymbol{\varepsilon}})^{2}/\sigma^{2}$
follows a chi-squared distribution with degrees of freedom 1. Hence, Markov’s
inequality shows that

${\mathbb{P}}\left\\{n^{-1}{\boldsymbol{\varepsilon}}^{\top}(\mathbf{P}_{\mathcal{S}}-\mathbf{P}_{\mathcal{I}_{0}}){\boldsymbol{\varepsilon}}>\frac{2\sigma^{2}}{n}\right\\}\leq\frac{1}{2}.$
(S1.18)

Recall that

$\left\|{\boldsymbol{\gamma}}_{\mathcal{D}}\right\|_{2}\geq\frac{\sigma}{\sqrt{n}},\quad\text{for
all $\mathcal{D}\in\mathscr{C}_{j_{0}}$}.$

Next, by Sudakov’s lower bound, we have

$\mathbb{E}\left(\max_{\mathcal{D}\in\mathscr{C}_{j_{0}}}\widehat{{\boldsymbol{\gamma}}}_{\mathcal{D}}^{\top}\widetilde{{\boldsymbol{\varepsilon}}}\right)\geq\sigma{\mathscr{E}^{*}}_{\mathcal{T}_{\mathcal{I}_{0}}}\sqrt{\log(p-s)}\geq\frac{{\mathscr{E}^{*}}_{\mathcal{T}_{\mathcal{I}_{0}}}}{\sqrt{2}}\sqrt{\log(ep)}.$

The last inequality uses the fact that $s<p/2$ and $p>16e^{3}$. Again, an
application of Borell-TIS inequality yields

${\mathbb{P}}\left\\{\max_{\mathcal{D}\in\mathscr{C}_{j_{0}}}\widehat{{\boldsymbol{\gamma}}}_{\mathcal{D}}^{\top}\widetilde{{\boldsymbol{\varepsilon}}}\geq\sigma{\mathscr{E}^{*}}_{\mathcal{T}_{\mathcal{I}_{0}}}\sqrt{\log(ep)}-c_{\mathcal{T}}\sigma\sqrt{\log(ep)}\right\\}\geq
1-(ep)^{-c_{\mathcal{T}}^{2}/2},$ (S1.19)

for any $c_{\mathcal{T}}>0$. Choosing
$c_{\mathcal{T}}={\mathscr{E}^{*}}_{\mathcal{T}_{\mathcal{I}_{0}}}(2^{-1/2}-2^{-1})$,
and using the fact that
${\mathscr{E}^{*}}_{\mathcal{T}_{\mathcal{I}_{0}}}^{2}\geq
16\\{\log(ep)\\}^{-1}$, we get

${\mathbb{P}}\left\\{\max_{\mathcal{D}\in\mathscr{C}_{j_{0}}}\widehat{{\boldsymbol{\gamma}}}_{\mathcal{D}}^{\top}\widetilde{{\boldsymbol{\varepsilon}}}\geq\sigma{\mathscr{E}^{*}}_{\mathcal{T}_{\mathcal{I}_{0}}}\sqrt{\log(ep)}/2\right\\}\geq
1-e^{-1}.$ (S1.20)

Let us define
$\widehat{\tau}_{j_{0}}=\max_{\mathcal{D}\in\mathscr{C}_{j_{0}}}\Delta(\mathcal{D})\beta_{j_{0}}^{2}.$,
and by construction we have
$\widehat{\tau}(s)=\max_{j_{0}\in\mathcal{S}}\widehat{\tau}_{j_{0}}$. If
$\widehat{\tau}_{j_{0}}^{1/2}\leq\sigma{\mathscr{E}^{*}}_{\mathcal{T}_{\mathcal{I}_{0}}}\sqrt{\log(ep)}/(2n^{1/2})$,
then we have

$\min_{\mathcal{D}\in\mathscr{C}_{j_{0}}}\left[\left\|{\boldsymbol{\gamma}}_{\mathcal{D}}\right\|_{2}\left\\{\left\|{\boldsymbol{\gamma}}_{\mathcal{D}}\right\|_{2}-2n^{-1/2}\widehat{{\boldsymbol{\gamma}}}_{\mathcal{D}}^{\top}\widetilde{{\boldsymbol{\varepsilon}}}\right\\}\right]\leq-\frac{\sigma^{2}{\mathscr{E}^{*}}_{\mathcal{T}_{\mathcal{I}_{0}}}\sqrt{\log(ep)}}{2n}$

Thus, using (S1.18) and the above display we have

${\mathbb{P}}(\mathcal{S}\notin{\arg\min}_{\mathcal{D}:\left|\mathcal{D}\right|=s}R_{\mathcal{D}})={\mathbb{P}}\left(\min_{\mathcal{D}\in\mathscr{C}_{j_{0}}}n^{-1}(R_{\mathcal{D}}-R_{\mathcal{S}})<0\right)\geq
1-\frac{1}{e}-\frac{1}{2}\geq\frac{1}{10},$

as we have
${\mathscr{E}^{*}}_{\mathcal{T}_{\mathcal{I}_{0}}}\sqrt{\log(ep)}>4$
(Assumption 2). Thus the necessary condition turns out to be

$\widehat{\tau}_{j_{0}}\geq\frac{{\mathscr{E}^{*}}_{\mathcal{T}_{\mathcal{I}_{0}}}^{2}}{4}\frac{\sigma^{2}\log(ep)}{n}.$
(S1.21)

#### Second decomposition:

We again start with the difference of RSS between a candidate model
$\mathcal{D}\in\mathscr{A}_{s,k}$ and the true model $\mathcal{S}$:

$\displaystyle n^{-1}(R_{\mathcal{D}}-R_{\mathcal{S}})$ (S1.22)
$\displaystyle=n^{-1}\\{\mathbf{y}^{\top}(\mathbb{I}_{n}-\mathbf{P}_{\mathcal{D}})\mathbf{y}-\mathbf{y}^{\top}(\mathbb{I}_{n}-\mathbf{P}_{\mathcal{S}})\mathbf{y}\\}$
$\displaystyle=n^{-1}\left\\{(\mathbf{X}_{\mathcal{S}\setminus\mathcal{D}}{\boldsymbol{\beta}}_{\mathcal{S}\setminus\mathcal{D}}+{\boldsymbol{\varepsilon}})^{\top}(\mathbb{I}_{n}-\mathbf{P}_{\mathcal{D}})(\mathbf{X}_{\mathcal{S}\setminus\mathcal{D}}{\boldsymbol{\beta}}_{\mathcal{S}\setminus\mathcal{D}}+{\boldsymbol{\varepsilon}})-{\boldsymbol{\varepsilon}}^{\top}(\mathbb{I}_{n}-\mathbf{P}_{\mathcal{S}}){\boldsymbol{\varepsilon}}\right\\}$
$\displaystyle={\boldsymbol{\beta}}_{\mathcal{S}\setminus\mathcal{D}}^{\top}\Gamma(\mathcal{D}){\boldsymbol{\beta}}_{\mathcal{S}\setminus\mathcal{D}}-2\left\\{n^{-1}(\mathbb{I}_{n}-\mathbf{P}_{\mathcal{D}})\mathbf{X}_{\mathcal{S}\setminus\mathcal{D}}{\boldsymbol{\beta}}_{\mathcal{S}\setminus\mathcal{D}}\right\\}^{\top}\widetilde{{\boldsymbol{\varepsilon}}}-n^{-1}{\boldsymbol{\varepsilon}}^{\top}(\mathbf{P}_{\mathcal{D}}-\mathbf{P}_{\mathcal{S}}){\boldsymbol{\varepsilon}}.$

First of all, in order achieve model consistency, the following is necessary
for any $k\in[s]$:

$\min_{\mathcal{D}\in\mathscr{C}_{j_{0}}}n^{-1}(R_{\mathcal{D}}-R_{\mathcal{S}})>0.$
(S1.23)

Recall that
$\widehat{\tau}_{j_{0}}=\max_{\mathcal{D}\in\mathscr{C}_{j_{0}}}\Gamma(\mathcal{D})\beta_{j_{0}}^{2}$.
Next we note that

$\displaystyle\min_{\mathcal{D}\in\mathcal{D}_{j_{0}}}n^{-1}(R_{\mathcal{D}}-R_{\mathcal{S}})$
$\displaystyle\leq\min_{\mathcal{D}\in\mathscr{C}_{j_{0}}}\\{{\boldsymbol{\beta}}_{\mathcal{S}\setminus\mathcal{D}}^{\top}\Gamma(\mathcal{D}){\boldsymbol{\beta}}_{\mathcal{S}\setminus\mathcal{D}}-2n^{-1/2}{\boldsymbol{\gamma}}_{\mathcal{D}}^{\top}\widetilde{{\boldsymbol{\varepsilon}}}-n^{-1}{\boldsymbol{\varepsilon}}^{\top}(\mathbf{P}_{\mathcal{D}}-\mathbf{P}_{\mathcal{S}}){\boldsymbol{\varepsilon}}\\}$
(S1.24)
$\displaystyle\leq\widehat{\tau}_{j_{0}}+2n^{-1/2}\max_{\mathcal{D}\in\mathscr{C}_{j_{0}}}\left|{\boldsymbol{\gamma}}_{\mathcal{D}}^{\top}\widetilde{{\boldsymbol{\varepsilon}}}\right|-n^{-1}\max_{\mathcal{D}\in\mathscr{C}_{j_{0}}}{\boldsymbol{\varepsilon}}^{\top}(\mathbf{P}_{\mathcal{D}}-\mathbf{P}_{\mathcal{S}}){\boldsymbol{\varepsilon}}$
$\displaystyle\leq\widehat{\tau}_{j_{0}}+2(\widehat{\tau}_{j_{0}}/n)^{1/2}\max_{\mathcal{D}\in\mathscr{C}_{j_{0}}}\left|\widehat{{\boldsymbol{\gamma}}}_{\mathcal{D}}^{\top}\widetilde{{\boldsymbol{\varepsilon}}}\right|-n^{-1}\max_{\mathcal{D}\in\mathscr{C}_{j_{0}}}{\boldsymbol{\varepsilon}}^{\top}(\mathbf{P}_{\mathcal{D}}-\mathbf{P}_{\mathcal{S}}){\boldsymbol{\varepsilon}}.$

Similar to the proof of Theorem 1, we define
$f_{\mathcal{D}}:=\widehat{{\boldsymbol{\gamma}}}_{\mathcal{D}}^{\top}\widetilde{{\boldsymbol{\varepsilon}}}$
and
$\left\|f\right\|:=\max_{\mathcal{D}\in\mathscr{C}_{j_{0}}}f_{\mathcal{D}}$.
Hence, we have

$\max_{\mathcal{D}\in\mathscr{C}_{j_{0}}}\left|\widehat{{\boldsymbol{\gamma}}}_{\mathcal{D}}^{\top}\widetilde{{\boldsymbol{\varepsilon}}}\right|=\max_{\mathcal{D}\in\mathscr{C}_{j_{0}}}f_{\mathcal{D}}\vee(-f_{\mathcal{D}})$
(S1.25)

By Borell-TIS inequality (Adler et al.,, 2007, Theorem 2.1.1), we have

${\mathbb{P}}\left\\{\left\|f\right\|-\mathbb{E}(\left\|f\right\|)\geq\sigma
u\right\\}\leq\exp\left(-\frac{u^{2}}{2}\right),$

for all $u>0$. Setting $u=c_{\mathcal{T}}\sqrt{\log(ep)}$ we get

${\mathbb{P}}\left\\{\left\|f\right\|-\mathbb{E}(\left\|f\right\|)\geq
c_{\mathcal{T}}\sigma\sqrt{\log(ep)}\right\\}\leq(ep)^{-c_{\mathcal{T}}^{2}/2}.$

$\mathbb{E}(\left\|f\right\|)\leq
4\sqrt{2}\mathscr{E}_{\mathcal{T}_{\mathcal{I}_{0}}}\sigma\sqrt{\log(ep)},$

which ultimately yields

${\mathbb{P}}\left\\{\left\|f\right\|\geq(4\sqrt{2}\mathscr{E}_{\mathcal{T}_{\mathcal{I}_{0}}}+c_{\mathcal{T}})\sigma\sqrt{\log(ep)}\right\\}\leq(ep)^{-c_{\mathcal{T}}^{2}/2}.$

Finally. using (S1.25) we have the following for any $c_{\mathcal{T}}>0$:

${\mathbb{P}}\left\\{\max_{\mathcal{D}\in\mathscr{C}_{j_{0}}}\left|\widehat{{\boldsymbol{\gamma}}}_{\mathcal{D}}^{\top}\widetilde{E}\right|\geq(4\sqrt{2}\mathscr{E}_{\mathcal{T}_{\mathcal{I}_{0}}}+c_{\mathcal{T}})\sigma\sqrt{\log(ep)}\right\\}\leq
2(ep)^{-c_{\mathcal{T}}^{2}/2}.$ (S1.26)

Next, we will lower bound the quadratic term in Equation (S1.24) with high
probability. similar to the proof of Theorem 1, we consider the decomposition

$\max_{\mathcal{D}\in\mathscr{C}_{j_{0}}}n^{-1}{\boldsymbol{\varepsilon}}^{\top}(\mathbf{P}_{\mathcal{D}}-\mathbf{P}_{\mathcal{S}}){\boldsymbol{\varepsilon}}=\max_{j\notin\mathcal{S}}n^{-1}{\boldsymbol{\varepsilon}}^{\top}(\widehat{\mathbf{u}}_{j}\widehat{\mathbf{u}}_{j}^{\top}-\widehat{\mathbf{u}}_{j_{0}}\widehat{\mathbf{u}}_{j_{0}}^{\top}){\boldsymbol{\varepsilon}}.$

For the maximal process we will use Theorem 2.10 of Adamczak, (2015). We begin
with the definition of concentration property.

######  Definition S1.1 (Adamczak, (2015)).

Let $Z$ be random vector in $\mathbb{R}^{n}$. We say that $Z$ has
concentration property with constant $K$ if for every 1-Lipschitz function
$\varphi:\mathbb{R}^{n}\to\mathbb{R}$, we have
$\mathbb{E}\left|\varphi(Z)\right|<\infty$ and for every $u>0$,

${\mathbb{P}}\left(\left|\varphi(Z)-\mathbb{E}(\varphi(Z))\right|\geq
u\right)\leq 2\exp(-u^{2}/K^{2}).$

Note that the Gaussian vector ${\boldsymbol{\varepsilon}}/\sigma$ enjoys
concentration property with $K=\sqrt{2}$ (Boucheron et al.,, 2013, Theorem
5.6). Let
$Q_{1}:=\max_{j\notin\mathcal{S}}{\boldsymbol{\varepsilon}}^{\top}(\widehat{\mathbf{u}}_{j}\widehat{\mathbf{u}}_{j}^{\top}-\widehat{\mathbf{u}}_{j_{0}}\widehat{\mathbf{u}}_{j_{0}}^{\top}){\boldsymbol{\varepsilon}}$.
By Theorem 2.10 of Adamczak, (2015) we conclude that

${\mathbb{P}}\left\\{n^{-1}\left|Q_{1}-\mathbb{E}(Q_{1})\right|\geq
t\sigma^{2}\right\\}\leq
2\exp\left\\{-\frac{1}{2}\min\left(\frac{n^{2}t^{2}}{16},\frac{nt}{2}\right)\right\\}.$

Setting $t=2\delta\log(ep)/n$ in the above equation we get

${\mathbb{P}}\left\\{n^{-1}\left|Q_{1}-\mathbb{E}(Q_{1})\right|\geq
2\delta\sigma^{2}\log(ep)/n\right\\}\leq 2(ep)^{-\frac{\delta}{2}}.$ (S1.27)

Next, we will lower bound the expected value of $Q_{1}$. First, note the
following:

$\displaystyle\mathbb{E}(Q_{1})$
$\displaystyle=\mathbb{E}\left\\{\max_{j\notin\mathcal{S}}{\boldsymbol{\varepsilon}}^{\top}(\widehat{\mathbf{u}}_{j}\widehat{\mathbf{u}}_{j}^{\top}-\widehat{\mathbf{u}}_{j_{0}}\widehat{\mathbf{u}}_{j_{0}}^{\top}){\boldsymbol{\varepsilon}}\right\\}$
(S1.28)
$\displaystyle=\mathbb{E}\left\\{\max_{j\notin\mathcal{S}}(\widehat{\mathbf{u}}_{j}^{\top}{\boldsymbol{\varepsilon}})^{2}\right\\}-\sigma^{2}$
$\displaystyle\geq\left\\{\mathbb{E}\max_{j\in\mathcal{S}^{c}}\;(\widehat{\mathbf{u}}_{j}^{\top}{\boldsymbol{\varepsilon}})\vee(-\widehat{\mathbf{u}}_{j}^{\top}{\boldsymbol{\varepsilon}})\right\\}^{2}-\sigma^{2}$

Define the set

$\mathcal{U}_{\rm
sym}:=\\{\widehat{\mathbf{u}}_{j}\;:j\in\mathcal{S}^{c}\\}\cup\\{-\widehat{\mathbf{u}}_{j}\;:j\in\mathcal{S}^{c}\\}.$

We denote by $\widetilde{\mathbf{u}}_{j}$ a generic element of
$\mathcal{U}_{\rm sym}$. Thus for any two elements
$\widetilde{\mathbf{u}}_{j},\widetilde{\mathbf{u}}_{k}$, we have

$\left\|\widetilde{\mathbf{u}}_{j}-\widetilde{\mathbf{u}}_{k}\right\|_{2}\geq\min\left\\{\left\|\widehat{\mathbf{u}}_{j}-\widehat{\mathbf{u}}_{k}\right\|_{2},\left\|\widehat{\mathbf{u}}_{j}+\widehat{\mathbf{u}}_{k}\right\|_{2}\right\\}\geq\left\|\widehat{\mathbf{u}}_{j}\widehat{\mathbf{u}}_{j}^{\top}-\widehat{\mathbf{u}}_{k}\widehat{\mathbf{u}}_{k}^{\top}\right\|_{\rm
op}.$

By Sudakov’s lower bound, we have

$\displaystyle\mathbb{E}\max_{j\in\mathcal{S}^{c}}\;(\widehat{\mathbf{u}}_{j}^{\top}{\boldsymbol{\varepsilon}})\vee(-\widehat{\mathbf{u}}_{j}^{\top}{\boldsymbol{\varepsilon}})$
$\displaystyle\geq\sup_{\delta>0}\frac{\sigma\delta}{2}\sqrt{\log\mathcal{M}(\delta,\mathcal{U}_{\rm
sym},\left\|\cdot\right\|_{2})}$
$\displaystyle\geq\sup_{\delta>0}\frac{\sigma\delta}{2}\sqrt{\log\mathcal{M}(\delta,\\{\widehat{\mathbf{u}}_{j}\widehat{\mathbf{u}}_{j}^{\top}\\}_{j\in\mathcal{S}^{c}},\left\|\cdot\right\|_{\rm
op})}$
$\displaystyle\geq\sigma\frac{{\mathscr{E}^{*}}_{\mathcal{G}_{\mathcal{I}_{0}}}}{\sqrt{4/3}}\sqrt{\log(ep)}.$

The last inequality uses the fact that $p>16e^{3}$. Finally, (S1.28) yields

$\mathbb{E}(Q_{1})\geq\frac{\sigma^{2}{\mathscr{E}^{*}}_{\mathcal{G}_{\mathcal{I}_{0}}}^{2}\log(ep)}{4/3}-\sigma^{2}.$

Thus, combined with (S1.27) we finally get

${\mathbb{P}}\left[Q_{1}\geq(3/4)\sigma^{2}{\mathscr{E}^{*}}_{\mathcal{G}_{\mathcal{I}_{0}}}^{2}\log(ep)-\sigma^{2}-2\delta\sigma^{2}\log(ep)\right]\geq
1-2(ep)^{-\delta/2}.$

Setting
$\delta=\frac{{\mathscr{E}^{*}}_{\mathcal{G}_{\mathcal{I}_{0}}}^{2}}{8}$, we
get

${\mathbb{P}}\left[Q_{1}\geq\frac{\sigma^{2}{\mathscr{E}^{*}}_{\mathcal{G}_{\mathcal{I}_{0}}}^{2}}{2}\log(ep)-\sigma^{2}\right]\geq
1-2(ep)^{-\frac{{\mathscr{E}^{*}}_{\mathcal{G}_{\mathcal{I}_{0}}}^{2}}{16}}.$

Thus, finally combining the above with (S1.24) and (S1.26) we get

$\displaystyle\min_{\mathcal{D}\in\mathscr{C}_{j_{0}}}n^{-1}(R_{\mathcal{D}}-R_{\mathcal{S}})$
(S1.29)
$\displaystyle\leq\widehat{\tau}_{j_{0}}+2\widehat{\tau}_{j_{0}}^{1/2}(4\sqrt{2}\mathscr{E}_{\mathcal{T}_{\mathcal{I}_{0}}}+c_{\mathcal{T}})\sigma\sqrt{\frac{\log(ep)}{n}}-n^{-1}\left(\frac{\sigma^{2}{\mathscr{E}^{*}}_{\mathcal{G}_{\mathcal{I}_{0}}}^{2}}{2}\log(ep)-\sigma^{2}\right)$

with probability at least
$1-2(ep)^{-c_{\mathcal{T}}^{2}/2}-2(ep)^{-\frac{{\mathscr{E}^{*}}_{\mathcal{G}_{\mathcal{I}_{0}}}^{2}}{16}}$.
Thus, for large $p$ we have

$\displaystyle\min_{\mathcal{D}\in\mathscr{C}_{j_{0}}}n^{-1}(R_{\mathcal{D}}-R_{\mathcal{D}})$
$\displaystyle\leq\widehat{\tau}_{j_{0}}+2\widehat{\tau}_{j_{0}}^{1/2}(4\sqrt{2}\mathscr{E}_{\mathcal{T}_{\mathcal{I}_{0}}}+c_{\mathcal{T}})\sigma\sqrt{\frac{\log(ep)}{n}}-\frac{\sigma^{2}{\mathscr{E}^{*}}_{\mathcal{G}_{\mathcal{I}_{0}}}^{2}}{4}\frac{\log(ep)}{n}$

with probability at least
$1-2(ep)^{c_{\mathcal{T}}^{2}/2}-2(ep)^{-\frac{{\mathscr{E}^{*}}_{\mathcal{G}_{\mathcal{I}_{0}}}^{2}}{16}}$.
Now choose $c_{\mathcal{T}}=4/\sqrt{\log(ep)}$ and use Assumption 2 to get

$\displaystyle\min_{\mathcal{D}\in\mathscr{C}_{j_{0}}}n^{-1}(R_{\mathcal{D}}-R_{\mathcal{D}})$
$\displaystyle\leq\widehat{\tau}_{j_{0}}+2\widehat{\tau}_{j_{0}}^{1/2}\left(4\sqrt{2}\mathscr{E}_{\mathcal{T}_{\mathcal{I}_{0}}}+\frac{4}{\sqrt{\log(ep)}}\right)\sigma\sqrt{\frac{\log(ep)}{n}}-\frac{\sigma^{2}{\mathscr{E}^{*}}_{\mathcal{G}_{\mathcal{I}_{0}}}^{2}}{4}\frac{\log(ep)}{n}$
$\displaystyle\leq\widehat{\tau}_{j_{0}}+8\sqrt{2}\widehat{\tau}_{j_{0}}^{1/2}\left(\mathscr{E}_{\mathcal{T}_{\mathcal{I}_{0}}}+\frac{1}{\sqrt{2\log(ep)}}\right)\sigma\sqrt{\frac{\log(ep)}{n}}-\frac{\sigma^{2}{\mathscr{E}^{*}}_{\mathcal{G}_{\mathcal{I}_{0}}}^{2}}{4}\frac{\log(ep)}{n}$
$\displaystyle\leq\widehat{\tau}_{j_{0}}+10\sqrt{2}\widehat{\tau}_{j_{0}}^{1/2}\mathscr{E}_{\mathcal{T}_{\mathcal{I}_{0}}}\sigma\sqrt{\frac{\log(ep)}{n}}-\frac{\sigma^{2}{\mathscr{E}^{*}}_{\mathcal{G}_{\mathcal{I}_{0}}}^{2}}{4}\frac{\log(ep)}{n}$

with probability at least $1/5$.

Thus, in light of (S1.23), the following is necessary:

$\widehat{\tau}_{j_{0}}\geq\left\\{\frac{\sqrt{200\mathscr{E}_{\mathcal{T}_{\mathcal{I}_{0}}}^{2}+{\mathscr{E}^{*}}_{\mathcal{G}_{\mathcal{I}_{0}}}^{2}}-10\sqrt{2}\mathscr{E}_{\mathcal{T}_{\mathcal{I}_{0}}}}{2}\right\\}^{2}\frac{\sigma^{2}\log(ep)}{n}.$
(S1.30)

##### Case 1:

If
$\mathscr{E}_{\mathcal{T}_{\mathcal{I}_{0}}}\leq{\mathscr{E}^{*}}_{\mathcal{G}_{\mathcal{I}_{0}}}$,
then the right hand side of (S1.30) is lower bounded by

$\frac{{\mathscr{E}^{*}}_{\mathcal{G}_{\mathcal{I}_{0}}}^{2}}{(\sqrt{201}+10\sqrt{2})^{2}}\frac{\sigma^{2}\log(ep)}{n}.$

Thus (S1.30) yields the necessary condition

$\widehat{\tau}_{j_{0}}\geq\frac{{\mathscr{E}^{*}}_{\mathcal{G}_{\mathcal{I}_{0}}}^{2}}{(\sqrt{201}+10\sqrt{2})^{2}}\frac{\sigma^{2}\log(ep)}{n}.$

Combining this with (S1.21) we have the necessary condition to be

$\widehat{\tau}_{j_{0}}\geq\tilde{C}_{1}\max\\{{\mathscr{E}^{*}}_{\mathcal{T}_{\mathcal{I}_{0}}}^{2},{\mathscr{E}^{*}}_{\mathcal{G}_{\mathcal{I}_{0}}}^{2}\\}\frac{\sigma^{2}\log(ep)}{n},$

for a universal constant $\tilde{C}_{1}$.

##### Case 2:

If
${\mathscr{E}^{*}}_{\mathcal{G}_{\mathcal{I}_{0}}}\leq{\mathscr{E}^{*}}_{\mathcal{T}_{\mathcal{I}_{0}}}$,
then using the inequality $\sqrt{1+t}-\sqrt{t}<1$ for all $t>0$, we can
conclude that the right hand side of (S1.30) is always smaller than

$\frac{{\mathscr{E}^{*}}_{\mathcal{G}_{\mathcal{I}_{0}}}^{2}}{4}\frac{\sigma^{2}\log(ep)}{n},$

which is further smaller than

$\frac{{\mathscr{E}^{*}}_{\mathcal{T}_{\mathcal{I}_{0}}}^{2}}{4}\frac{\sigma^{2}\log(ep)}{n}.$

Thus, combining this with (S1.21) we have the necessary condition to be

$\widehat{\tau}_{j_{0}}\geq\frac{{\mathscr{E}^{*}}_{\mathcal{T}_{\mathcal{I}_{0}}}^{2}}{4}\frac{\sigma^{2}\log(ep)}{n}=\max\\{{\mathscr{E}^{*}}^{2}_{\mathcal{T}_{\mathcal{I}_{0}}},{\mathscr{E}^{*}}_{\mathcal{G}_{\mathcal{I}_{0}}}^{2}\\}\frac{\sigma^{2}\log(ep)}{4n}.$

Combining all these cases we finally have the following necessary condition
for consistent model selection with $C_{1},C_{2}>0$ being some absolute
constants:

$\widehat{\tau}_{j_{0}}\geq
C_{1}\max\\{{\mathscr{E}^{*}}_{\mathcal{T}_{\mathcal{I}_{0}}}^{2},{\mathscr{E}^{*}}_{\mathcal{G}_{\mathcal{I}_{0}}}^{2}\\}\frac{\sigma^{2}\log(ep)}{n},\quad\text{if
${\mathscr{E}^{*}}_{\mathcal{G}_{\mathcal{I}_{0}}}\notin({\mathscr{E}^{*}}_{\mathcal{T}_{\mathcal{I}_{0}}},\mathscr{E}_{\mathcal{T}_{\mathcal{I}_{0}}})$},$
(S1.31)

or,

$\widehat{\tau}_{j_{0}}\geq
C_{2}\max\left\\{{\mathscr{E}^{*}}_{\mathcal{T}_{\mathcal{I}_{0}}}^{2},\left(\sqrt{200\mathscr{E}_{\mathcal{T}_{\mathcal{I}_{0}}}^{2}+{\mathscr{E}^{*}}_{\mathcal{G}_{\mathcal{I}_{0}}}^{2}}-10\sqrt{2}\mathscr{E}_{\mathcal{T}_{\mathcal{I}_{0}}}\right)^{2}\right\\}\frac{\sigma^{2}\log(ep)}{n},$

if
${\mathscr{E}^{*}}_{\mathcal{G}_{\mathcal{I}_{0}}}\in({\mathscr{E}^{*}}_{\mathcal{T}_{\mathcal{I}_{0}}},\mathscr{E}_{\mathcal{T}_{\mathcal{I}_{0}}})$.

If there exists $\mathcal{I}_{0}\in\mathcal{J}$ such that
${\mathscr{E}^{*}}_{\mathcal{T}_{\mathcal{I}_{0}}}/\mathscr{E}_{\mathcal{T}_{\mathcal{I}_{0}}}\in(\alpha,1)$,
then in the preceding case it follows that the last display can be simplified
in the following form:

$\widehat{\tau}_{j_{0}}\geq
C_{2}\max\left\\{{\mathscr{E}^{*}}_{\mathcal{T}_{\mathcal{I}_{0}}}^{2},A_{\alpha}{\mathscr{E}^{*}}_{\mathcal{G}_{\mathcal{I}_{0}}}^{2}\right\\}\frac{\sigma^{2}\log(ep)}{n},\quad\text{if
${\mathscr{E}^{*}}_{\mathcal{G}_{\mathcal{I}_{0}}}\in({\mathscr{E}^{*}}_{\mathcal{T}_{\mathcal{I}_{0}}},\mathscr{E}_{\mathcal{T}_{\mathcal{I}_{0}}})$},$

where

$A_{\alpha}=\left(\sqrt{\frac{200}{\alpha^{2}}+1}+\frac{10\sqrt{2}}{\alpha}\right)^{-1}.$

Thus, using the fact that $A_{\alpha}<1$ and combining the previous three
displays we have the necessary condition to be

$\widehat{\tau}_{j_{0}}\geq
C_{2}A_{\alpha}\max\left\\{{\mathscr{E}^{*}}_{\mathcal{T}_{\mathcal{I}_{0}}}^{2},{\mathscr{E}^{*}}_{\mathcal{G}_{\mathcal{I}_{0}}}^{2}\right\\}\frac{\sigma^{2}\log(ep)}{n}.$
(S1.32)

Since, (S1.31) and (S1.32) hold for all choices of $j_{0}$ and
$\mathcal{I}_{0}$ depending on whether the
${\mathscr{E}^{*}}_{\mathcal{G}_{\mathcal{I}_{0}}}\in({\mathscr{E}^{*}}_{\mathcal{T}_{\mathcal{I}_{0}}},\mathscr{E}_{\mathcal{T}_{\mathcal{I}_{0}}})$
is satisfied or not, the claim follows. ∎

###  S1.5 Correlated random feature model example

Consider the model (1). We assume that the rows of $\mathbf{X}$ are
independently generated from $p$-dimensional multivariate Gaussian
distribution with mean-zero and variance-covariance matrix

${\boldsymbol{\Sigma}}=\begin{pmatrix}1&c{\mathbbm{1}}_{p-1}^{\top}\\\
c{\mathbbm{1}}_{p-1}&(1-r)\mathbb{I}_{p-1}+r{\mathbbm{1}}_{p-1}{\mathbbm{1}}_{p-1}^{\top}\end{pmatrix},$

where $c\in[0,0.997],r\in[0,1)$ and true model is $\mathcal{S}=\\{1\\}$. We
need to further impose the restriction

$c^{2}<r+\frac{1-r}{p-1}$

to ensure positive-definiteness of $\Sigma$. In this case

$\widehat{\tau}=\beta_{1}^{2}\min_{j\neq
1}\left\\{\frac{\left\|\mathbf{X}_{1}\right\|_{2}^{2}}{n}-\frac{(\mathbf{X}_{1}^{\top}\mathbf{X}_{j}/n)^{2}}{\left\|\mathbf{X}_{j}\right\|^{2}/n}\right\\}=\frac{\beta_{1}^{2}\left\|\mathbf{X}_{1}\right\|_{2}^{2}}{n}-\beta_{1}^{2}\max_{j\neq
1}\frac{(\mathbf{X}_{1}^{\top}\mathbf{X}_{j}/n)^{2}}{\left\|\mathbf{X}_{j}\right\|_{2}^{2}/n}.$

We start with providing an upper bound on the margin quantity
$\widehat{\tau}$. Using Equation (S5.3) and (S5.4), for any
$\varepsilon_{n,p}\in(0,1)$ we get

${\mathbb{P}}\left(\frac{\left\|\mathbf{X}_{j}\right\|_{2}^{2}}{n}\geq
1+\varepsilon_{n,p}\right)\leq\exp(-n\varepsilon_{n,p}^{2}/16),\quad\forall
j\in[p].$ (S1.33)

${\mathbb{P}}\left(\frac{\left\|\mathbf{X}_{j}\right\|_{2}^{2}}{n}\leq
1-\varepsilon_{n,p}\right)\leq\exp(-n\varepsilon_{n,p}^{2}/4),\quad\forall
j\in[p].$ (S1.34)

Using Equation (S1.33) and Equation (S1.34), we also get the following:

$\displaystyle{\mathbb{P}}\left(\max_{j\neq
1}\left|\frac{\left\|\mathbf{X}_{j}\right\|_{2}^{2}}{n}-1\right|\geq\varepsilon_{n,p}\right)\leq
2p\exp(-n\varepsilon_{n,p}^{2}/16).$ (S1.35)

Let $\mathbf{X}_{j}=(x_{1,j},\ldots,x_{n,j})^{\top}$ for all $j\in[p]$. Now
note that
$\left\|x_{u,1}x_{u,j}\right\|_{\psi_{1}}\leq\left\|x_{u,1}\right\|_{\psi_{2}}\left\|x_{u,2}\right\|_{\psi_{2}}\leq
4.$ Thus, by Berstein’s inequality, we have

${\mathbb{P}}\left(\left|\frac{\mathbf{X}_{1}^{\top}\mathbf{X}_{j}}{n}-c\right|>\varepsilon_{n,p}\right)\leq\exp\left(-Cn\min\\{\varepsilon_{n,p},\varepsilon_{n,p}^{2}\\}\right),$

where $C>0$ is a universal constant. Thus, we have

${\mathbb{P}}\left(\max_{j\neq
1}\left|\frac{\mathbf{X}_{1}^{\top}\mathbf{X}_{j}}{n}-c\right|>\varepsilon_{n,p}\right)\leq
p\exp(-Cn\min\\{\varepsilon_{n,p},\varepsilon_{n,p}^{2}\\}).$ (S1.36)

Combining (S1.34), (S1.35) and (S1.36) we have

$\displaystyle{\mathbb{P}}\left[\left\\{1+\varepsilon_{n,p}-\frac{(c-\varepsilon_{n,p})^{2}}{1+\varepsilon_{n,p}}\right\\}\geq\frac{\widehat{\tau}}{\beta_{1}^{2}}\geq\left\\{1-\varepsilon_{n,p}-\frac{(c+\varepsilon_{n,p})^{2}}{1-\varepsilon_{n,p}}\right\\}\right]$
(S1.37) $\displaystyle\geq
1-\exp(-n\varepsilon_{n,p}^{2}/16)-2p\exp(-n\varepsilon_{n,p}^{2}/4)-p\exp(-C\varepsilon_{n,p}^{2}n)$
$\displaystyle=1+o(1/p),$

if $\varepsilon_{n,p}\asymp\\{(\log p)/n\\}^{1/2}$ and $(\log p)/n$ is small
enough. Similarly, due to Bernstein’s inequality, it can also be shown that

${\mathbb{P}}\left(\max_{j,k\neq
1}\left|\frac{\mathbf{X}_{k}^{\top}\mathbf{X}_{j}}{n}-r\right|\leq\varepsilon_{n,p}\right)\geq
1-p^{2}\exp(-Cn\varepsilon_{n,p}^{2})=1+o(1/p),$ (S1.38)

where $r\in[0,1)$ and with the same conditions on $\varepsilon_{n,p}$.

Next, we will analyze the geometric quantities. In this case, we have

$\widehat{{\boldsymbol{\gamma}}}_{j}=\frac{\mathbf{X}_{1}-\frac{\mathbf{X}_{j}^{\top}\mathbf{X}_{1}}{\left\|\mathbf{X}_{j}\right\|^{2}}.\mathbf{X}_{j}}{\sqrt{\left\|\mathbf{X}_{1}\right\|^{2}-\frac{(\mathbf{X}_{1}^{\top}\mathbf{X}_{j})^{2}}{\left\|\mathbf{X}_{j}\right\|^{2}}}}.$

Note that

$\displaystyle\left\|\widehat{{\boldsymbol{\gamma}}}_{j}-\widehat{{\boldsymbol{\gamma}}}_{k}\right\|_{2}^{2}$
$\displaystyle=2(1-\widehat{{\boldsymbol{\gamma}}}_{j}^{\top}\widehat{{\boldsymbol{\gamma}}}_{k})$

and

$\displaystyle\widehat{{\boldsymbol{\gamma}}}_{j}^{\top}\widehat{{\boldsymbol{\gamma}}}_{k}$
$\displaystyle=\dfrac{\left\|\mathbf{X}_{1}\right\|^{2}/n-\frac{(\mathbf{X}_{j}^{\top}\mathbf{X}_{1}/n)^{2}}{\left\|\mathbf{X}_{j}\right\|^{2}/n}-\frac{(\mathbf{X}_{k}^{\top}\mathbf{X}_{1}/n)^{2}}{\left\|\mathbf{X}_{k}\right\|^{2}/n}+\frac{(\mathbf{X}_{j}^{\top}\mathbf{X}_{1}/n)(\mathbf{X}_{k}^{\top}\mathbf{X}_{1}/n)(\mathbf{X}_{j}^{\top}\mathbf{X}_{k}/n)}{(\left\|\mathbf{X}_{j}\right\|^{2}/n)(\left\|\mathbf{X}_{k}\right\|^{2}/n)}}{\sqrt{\left\|\mathbf{X}_{1}\right\|^{2}/n-\frac{(\mathbf{X}_{1}^{\top}\mathbf{X}_{j}/n)^{2}}{\left\|\mathbf{X}_{j}\right\|^{2}/n}}\sqrt{\left\|\mathbf{X}_{1}\right\|^{2}/n-\frac{(\mathbf{X}_{1}^{\top}\mathbf{X}_{k}/n)^{2}}{\left\|\mathbf{X}_{k}\right\|^{2}/n}}}.$

Next, we consider the event

$\mathcal{G}_{n}:=\left\\{\max_{j\in[p]}\left|\frac{\left\|\mathbf{X}_{j}\right\|_{2}^{2}}{n}-1\right|\leq\varepsilon_{n,p},\max_{j\neq
1}\left|\frac{\mathbf{X}_{1}^{\top}\mathbf{X}_{j}}{n}-c\right|\leq\varepsilon_{n,p},\max_{j,k\neq
1}\left|\frac{\mathbf{X}_{k}^{\top}\mathbf{X}_{j}}{n}-r\right|\leq\varepsilon_{n,p}\right\\}.$

Due to (S1.33), (S1.34), (S1.36) and (S1.38) we have
${\mathbb{P}}(\mathcal{G}_{n})=1+o(1/p)$. Also, for large $n,p$, the value of
$\varepsilon_{n,p}$ can be chosen such that $\varepsilon_{n,p}<0.001$ so that
$c+\varepsilon_{n,p}<0.998$ for all $c\in[0,0.997]$.

##### Complexity of unexplained signals:

Let $\mathbf{u}:=(u_{1},u_{2},u_{3})\in\mathbb{R}^{3}$ and
$\mathbf{t}:=(t_{1},t_{2},t_{2})\in\mathbb{R}^{3}$. Define the function

$\Phi(\mathbf{u},\mathbf{t}):=\frac{u_{1}-(t_{1}^{2}/u_{2})-(t_{2}^{2}/u_{3})+(t_{1}t_{2}t_{3})/(u_{2}u_{3})}{\sqrt{u_{1}-t_{1}^{2}/u_{2}}\sqrt{u_{1}-t_{2}^{2}/u_{3}}},$

where

$(u_{1},u_{2},u_{3},t_{1},t_{2},t_{3})\in\underbrace{[0.999,1.001]\times[0.999,1.001]\times[0.999,1.001]\times[0,0.998]\times[0,0.998]\times[0,1]}_{:=\mathcal{K}}.$

It is easy to see that the function $\Phi$ is continuously differentiable on
the compact set $\mathcal{K}$. Hence, there exists a universal constant $L>0$
such that

$\left|\Phi(\mathbf{u},\mathbf{t})-\Phi(\mathbf{u}^{\prime},\mathbf{t}^{\prime})\right|\leq
L(\left\|\mathbf{u}-\mathbf{u}^{\prime}\right\|_{1}+\left\|\mathbf{t}-\mathbf{t}^{\prime}\right\|_{1}).$

Noting the fact that

$\widehat{{\boldsymbol{\gamma}}}_{j}^{\top}{\boldsymbol{\gamma}}_{k}=\Phi\left(\frac{\left\|\mathbf{X}_{1}\right\|_{2}^{2}}{n},\frac{\left\|\mathbf{X}_{j}\right\|_{2}^{2}}{n},\frac{\left\|\mathbf{X}_{K}\right\|_{2}^{2}}{n},\frac{\mathbf{X}_{1}^{\top}\mathbf{X}_{j}}{n},\frac{\mathbf{X}_{1}^{\top}\mathbf{X}_{k}}{n},\frac{\mathbf{X}_{j}^{\top}\mathbf{X}_{k}}{n}\right),$

it follows that on the event $\mathcal{G}_{n}$, the following holds for all
$j,k\in[p]\setminus\\{1\\}:$

$\displaystyle\left|\left\|\widehat{{\boldsymbol{\gamma}}}_{j}-\widehat{{\boldsymbol{\gamma}}}_{k}\right\|_{2}^{2}-\frac{2c^{2}(1-r)}{1-c^{2}}\right|\leq
12L\varepsilon_{n,p},$
$\displaystyle\Rightarrow\sqrt{\max\left\\{\frac{2c^{2}(1-r)}{1-c^{2}}-12L\varepsilon_{n,p},0\right\\}}\leq\left\|\widehat{{\boldsymbol{\gamma}}}_{j}-\widehat{{\boldsymbol{\gamma}}}_{k}\right\|_{2}\leq\sqrt{\frac{2c^{2}(1-r)}{1-c^{2}}+12L\varepsilon_{n,p}}.$

Hence we have

$\displaystyle\mathscr{E}_{\mathcal{T}_{\varnothing}}$
$\displaystyle=\\{\log(ep)\\}^{-1/2}\Big{[}\int_{0}^{\sqrt{\left(\frac{2c^{2}(1-r)}{1-c^{2}}-12L\varepsilon_{n,p}\right)\vee
0}}\sqrt{\log\mathcal{N}(\varepsilon,\mathcal{T}_{\varnothing},\left\|\cdot\right\|_{2})}\;d\varepsilon$
$\displaystyle+\int_{\sqrt{\left(\frac{2c^{2}(1-r)}{1-c^{2}}-12L\varepsilon_{n,p}\right)\vee
0}}^{\sqrt{\frac{2c^{2}(1-r)}{1-c^{2}}+12L\varepsilon_{n,p}}}\sqrt{\log\mathcal{N}(\varepsilon,\mathcal{T}_{\varnothing},\left\|\cdot\right\|_{2})}\;d\varepsilon\Big{]}.$

Applying Lemma S5.4 on the second integral, it follows that

$\omega_{n,p}\sqrt{\frac{\log
p}{\log(ep)}}\leq\mathscr{E}_{\mathcal{T}_{\varnothing}}\leq\left(\omega_{n,p}+\sqrt{24L\varepsilon_{n,p}}\right)\sqrt{\frac{\log
p}{\log(ep)}},$

where
$\omega_{n,p}=\sqrt{\left(\frac{2c^{2}(1-r)}{1-c^{2}}-12L\varepsilon_{n,p}\right)\vee
0}$. Thus, for $c=0$ we have
$0\leq\mathscr{E}_{\mathcal{T}_{\varnothing}}\leq\sqrt{24L\varepsilon_{n,p}}$.
For any fixed $c>0$ and $r\in[0,1)$ we have

$\mathscr{E}_{\mathcal{T}_{\varnothing}}\sim\left\\{\frac{2c^{2}(1-r)}{1-c^{2}}\right\\}^{1/2}\quad\text{for
large $n,p$}.$ (S1.39)

##### Complexity of spurious projections:

For $j,k\neq 1$, let $\theta_{j,k}$ denote the angle between $\mathbf{X}_{j}$
and $\mathbf{X}_{k}$.

$\left\|\mathbf{P}_{j}-\mathbf{P}_{k}\right\|_{{\rm
op}}=\sin(\theta_{j,k})=\sqrt{1-\cos^{2}(\theta_{j,k})}=\sqrt{1-\left(\frac{\mathbf{X}_{j}^{\top}\mathbf{X}_{k}}{\left\|\mathbf{X}_{j}\right\|_{2}\left\|\mathbf{X}_{k}\right\|_{2}}\right)^{2}}.$

By a similar argument as above, we can conclude that there exists a universal
constant $M>0$ such that on the event $\mathcal{G}_{n}$ we have,

$\left|\left\|\mathbf{P}_{j}-\mathbf{P}_{k}\right\|_{\rm
op}^{2}-(1-r^{2})\right|\leq M\varepsilon_{n,p},\quad\text{for all
$j,k\in[p]\setminus\\{1\\}$}.$

Thus, for any fixed $r\in[0,1)$ we have

$\mathscr{E}_{\mathcal{G}_{\varnothing}}\sim(1-r^{2})^{1/2}.$ (S1.40)

##  S2 Generalized linear model

In this section, we will focus on the best subset selection problem under
generalized linear models (GLM). Similar to the linear regression setup, we
will also adopt the fixed design setup in this case. In particular, given the
data matrix
$\mathbf{X}:=(\mathbf{x}_{1},\ldots,\mathbf{x}_{n})^{\top}\in\mathbb{R}^{n\times
p}$ we observe the responses $\mathbf{y}:=(y_{1},\ldots,y_{n})^{\top}$ coming
from the distribution

$f_{\mathbf{x},{\boldsymbol{\beta}}^{*}}(y):=h(y)\exp\left\\{\frac{y(\mathbf{x}^{\top}{\boldsymbol{\beta}}^{*})-b(\mathbf{x}^{\top}{\boldsymbol{\beta}}^{*})}{\phi}\right\\}=h(y)\exp\left\\{\frac{y\eta-b(\eta)}{\phi}\right\\}.$
(S2.1)

Here $\eta=\mathbf{x}^{\top}{\boldsymbol{\beta}}^{*}$ is linear predictor and
${\boldsymbol{\beta}}^{*}$ is true parameter with
$\left\|{\boldsymbol{\beta}}^{*}\right\|_{0}=s$ and support $\mathcal{S}$. The
functions $b:\mathbb{R}\to\mathbb{R}$ and $h:\mathbb{R}\to\mathbb{R}$ are
known and specific to modeling assumptions. Examples include several well-
known models such as

  1. 1.

Linear regression: Consider the linear regression model
$y=\mathbf{x}^{\top}{\boldsymbol{\beta}}^{*}+\varepsilon$, where
$\varepsilon\sim{\sf N}(0,\sigma^{2})$. In this case
$h(y)=\exp\\{-y^{2}/(2\sigma^{2})\\}$ and $b(u)=u^{2}/2$.

  2. 2.

Logistic regression: In this model
$y\sim{\sf{Ber}}(1/(1+\exp(-\mathbf{x}^{\top}{\boldsymbol{\beta}}^{*})))$.
Standard calculations show that $h(y)=1$ and $b(u)=\log(1+e^{u})$.

For the purpose of model selection, we choose the loss function to be the
scaled negative log-likelihood function

$\mathcal{L}({\boldsymbol{\beta}};\\{(\mathbf{x}_{i},y_{i})\\}_{i\in[n]})=\frac{2}{n}\sum_{i\in[n]}\ell({\boldsymbol{\beta}};(\mathbf{x}_{i},y_{i})),$

where
$\ell({\boldsymbol{\beta}};(\mathbf{x},y))=-y(\mathbf{x}^{\top}{\boldsymbol{\beta}})+b(\mathbf{x}^{\top}{\boldsymbol{\beta}}).$
Furthermore, for a candidate model $\mathcal{D}\in\mathscr{A}_{s}$ and
$\tilde{{\boldsymbol{\beta}}}\in\mathbb{R}^{s}$, define the restricted version
of the scaled negative log-likelihood function as
$\mathcal{L}_{\mathcal{D}}(\tilde{{\boldsymbol{\beta}}};\\{(\mathbf{x}_{i,\mathcal{D}},y_{i})\\}_{i\in[n]}):=(n/2)^{-1}\sum_{i\in[n]}\ell_{\mathcal{D}}(\tilde{{\boldsymbol{\beta}}};(\mathbf{x}_{i,\mathcal{D}},y_{i}))$
with
$\ell_{\mathcal{D}}(\tilde{{\boldsymbol{\beta}}};(\mathbf{x}_{\mathcal{D}},y)):=-y(\mathbf{x}_{\mathcal{D}}^{\top}\tilde{{\boldsymbol{\beta}}})+b(\mathbf{x}_{\mathcal{D}}^{\top}\tilde{{\boldsymbol{\beta}}})$.
Let $\widehat{{\boldsymbol{\beta}}}_{\mathcal{D}}$ be the unique minimizer of
$\mathcal{L}_{\mathcal{D}}(\tilde{{\boldsymbol{\beta}}};\\{(\mathbf{x}_{i,\mathcal{D}},y_{i})\\}_{i\in[n]})$.
Under the oracle knowledge of sparsity $s$, BSS solves for

$\widehat{\mathcal{S}}_{\rm
best}=n^{-1}{\arg\min}_{\mathcal{D}:\left|\mathcal{D}\right|=s}\mathcal{L}_{\mathcal{D}}(\widehat{{\boldsymbol{\beta}}}_{\mathcal{D}};\\{(\mathbf{x}_{i,\mathcal{D}},y_{i})\\}_{i\in[n]}).$

Next, we will introduce the quantities that capture the degree of separation
between the true model $\mathcal{S}$ and a candidate model
$\mathcal{D}\in\mathscr{A}_{s}$ and characterize the identifiability margin
for model selection consistency. Let
$\mathcal{P}_{\mathcal{D},\tilde{{\boldsymbol{\beta}}}}$ be probability
measure corresponding to the joint density
$\prod_{i\in[n]}f_{\mathbf{x}_{i,\mathcal{D}},\tilde{{\boldsymbol{\beta}}}}(y_{i})$
and define

$\displaystyle\Delta_{\rm kl}(\mathcal{D})$
$\displaystyle:=\frac{2\phi}{n}\min_{\tilde{{\boldsymbol{\beta}}}\in\mathbb{R}^{s}}{\sf
KL}\left(\mathcal{P}_{\mathcal{S},{\boldsymbol{\beta}}^{*}_{\mathcal{S}}}\;\big{\|}\;\mathcal{P}_{\mathcal{D},\tilde{{\boldsymbol{\beta}}}}\right)$
$\displaystyle=\frac{2}{n}\sum_{i=1}^{n}\left\\{(\mathbf{x}_{i,\mathcal{S}}^{\top}{\boldsymbol{\beta}}^{*}_{\mathcal{S}})b^{\prime}(\mathbf{x}_{i,\mathcal{S}}^{\top}{\boldsymbol{\beta}}^{*}_{\mathcal{S}})-b(\mathbf{x}_{i,\mathcal{S}}^{\top}{\boldsymbol{\beta}}^{*}_{\mathcal{S}})\right\\}$
$\displaystyle\quad-\max_{\tilde{{\boldsymbol{\beta}}}\in\mathbb{R}^{s}}\frac{2}{n}\sum_{i=1}^{n}\left\\{(\mathbf{x}_{i,\mathcal{D}}^{\top}\tilde{{\boldsymbol{\beta}}})b^{\prime}(\mathbf{x}_{i,\mathcal{S}}^{\top}{\boldsymbol{\beta}}^{*}_{\mathcal{S}})-b(\mathbf{x}_{i,\mathcal{D}}^{\top}\tilde{{\boldsymbol{\beta}}})\right\\}.$

The above quantity can be thought of as the degree of model separation as it
measures the minimum KL-distance between the likelihood generated by the data
under $(\mathcal{S},{\boldsymbol{\beta}}^{*}_{\mathcal{S}})$ and the
likelihood generated by $\mathcal{D}$ and all possible choices of
${\boldsymbol{\beta}}\in\mathbb{R}^{p}$ with the support in $\mathcal{D}$. Let
$\bar{{\boldsymbol{\beta}}}_{\mathcal{D}}$ be the minimizer of the
optimization problem in the above display, i.e.,

$\bar{{\boldsymbol{\beta}}}_{\mathcal{D}}:={\arg\min}_{\tilde{{\boldsymbol{\beta}}}\in\mathbb{R}^{s}}{\sf
KL}\left(\mathcal{P}_{\mathcal{S},{\boldsymbol{\beta}}^{*}_{\mathcal{S}}}\;\big{\|}\;\mathcal{P}_{\mathcal{D},\tilde{{\boldsymbol{\beta}}}}\right).$

By definition it follows that
$\bar{{\boldsymbol{\beta}}}_{\mathcal{S}}={\boldsymbol{\beta}}^{*}_{\mathcal{S}}$.
Also, note that for
$\mathcal{P}_{\mathcal{D},\bar{{\boldsymbol{\beta}}}_{\mathcal{D}}}$, the
natural parameter of the density function is
$\mathbf{X}_{\mathcal{D}}\bar{{\boldsymbol{\beta}}}_{\mathcal{D}}$. Thus, one
can also measure the separation between two models through the mutual distance
between the corresponding natural parameters. This motivates the definition of
the second measure of separability between the true model $\mathcal{S}$ and
candidate model $\mathcal{D}$:

$\Delta_{\rm
par}(\mathcal{D}):=\frac{\left\|\mathbf{X}_{\mathcal{S}}{\boldsymbol{\beta}}_{\mathcal{S}}^{*}-\mathbf{X}_{\mathcal{D}}\bar{{\boldsymbol{\beta}}}_{\mathcal{D}}\right\|_{2}^{2}}{n}.$

Note that, under the linear regression model with isotropic Gaussian error,
both $\Delta_{\rm kl}(\mathcal{D})$ and $\Delta_{\rm par}(\mathcal{D})$
becomes equal to the quantity
${\boldsymbol{\beta}}_{\mathcal{S}\setminus\mathcal{D}}^{\top}\Gamma(\mathcal{D}){\boldsymbol{\beta}}_{\mathcal{S}\setminus\mathcal{D}}$.
To see this, recall that for linear regression model $b(u)=u^{2}/2$ and the
KL-divergence ${\sf
KL}(\mathcal{P}_{\mathcal{S},{\boldsymbol{\beta}}^{*}_{\mathcal{S}}}\;\big{\|}\;\mathcal{P}_{\mathcal{D},\bar{{\boldsymbol{\beta}}}_{\mathcal{D}}})=\left\|\mathbf{X}_{\mathcal{S}}{\boldsymbol{\beta}}_{\mathcal{S}}^{*}-\mathbf{X}_{\mathcal{D}}\bar{{\boldsymbol{\beta}}}_{\mathcal{D}}\right\|_{2}^{2}/(2\sigma^{2})$.
Thus, from the definition of $\bar{{\boldsymbol{\beta}}}_{\mathcal{D}}$, it
immediately follows that
$\mathbf{X}_{\mathcal{D}}\bar{{\boldsymbol{\beta}}}_{\mathcal{D}}=\mathbf{P}_{\mathcal{D}}\mathbf{X}_{\mathcal{S}}{\boldsymbol{\beta}}^{*}_{\mathcal{S}}$.
Later, we will see that these two notions of distances are equivalent under
certain regularity conditions on the link function $b(\cdot)$.

###  S2.1 Identifiability margin and two complexities

In this section we will introduce the identifiability margin and the two
complexities similar the case of linear model. We consider the following
identifiability margin:

$\tilde{\tau}_{*}(s):=\min_{\mathcal{D}\in\mathscr{A}_{s}}\frac{\Delta_{\rm
kl}(\mathcal{D})}{\left|\mathcal{D}\setminus\mathcal{S}\right|}.$

We assume that $\tilde{\tau}_{*}(s)>0$ to avoid non-identifiability issue.
Next, we consider the transformed features as follows:

$\widetilde{\mathbf{X}}_{\mathcal{D}}={\boldsymbol{\Lambda}}_{\mathcal{D}}^{1/2}\mathbf{X}_{\mathcal{D}},$

where
${\boldsymbol{\Lambda}}_{\mathcal{D}}=\mathbf{diag}(b^{\prime\prime}(\mathbf{x}_{1,\mathcal{D}}^{\top}\bar{{\boldsymbol{\beta}}}_{\mathcal{D}}),\ldots,b^{\prime\prime}(\mathbf{x}_{n,\mathcal{D}}^{\top}\bar{{\boldsymbol{\beta}}}_{\mathcal{D}}))$.
Let $\widetilde{\mathbf{P}}_{\mathcal{D}}$ be orthogonal projection matrices
onto the columnspace of $\widetilde{\mathbf{X}}_{\mathcal{D}}$. Let
$\widetilde{\mathbf{P}}_{\mathcal{I}\mid\mathcal{D}}$ be the orthogonal
projector onto the columnspace of
$[\widetilde{\mathbf{X}}_{\mathcal{D}}]_{\mathcal{I}}$. Now we define the
following sets of residualized signals and spurious projections:

${\widetilde{\mathcal{T}}}_{\mathcal{I}}^{(s)}=\left\\{\frac{\mathbf{X}_{\mathcal{D}}\bar{{\boldsymbol{\beta}}}_{\mathcal{D}}-\mathbf{X}_{\mathcal{S}}{\boldsymbol{\beta}}_{\mathcal{S}}^{*}}{\left\|\mathbf{X}_{\mathcal{D}}\bar{{\boldsymbol{\beta}}}_{\mathcal{D}}-\mathbf{X}_{\mathcal{S}}{\boldsymbol{\beta}}_{\mathcal{S}}^{*}\right\|_{2}}:\mathcal{D}\in\mathscr{A}_{\mathcal{I}}\right\\},$

${\widetilde{\mathcal{G}}}_{\mathcal{I}}^{(s)}=\left\\{\widetilde{\mathbf{P}}_{\mathcal{D}}-\widetilde{\mathbf{P}}_{\mathcal{I}\mid\mathcal{D}}:\mathcal{D}\in\mathscr{A}_{\mathcal{I}}\right\\}.$

The complexity measures for these two sets are
$\mathscr{E}_{{\widetilde{\mathcal{T}}}_{\mathcal{I}}^{(s)}}$ and
$\mathscr{E}_{{\widetilde{\mathcal{G}}}_{\mathcal{I}^{(s)}}}$ respectively,
which are defined in the same way as the complexity measures in Section 3.2
and Section 3.3 of the main paper.

###  S2.2 Main results

In this section we will state the main result analogous to the Theorem 1 in
the main paper. We begin with some standard assumption necessary for the
theoretical analysis for GLM models.

######  Assumption S2.1 (Features and parameters).

We assume the following conditions:

  1. (a)

There exists positive constants $x_{0}$ and $R_{0}$ such that
$\max_{i\in[n]}\left\|\mathbf{x}_{i}\right\|_{\infty}\leq x_{0}$ and
$\left\|{\boldsymbol{\beta}}^{*}\right\|_{1}\leq R_{0}$.

  2. (b)

There exists a constant $\kappa_{0}>0$ such that

$\min_{\mathcal{D}\subset[p]:\left|\mathcal{D}\right|=s}\lambda_{\min}\left(\mathbf{X}_{\mathcal{D}}^{\top}\mathbf{X}_{\mathcal{D}}/n\right)\geq\kappa_{0}.$

  3. (c)

There exists constant $M>0$ such that

$\max_{\mathcal{D}\subset[p]:\left|\mathcal{D}\right|=s}\left\|\frac{1}{n}\sum_{i\in[n]}\mathbf{x}_{i,\mathcal{D}}\otimes\mathbf{x}_{i,\mathcal{D}}\otimes\mathbf{x}_{i,\mathcal{D}}\right\|_{{\rm
op}}\leq M.$

  4. (d)

There exists a constant $R>0$ such that
$\max_{\mathcal{D}\in\mathscr{A}_{s}}\max_{i\in[n]}\left|\mathbf{x}_{i,\mathcal{D}}^{\top}\bar{{\boldsymbol{\beta}}}_{\mathcal{D}}\right|\leq
x_{0}R$.

  5. (e)

The design matrix $\mathbf{X}$ enjoys the following property: