# Optimality of Graphlet Screening in High Dimensional Variable Selection

Jiashun Jin<EMAIL_ADDRESS>
Department of Statistics
Carnegie Mellon University
Pittsburgh, PA 15213, USA Cun-Hui Zhang<EMAIL_ADDRESS>
Department of Statistics
Rutgers University
Piscataway, NJ 08854, USA Qi Zhang<EMAIL_ADDRESS>
Department of Biostatistics $\&$ Medical Informatics
University of Wisconsin-Madison
Madison, WI 53705, USA

###### Abstract

Consider a linear model $Y=X\beta+\sigma z$, where $X$ has $n$ rows and $p$
columns and $z\sim N(0,I_{n})$. We assume both $p$ and $n$ are large,
including the case of $p\gg n$. The unknown signal vector $\beta$ is assumed
to be sparse in the sense that only a small fraction of its components is
nonzero. The goal is to identify such nonzero coordinates (i.e., variable
selection).

We are primarily interested in the regime where signals are both rare and weak
so that successful variable selection is challenging but is still possible. We
assume the Gram matrix $G=X^{\prime}X$ is sparse in the sense that each row
has relatively few large entries (diagonals of $G$ are normalized to $1$). The
sparsity of $G$ naturally induces the sparsity of the so-called Graph of
Strong Dependence (GOSD). The key insight is that there is an interesting
interplay between the signal sparsity and graph sparsity: in a broad context,
the signals decompose into many small-size components of GOSD that are
disconnected to each other.

We propose Graphlet Screening for variable selection. This is a two-step
Screen and Clean procedure, where in the first step, we screen subgraphs of
GOSD with sequential $\chi^{2}$-tests, and in the second step, we clean with
penalized MLE. The main methodological innovation is to use GOSD to guide both
the screening and cleaning processes.

For any variable selection procedure $\hat{\beta}$, we measure its performance
by the Hamming distance between the sign vectors of $\hat{\beta}$ and $\beta$,
and assess the optimality by the minimax Hamming distance. Compared with more
stringent criterions such as exact support recovery or oracle property, which
demand strong signals, the Hamming distance criterion is more appropriate for
weak signals since it naturally allows a small fraction of errors.

We show that in a broad class of situations, Graphlet Screening achieves the
optimal rate of convergence in terms of the Hamming distance. Unlike Graphlet
Screening, well-known procedures such as the $L^{0}/L^{1}$-penalization
methods do not utilize local graphic structure for variable selection, so they
generally do not achieve the optimal rate of convergence, even in very simple
settings and even if the tuning parameters are ideally set.

The the presented algorithm is implemented as R-CRAN package ScreenClean and
in _matlab_ (available at
http://www.stat.cmu.edu/$\sim$jiashun/Research/software/GS-matlab/).

> Keywords: Asymptotic minimaxity, Graph of Least Favorables (GOLF), Graph of
> Strong Dependence (GOSD), Graphlet Screening (GS), Hamming distance, phase
> diagram, Rare and Weak signal model, Screen and Clean, sparsity.

##  1 Introduction

Consider a linear regression model

$Y=X\beta+\sigma z,\qquad X=X_{n,p},\qquad z\sim N(0,I_{n}).$ (1.1)

We write

$X=[x_{1},x_{2},\ldots,x_{p}],\qquad\mbox{and}\qquad
X^{\prime}=[X_{1},X_{2},\ldots,X_{n}],$ (1.2)

so that $x_{j}$ is the $j$-th design vector and $X_{i}$ is the $i$-th sample.
Motivated by the recent interest in ‘Big Data’, we assume both $p$ and $n$ are
large but $p\geq n$ (though this should not be taken as a restriction). The
vector $\beta$ is unknown to us, but is presumably sparse in the sense that
only a small proportion of its entries is nonzero. Calling a nonzero entry of
$\beta$ a signal, the main interest of this paper is to identify all signals
(i.e., variable selection).

Variable selection is one of the most studied problem in statistics. However,
there are important regimes where our understanding is very limited.

One of such regimes is the rare and weak regime, where the signals are both
rare (or sparse) and individually weak. Rare and weak signals are frequently
found in research areas such as Genome-wide Association Study (GWAS) or next
generation sequencing. Unfortunately, despite urgent demand in applications,
the literature of variable selection has been focused on the regime where the
signals are rare but individually strong. This motivates a revisit to variable
selection, focusing on the rare and weak regime.

For variable selection in this regime, we need new methods and new theoretical
frameworks. In particular, we need a loss function that is appropriate for
rare and weak signals to evaluate the optimality. In the literature, given a
variable selection procedure $\hat{\beta}$, we usually use the probability of
exact recovery $P(\hbox{\rm sgn}(\hat{\beta})\neq\hbox{\rm sgn}(\beta))$ as
the measure of loss (Fan and Li, 2001); $\hbox{\rm sgn}(\hat{\beta})$ and
$\hbox{\rm sgn}(\beta)$ are the sign vectors of $\hat{\beta}$ and $\beta$
respectively. In the rare and weak regime, the signals are so rare and weak
that exact recovery is impossible, and the Hamming distance between $\hbox{\rm
sgn}(\hat{\beta})$ and $\hbox{\rm sgn}(\beta)$ is a more appropriate measure
of loss.

Our focus on the rare and weak regime and the Hamming distance loss provides
new perspectives to variable selection, in methods and in theory.

Throughout this paper, we assume the diagonals of the Gram matrix

$G=X^{\prime}X$ (1.3)

are normalized to $1$ (and approximately $1$ in the random design model),
instead of $n$ as often used in the literature. The difference between two
normalizations is non-essential, but the signal vector $\beta$ are different
by a factor of $n^{1/2}$.

We also assume the Gram matrix $G$ is ‘sparse’ (aka. graph sparsity) in the
sense that each of its rows has relatively few large entries. Signal sparsity
and graph sparsity can be simultaneously found in the following application
areas.

  * •

Compressive sensing. We are interested in a very high dimensional sparse
vector $\beta$. The goal is to store or transmit $n$ linear functionals of
$\beta$ and then reconstruct it. For $1\leq i\leq n$, we choose a
$p$-dimensional coefficient vector $X_{i}$ and observe
$Y_{i}=X_{i}^{\prime}\beta+\sigma z_{i}$ with an error $\sigma z_{i}$. The so-
called Gaussian design is often considered (Donoho, 2006a, b; Bajwa et al,
2007), where $X_{i}\stackrel{{\scriptstyle iid}}{{\sim}}N(0,\Omega/n)$ and
$\Omega$ is sparse; the sparsity of $\Omega$ induces the sparsity of
$G=X^{\prime}X$.

  * •

Genetic Regulatory Network (GRN). For $1\leq i\leq n$,
$W_{i}=(W_{i}(1),\ldots,W_{i}(p))^{\prime}$ represents the expression level of
$p$ different genes of the $i$-th patient. Approximately,
$W_{i}\stackrel{{\scriptstyle iid}}{{\sim}}N(\alpha,\Sigma)$, where the
contrast mean vector $\alpha$ is sparse reflecting that only few genes are
differentially expressed between a normal patient and a diseased one (Peng et
al, 2009). Frequently, the concentration matrix $\Omega=\Sigma^{-1}$ is
believed to be sparse, and can be effectively estimated in some cases (e.g.,
Bickel and Levina (2008) and Cai et al (2010)), or can be assumed as known in
others, with the so-called “data about data” available (Li and Li, 2011). Let
$\hat{\Omega}$ be a positive-definite estimate of $\Omega$, the setting can be
re-formulated as the linear model
$(\hat{\Omega})^{1/2}Y\approx\Omega^{1/2}Y\sim N(\Omega^{1/2}\beta,I_{p})$,
where $\beta=\sqrt{n}\alpha$ and the Gram matrix $G\approx\Omega$, and both
are sparse.

Other examples can be found in Computer Security (Ji and Jin, 2011) and Factor
Analysis (Fan et al, 2011).

The sparse Gram matrix $G$ induces a sparse graph which we call the Graph of
Strong Dependence (GOSD), denoted by ${\cal G}=(V,E)$, where
$V=\\{1,2,\ldots,p\\}$ and there is an edge between nodes $i$ and $j$ if and
only the design vectors $x_{i}$ and $x_{j}$ are strongly correlated. Let

$S=S(\beta)=\\{1\leq j\leq p:\beta_{j}\neq 0\\}$ (1.4)

be the support of $\beta$ and ${\cal G}_{S}$ be the subgraph of ${\cal G}$
formed by all nodes in $S$. The key insight is that, there is an interesting
interaction between signal sparsity and graph sparsity, which yields the
subgraph ${\cal G}_{S}$ decomposable: ${\cal G}_{S}$ splits into many
“graphlet”; each “graphlet” is a small-size component and different components
are not connected (in ${\cal G}_{S}$).

While we can always decompose ${\cal G}_{S}$ in this way, our emphasis in this
paper is that, in many cases, the maximum size of the graphlets is small; see
Lemma 1 and related discussions.

The decomposability of ${\cal G}_{S}$ motivates a new approach to variable
selection, which we call Graphlet Screening (GS). GS is a Screen and Clean
method (Wasserman and Roeder, 2009). In the screening stage, we use
multivariate screening to identify candidates for all the graphlets. Let
$\hat{S}$ be all the nodes that survived the screening, and let ${\cal
G}_{\hat{S}}$ be the subgraph of GOSD formed by all nodes in $\hat{S}$.
Although $\hat{S}$ is expected to be somewhat larger than $S$, the subgraph
${\cal G}_{\hat{S}}$ is still likely to resemble ${\cal G}_{S}$ in structure
in the sense that it, too, splits into many small-size disconnected
components. We then clean each component separately to remove false positives.

The objective of the paper is two-fold.

  * •

To propose a “fundamentally correct” solution in the rare and weak paradigm
along with a computationally fast algorithm for the solution.

  * •

To show that GS achieves the optimal rate of convergence in terms of the
Hamming distance, and achieves the optimal phase diagram for variable
selection.

Phase diagram can be viewed as an optimality criterion which is especially
appropriate for rare and weak signals. See Donoho and Jin (2004) and Jin
(2009) for example.

In the settings we consider, most popular approaches are not rate optimal; we
explain this in Sections 1.1-1.3. In Section 1.4, we explain the basic idea of
GS and why it works.

###  1.1 Non-optimality of the $L^{0}$-penalization method for rare and weak
signals

When $\sigma=0$, Model (1.1) reduces to the “noiseless” model $Y=X\beta$. In
this model, Donoho and Stark (1989) (see also Donoho and Huo (2001)) reveals a
fundamental phenomenon on sparse representation. Fix $(X,Y)$ and consider the
equation $Y=X\beta$. Since $p>n$, the equation has infinitely many solutions.
However, a very sparse solution, if exists, is unique under mild conditions on
the design $X$, with all other solutions being much denser. In fact, if the
sparsest solution $\beta_{0}$ has $k$ elements, then all other solutions of
the equation $Y=X\beta$ must have at least $(\mathrm{rank}(X)-k+1)$ nonzero
elements, and $\mathrm{rank}(X)=n$ when $X$ is in a “general position”.

From a practical viewpoint, we frequently believe that this unique sparse
solution is the truth (i.e., Occam’s razor). Therefore, the problem of
variable selection can be solved by some global methods designed for finding
the sparsest solution to the equation $Y=X\beta$.

Since the $L^{0}$-norm is (arguably) the most natural way to measure the
sparsity of a vector, the above idea suggests that the $L^{0}$-penalization
method is a “fundamentally correct” (but computationally intractable) method
for variable selection, provided that some mild conditions hold (noiseless,
Signal-to-Noise Ratio (SNR) is high, signals are sufficiently sparse (Donoho
and Stark, 1989; Donoho and Huo, 2001)).

Motivated by this, in the past two decades, a long list of computationally
tractable algorithms have been proposed that approximate the solution of the
$L^{0}$-penalization method, including the lasso, SCAD, MC+, and many more
(Akaike, 1974; Candes and Tao, 2007; Efron et al, 2004; Fan and Li, 2001;
Schwarz, 1978; Tibshirani, 1996; Zhang, 2010, 2011; Zhao and Yu, 2006; Zou,
2006).

With that being said, we must note that these methodologies were built upon a
framework with four tightly woven core components: “signals are rare but
strong”, “the truth is also the sparsest solution to $Y=X\beta$”, “probability
of exact recovery is an appropriate loss function”, and “$L^{0}$-penalization
method is a fundamentally correct approach”. Unfortunately, when signals are
rare and weak, such a framework is no longer suitable.

  * •

When signals are “rare and weak”, the fundamental uniqueness property of the
sparse solution in the noiseless case is no longer valid in the noisy case.
Consider the model $Y=X\beta+\sigma z$ and suppose that a sparse $\beta_{0}$
is the true signal vector. There are many vectors $\beta$ that are small
perturbations of $\beta_{0}$ such that the two models $Y=X\beta+\sigma z$ and
$Y=X\beta_{0}+\sigma z$ are indistinguishable (i.e., all tests are
asymptotically powerless). In the “rare and strong” regime, $\beta_{0}$ is the
sparsest solution among all such “eligible” solutions of $Y=X\beta+\sigma z$.
However, this claim no longer holds in the “rare and weak” regime and the
principle of Occam’s razor may not be as relevant as before.

  * •

The $L^{0}$-penalization method is originally designed for “rare and strong”
signals where “exact recovery” is used to measure its performance (Donoho and
Stark, 1989; Donoho and Huo, 2001; Donoho, 2006a). When we must consider “rare
and weak” signals and when we use the Hamming distance as the loss function,
it is unclear whether the $L^{0}$-penalization method is still “fundamentally
correct”.

In fact, in Section 2.8 (see also Ji and Jin (2011)), we show that the
$L^{0}$-penalization method is not optimal in Hamming distance when signals
are rare and weak, even with very simple designs (i.e., Gram matrix is
tridiagonal or block-wise) and even when the tuning parameter is ideally set.
Since the $L^{0}$-penalization method is used as the benchmark in the
development of many other penalization methods, its sub-optimality is expected
to imply the sub-optimality of other methods designed to match its performance
(e.g., lasso, SCAD, MC+).

###  1.2 Limitation of Univariate Screening and UPS

Univariate Screening (also called marginal regression or Sure Screening (Fan
and Lv, 2008; Genovese et al, 2012)) is a well-known variable selection
method. For $1\leq j\leq p$, recall that $x_{j}$ is the $j$-th column of $X$.
Univariate Screening selects variables with large marginal correlations:
$|(x_{j},Y)|$, where $(\cdot,\cdot)$ denotes the inner product. The method is
computationally fast, but it can be seriously corrupted by the so-called
phenomenon of “signal cancellation” (Wasserman and Roeder, 2009). In our model
(1.1)-(1.3), the SNR associated with $(x_{j},Y)$ is

$\frac{1}{\sigma}\sum_{\ell=1}^{p}(x_{j},x_{\ell})\beta_{\ell}=\frac{\beta_{j}}{\sigma}+\frac{1}{\sigma}\sum_{\ell\neq
j}(x_{j},x_{\ell})\beta_{\ell}.$ (1.5)

“Signal cancellation” happens if SNR is significantly smaller than
$\beta_{j}/\sigma$. For this reason, the success of Univariate Screening needs
relatively strong conditions (e.g., Faithfulness Condition (Genovese et al,
2012)), under which signal cancellation does not have a major effect.

In Ji and Jin (2011), Ji and Jin proposed Univariate Penalized Screening (UPS)
as a refinement of Univariate Screening, where it was showed to be optimal in
the rare and weak paradigm, for the following two scenarios. The first
scenario is where the nonzero effects of variables are all positively
correlated: $(x_{j}\beta_{j})^{\prime}(x_{k}\beta_{k})\geq 0$ for all
$\\{j,k\\}$. This guarantees the faithfulness of the univariate association
test. The second scenario is a Bernoulli model where the “signal cancellation”
only has negligible effects over the Hamming distance of UPS.

With that being said, UPS attributes its success mostly to the cleaning stage;
the screening stage of UPS uses nothing but Univariate Screening, so UPS does
not adequately address the challenge of “signal cancellation”. For this
reason, we should not expect UPS to be optimal in much more general settings.

###  1.3 Limitations of Brute-force Multivariate Screening

One may attempt to overcome “signal cancellation” by multivariate screening,
with Brute-force Multivariate Screening (BMS) being the most straightforward
version. Fix an integer $1\leq m_{0}\ll p$. BMS consists of a series of
screening phases, indexed by $m$, $1\leq m\leq m_{0}$, that are increasingly
more ambitious. In Phase-$m$ BMS, we test the significance of the association
between $Y$ and any set of $m$ different design variables
$\\{x_{j_{1}},x_{j_{2}},\ldots,x_{j_{m}}\\}$, $j_{1}<j_{2}<\ldots<j_{m}$, and
retain all such design variables if the test is significant. The problem of
BMS is, it enrolls too many candidates for screening, which is both
unnecessary and unwise.

  * •

(Screening inefficiency). In Phase-$m$ of BMS, we test about $\binom{p}{m}$
hypotheses involving different subsets of $m$ design variables. The larger the
number of hypotheses we consider, the higher the threshold we need to set for
the tests, in order to control the false positives. When we enroll too many
candidates for hypothesis testing, we need signals that are stronger than
necessary in order for them to survive the screening.

  * •

(Computational challenge). Testing $\binom{p}{m}$ hypotheses is
computationally infeasible when $p$ is large, even when $m$ is very small
(say, $(p,m)=(10^{4},3)$).

###  1.4 Graphlet Screening: how it is different and how it works

Graphlet Screening (GS) uses a similar screening strategy as BMS does, except
for a major difference. When it comes to the test of significance between $Y$
and design variables $\\{x_{j_{1}},x_{j_{2}},\ldots,x_{j_{m}}\\}$,
$j_{1}<j_{2}<\ldots<j_{m}$, GS only carries out such a test if
$\\{j_{1},j_{2},\ldots,j_{m}\\}$ is a connected subgraph of the GOSD.
Otherwise, the test is safely skipped!

Fixing an appropriate threshold $\delta>0$, we let $\Omega^{*,\delta}$ be the
regularized Gram matrix:

$\Omega^{*,\delta}(i,j)=G(i,j)1\\{|G(i,j)|\geq\delta\\},\qquad 1\leq i,j\leq
p.$ (1.6)

The GOSD ${\cal G}\equiv{\cal G}^{*,\delta}=(V,E)$ is the graph where
$V=\\{1,2,\ldots,p\\}$ and there is an edge between nodes $i$ and $j$ if and
only if $\Omega^{*,\delta}(i,j)\neq 0$. See Section 2.6 for the choice of
$\delta$.

Remark. GOSD and ${\cal G}$ are generic terms which vary from case to case,
depending on $G$ and $\delta$. GOSD is very different from the Bayesian
conditional independence graphs (Pearl, 2000).

Fixing $m_{0}\geq 1$ as in BMS, we define

${\cal A}(m_{0})={\cal A}(m_{0};G,\delta)=\\{\mbox{all connected subgraphs of
${\cal G}^{*,\delta}$ with size $\leq m_{0}$}\\}.$ (1.7)

GS is a Screen and Clean method, consisting of a graphical screening step (GS-
step) and a graphical cleaning step (GC-step).

  * •

$GS$-step. We test the significance of association between $Y$ and
$\\{x_{j_{1}},x_{j_{2}},\ldots,x_{j_{m}}\\}$ if and only if
$\\{j_{1},j_{2},\ldots,j_{m}\\}\in{\cal A}(m_{0})$ (i.e., graph guided
multivariate screening). Once $\\{j_{1},\ldots,j_{m}\\}$ is retained, it
remains there until the end of the GS-step.

  * •

$GC$-step. The set of surviving nodes decompose into many small-size
components, which we fit separately using an efficient low-dimensional test
for small graphs.

GS is similar to Wasserman and Roeder (2009) for both of them have a screening
and a cleaning stage, but is more sophisticated. For clarification, note that
Univariate Screening or BMS introduced earlier does not contain a cleaning
stage and can be viewed as a counterpart of the GS-step.

We briefly explain why GS works. We discuss the GS-step and GC-step
separately.

Consider the GS-step first. Compared with BMS, the GS-step recruits far fewer
candidates for screening, so it is able to overcome the two major shortcomings
of BMS aforementioned: high computational cost and low statistical efficiency.
In fact, fix $K\geq 1$ and suppose ${\cal G}^{*,\delta}$ is $K$-sparse (see
Section 1.5 for the definition). By a well-known result in graph theory
(Frieze and Molloy, 1999),

$|{\cal A}(m_{0})|\leq Cm_{0}p(eK)^{m_{0}}.$ (1.8)

The right hand side is much smaller than the term $\binom{p}{m_{0}}$ as we
encounter in BMS.

At the same time, recall that $S=S(\beta)$ is the support of $\beta$. Let
${\cal G}_{S}\equiv{\cal G}_{S}^{*,\delta}$ be the subgraph of ${\cal
G}^{*,\delta}$ consisting all signal nodes. We can always split ${\cal
G}_{S}^{*,\delta}$ into “graphlets” (arranged lexicographically) as follows:

${\cal G}_{S}^{*,\delta}={\cal G}_{S,1}^{*,\delta}\cup{\cal
G}_{S,2}^{*,\delta}\ldots\cup{\cal G}_{S,M}^{*,\delta},$ (1.9)

where each ${\cal G}_{S,i}^{*,\delta}$ is a component (i.e., a maximal
connected subgraph) of ${\cal G}_{S}^{*,\delta}$, and different ${\cal
G}_{S,i}^{*,\delta}$ are not connected in ${\cal G}_{S}^{*,\delta}$. Let

$m_{0}^{*}=m_{0}^{*}(S(\beta),G,\delta)=\max_{1\leq i\leq M}|{\cal
G}_{S,i}^{*,\delta}|$ (1.10)

be the maximum size of such graphlets (note that $M$ also depends on
($S(\beta)$, $G$, $\delta$)).

In many cases, $m_{0}^{*}$ is small. One such case is when we have a Bernoulli
signal model.

###### Lemma 1

Fix $K\geq 1$ and $\epsilon>0$. If ${\cal G}^{*,\delta}$ is $K$-sparse and
$\hbox{\rm sgn}(|\beta_{1}|),\hbox{\rm sgn}(|\beta_{2}|),\ldots,\hbox{\rm
sgn}(|\beta_{p}|)$ are $iid$ from $\mathrm{Bernoulli}(\epsilon)$, then except
for a probability $p(e\epsilon K)^{m_{0}+1}$,
$m_{0}^{*}(S(\beta),G,\delta)\leq m_{0}$.

Lemma 1 is not tied to the Bernoulli model and holds more generally. For
example, it holds when $\\{\hbox{\rm sgn}(|\beta_{i}|)\\}_{i=1}^{p}$ are
generated according to certain Ising models (Ising, 1925).

We recognize that in order for the GS-step to be efficient both in screening
and in computation, it is sufficient that

$m_{0}\geq m_{0}^{*}.$ (1.11)

In fact, first, if (1.11) holds, then for each $1\leq\ell\leq M$, ${\cal
G}_{S,\ell}^{*,\delta}\in{\cal A}(m_{0})$. Therefore, at some point of the
screening process of the GS-step, we must have considered a significance test
between $Y$ and the set of design variables $\\{x_{j}:j\in{\cal
G}_{S,\ell}^{*,\delta}\\}$. Consequently, the GS-step is able to overcome the
“signal cancellations” (the explanation is a little bit long, and we slightly
defer it). Second, since $m_{0}^{*}$ is small in many situations, we could
choose a relatively small $m_{0}$ such that (1.11) holds. When $m_{0}$ is
small, as long as $K$ is small or moderately large, the GS-step is
computationally feasible. In fact, the right hand side of (1.8) is only larger
than $p$ by a moderate factor. See Section 2.2 for more discussion on the
computation complexity.

We now explain the first point above. The notations below are frequently used.

###### Definition 2

For $X$ in Models (1.1)-(1.2) and any subset ${\cal I}\subset\\{1,2,...,p\\}$,
let $P^{{\cal I}}=P^{{\cal I}}(X)$ be the projection from $\mathbb{R}^{n}$ to
the subspace spanned by $\\{x_{j}:j\in{\cal I}\\}$.

###### Definition 3

For an $n\times p$ matrix $A$ and sets ${\cal I}\subset\\{1,\ldots,n\\}$ and
${\cal J}\subset\\{1,\ldots,p\\}$, $A^{{\cal I},{\cal J}}$ is the $|{\cal
I}|\times|{\cal J}|$ sub-matrix formed by restricting the rows of $A$ to
${\cal I}$ and columns to ${\cal J}$.

When $p=1$, $A$ is a vector, and $A^{{\cal I}}$ is the sub-vector of $A$
formed by restricting the rows of $A$ to ${\cal I}$ . When ${\cal
I}=\\{1,2,\ldots,n\\}$ (or ${\cal J}=\\{1,2,\ldots,p\\}$), we write $A^{{\cal
I},{\cal J}}$ as $A^{\otimes,{\cal J}}$ (or $A^{{\cal I},\otimes}$). Note that
indices in ${\cal I}$ or ${\cal J}$ are not necessarily sorted ascendingly.

Recall that for each $1\leq\ell\leq M$, at some point of the GS-step, we must
have considered a significance test between $Y$ and the set of design
variables $\\{x_{j}:j\in{\cal G}_{S,\ell}^{*,\delta}\\}$. By (1.9), we rewrite
Model (1.1) as

$Y=\sum_{\ell=1}^{M}X^{\otimes,{\cal G}_{S,\ell}^{*,\delta}}\beta^{{\cal
G}_{S,\ell}^{*,\delta}}+\sigma z,\qquad z\sim N(0,I_{n}).$

The key is the set of matrices $\\{X^{\otimes,{\cal
G}_{S,\ell}^{*,\delta}}:1\leq\ell\leq M\\}$ are nearly orthogonal (i.e., for
any column $\xi$ of $X^{\otimes,{\cal G}_{S,k}^{*,\delta}}$ and any column
$\eta$ of $X^{\otimes,{\cal G}_{S,\ell}^{*,\delta}}$, $|(\xi,\eta)|$ is small
when $k\neq\ell$).

When we test the significance between $Y$ and $\\{x_{j},j\in{\cal
G}_{S,\ell}^{*,\delta}\\}$, we are testing the null hypothesis $\beta^{{\cal
G}_{S,\ell}^{*,\delta}}=0$ against the alternative $\beta^{{\cal
G}_{S,\ell}^{*,\delta}}\neq 0$. By the near orthogonality aforementioned,
approximately, $(X^{\otimes,{\cal G}_{S,\ell}^{*,\delta}})^{\prime}Y$ is a
sufficient statistic for $\beta^{{\cal G}_{S,\ell}^{*,\delta}}$, and the
optimal test is based on the $\chi^{2}$-test statistic $\|P^{{\cal
G}_{S,\ell}^{*,\delta}}Y\|^{2}$.

The near orthogonality also implies that significant “signal cancellation”
only happens among signals within the same graphlet. When we screen each
graphlet as a whole using the $\chi^{2}$-statistic above, “signal
cancellation” between different graphlets only has negligible effects. In this
way, GS-step is able to retain all nodes in ${\cal G}_{S,\ell}^{*,\delta}$ in
a nearly optimal way, and so overcome the challenge of “signal cancellation”.
This explains the first point.

Note that the GS-step consists of a sequence of sub-steps, each sub-step is
associated with an element of ${\cal A}(m_{0})$. When we screen ${\cal
G}_{S,\ell}^{*,\delta}$ as a whole, it is possible some of the nodes have
already been retained in the previous sub-steps. In this case, we implement
the $\chi^{2}$-test slightly differently, but the insight is similar. See
Section 2.1 for details.

We now discuss the GC-step. Let $\hat{S}$ be all the surviving nodes of the
GS-step, and let ${\cal G}_{\hat{S}}^{*,\delta}$ be the subgraph of ${\cal
G}^{*,\delta}$ formed by confining all nodes to $\hat{S}$. Similarly, we have
(a) the decomposition ${\cal G}_{\hat{S}}^{*,\delta}={\cal
G}_{\hat{S},1}^{*,\delta}\cup{\cal G}_{\hat{S},2}^{*,\delta}\ldots\cup{\cal
G}_{\hat{S},\hat{M}}^{*,\delta}$, (b) the near orthogonality between the
$\hat{M}$ different matrices, each is formed by $\\{x_{j}:j\in{\cal
G}_{\hat{S},\ell}^{*,\delta}\\}$. Moreover, a carefully tuned screening stage
of the GS ensures that most of the components ${\cal
G}_{\hat{S},\ell}^{*,\delta}$ are only small perturbations of their
counterparts in the decomposition of ${\cal G}_{S}^{*,\delta}={\cal
G}_{S,1}^{*,\delta}\cup{\cal G}_{S,2}^{*,\delta}\ldots\cup{\cal
G}_{S,M}^{*,\delta}$ as in (1.9), and the maximum size of ${\cal
G}_{\hat{S},\ell}^{*,\delta}$ is not too much larger than
$m_{0}^{*}=m_{0}^{*}(S(\beta),G,\delta)$. Together, these allow us to clean
${\cal G}_{\hat{S},\ell}^{*,\delta}$ separately, without much loss of
efficiency. Since the maximum size of ${\cal G}_{\hat{S},\ell}^{*,\delta}$ is
small, the computational complexity in the cleaning stage is moderate.

###  1.5 Content

The remaining part of the paper is organized as follows. In Section 2, we show
that GS achieves the minimax Hamming distance in the Asymptotic Rare and Weak
(ARW) model, and use the phase diagram to visualize the optimality of GS, and
to illustrate the advantage of GS over the $L^{0}/L^{1}$-penalization methods.
In Section 3, we explain that GS attributes its optimality to the so-called
Sure Screening property and the Separable After Screening property, and use
these two properties to prove our main result, Theorem 8. Section 4 contains
numeric results, Section 5 discusses more connections to existing literature
and possible extensions of GS, and Section 6 contains technical proofs.

Below are some notations we use in this paper. $L_{p}$ denotes a generic
multi-$\log(p)$ term that may vary from occurrence to occurrence; see
Definition 5. For a vector $\beta\in R^{p}$, $\|\beta\|_{q}$ denotes the
$L^{q}$-norm, and when $q=2$, we drop $q$ for simplicity. For two vectors
$\alpha$ and $\beta$ in $R^{p}$, $\alpha\circ\beta\in R^{p}$ denotes the
vector in $R^{p}$ that satisfies $(\alpha\circ\beta)_{i}=\alpha_{i}\beta_{i}$,
$1\leq i\leq p$; “$\circ"$ is known as the Hadamard product.

For an $n\times p$ matrix $A$, $\|A\|_{\infty}$ denotes the matrix
$L^{\infty}$-norm, and $\|A\|$ denotes the spectral norm (Horn and Johnson,
1990). Recall that for two sets ${\cal I}$ and ${\cal J}$ such that ${\cal
I}\subset\\{1,2,\ldots,n\\}$ and ${\cal J}\subset\\{1,2,\ldots,p\\}$,
$A^{{\cal I},{\cal J}}$ denotes the submatrix of $A$ formed by restricting the
rows and columns of $A$ to ${\cal I}$ and ${\cal J}$, respectively. Note that
the indices in ${\cal I}$ and ${\cal J}$ are not necessarily sorted in the
ascending order. In the special case where ${\cal I}=\\{1,2,\ldots,n\\}$ (or
${\cal J}=\\{1,2,\ldots,p\\}$), we write $A^{{\cal I},{\cal J}}$ as
$A^{\otimes,{\cal J}}$ (or $A^{{\cal I},\otimes}$). In the special case where
$n=p$ and $A$ is positive definite, $\lambda_{k}^{*}(A)$ denotes the minimum
eigenvalue of all the size $k$ principal submatrices of $A$, $1\leq k\leq p$.
For $X$ in (1.1), $P^{{\cal I}}$ denotes the projection to the column space of
$X^{\otimes,{\cal I}}$.

Recall that in Model (1.1), $Y=X\beta+\sigma z$. Fixing a threshold
$\delta>0$. Let $G=X^{\prime}X$ be the Gram matrix, and let
$\Omega^{*,\delta}$ be the regularized Gram matrix defined by
$\Omega^{*,\delta}(i,j)=G(i,j)1\\{|G(i,j)|\geq\delta\\}$, $1\leq i,j\leq p$.
Let ${\cal G}^{*,\delta}$ be the graph where each index in
$\\{1,2,\ldots,p\\}$ is a node, and there is an edge between node $i$ and node
$j$ if and only if $\Omega^{*,\delta}(i,j)\neq 0$. We let $S(\beta)$ be the
support of $\beta$, and denote ${\cal G}_{S}^{*,\delta}$ by the subgraph of
${\cal G}^{*,\delta}$ formed by all nodes in $S(\beta)$. We call ${\cal
G}^{*,\delta}$ the Graph of Strong Dependence (GOSD) and sometimes write it by
${\cal G}$ for short. The GOSD and ${\cal G}$ are generic notations which
depend on $(G,\delta)$ and may vary from occurrence to occurrence. We also
denote ${\cal G}^{\diamond}$ by the Graph of Least Favorable (GOLF). GOLF only
involves the study of the information lower bound. For an integer $K\geq 1$, a
graph ${\cal G}$, and one of its subgraph ${\cal I}_{0}$, we write ${\cal
I}_{0}\triangleleft{\cal G}$ if and only if ${\cal I}_{0}$ is a component of
${\cal G}$ (i.e., a maximal connected subgraph of ${\cal G}$), and we call
${\cal G}$ $K$-sparse if its maximum degree is no greater than $K$.

##  2 Main results

In Section 2.1, we formally introduce GS. In Section 2.2, we discuss the
computational complexity of GS. In Sections 2.3-2.6, we show that GS achieves
the optimal rate of convergence in the Asymptotic Rare and Weak model. In
Sections 2.7-2.8, we introduce the notion of phase diagram and use it to
compare GS with the $L^{0}/L^{1}$-penalization methods. We conclude the
section with a summary in Section 2.9.

###  2.1 Graphlet Screening: the procedure

GS consists of a GS-step and a GC-step. We describe two steps separately.
Consider the $GS$-step first. Fix $m_{0}\geq 1$ and $\delta>0$, recall that
${\cal G}^{*,\delta}$ denotes the GOSD and ${\cal A}(m_{0})$ consists of all
connected subgraphs of ${\cal G}^{*,\delta}$ with size $\leq m_{0}$.

  * •

Initial sub-step. Let ${\cal U}_{p}^{*}=\emptyset$. List all elements in
${\cal A}(m_{0})$ in the ascending order of the number of nodes it contains,
with ties broken lexicographically. Since a node is thought of as connected to
itself, the first $p$ connected subgraphs on the list are simply the nodes
$1,2,\ldots,p$. We screen all connected subgraphs in the order they are
listed.

  * •

Updating sub-step. Let ${\cal I}_{0}$ be the connected subgraph under
consideration, and let ${\cal U}_{p}^{*}$ be the current set of retained
indices. We update ${\cal U}_{p}^{*}$ with a $\chi^{2}$ test as follows. Let
$\hat{F}={\cal I}_{0}\cap{\cal U}_{p}^{*}$ and $\hat{D}={\cal
I}_{0}\setminus{\cal U}_{p}^{*}$, so that $\hat{F}$ is the set of nodes in
${\cal I}_{0}$ that have already been accepted, and $\hat{D}$ is the set of
nodes in ${\cal I}_{0}$ that is currently under investigation. Note that no
action is needed if $\hat{D}=\emptyset$. For a threshold
$t(\hat{D},\hat{F})>0$ to be determined, we update ${\cal U}_{p}^{*}$ by
adding all nodes in $\hat{D}$ to it if

$T(Y,\hat{D},\hat{F})=\|P^{{\cal
I}_{0}}Y\|^{2}-\|P^{\hat{F}}Y\|^{2}>t(\hat{D},\hat{F}),$ (2.12)

and we keep ${\cal U}_{p}^{*}$ the same otherwise (by default,
$\|P^{\hat{F}}Y\|=0$ if $\hat{F}=\emptyset$). We continue this process until
we finish screening all connected subgraphs on the list. The final set of
retained indices is denoted by ${\cal U}_{p}^{*}$.

See Table $1$ for a recap of the procedure. In the $GS$-step, once a node is
kept in any sub-stage of the screening process, it remains there until the end
of the $GS$-step (however, it may be killed in the $GC$-step). This has a
similar flavor to that of the Forward regression.

In principle, the procedure depends on how the connected subgraphs of the same
size are initially ordered, and different ordering could give different
numeric results. However, such differences are usually negligibly small.
Alternatively, one could revise the procedure so that it does not depend on
the ordering. For example, in the updating sub-step, we could choose to update
${\cal U}_{p}^{*}$ only when we finish screening all connected sub-graphs of
size $k$, $1\leq k\leq m_{0}$. While the theoretic results below continue to
hold if we revise GS in this way, we must note that from a numeric
perspective, the revision would not produce a very different result. For
reasons of space, we skip discussions along this line.

The GS-step uses a set of tuning parameters:

${\cal Q}\equiv\\{t(\hat{D},\hat{F}):\mbox{$(\hat{D},\hat{F})$ are as defined
in (\ref{DefineDF})}\\}.$

A convenient way to set these parameters is to let
$t(\hat{D},\hat{F})=2\sigma^{2}q\log p$ for a fixed $q>0$ and all
$(\hat{D},\hat{F})$. More sophisticated choices are given in Section 2.6.

The $GS$-step has two important properties: Sure Screening and Separable After
Screening (SAS). With tuning parameters ${\cal Q}$ properly set, the Sure
Screening property says that ${\cal U}_{p}^{*}$ retains all but a negligible
fraction of the signals. The SAS property says that as a subgraph of ${\cal
G}^{*,\delta}$, ${\cal U}_{p}^{*}$ decomposes into many disconnected
components, each has a size $\leq\ell_{0}$ for a fixed small integer
$\ell_{0}$. Together, these two properties enable us to reduce the original
large-scale regression problem to many small-size regression problems that can
be solved parallelly in the $GC$-step. See Section 3 for elaboration on these
ideas.

We now discuss the $GC$-step. For any $1\leq j\leq p$, we have either
$j\notin{\cal U}_{p}^{*}$, or that there is a unique connected subgraph ${\cal
I}_{0}$ such that $j\in{\cal I}_{0}\lhd{\cal U}_{p}^{*}$. In the first case,
we estimate $\beta_{j}$ as $0$. In the second case, for two tuning parameters
$u^{gs}>0$ and $v^{gs}>0$, we estimate the whole set of variables
$\beta^{{\cal I}_{0}}$ by minimizing the functional

$\|P^{{\cal I}_{0}}(Y-X^{\otimes,{\cal
I}_{0}}\xi)\|^{2}+(u^{gs})^{2}\|\xi\|_{0}$ (2.13)

over all $|{\cal I}_{0}|\times 1$ vectors $\xi$, each nonzero coordinate of
which $\geq v^{gs}$ in magnitude. The resultant estimator is the final
estimate of GS which we denote by
$\hat{\beta}^{gs}=\hat{\beta}^{gs}(Y;\delta,{\cal Q},u^{gs},v^{gs},X,p,n)$.
See Section 1.5 for notations used in this paragraph.

Table 1: Graphlet Screening Algorithm. $GS$-step: | List ${\cal G}^{*,\delta}$-connected submodels ${\cal I}_{0,k}$ with $|{\cal I}_{0,1}|\leq|{\cal I}_{0,2}|\leq\cdots\leq m_{0}$
---|---
| Initialization: ${\cal U}_{p}^{*}=\emptyset$ and $k=1$
| Test $H_{0}:{\cal I}_{0,k}\cap{\cal U}_{p}^{*}$ against $H_{1}:{\cal
I}_{0,k}$ with $\chi^{2}$ test (2.12)
| Update: ${\cal U}_{p}^{*}\leftarrow{\cal U}_{p}^{*}\cup{\cal I}_{0,k}$ if
$H_{0}$ rejected, $k\leftarrow k+1$
$GC$-step: | As a subgraph of ${\cal G}^{*,\delta}$, ${\cal U}_{p}^{*}$ decomposes into many components ${\cal I}_{0}$
| Use the $L^{0}$-penalized test (2.13) to select a subset ${\hat{\cal
I}}_{0}$ of each ${\cal I}_{0}$
| Return the union of ${\hat{\cal I}}_{0}$ as the selected model

Sometimes for linear models with random designs, the Gram matrix $G$ is very
noisy, and GS is more effective if we use it iteratively for a few times
($\leq 5$). This can be implemented in a similar way as that in Ji and Jin
(2011, Section 3). Here, the main purpose of iteration is to denoise $G$, not
for variable selection. See Ji and Jin (2011, Section 3) and Section 4 for
more discussion.

###  2.2 Computational complexity

If we exclude the overhead of obtaining ${\cal G}^{*,\delta}$, then the
computation cost of GS contains two parts, that of the GS-step and that of the
GC-step. In each part, the computation cost hinges on the sparsity of ${\cal
G}^{*,\delta}$. In Section 2.3, we show that with a properly chosen $\delta$,
for a wide class of design matrices, ${\cal G}^{*,\delta}$ is $K$-sparse for
some $K=K_{p}\leq C\log^{\alpha}(p)$ as $p\rightarrow\infty$, where $\alpha>0$
is a constant. As a result (Frieze and Molloy, 1999),

$|{\cal A}(m_{0})|\leq pm_{0}(eK_{p})^{m_{0}}\leq
Cm_{0}p\log^{m_{0}\alpha}(p).$ (2.14)

We now discuss two parts separately.

In the GS-step, the computation cost comes from that of listing all elements
in ${\cal A}(m_{0})$, and that of screening all connected-subgraphs in ${\cal
A}(m_{0})$. Fix $1\leq k\leq m_{0}$. By (2.14) and the fact that every size
$k$ ($k>1$) connected subgraph at least contains one size $k-1$ connected
subgraph, greedy algorithm can be used to list all sub-graphs with size $k$
with computational complexity $\leq Cp(K_{p}k)^{k}\leq Cp\log^{k\alpha}(p)$,
and screening all connected subgraphs of size $k$ has computational complexity
$\leq Cnp\log^{k\alpha}(p)$. Therefore, the computational complexity of the
GS-step $\leq Cnp(\log(p))^{(m_{0}+1)\alpha}$.

The computation cost of the $GC$-step contains the part of breaking ${\cal
U}_{p}^{*}$ into disconnected components, and that of cleaning each component
by minimizing (2.13). As a well-known application of the breadth-first search
(Hopcroft and Tarjan, 1973), the first part $\leq|{\cal U}_{p}^{*}|(K_{p}+1)$.
For the second part, by the SAS property of the $GS$-step (i.e., Lemma 16),
for a broad class of design matrices, with the tuning parameters chosen
properly, there is a fixed integer $\ell_{0}$ such that with overwhelming
probability, $|{\cal I}_{0}|\leq\ell_{0}$ for any ${\cal I}_{0}\lhd{\cal
U}_{p}^{*}$. As a result, the total computational cost of the $GC$-step is no
greater than $C(2^{\ell_{0}}\log^{\alpha}(p))|{\cal U}_{p}^{*}|n$, which is
moderate.

The computational complexity of GS is only moderately larger than that of
Univariate Screening or UPS (Ji and Jin, 2011). UPS uses univariate
thresholding for screening which has a computational complexity of $O(np)$,
and GS implements multivariate screening for all connected subgraphs in ${\cal
A}(m_{0})$, which has a computational complexity $\leq
Cnp(\log(p))^{(m_{0}+1)\alpha}$. The latter is only larger by a
multi-$\log(p)$ term.

###  2.3 Asymptotic Rare and Weak model and Random Design model

To analyze GS, we consider the regression model $Y=X\beta+\sigma z$ as in
(1.1), and use an Asymptotic Rare and Weak (ARW) model for $\beta$ and a
random design model for $X$.

We introduce the ARW first. Fix parameters $\epsilon\in(0,1)$, $\tau>0$, and
$a\geq 1$. Let $b=(b_{1},\ldots,b_{p})^{\prime}$ be the $p\times 1$ random
vector where

$b_{i}\stackrel{{\scriptstyle iid}}{{\sim}}\mathrm{Bernoulli}(\epsilon).$
(2.15)

We model the signal vector $\beta$ in Model (1.1) by

$\beta=b\circ\mu,$ (2.16)

where “$\circ$” denotes the Hadamard product (see Section 1.5) and
$\mu\in\Theta_{p}^{*}(\tau,a)$, with

$\Theta_{p}^{*}(\tau,a)=\\{\mu\in\Theta_{p}(\tau),\|\mu\|_{\infty}\leq
a\tau\\},\qquad\Theta_{p}(\tau)=\\{\mu\in\mathbb{R}^{p}:\,|\mu_{i}|\geq\tau,1\leq
i\leq p\\}.$ (2.17)

In this model, $\epsilon$ calibrates the sparsity level and $\tau$ calibrates
the minimum signal strength. We are primarily interested in the case where
$\epsilon$ is small and $\tau$ is smaller than the required signal strength
for the exact recovery of the support of $\beta$, so the signals are both rare
and weak. The constraint of $\|\mu\|_{\infty}\leq a\tau_{p}$ is mainly for
technical reasons (only needed for Lemma 16); see Section 2.6 for more
discussions.

We let $p$ be the driving asymptotic parameter, and tie $(\epsilon,\tau)$ to
$p$ through some fixed parameters. In detail, fixing $0<\vartheta<1$, we model

$\epsilon=\epsilon_{p}=p^{-\vartheta}.$ (2.18)

For any fixed $\vartheta$, the signals become increasingly sparser as
$p\rightarrow\infty$. Also, as $\vartheta$ ranges, the sparsity level ranges
from very dense to very sparse, and covers all interesting cases.

It turns out that the most interesting range for $\tau$ is
$\tau=\tau_{p}=O(\sqrt{\log(p)})$. In fact, when
$\tau_{p}\ll\sigma\sqrt{\log(p)}$, the signals are simply too rare and weak so
that successful variable selection is impossible. On the other hand, exact
support recovery requires $\tau\gtrsim\sigma\sqrt{2\log p}$ for orthogonal
designs and possibly even larger $\tau$ for correlated designs. In light of
this, we fix $r>0$ and calibrate $\tau$ by

$\tau=\tau_{p}=\sigma\sqrt{2r\log(p)}.$ (2.19)

Next, consider the random design model. The use of random design model is
mainly for simplicity in presentation. The main results in the paper can be
translated to fixed design models with a careful modification of the
notations; see Corollary 7 and Section 5.

For any positive definite matrix $A$, let $\lambda(A)$ be the smallest
eigenvalue, and let

$\lambda_{k}^{*}(\Omega)=\min\\{\lambda(A):\mbox{$A$ is a $k\times k$
principle submatrix of $\Omega$}\\}.$ (2.20)

For $m_{0}$ as in the GS-step, let $g=g(m_{0},\vartheta,r)$ be the smallest
integer such that

$g\geq\max\\{m_{0},(\vartheta+r)^{2}/(2\vartheta r)\\}.$ (2.21)

Fixing a constant $c_{0}>0$, introduce

${\cal M}_{p}(c_{0},g)=\\{\Omega:\mbox{$p\times p$ correlation matrix,
$\lambda_{g}^{*}(\Omega)\geq c_{0}$}\\}.$ (2.22)

Recall $X_{i}$ is the $i$-th row of $X$; see (1.2). In the random design
model, we fix an $\Omega\in{\cal M}(c_{0},g)$ ($\Omega$ is unknown to us), and
assume

$X_{i}\stackrel{{\scriptstyle iid}}{{\sim}}N(0,\frac{1}{n}\Omega),\qquad 1\leq
i\leq n.$ (2.23)

In the literature, this is called the Gaussian design, which can be found in
Compressive Sensing (Bajwa et al, 2007), Computer Security (Dinur and Nissim,
2003), and other application areas.

At the same time, fixing $\kappa\in(0,1)$, we model the sample size $n$ by

$n=n_{p}=p^{\kappa}.$ (2.24)

As $p\rightarrow\infty$, $n_{p}$ becomes increasingly large but is still much
smaller than $p$. We assume

$\kappa>(1-\vartheta),$ (2.25)

so that $n_{p}\gg p\epsilon_{p}$. Note $p\epsilon_{p}$ is approximately the
total number of signals. Condition (2.25) is almost necessary for successful
variable selection (Donoho, 2006a, b).

###### Definition 4

We call model (2.15)-(2.19) for $\beta$ the Asymptotic Rare Weak model
$\mathrm{ARW}(\vartheta,r,a,\mu)$, and call Model (2.23)-(2.25) for $X$ the
Random Design model $\mathrm{RD}(\vartheta,\kappa,\Omega)$.

###  2.4 Minimax Hamming distance

In many works on variables selection, one assesses the optimality by the
‘oracle property’, where the probability of non-exact recovery $P(\hbox{\rm
sgn}(\hat{\beta})\neq\hbox{\rm sgn}(\beta))$ is the loss function. When
signals are rare and weak, $P(\hbox{\rm sgn}(\hat{\beta})\neq\hbox{\rm
sgn}(\beta))\approx 1$ and ‘exact recovery’ is usually impossible. A more
appropriate loss function is the Hamming distance between $\hbox{\rm
sgn}(\hat{\beta})$ and $\hbox{\rm sgn}(\beta)$.

For any fixed $\beta$ and any variable selection procedure $\hat{\beta}$, we
measure the performance by the Hamming distance:

$h_{p}(\hat{\beta},\beta\bigl{|}X)=E\Bigl{[}\sum_{j=1}^{p}1\bigl{\\{}\hbox{\rm
sgn}(\hat{\beta}_{j})\neq\hbox{\rm
sgn}(\beta_{j})\bigr{\\}}\bigr{|}X\Bigr{]}.$

In the Asymptotic Rare Weak model, $\beta=b\circ\mu$, and
$(\epsilon_{p},\tau_{p})$ depend on $p$ through $(\vartheta,r)$, so the
overall Hamming distance for $\hat{\beta}$ is

$H_{p}(\hat{\beta};\epsilon_{p},n_{p},\mu,\Omega)=E_{\epsilon_{p}}E_{\Omega}\bigl{[}h_{p}(\hat{\beta},\beta\bigr{|}X)\bigr{]}\equiv
E_{\epsilon_{p}}E_{\Omega}\bigl{[}h_{p}(\hat{\beta},b\circ\mu\bigr{|}X)\bigr{]},$

where $E_{\epsilon_{p}}$ is the expectation with respect to the law of $b$,
and $E_{\Omega}$ is the expectation with respect to the law of $X$; see (2.15)
and (2.23). Finally, the minimax Hamming distance is

$\mathrm{Hamm}_{p}^{*}(\vartheta,\kappa,r,a,\Omega)=\inf_{\hat{\beta}}\sup_{\mu\in\Theta_{p}^{*}(\tau_{p},a)}\bigl{\\{}H_{p}(\hat{\beta};\epsilon_{p},n_{p},\mu,\Omega)\bigr{\\}}.$
(2.26)

The Hamming distance is no smaller than the sum of the expected number of
signal components that are misclassified as noise and the expected number of
noise components that are misclassified as signal.

###  2.5 Lower bound for the minimax Hamming distance, and GOLF

We first construct lower bounds for “local risk” at different $j$, $1\leq
j\leq p$, and then aggregate them to construct a lower bound for the global
risk. One challenge we face is the least favorable configurations for
different $j$ overlap with each other. We resolve this by exploiting the
sparsity of a new graph to be introduced: Graph of Least Favorable (GOLF).

To recap, the model we consider is Model (1.1), where

$\beta\text{ is modeled by
}\mathrm{ARW}(\vartheta,r,a,\mu),\qquad\mbox{and}\qquad X\text{ is modeled by
}\mathrm{RD}(\vartheta,\kappa,\Omega).$ (2.27)

Fix $1\leq j\leq p$. The “local risk” at an index $j$ is the risk of
estimating the set of variables $\\{\beta_{k}:d(k,j)\leq g\\}$, where $g$ is
defined in (2.21) and $d(j,k)$ denotes the geodesic distance between $j$ and
$k$ in the graph ${\cal G}^{*,\delta}$. The goal is to construct two subsets
$V_{0}$ and $V_{1}$ and two realizations of $\beta$, $\beta^{(0)}$ and
$\beta^{(1)}$ (in the special case of $V_{0}=V_{1}$, we require $\hbox{\rm
sgn}(\beta^{(0)})\neq\hbox{\rm sgn}(\beta^{(1)})$), such that $j\in V_{0}\cup
V_{1}$ and

$\mbox{If $k\notin V_{0}\cup V_{1}$, $\beta^{(0)}_{k}=\beta_{k}^{(1)}$;
otherwise, $\beta^{(i)}_{k}\neq 0$ if and only if $k\in V_{i}$, $i=0,1$}.$

In the literature, it is known that how well we can estimate
$\\{\beta_{k}:d(k,j)\leq g\\}$ depends on how well we can separate two
hypotheses (where $\beta^{(0)}$ and $\beta^{(1)}$ are assumed as known):

$H_{0}^{(j)}:Y=X\beta^{(0)}+\sigma z\qquad\mbox{vs}.\qquad
H_{1}^{(j)}:Y=X\beta^{(1)}+\sigma z,\qquad z\sim N(0,I_{n}).$ (2.28)

The least favorable configuration for the local risk at index $j$ is the
quadruple $(V_{0},V_{1},\beta^{(0)},\beta^{(1)})$ for which two hypotheses are
the most difficult to separate.

For any $V\subset\\{1,2,\ldots,p\\}$, let $I_{V}$ be the indicator vector of
$V$ such that for any $1\leq k\leq p$, the $k$-th coordinate of $I_{V}$ is $1$
if $k\in V$ and is $0$ otherwise. Define

$B_{V}=\\{I_{V}\circ\mu:\;\mu\in\Theta_{p}^{*}(\tau_{p},a)\\},$ (2.29)

where we recall “$\circ$” denotes the Hadamard product (see Section 1.5).
Denote for short $\theta^{(i)}=I_{V_{0}\cup V_{1}}\circ\beta^{(i)}$, and so
$\beta^{(1)}-\beta^{(0)}=\theta^{(1)}-\theta^{(0)}$ and $\theta^{(i)}\in
B_{V_{i}}$, $i=0,1$. Introduce

$\alpha(\theta^{(0)},\theta^{(1)})=\alpha(\theta^{(0)},\theta^{(1)};V_{0},V_{1},\Omega,a)=\tau_{p}^{-2}(\theta^{(0)}-\theta^{(1)})^{\prime}\Omega(\theta^{(0)}-\theta^{(1)}).$

For the testing problem in (2.28), the optimal test is to reject $H_{0}^{(j)}$
if and only if
$(\theta^{(1)}-\theta^{(0)})^{\prime}X^{\prime}(Y-X\beta^{(0)})\geq
t\sigma\tau_{p}\sqrt{\alpha(\theta^{(0)},\theta^{(1)})}$ for some threshold
$t>0$ to be determined. In the ARW and RD models, $P(\beta_{k}\neq 0,\forall
k\in V_{i})\sim\epsilon_{p}^{|V_{i}|}$, $i=0,1$, and
$(\theta^{(0)}-\theta^{(1)})^{\prime}G(\theta^{(0)}-\theta^{(1)})\approx(\theta^{(0)}-\theta^{(1)})^{\prime}\Omega(\theta^{(0)}-\theta^{(1)})$,
since the support of $\theta^{(0)}-\theta^{(1)}$ is contained in a small-size
set $V_{0}\cup V_{1}$. Therefore the sum of Type I and Type II error of any
test associated with (2.28) is no smaller than (up to some negligible
differences)

$\epsilon_{p}^{|V_{0}|}\bar{\Phi}(t)+\epsilon_{p}^{|V_{1}|}\Phi\bigl{(}t-(\tau_{p}/\sigma)[\alpha(\theta^{(0)},\theta^{(1)})]^{1/2}\bigr{)},$
(2.30)

where $\bar{\Phi}=1-\Phi$ is the survival function of $N(0,1)$.

For a lower bound for the “local risk” at $j$, we first optimize the quantity
in (2.30) over all $\theta^{(0)}\in B_{V_{0}}$ and $\theta^{(1)}\in
B_{V_{1}}$, and then optimize over all $(V_{0},V_{1})$ subject to $j\in
V_{0}\cup V_{1}$. To this end, define
$\alpha^{*}(V_{0},V_{1})=\alpha^{*}(V_{0},V_{1};a,\Omega)$,
$\eta(V_{0},V_{1})=\eta(V_{0},V_{1};\vartheta,r,a,\Omega)$, and
$\rho_{j}^{*}=\rho_{j}^{*}(\vartheta,r,a,\Omega)$ by

$\alpha^{*}(V_{0},V_{1})=\min\bigl{\\{}\alpha(\theta^{(0)},\theta^{(1)};V_{0},V_{1},\Omega,a):\;\theta^{(i)}\in
B_{V_{i}},i=0,1,\hbox{\rm sgn}(\theta^{(0)})\neq\hbox{\rm
sgn}(\theta^{(1)})\\},$ (2.31)
$\eta(V_{0},V_{1})=\max\\{|V_{0}|,|V_{1}|\\}\vartheta+\frac{1}{4}\left[\left(\sqrt{\alpha^{*}(V_{0},V_{1})r}-\frac{\bigl{|}(|V_{1}|-|V_{0}|)\bigr{|}\vartheta}{\sqrt{\alpha^{*}(V_{0},V_{1})r}}\right)_{+}\right]^{2},$

and

$\rho_{j}^{*}(\vartheta,r,a,\Omega)=\min_{\\{(V_{0},V_{1}):j\in V_{1}\cup
V_{0}\\}}\eta(V_{0},V_{1}).$ (2.32)

The following shorthand notation is frequently used in this paper, which
stands for a generic multi-$\log(p)$ term that may vary from one occurrence to
another.

###### Definition 5

$L_{p}>0$ denotes a multi-$\log(p)$ term such that when $p\rightarrow\infty$,
for any $\delta>0$, $L_{p}p^{\delta}\rightarrow\infty$ and
$L_{p}p^{-\delta}\rightarrow 0$.

By (2.30) and Mills’ ratio (Wasserman and Roeder, 2009), a lower bound for the
“local risk” at $j$ is

$\displaystyle\qquad\qquad\sup_{\\{(V_{0},V_{1}):\;j\in V_{0}\cup
V_{1}\\}}\bigl{\\{}\inf_{t}\bigl{[}\epsilon_{p}^{|V_{0}|}\bar{\Phi}(t)+\epsilon_{p}^{|V_{1}|}\Phi\bigl{(}t-(\tau_{p}/\sigma)[\alpha^{*}(V_{0},V_{1})]^{1/2}\bigr{)}\bigr{]}\bigr{\\}}$
(2.33) $\displaystyle=$ $\displaystyle\sup_{\\{(V_{0},V_{1}):\;j\in V_{0}\cup
V_{1}\\}}\big{\\{}L_{p}\exp\bigl{(}-\eta(V_{0},V_{1})\cdot\log(p)\bigr{)}\bigr{\\}}=L_{p}\exp(-\rho_{j}^{*}(\vartheta,r,a,\Omega)\log(p)).$
(2.34)

We now aggregate such lower bounds for “local risk” for a global lower bound.
Since the “least favorable” configurations of $(V_{0},V_{1})$ for different
$j$ may overlap with each other, we need to consider a graph as follows.
Revisit the optimization problem in (2.31) and let

$(V_{0j}^{*},V_{1j}^{*})=\mathrm{argmin}_{\\{(V_{0},V_{1}):j\in V_{1}\cup
V_{0}\\}}\eta(V_{0},V_{1};\vartheta,r,a,\Omega).$ (2.35)

When there is a tie, pick the pair that appears first lexicographically.
Therefore, for any $1\leq j\leq p$, $V_{0j}^{*}\cup V_{1j}^{*}$ is uniquely
defined. In Lemma 22 of the appendix, we show that $|V_{0j}^{*}\cup
V_{1j}^{*}|\leq(\vartheta+r)^{2}/(2\vartheta r)$ for all $1\leq j\leq p$.

We now define a new graph, Graph of Least Favorable (GOLF), ${\cal
G}^{\diamond}=(V,E)$, where $V=\\{1,2,\ldots,p\\}$ and there is an edge
between $j$ and $k$ if and only if $(V_{0j}^{*}\cup V_{1j}^{*})$ and
$(V_{0k}^{*}\cup V_{1k}^{*})$ have non-empty intersections. Denote the maximum
degree of GOLF by $d_{p}({\cal G}^{\diamond})$.

###### Theorem 6

Fix $(\vartheta,\kappa)\in(0,1)^{2}$, $r>0$, and $a\geq 1$ such that
$\kappa>(1-\vartheta)$, and let ${\cal M}_{p}(c_{0},g)$ be as in (2.22).
Consider Model (1.1) where $\beta$ is modeled by $ARW(\vartheta,r,a,\mu)$ and
$X$ is modeled by $RD(\vartheta,\kappa,\Omega)$ and $\Omega\in{\cal
M}_{p}(c_{0},g)$ for sufficiently large $p$. Then as $p\rightarrow\infty$,
$\mathrm{Hamm}_{p}^{*}(\vartheta,\kappa,r,a,\Omega)\geq
L_{p}\bigl{[}d_{p}({\cal
G}^{\diamond})\bigr{]}^{-1}\sum_{j=1}^{p}p^{-\rho_{j}^{*}(\vartheta,r,a,\Omega)}$.

A similar claim holds for deterministic design models; the proof is similar so
we omit it.

###### Corollary 7

For deterministic design models, the parallel lower bound holds for the
minimax Hamming distance with $\Omega$ replaced by $G$ in the calculation of
$\rho_{j}^{*}(\vartheta,r,a,\Omega)$ and $d_{p}({\cal G}^{\diamond})$.

Remark. The lower bounds contain a factor of $\bigl{[}d_{p}({\cal
G}^{\diamond})\bigr{]}^{-1}$. In many cases including that considered in our
main theorem (Theorem 8), this factor is a multi-$\log(p)$ term so it does not
have a major effect. In some other cases, the factor $\bigl{[}d_{p}({\cal
G}^{\diamond})\bigr{]}^{-1}$ could be much smaller, say, when the GOSD has one
or a few hubs, the degrees of which grow algebraically fast as $p$ grows. In
these cases, the associated GOLF may (or may not) have large-degree hubs. As a
result, the lower bounds we derive could be very conservative, and can be
substantially improved if we treat the hubs, neighboring nodes of the hubs,
and other nodes separately. For the sake of space, we leave such discussion to
future work.

Remark. A similar lower bound holds if the condition
$\mu\in\Theta^{*}_{p}(\tau_{p},a)$ of ARW is replaced by
$\mu\in\Theta_{p}(\tau_{p})$. In (2.31), suppose we replace
$\Theta_{p}^{*}(\tau_{p},a)$ by $\Theta_{p}(\tau_{p})$, and the minimum is
achieved at
$(\theta^{(0)},\theta^{(1)})=(\theta_{*}^{(0)}(V_{0},V_{1};\Omega),\theta_{*}^{(1)}(V_{0},V_{1};\Omega))$.
Let $g=g(m_{0},\vartheta,r)$ be as in (2.21) and define

$a_{g}^{*}(\Omega)=\max_{\\{(V_{0},V_{1}):|V_{0}\cup V_{1}|\leq
g\\}}\\{\|\theta_{*}^{(0)}(V_{0},V_{1};\Omega)\|_{\infty},\|\theta_{*}^{(1)}(V_{0},V_{1};\Omega)\|_{\infty}\\}.$

By elementary calculus, it is seen that for $\Omega\in{\cal M}_{p}(c_{0},g)$,
there is a a constant $C=C(c_{0},g)$ such that $a_{g}^{*}(\Omega)\leq C$. If
additionally we assume

$a>a_{g}^{*}(\Omega),$ (2.36)

then $\alpha^{*}(V_{0},V_{1})=\alpha^{*}(V_{0},V_{1};\Omega,a)$,
$\eta(V_{0},V_{1};\Omega,a,\vartheta,r)$, and
$\rho_{j}^{*}(\vartheta,r,a,\Omega)$ do not depend on $a$. Especially, we can
derive an alternative formula for $\rho_{j}^{*}(\vartheta,r,a,\Omega)$; see
Lemma 18 for details.

When (2.36) holds, $\Theta_{p}^{*}(\tau_{p},a)$ is broad enough in the sense
that the least favorable configurations
$(V_{0},V_{1},\beta^{(0)},\beta^{(1)})$ for all $j$ satisfy
$\|\beta^{(i)}\|_{\infty}\leq a\tau_{p}$, $i=0,1$. Consequently, neither the
minimax rate nor GS needs to adapt to $a$. In Section 2.6, we assume (2.36)
holds; (2.36) is a mild condition for it only involves small-size sub-matrices
of $\Omega$.

###  2.6 Upper bound and optimality of Graphlet Screening

Fix constants $\gamma\in(0,1)$ and $A>0$. Let ${\cal M}_{p}(c_{0},g)$ be as in
(2.22). In this section, we further restrict $\Omega$ to the following set:

${\cal M}_{p}^{*}(\gamma,c_{0},g,A)=\Bigl{\\{}\Omega\in{\cal
M}_{p}(c_{0},g):\,\sum_{j=1}^{p}|\Omega(i,j)|^{\gamma}\leq A,\,1\leq i\leq
p\Bigr{\\}}.$ (2.37)

Note that any $\Omega\in{\cal M}_{p}^{*}(\gamma,c_{0},g,A)$ is sparse in the
sense that each row of $\Omega$ has relatively few large coordinates. The
sparsity of $\Omega$ implies the sparsity of the Gram matrix $G$, since small-
size sub-matrices of $G$ approximately equal to their counterparts of
$\Omega$.

In GS, when we regularize GOSD (see (1.6)), we set the threshold $\delta$ by

$\delta=\delta_{p}=1/\log(p).$ (2.38)

Such a choice for threshold is mainly for convenience, and can be replaced by
any term that tends to $0$ logarithmically fast as $p\rightarrow\infty$.

For any subsets $D$ and $F$ of $\\{1,2,\ldots,p\\}$, define
$\omega(D,F;\Omega)=\omega(D,F;\vartheta,r,a,\Omega,p)$ by

$\omega(D,F;\Omega)=\min_{\xi\in\mathbb{R}^{|D|},\min_{i\in D}|\xi_{i}|\geq
1}\Big{\\{}\xi^{\prime}\bigl{(}\Omega^{D,D}-\Omega^{D,F}(\Omega^{F,F})^{-1}\Omega^{F,D}\bigr{)}\xi\Big{\\}},$
(2.39)

We choose the tuning parameters in the $GS$-step in a way such that

$t(\hat{D},\hat{F})=2\sigma^{2}q(\hat{D},\hat{F})\log p,$ (2.40)

where $q=q(\hat{D},\hat{F})>0$ satisfies (for short,
$\omega=\omega({\hat{D}},{\hat{F}};\Omega)$)

$\left\\{\begin{array}[]{ll}\sqrt{q_{0}}\leq\sqrt{q}\leq\sqrt{\omega
r}-\sqrt{\frac{(\vartheta+\omega r)^{2}}{4\omega
r}-\frac{|\hat{D}|+1}{2}\vartheta},&\ \mbox{$|\hat{D}|$ is odd \& $\omega
r/\vartheta>|\hat{D}|+(|\hat{D}|^{2}-1)^{1/2}$},\\\
\sqrt{q_{0}}\leq\sqrt{q}\leq\sqrt{\omega r}-\sqrt{\frac{1}{4}\omega
r-\frac{1}{2}|\hat{D}|\vartheta},&\ \mbox{$|\hat{D}|$ is even \& $\omega
r/\vartheta\geq 2|\hat{D}|$},\\\ \mbox{$q$ is a constant such that $q\geq
q_{0}$},&\ \mbox{otherwise}.\end{array}\right.$ (2.41)

We set the $GC$-step tuning parameters by

$u^{gs}=\sigma\sqrt{2\vartheta\log p},\qquad
v^{gs}=\tau_{p}=\sigma\sqrt{2r\log p}.$ (2.42)

The main theorem of this paper is the following theorem.

###### Theorem 8

Fix $m_{0}\geq 1$, $(\vartheta,\gamma,\kappa)\in(0,1)^{3}$, $r>0$, $c_{0}>0$,
$g>0$, $a>1$, $A>0$ such that $\kappa>1-\vartheta$ and (2.21) is satisfied.
Consider Model (1.1) where $\beta$ is modeled by $ARW(\vartheta,r,a,\mu)$, $X$
is modeled by $RD(\vartheta,\kappa,\Omega)$, and where $\Omega\in{\cal
M}_{p}^{*}(\gamma,c_{0},g,A)$ and $a>a_{g}^{*}(\Omega)$ for sufficiently large
$p$. Let $\hat{\beta}^{gs}=\hat{\beta}^{gs}(Y;\delta,{\cal
Q},u^{gs},v^{gs},X,p,n)$ be the Graphlet Screening procedure defined as in
Section 2.1, where the tuning parameters $(\delta,{\cal Q},u^{gs},v^{gs})$ are
set as in (2.38)-(2.42). Then as $p\rightarrow\infty$,
$\sup_{\mu\in\Theta_{p}^{*}(\tau_{p},a)}H_{p}(\hat{\beta}^{gs};\epsilon_{p},n_{p},\mu,\Omega)\leq
L_{p}\Big{[}p^{1-(m_{0}+1)\vartheta}+\sum_{j=1}^{p}p^{-\rho_{j}^{*}(\vartheta,r,a,\Omega)}\Big{]}+o(1)$.

Note that $\rho_{j}^{*}=\rho_{j}^{*}(\vartheta,r,a,\Omega)$ does not depend on
$a$. Also, note that in the most interesting range,
$\sum_{j=1}^{p}p^{-\rho_{j}^{*}}\gg 1$. So if we choose $m_{0}$ properly large
(e.g., $(m_{0}+1)\vartheta>1$), then
$\sup_{\mu\in\Theta_{p}^{*}(\tau_{p},a)}H_{p}(\hat{\beta}^{gs};\epsilon_{p},n_{p},\mu,\Omega)\leq
L_{p}\sum_{j=1}^{p}p^{-\rho_{j}^{*}(\vartheta,r,a,\Omega)}$. Together with
Theorem 6, this says that GS achieves the optimal rate of convergence,
adaptively to all $\Omega$ in ${\cal M}_{p}^{*}(\gamma,c_{0},g,A)$ and
$\beta\in\Theta_{p}^{*}(\tau_{p},a)$. We call this property optimal
adaptivity. Note that since the diagonals of $\Omega$ are scaled to $1$
approximately, $\kappa\equiv\log(n_{p})/\log(p)$ does not have a major
influence over the convergence rate, as long as (2.25) holds.

Remark. Theorem 8 addresses the case where (2.36) holds so
$a>a_{g}^{*}(\Omega)$. We now briefly discuss the case where
$a<a_{g}^{*}(\Omega)$. In this case, the set $\Theta_{p}^{*}(\tau_{p},a)$
becomes sufficiently narrow and $a$ starts to have some influence over the
optimal rate of convergence, at least for some choices of $(\vartheta,r)$. To
reflect the role of $a$, we modify GS as follows: (a) in the $GC$-step (2.13),
limit $\xi$ to the class where either $\xi_{i}=0$ or
$\tau_{p}\leq|\xi_{i}|\leq a\tau_{p}$, and (b) in the $GS$-step, replacing the
$\chi^{2}$-screening by the likelihood based screening procedure; that is,
when we screen ${\cal I}_{0}=\hat{D}\cup\hat{F}$, we accept nodes in $\hat{D}$
only when $h(\hat{F})>h({\cal I}_{0})$, where for any subset
$D\subset\\{1,2,\ldots,p\\}$,
$h(D)=\min\bigl{\\{}\frac{1}{2}\|P^{D}(Y-X^{\otimes,D}\xi)\|^{2}+\vartheta\sigma^{2}\log(p)|D|\big{\\}}$,
where the minimum is computed over all $|D|\times 1$ vectors $\xi$ whose
nonzero elements all have magnitudes in $[\tau_{p},a\tau_{p}]$. From a
practical point of view, this modified procedure depends more on the
underlying parameters and is harder to implement than is GS. However, this is
the price we need to pay when $a$ is small. Since we are primarily interested
in the case of relatively larger $a$ (so that $a>a_{g}^{*}(\Omega)$ holds), we
skip further discussion along this line.

###  2.7 Phase diagram and examples where $\rho_{j}^{*}(\vartheta,r,a,\Omega)$
have simple forms

In general, the exponents $\rho_{j}^{*}(\vartheta,r,a,\Omega)$ may depend on
$\Omega$ in a complicated way. Still, from time to time, one may want to find
a simple expression for $\rho_{j}^{*}(\vartheta,r,a,\Omega)$. It turns out
that in a wide class of situations, simple forms for
$\rho_{j}^{*}(\vartheta,r,a,\Omega)$ are possible. The surprise is that, in
many examples, $\rho_{j}^{*}(\vartheta,r,a,\Omega)$ depends more on the trade-
off between the parameters $\vartheta$ and $r$ (calibrating the signal
sparsity and signal strength, respectively), rather than on the large
coordinates of $\Omega$.

We begin with the following theorem, which is proved in Ji and Jin (2011,
Theorem 1.1).

###### Theorem 9

Fix $(\vartheta,\kappa)\in(0,1)$, $r>0$, and $a>1$ such that
$\kappa>(1-\vartheta)$. Consider Model (1.1) where $\beta$ is modeled by
$ARW(\vartheta,r,a,\mu)$ and $X$ is modeled by $RD(\vartheta,\kappa,\Omega)$.
Then as $p\rightarrow\infty$,

$\frac{\mathrm{Hamm}_{p}^{*}(\vartheta,\kappa,r,a,\Omega)}{p^{1-\vartheta}}\gtrsim\left\\{\begin{array}[]{ll}1,&\qquad
0<r<\vartheta,\\\ L_{p}p^{-(r-\vartheta)^{2}/(4r)},&\qquad
r>\vartheta.\end{array}\right.$

Note that $p^{1-\vartheta}$ is approximately the number of signals. Therefore,
when $r<\vartheta$, the number of selection errors can not get substantially
smaller than the number of signals. This is the most difficult case where no
variable selection method can be successful.

In this section, we focus on the case $r>\vartheta$, so that successful
variable selection is possible. In this case, Theorem 9 says that a universal
lower bound for the Hamming distance is $L_{p}p^{1-(\vartheta+r)^{2}/(4r)}$.
An interesting question is, to what extend, this lower bound is tight.

Recall that $\lambda_{k}^{*}(\Omega)$ denotes the minimum of smallest
eigenvalues across all $k\times k$ principle submatrices of $\Omega$, as
defined in (2.20). The following corollaries are proved in Section 6.

###### Corollary 10

Suppose the conditions of Theorem 8 hold, and that additionally,
$1<r/\vartheta<3+2\sqrt{2}\approx 5.828$, and $|\Omega(i,j)|\leq
4\sqrt{2}-5\approx 0.6569$ for all $1\leq i,j\leq p$, $i\neq j$. Then as
$p\rightarrow\infty$,
$\mathrm{Hamm}_{p}^{*}(\vartheta,\kappa,r,a,\Omega)=L_{p}p^{1-(\vartheta+r)^{2}/(4r)}$.

###### Corollary 11

Suppose the conditions of Theorem 8 hold. Also, suppose that
$1<r/\vartheta<5+2\sqrt{6}\approx 9.898$, and that
$\lambda_{3}^{*}(\Omega)\geq 2(5-2\sqrt{6})\approx 0.2021$,
$\lambda_{4}^{*}(\Omega)\geq 5-2\sqrt{6}\approx 0.1011$, and
$|\Omega(i,j)|\leq 8\sqrt{6}-19\approx 0.5959$ for all $1\leq i,j\leq p$,
$i\neq j$. Then as $p\rightarrow\infty$,
$\mathrm{Hamm}_{p}^{*}(\vartheta,\kappa,r,a,\Omega)=L_{p}p^{1-(\vartheta+r)^{2}/(4r)}$.

In these corollaries, the conditions on $\Omega$ are rather relaxed. Somewhat
surprisingly, the off-diagonals of $\Omega$ do not necessarily have a major
influence on the optimal rate of convergence, as one might have expected.

Figure 1: Phase diagram for $\Omega=I_{p}$ (left), for $\Omega$ satisfying
conditions of Corollary 10 (middle), and for $\Omega$ satisfying conditions of
Corollary 11 (right). Red line: $r=\vartheta$. Solid red curve:
$r=\rho(\vartheta,\Omega)$. In each of the last two panels, the blue line
intersects with the red curve at $(\vartheta,r)=(1/2,[3+2\sqrt{2}]/2)$
(middle) and $(\vartheta,r)=(1/3,[5+2\sqrt{6}]/3)$ (right), which splits the
red solid curve into two parts; the part to the left is illustrative for it
depends on $\Omega$ in a complicated way; the part to the right, together with
the dashed red curve, represent $r=(1+\sqrt{1-\vartheta})^{2}$ (in the left
panel, this is illustrated by the red curve).

Note also that by Theorem 8, under the condition of either Corollaries 10 or
Corollary 11, GS achieves the optimal rate in that

$\sup_{\mu\in\Theta_{p}^{*}(\tau_{p},a)}H_{p}(\hat{\beta}^{gs};\epsilon_{p},n_{p},\mu,\Omega)\leq
L_{p}p^{1-(\vartheta+r)^{2}/(4r)}.$ (2.43)

Together, Theorem 9, Corollaries 10-11, and (2.43) have an interesting
implication on the so-called phase diagram. Call the two-dimensional parameter
space $\\{(\vartheta,r):0<\vartheta<1,r>0\\}$ the phase space. There are two
curves $r=\vartheta$ and $r=\rho(\vartheta,\Omega)$ (the latter can be thought
of as the solution of
$\sum_{j=1}^{p}p^{-\rho^{*}_{j}(\vartheta,r,a,\Omega)}=1$; recall that
$\rho_{j}^{*}(\vartheta,r,a,\Omega)$ does not depend on $a$) that partition
the whole phase space into three different regions:

  * •

Region of No Recovery. $\\{(\vartheta,r):0<r<\vartheta,0<\vartheta<1\\}$. In
this region, as $p\rightarrow\infty$, for any $\Omega$ and any procedures, the
minimax Hamming error equals approximately to the total expected number of
signals. This is the most difficult region, in which no procedure can be
successful in the minimax sense.

  * •

Region of Almost Full Recovery.
$\\{(\vartheta,r):\vartheta<r<\rho(\vartheta,\Omega)\\}$. In this region, as
$p\rightarrow\infty$, the minimax Hamming distance satisfies
$1\ll\mathrm{Hamm}_{p}^{*}(\vartheta,\kappa,r,a,\Omega)\ll p^{1-\vartheta}$,
and it is possible to recover most of the signals, but it is impossible to
recover all of them.

  * •

Region of Exact Recovery. In this region, as $p\rightarrow\infty$, the minimax
Hamming distance $\mathrm{Hamm}_{p}^{*}(\vartheta,\kappa,r,a,\Omega)=o(1)$,
and it is possible to exactly recover all signals with overwhelming
probability.

In general, the function $\rho(\vartheta,\Omega)$ depends on $\Omega$ in a
complicated way. However, by Theorem 9 and Corollaries 10-11, we have the
following conclusions. First, for all $\Omega$ and $a>1$,
$\rho(\vartheta,\Omega)\geq(1+\sqrt{1-\vartheta})^{2}$ for all
$0<\vartheta<1$. Second, in the simplest case where $\Omega=I_{p}$,
$\mathrm{Hamm}_{p}^{*}(\vartheta,\kappa,r,a,\Omega)=L_{p}p^{1-(\vartheta+r)^{2}/(4r)}$,
and $\rho(\vartheta,\Omega)=(1+\sqrt{1-\vartheta})^{2}$ for all
$0<\vartheta<1$. Third, under the conditions of Corollary 10,
$\rho(\vartheta,\Omega)=(1+\sqrt{1-\vartheta})^{2}$ if $1/2<\vartheta<1$.
Last, under the conditions of Corollary 11,
$\rho(\vartheta,\Omega)=(1+\sqrt{1-\vartheta})^{2}$ if $1/3<\vartheta<1$. The
phase diagram for the last three cases are illustrated in Figure 1. The blue
lines are $r/\vartheta=3+2\sqrt{2}$ (middle) and $r/\vartheta=5+2\sqrt{6}$
(right).

Corollaries 10-11 can be extended to more general situations, where
$r/\vartheta$ may get arbitrary large, but consequently, we need stronger
conditions on $\Omega$. Towards this end, we note that for any $(\vartheta,r)$
such that $r>\vartheta$, we can find a unique integer $N=N(\vartheta,r)$ such
that $2N-1\leq(\vartheta/r+r/\vartheta)/2<2N+1$. Suppose that for any $2\leq
k\leq 2N-1$,

$\lambda_{k}^{*}(\Omega)\geq\max_{\\{(k+1)/2\leq
j\leq\min\\{k,N\\}\\}}\Big{\\{}\frac{(r/\vartheta+\vartheta/r)/2-2j+2+\sqrt{[(r/\vartheta+\vartheta/r)/2-2j+2]^{2}-1}}{(2k-2j+1)(r/\vartheta)}\Big{\\}},$
(2.44)

and that for any $2\leq k\leq 2N$,

$\lambda_{k}^{*}(\Omega)\geq\max_{\\{k/2\leq
j\leq\min\\{k-1,N\\}\\}}\Big{\\{}\frac{(r/\vartheta+\vartheta/r)/2+1-2j}{(k-j)(r/\vartheta)}\Big{\\}}.$
(2.45)

Then we have the following corollary.

###### Corollary 12

Suppose the conditions in Theorem 8 and that in (2.44)-(2.45) hold. Then as
$p\rightarrow\infty$,
$\mathrm{Hamm}_{p}^{*}(\vartheta,\kappa,r,a,\Omega)=L_{p}p^{1-(\vartheta+r)^{2}/(4r)}$.

The right hand sides of (2.44)-(2.45) decrease with $(r/\vartheta)$. For a
constant $s_{0}>1$, (2.44)-(2.45) hold for all $1<r/\vartheta\leq s_{0}$ as
long as they hold for $r/\vartheta=s_{0}$. Hence Corollary 12 implies a
similar partition of the phase diagram as do Corollaries 10-11.

Remark. Phase diagram can be viewed as a new criterion for assessing the
optimality, which is especially appropriate for rare and weak signals. The
phase diagram is a partition of the phase space
$\\{(\vartheta,r):0<\vartheta<1,r>0\\}$ into different regions where
statistical inferences are distinctly different. In general, a phase diagram
has the following four regions:

  * •

An “exact recovery” region corresponding to the “rare and strong” regime in
which high probability of completely correct variable selection is feasible.

  * •

An “almost full recovery” region as a part of the “rare and weak” regime in
which completely correct variable selection is not achievable with high
probability but variable selection is still feasible in the sense that with
high probability, the number of incorrectly selected variables is a small
fraction of the total number of signals.

  * •

A “detectable” region in which variable selection is infeasible but the
detection of the existence of a signal (somewhere) is feasible (e.g., by the
Higher Criticism method).

  * •

An “undetectable” region where signals are so rare and weak that nothing can
be sensibly done.

In the sparse signal detection (Donoho and Jin, 2004) and classification (Jin,
2009) problems, the main interest is to find the detectable region, so that
the exact recovery and almost full recovery regions were lumped into a single
“estimable” region (e.g., Donoho and Jin (2004, Figure 1)). For variable
selection, the main interest is to find the boundaries of the almost full
discovery region so that the detectable and non-detectable regions are lumped
into a single “no recovery” region as in Ji and Jin (2011) and Figure 1 of
this paper.

Variable selection in the “almost full recovery” region is a new and
challenging problem. It was studied in Ji and Jin (2011) when the effect of
signal cancellation is negligible, but the hardest part of the problem was
unsolved in Ji and Jin (2011). This paper (the second in this area) deals with
the important issue of signal cancellation, in hopes of gaining a much deeper
insight on variable selection in much broader context.

###  2.8 Non-optimaility of subset selection and the lasso

Subset selection (also called the $L^{0}$-penalization method) is a well-known
method for variable selection, which selects variables by minimizing the
following functional:

$\frac{1}{2}\|Y-X\beta\|^{2}+\frac{1}{2}(\lambda_{ss})^{2}\|\beta\|_{0},$
(2.46)

where $\|\beta\|_{q}$ denotes the $L^{q}$-norm, $q\geq 0$, and
$\lambda_{ss}>0$ is a tuning parameter. The AIC, BIC, and RIC are methods of
this type (Akaike, 1974; Schwarz, 1978; Foster and George, 1994). Subset
selection is believed to have good “theoretic property”, but the main drawback
of this method is that it is computationally NP hard. To overcome the
computational challenge, many relaxation methods are proposed, including but
are not limited to the lasso (Chen et al, 1998; Tibshirani, 1996), SCAD (Fan
and Li, 2001), MC+ (Zhang, 2010), and Dantzig selector (Candes and Tao, 2007).
Take the lasso for example. The method selects variables by minimizing

$\frac{1}{2}\|Y-X\beta\|^{2}+\lambda_{lasso}\|\beta\|_{1},$ (2.47)

where the $L^{0}$-penalization is replaced by the $L^{1}$-penalization, so the
functional is convex and the optimization problem is solvable in polynomial
time under proper conditions.

Somewhat surprisingly, subset selection is generally rate non-optimal in terms
of selection errors. This sub-optimality of subset selection is due to its
lack of flexibility in adapting to the “local” graphic structure of the design
variables. Similarly, other global relaxation methods are sub-optimal as well,
as the subset selection is the “idol” these methods try to mimic. To save
space, we only discuss subset selection and the lasso, but a similar
conclusion can be drawn for SCAD, MC+, and Dantzig selector.

For mathematical simplicity, we illustrate the point with an idealized
regression model where the Gram matrix $G=X^{\prime}X$ is diagonal block-wise
and has $2\times 2$ blocks

$G(i,j)=1\\{i=j\\}+h_{0}\cdot 1\\{\mbox{$|j-i|=1$, $max(i,j)$ is
even}\\},\;\;|h_{0}|<1,\;1\leq i,j\leq p.$ (2.48)

Using an idealized model is mostly for technical convenience, but the non-
optimality of subset selection or the lasso holds much more broadly than what
is considered here. On the other hand, using a simple model is sufficient
here: if a procedure is non-optimal in an idealized case, we can not expect it
to be optimal in a more general context.

At the same time, we continue to model $\beta$ with the Asymptotic Rare and
Weak model ARW$(\vartheta,r,a,\mu)$, but where we relax the assumption of
$\mu\in\Theta_{p}^{*}(\tau_{p},a)$ to that of $\mu\in\Theta_{p}(\tau_{p})$ so
that the strength of each signal $\geq\tau_{p}$ (but there is no upper bound
on the strength). Consider a variable selection procedure
$\hat{\beta}^{\star}$, where $\star=gs,ss,lasso$, representing GS, subset
selection, and the lasso (where the tuning parameters for each method are
ideally set; for the worst-case risk considered below, the ideal tuning
parameters depend on $(\vartheta,r,p,h_{0})$ but do not depend on $\mu$).
Since the index groups $\\{2j-1,2j\\}$ are exchangeable in (2.48) and the ARW
models, the Hamming error of $\beta^{\star}$ in its worst case scenario has
the form of
$\sup_{\\{\mu\in\Theta_{p}(\tau_{p})\\}}H_{p}(\hat{\beta}^{\star};\epsilon_{p},\mu,G)=L_{p}p^{1-\rho_{\star}(\vartheta,r,h_{0})}$.

We now study $\rho_{\star}(\vartheta,r,h_{0})$. Towards this end, we first
introduce
$\rho_{lasso}^{(3)}(\vartheta,r,h_{0})=\bigl{\\{}(2|h_{0}|)^{-1}[(1-h_{0}^{2})\sqrt{r}-\sqrt{(1-h_{0}^{2})(1-|h_{0}|)^{2}r-4|h_{0}|(1-|h_{0}|)\vartheta}]\bigr{\\}}^{2}$
and
$\rho_{lasso}^{(4)}(\vartheta,r,h_{0})=\vartheta+\frac{(1-|h_{0}|)^{3}(1+|h_{0}|)}{16h_{0}^{2}}\bigl{[}(1+|h_{0}|)\sqrt{r}-\sqrt{(1-|h_{0}|)^{2}r-4|h_{0}|\vartheta/(1-h_{0}^{2})}\bigr{]}^{2}$.
We then let

$\rho_{ss}^{(1)}(\vartheta,r,h_{0})=\left\\{\begin{array}[]{ll}2\vartheta,&\qquad\quad
r/\vartheta\leq 2/(1-h_{0}^{2})\\\
\,[2\vartheta+(1-h_{0}^{2})r]^{2}/[4(1-h_{0}^{2})r],&\qquad\quad
r/\vartheta>2/(1-h_{0}^{2})\end{array},\right.$
$\rho_{ss}^{(2)}(\vartheta,r,h_{0})=\left\\{\begin{array}[]{ll}2\vartheta,&\quad
r/\vartheta\leq 2/(1-|h_{0}|)\\\
2[\sqrt{2(1-|h_{0}|)r}-\sqrt{(1-|h_{0}|)r-\vartheta}]^{2},&\quad
r/\vartheta>2/(1-|h_{0}|)\end{array},\right.$
$\rho_{lasso}^{(1)}(\vartheta,r,h_{0})=\left\\{\begin{array}[]{ll}2\vartheta,&\qquad\qquad
r/\vartheta\leq 2/(1-|h_{0}|)^{2}\\\
\rho_{lasso}^{(3)}(\vartheta,r,h_{0}),&\qquad\qquad
r/\vartheta>2/(1-|h_{0}|)^{2}\end{array},\right.$

and

$\rho_{lasso}^{(2)}(\vartheta,r,h_{0})=\left\\{\begin{array}[]{ll}2\vartheta,&\qquad\qquad
r/\vartheta\leq(1+|h_{0}|)/(1-|h_{0}|)^{3}\\\
\rho_{lasso}^{(4)}(\vartheta,r,h_{0}),&\qquad\qquad
r/\vartheta>(1+|h_{0}|)/(1-|h_{0}|)^{3}\end{array}.\right.$

The following theorem is proved in Section 6.

###### Theorem 13

Fix $\vartheta\in(0,1)$ and $r>0$ such that $r>\vartheta$. Consider Model
(1.1) where $\beta$ is modeled by $ARW(\vartheta,r,a,\mu)$ and $X$ satisfies
(2.48). For GS, we set the tuning parameters $(\delta,m_{0})=(0,2)$, and set
$({\cal Q},u^{gs},v^{gs})$ as in (2.40)-(2.42). For subset selection as in
(2.46) and the lasso as in (2.47), we set their tuning parameters ideally
given that $(\vartheta,r)$ are known. Then as $p\rightarrow\infty$,

$\rho_{gs}(\vartheta,r,h_{0})=\min\bigl{\\{}\frac{(\vartheta+r)^{2}}{4r},\,\vartheta+\frac{(1-|h_{0}|)}{2}r,\,2\vartheta+\frac{\\{[(1-h_{0}^{2})r-\vartheta]_{+}\\}^{2}}{4(1-h_{0}^{2})r}\bigr{\\}},$
(2.49)
$\rho_{ss}(\vartheta,r,h_{0})=\min\bigl{\\{}\frac{(\vartheta+r)^{2}}{4r},\,\vartheta+\frac{(1-|h_{0}|)}{2}r,\,\rho_{ss}^{(1)}(\vartheta,r,h_{0}),\rho_{ss}^{(2)}(\vartheta,r,h_{0})\bigr{\\}},$
(2.50)

and

$\rho_{lasso}(\vartheta,r,h_{0})=\min\\{\frac{(\vartheta+r)^{2}}{4r},\,\vartheta+\frac{(1-|h_{0}|)r}{2(1+\sqrt{1-h_{0}^{2}})},\,\rho_{lasso}^{(1)}(\vartheta,r,h_{0}),\,\rho_{lasso}^{(2)}(\vartheta,r,h_{0})\bigr{\\}}.$
(2.51)

It can be shown that
$\rho_{gs}(\vartheta,r,h_{0})\geq\rho_{ss}(\vartheta,r,h_{0})\geq\rho_{lasso}(\vartheta,r,h_{0})$,
where depending on the choices of $(\vartheta,r,h_{0})$, we may have equality
or strict inequality (note that a larger exponent means a better error rate).
This fits well with our expectation, where as far as the convergence rate is
concerned, GS is optimal for all $(\vartheta,r,h_{0})$, so it outperforms the
subset selection, which in turn outperforms the lasso. Table 2 summarizes the
exponents for some representative $(\vartheta,r,h_{0})$. It is seen that
differences between these exponents become increasingly prominent when $h_{0}$
increase and $\vartheta$ decrease.

$\vartheta/r/h_{0}$ | .1/11/.8 | .3/9/.8 | .5/4/.8 | .1/4/.4 | .3/4/.4 | .5/4/.4 | .1/3/.2 | .3/3/.2
---|---|---|---|---|---|---|---|---
$\star=gs$ | 1.1406 | 1.2000 | 0.9000 | 0.9907 | 1.1556 | 1.2656 | 0.8008 | 0.9075
$\star=ss$ | 0.8409 | 0.9047 | 0.9000 | 0.9093 | 1.1003 | 1.2655 | 0.8007 | 0.9075
$\star=lasso$ | 0.2000 | 0.6000 | 0.7500 | 0.4342 | 0.7121 | 1.0218 | 0.6021 | 0.8919
Table 2: The exponents $\rho_{\star}(\vartheta,r,h_{0})$ in Theorem 13, where
$\star=gs,ss,lasso$.

As in Section 2.7, each of these methods has a phase diagram plotted in Figure
2, where the phase space partitions into three regions: Region of Exact
Recovery, Region of Almost Full Recovery, and Region of No Recovery.
Interestingly, the separating boundary for the last two regions are the same
for three methods, which is the line $r=\vartheta$. The boundary that
separates the first two regions, however, vary significantly for different
methods. For any $h_{0}\in(-1,1)$ and $\star=gs,ss,lasso$, the equation for
this boundary can be obtained by setting $\rho_{\star}(\vartheta,r,h_{0})=1$
(the calculations are elementary so we omit them). Note that the lower the
boundary is, the better the method is, and that the boundary corresponding to
the lasso is discontinuous at $\vartheta=1/2$. In the non-optimal region of
either subset selection or the lasso, the Hamming errors of the procedure are
much smaller than $p\epsilon_{p}$, so the procedure gives “almost full
recovery”; however, the rate of Hamming errors is slower than that of the
optimal procedure, so subset selection or the lasso is non-optimal in such
regions.

Subset selection and the lasso are rate non-optimal for they are so-called
one-step or non-adaptive methods (Ji and Jin, 2011), which use only one tuning
parameter, and which do not adapt to the local graphic structure. The non-
optimality can be best illustrated with the diagonal block-wise model
presented here, where each block is a $2\times 2$ matrix. Correspondingly, we
can partition the vector $\beta$ into many size $2$ blocks, each of which is
of the following three types (i) those have no signal, (ii) those have exactly
one signal, and (iii) those have two signals. Take the subset selection for
example. To best separate (i) from (ii), we need to set the tuning parameter
ideally. But such a tuning parameter may not be the “best” for separating (i)
from (iii). This explains the non-optimality of subset selection.

Seemingly, more complicated penalization methods that use multiple tuning
parameters may have better performance than the subset selection and the
lasso. However, it remains open how to design such extensions to achieve the
optimal rate for general cases. To save space, we leave the study along this
line to the future.

Figure 2: Phase diagrams for GS (top left), subset selection (top right), and
the lasso (bottom; zoom-in on the left and zoom-out on the right), where
$h_{0}=0.5$.

###  2.9 Summary

We propose GS as a new approach to variable selection. The key methodological
innovation is to use the GOSD to guide the multivariate screening. While a
brute-force $m$-variate screening has a computation cost of $O(p^{m}+np)$, GS
only has a computation cost of $L_{p}np$ (excluding the overhead of obtaining
the GOSD), by utilizing graph sparsity. Note that when the design matrix $G$
is approximately banded, say, all its large entries are confined to a diagonal
band with bandwidth $\leq K$, the overhead of GS can be reduced to $O(npK)$.
One such example is in Genome-Wide Association Study (GWAS), where $G$ is the
empirical Linkage Disequilibrium (LD) matrix, and $K$ can be as small as a few
tens. We remark that the lasso has a computational complexity of $O(npk)$,
where $k$, dominated by the number steps requiring re-evaluation of the
correlation between design vectors and updated residuals, could be smaller
than the $L_{p}$ term for GS (Wang et al, 2013).

We use asymptotic minimaxity of the Hamming distance as the criterion for
assessing optimality. Compared with existing literature on variable selection
where we use the oracle property or probability of exact support recovery to
assess optimality, our approach is mathematically more demanding, yet
scientifically more relevant in the rare/weak paradigm.

We have proved that GS achieves the optimal rate of convergence of Hamming
errors, especially when signals are rare and weak, provided that the Gram
matrix is sparse. Subset selection and the lasso are not rate optimal, even
with very simple Gram matrix $G$ and even when the tuning parameters are
ideally set. The sub-optimality of these methods is due to that they do not
take advantage of the ‘local’ graphical structure as GS does.

GS has three key tuning parameters: $q$ for the threshold level
$t(\hat{D},\hat{F})=2\sigma^{2}q\log p$ in the $GS$-step, and
$(u^{gs},v^{gs})=(\sigma\sqrt{2\vartheta\log p},\sigma\sqrt{2r\log p})$ in the
$GC$-step. While the choice of $q$ is reasonably flexible and a sufficiently
small fixed $q>0$ is usually adequate, the choice of $u^{gs}$ and $v^{gs}$ are
more directly tied to the signal sparsity and signal strength. Adaptive choice
of these tuning parameters is a challenging direction of further research. One
of our ideas to be developed in this direction is a subsampling scheme similar
to the Stability Selection (Meinsausen and Buhlmann, 2010). On the other hand,
as shown in our numeric results in Section 4, the performance of GS is
relatively insensitive to mis-specification of $(\epsilon_{p},\tau_{p})$; see
details therein.

##  3 Properties of Graphlet Screening, proof of Theorem 8

GS attributes the success to two important properties: the Sure Screening
property and the Separable After Screening (SAS) property.

The Sure Screening property means that in the $m_{0}$-stage $\chi^{2}$
screening, by picking an appropriate threshold, the set ${\cal U}_{p}^{*}$
(which is the set of retained indices after the GS-step) contains all but a
small fraction of true signals. Asymptotically, this fraction is comparably
smaller than the minimax Hamming errors, and so negligible. The SAS property
means that except for a negligible probability, as a subgraph of the GOSD,
${\cal U}_{p}^{*}$ decomposes into many disconnected components of the GOSD,
where the size of each component does not exceed a fixed integer. These two
properties ensure that the original regression problem reduces to many small-
size regression problems, and thus pave the way for the $GC$-step.

Below, we explain these ideas in detail, and conclude the section by the proof
of Theorem 8. Since the only place we need the knowledge of $\sigma$ is in
setting the tuning parameters, so without loss of generality, we assume
$\sigma=1$ throughout this section.

First, we discuss the $GS$-step. For short, write
$\hat{\beta}=\hat{\beta}^{gs}(Y;\delta,{\cal Q},u^{gs},v^{gs},X,p,n)$
throughout this section. We first discuss the computation cost of the
$GS$-step. As in Theorem 8, we take the threshold $\delta$ in ${\cal
G}^{*,\delta}$ to be $\delta=\delta_{p}=1/\log(p)$. The proof of the following
lemma is similar to that of Ji and Jin (2011, Lemma 2.2), so we omit it.

###### Lemma 14

Suppose the conditions of Theorem 8 hold, where we recall $\delta=1/\log(p)$,
and $\Omega^{*,\delta}$ is defined as in (1.6). As $p\rightarrow\infty$, with
probability $1-o(1/p^{2})$, $\|\Omega-\Omega^{*,\delta}\|_{\infty}\leq
C(\log(p))^{-(1-\gamma)}$, and ${\cal G}^{*,\delta}$ is $K$-sparse, where
$K\leq C(\log(p))^{1/\gamma}$.

Combining Lemma 14 and Frieze and Molloy (1999), it follows that with
probability $1-o(1/p^{2})$, ${\cal G}^{*,\delta}$ has at most
$p(Ce(\log(p))^{1/\gamma})^{m_{0}}$ connected subgraphs of size $\leq m_{0}$.
Note that the second factor is at most logarithmically large, so the
computation cost in the $GS$-step is at most $L_{p}p$ flops.

Consider the performance of the $GS$-step. The goal of this step is two-fold:
on one hand, it tries to retain as many signals as possible during the
screening; on the other hand, it tries to minimize the computation cost of the
$GC$-step by controlling the maximum size of all components of ${\cal
U}_{p}^{*}$. The key in the $GS$-step is to set the collection of thresholds
${\cal Q}$. The tradeoff is that, setting the thresholds too high may miss too
many signals during the screening, and setting the threshold too low may
increase the maximum size of the components in ${\cal U}_{p}^{*}$, and so
increase the computational burden of the $GC$-step. The following lemma
characterizes the Sure Screening property of GS, and is proved in Section 6.

###### Lemma 15

(Sure Screening). Suppose the settings and conditions are as in Theorem 8. In
the $m_{0}$-stage $\chi^{2}$ screening of the GS-step, if we set the
thresholds $t(\hat{D},\hat{F})$ as in (2.40), then as $p\rightarrow\infty$,
for any $\Omega\in{\cal M}_{p}^{*}(\gamma,c_{0},g,A)$,
$\sum_{j=1}^{p}P(\beta_{j}\neq 0,j\notin{\cal U}_{p}^{*})\leq
L_{p}[p^{1-(m_{0}+1)\vartheta}+\sum_{j=1}^{p}p^{-\rho_{j}^{*}(\vartheta,r,a,\Omega)}]+o(1)$.

Next, we formally state the SAS property. Viewing it as a subgraph of ${\cal
G}^{*,\delta}$, ${\cal U}_{p}^{*}$ decomposes into many disconnected
components ${\cal I}^{(k)}$, $1\leq k\leq N$, where $N$ is an integer that may
depend on the data.

###### Lemma 16

(SAS). Suppose the settings and conditions are as in Theorem 8. In the
$m_{0}$-stage $\chi^{2}$ screening in the GS-step, suppose we set the
thresholds $t(\hat{D},\hat{F})$ as in (2.40) such that $q(\hat{D},\hat{F})\geq
q_{0}$ for some constant $q_{0}=q_{0}(\vartheta,r)>0$. As
$p\rightarrow\infty$, under the conditions of Theorem 8, for any
$\Omega\in{\cal M}_{p}^{*}(\gamma,c_{0},g,A)$, there is a constant
$\ell_{0}=\ell_{0}(\vartheta,r,\kappa,\gamma,A,c_{0},g)>0$ such that with
probability at least $1-o(1/p)$, $|{\cal I}^{(k)}|\leq\ell_{0}$, $1\leq k\leq
N$.

We remark that a more convenient way of picking $q$ is to let

$\left\\{\begin{array}[]{ll}q_{0}\leq q\leq(\frac{\omega r+\vartheta}{2\omega
r})^{2}\omega r,&\qquad\mbox{$|\hat{D}|$ is odd \& $\omega
r/\vartheta>|\hat{D}|+(|\hat{D}|^{2}-1)^{1/2}$},\\\ q_{0}\leq
q\leq\frac{1}{4}\omega r,&\qquad\mbox{$|\hat{D}|$ is even \& $\omega
r/\vartheta\geq 2|\hat{D}|$},\end{array}\right.$ (3.52)

and let $q$ be any other number otherwise, with which both lemmas continue to
hold with this choice of $q$. Here, for short,
$\omega=\omega(\hat{D},\hat{F};\Omega)$. Note that numerically this choice is
comparably more conservative.

Together, the above two lemmas say that the $GS$-step makes only negligible
false non-discoveries, and decomposes ${\cal U}_{p}^{*}$ into many
disconnected components, each has a size not exceeding a fixed integer. As a
result, the computation cost of the following $GC$-step is moderate, at least
in theory.

We now discuss the $GC$-step. The key to understanding the $GC$-step is that
the original regression problem reduces to many disconnected small-size
regression problems. To see the point, define $\tilde{Y}=X^{\prime}Y$ and
recall that $G=X^{\prime}X$. Let ${\cal I}_{0}\lhd{\cal U}_{p}^{*}$ be a
component, we limit our attention to ${\cal I}_{0}$ by considering the
following regression problem:

$\tilde{Y}^{{\cal I}_{0}}=G^{{\cal I}_{0},\otimes}\beta+(X^{\prime}z)^{{\cal
I}_{0}},$ (3.53)

where $(X^{\prime}z)^{{\cal I}_{0}}\sim N(0,G^{{\cal I}_{0},{\cal
I}_{0}})\approx N(0,\Omega^{{\cal I}_{0},{\cal I}_{0}})$, and $G^{{\cal
I}_{0},\otimes}$ is a $|{\cal I}_{0}|\times p$ matrix according to our
notation. What is non-obvious here is that, the regression problem still
involves the whole vector $\beta$, and is still high-dimensional. To see the
point, letting $V=\\{1,2,\ldots,p\\}\setminus{\cal U}_{p}^{*}$, we write
$G^{{\cal I}_{0},\otimes}\beta=G^{{\cal I}_{0},{\cal I}_{0}}\beta^{{\cal
I}_{0}}+I+II$, where $I=\sum_{{\cal J}_{0}:{\cal J}_{0}\lhd{\cal
U}_{p}^{*},{\cal J}_{0}\neq{\cal I}_{0}}G^{{\cal I}_{0},{\cal
J}_{0}}\beta^{{\cal J}_{0}}$ and $II=G^{{\cal I}_{0},V}\beta^{V}$. First, by
Sure Screening property, $\beta^{V}$ contains only a negligible number of
signals, so we can think $II$ as negligible. Second, for any ${\cal
J}_{0}\neq{\cal I}_{0}$ and ${\cal J}_{0}\lhd{\cal U}_{p}^{*}$, by the SAS
property, ${\cal I}_{0}$ and ${\cal J}_{0}$ are disconnected and so the matrix
$G^{{\cal I}_{0},{\cal J}_{0}}$ is a small size matrix whose coordinates are
uniformly small. This heuristic is made precise in the proof of Theorem 8. It
is now seen that the regression problem in (3.53) is indeed low-dimensional:

$\tilde{Y}^{{\cal I}_{0}}\approx G^{{\cal I}_{0},{\cal I}_{0}}\beta^{{\cal
I}_{0}}+(X^{\prime}z)^{{\cal I}_{0}}\approx N(\Omega^{{\cal I}_{0},{\cal
I}_{0}}\beta^{{\cal I}_{0}},\Omega^{{\cal I}_{0},{\cal I}_{0}}),$ (3.54)

The above argument is made precise in Lemma 17, see details therein. Finally,
approximately, the $GC$-step is to minimize $\frac{1}{2}(\tilde{Y}^{{\cal
I}_{0}}-\Omega^{{\cal I}_{0},{\cal I}_{0}}\xi)^{\prime}(\Omega^{{\cal
I}_{0},{\cal I}_{0}})^{-1}(\tilde{Y}^{{\cal I}_{0}}-\Omega^{{\cal I}_{0},{\cal
I}_{0}}\xi)+\frac{1}{2}(u^{gs})^{2}\|\xi\|_{0}$, where each coordinate of
$\xi$ is either $0$ or $\geq v^{gs}$ in magnitude. Comparing this with (3.54),
the procedure is nothing but the penalized MLE of a low dimensional normal
model, and the main result follows by exercising basic statistical inferences.

We remark that in the $GC$-step, removing the constraints on the coordinates
of $\xi$ will not give the optimal rate of convergence. This is one of the
reasons why the classical subset selection procedure is rate non-optimal.
Another reason why the subset selection is non-optimal is that, the procedure
has only one tuning parameter, but GS has the flexibility of using different
tuning parameters in the $GS$-step and the $GC$-step. See Section 2.8 for more
discussion.

We are now ready for the proof of Theorem 8.

###  3.1 Proof of Theorem 8

For notational simplicity, we write
$\rho_{j}^{*}=\rho_{j}^{*}(\vartheta,r,a,\Omega)$. By Lemma 15,

$\sum_{j=1}^{p}P(\beta_{j}\neq 0,j\notin{\cal U}_{p}^{*})\leq
L_{p}[p^{1-(m_{0}+1)\vartheta}+\sum_{j=1}^{p}p^{-\rho_{j}^{*}}]+o(1).$ (3.55)

So to show the claim, it is sufficient to show

$\sum_{j=1}^{p}P(j\in{\cal U}_{p}^{*},\hbox{\rm sgn}(\beta_{j})\neq\hbox{\rm
sgn}(\hat{\beta}_{j}))\leq
L_{p}[\sum_{j=1}^{p}p^{-\rho_{j}^{*}}+p^{1-(m_{0}+1)\vartheta}]+o(1).$ (3.56)

Towards this end, let $S(\beta)$ be the support of $\beta$,
$\Omega^{*,\delta}$ be as in (1.6), and ${\cal G}^{*,\delta}$ be the GOSD. Let
${\cal U}_{p}^{*}$ be the set of retained indices after the GS-step. Note that
when $\hbox{\rm sgn}(\hat{\beta}_{j})\neq 0$, there is a unique component
${\cal I}_{0}$ such that $j\in{\cal I}_{0}\lhd{\cal U}_{p}^{*}$. For any
connected subgraph ${\cal I}_{0}$ of ${\cal G}^{*,\delta}$, let $B({\cal
I}_{0})=\\{\mbox{$k$: $k\notin{\cal I}_{0}$, $\Omega^{*,\delta}(k,\ell)\neq 0$
for some $\ell\in{\cal I}_{0}$, $1\leq k\leq p$}\\}$. Note that when ${\cal
I}_{0}$ is a component of ${\cal U}_{p}^{*}$, we must have $B({\cal
I}_{0})\cap{\cal U}_{p}^{*}=\emptyset$ as for any node in $B({\cal I}_{0})$,
there is at least one edge between it and some nodes in the component ${\cal
I}_{0}$. As a result,

$P(j\in{\cal I}_{0}\lhd{\cal U}_{p}^{*},B({\cal I}_{0})\cap
S(\beta)\neq\emptyset)\leq\sum_{{\cal I}_{0}:j\in{\cal I}_{0}}\sum_{k\in
B({\cal I}_{0})}P(k\notin{\cal U}_{p}^{*},\beta_{k}\neq 0),$ (3.57)

where the first summation is over all connected subgraphs that contains node
$j$. By Lemma 16, with probability at least $1-o(1/p)$, ${\cal G}^{*,\delta}$
is $K$-sparse with $K=C(\log(p))^{1/\gamma}$, and there is a finite integer
$\ell_{0}$ such that $|{\cal I}_{0}|\leq\ell_{0}$. As a result, there are at
most finite ${\cal I}_{0}$ such that the event $\\{j\in{\cal I}_{0}\lhd{\cal
U}_{p}^{*}\\}$ is non-empty, and for each of such ${\cal I}_{0}$, $B({\cal
I}_{0})$ contains at most $L_{p}$ nodes. Using (3.57) and Lemma 15, a direct
result is

$\sum_{j=1}^{p}P(j\in{\cal I}_{0}\lhd{\cal U}_{p}^{*},B({\cal I}_{0})\cap
S(\beta)\neq\emptyset)\leq
L_{p}[\sum_{j=1}^{p}p^{-\rho_{j}^{*}}+p^{1-(m_{0}+1)\vartheta}]+o(1).$ (3.58)

Comparing (3.58) with (3.56), to show the claim, it is sufficient to show that

$\sum_{j=1}^{p}P(\hbox{\rm sgn}(\beta_{j})\neq\hbox{\rm
sgn}(\hat{\beta}_{j}),j\in{\cal I}_{0}\lhd{\cal U}_{p}^{*},B({\cal I}_{0})\cap
S(\beta)=\emptyset)\leq
L_{p}[\sum_{j=1}^{p}p^{-\rho_{j}^{*}}+p^{1-(m_{0}+1)\vartheta}]+o(1).$ (3.59)

Fix $1\leq j\leq p$ and a connected subgraph ${\cal I}_{0}$ such that
$j\in{\cal I}_{0}$. For short, let $S$ be the support of $\beta^{{\cal
I}_{0}}$ and $\hat{S}$ be the support of $\hat{\beta}^{{\cal I}_{0}}$. The
event $\\{\hbox{\rm sgn}(\beta_{j})\neq\hbox{\rm
sgn}(\hat{\beta}_{j}),j\in{\cal I}_{0}\lhd{\cal U}_{p}^{*}\\}$ is identical to
the event of $\\{\hbox{\rm sgn}(\beta_{j})\neq\hbox{\rm
sgn}(\hat{\beta}_{j}),j\in S\cup\hat{S}\\}$. Moreover, Since ${\cal I}_{0}$
has a finite size, both $S$ and $\hat{S}$ have finite possibilities. So to
show (3.59), it is sufficient to show that for any fixed $1\leq j\leq p$,
connected subgraph ${\cal I}_{0}$, and subsets $S_{0},S_{1}\subset{\cal
I}_{0}$ such that $j\in S_{0}\cup S_{1}$,

$P(\hbox{\rm sgn}(\beta_{j})\neq\hbox{\rm
sgn}(\hat{\beta}_{j}),S=S_{0},\hat{S}=S_{1},j\in{\cal I}_{0}\lhd{\cal
U}_{p}^{*},B({\cal I}_{0})\cap S(\beta)=\emptyset)\leq
L_{p}[p^{-\rho_{j}^{*}}+p^{-(m_{0}+1)\vartheta}].$ (3.60)

We now show (3.60). The following lemma is proved in Ji and Jin (2011, A.4).

###### Lemma 17

Suppose the conditions of Theorem 8 hold. Over the event $\\{j\in{\cal
I}_{0}\lhd{\cal U}_{p}^{*}\\}\cap\\{B({\cal I}_{0})\cap
S(\beta)=\emptyset\\}$, $\|(\Omega\beta)^{{\cal I}_{0}}-\Omega^{{\cal
I}_{0},{\cal I}_{0}}\beta^{{\cal I}_{0}}\|_{\infty}\leq
C\tau_{p}(\log(p))^{-(1-\gamma)}$.

Write for short $\hat{M}=G^{{\cal I}_{0},{\cal I}_{0}}$ and $M=\Omega^{{\cal
I}_{0},{\cal I}_{0}}$. By definitions, $\hat{\beta}^{{\cal I}_{0}}$ is the
minimizer of the following functional
$Q(\xi)\equiv\frac{1}{2}(\tilde{Y}^{{\cal
I}_{0}}-\hat{M}\xi)^{\prime}\hat{M}^{-1}(\tilde{Y}^{{\cal
I}_{0}}-\hat{M}\xi)+\frac{1}{2}(u^{gs})^{2}\|\xi\|_{0}$, where $\xi$ is an
$|{\cal I}_{0}|\times 1$ vector whose coordinates are either $0$ or $\geq
v^{gs}$ in magnitude, $u^{gs}=\sqrt{2\vartheta\log(p)}$, and
$v^{gs}=\sqrt{2r\log(p)}$. In particular, $Q(\beta^{{\cal I}_{0}})\geq
Q(\hat{\beta}^{{\cal I}_{0}})$, or equivalently

$(\hat{\beta}^{{\cal I}_{0}}-\beta^{{\cal I}_{0}})^{\prime}(\tilde{Y}^{{\cal
I}_{0}}-\hat{M}\beta^{{\cal I}_{0}})\geq\frac{1}{2}(\hat{\beta}^{{\cal
I}_{0}}-\beta^{{\cal I}_{0}})^{\prime}\hat{M}(\hat{\beta}^{{\cal
I}_{0}}-\beta^{{\cal I}_{0}})+(|S_{1}|-|S_{0}|)\vartheta\log(p).$ (3.61)

Now, write for short $\delta=\tau_{p}^{-2}(\hat{\beta}^{{\cal
I}_{0}}-\beta^{{\cal I}_{0}})^{\prime}M(\hat{\beta}^{{\cal
I}_{0}}-\beta^{{\cal I}_{0}})$. First, by Schwartz inequality,
$[(\hat{\beta}^{{\cal I}_{0}}-\beta^{{\cal I}_{0}})^{\prime}(\tilde{Y}^{{\cal
I}_{0}}-\hat{M}\beta^{{\cal
I}_{0}})]^{2}\leq\delta\tau_{p}^{2}(\tilde{Y}^{{\cal
I}_{0}}-\hat{M}\beta^{{\cal I}_{0}})^{\prime}M^{-1}(\tilde{Y}^{{\cal
I}_{0}}-\hat{M}\beta^{{\cal I}_{0}})$. Second, by Lemma 17, $\tilde{Y}^{{\cal
I}_{0}}=w+M\beta^{{\cal I}_{0}}+rem$, where $w\sim N(0,M)$ and with
probability $1-o(1/p)$, $|rem|\leq C(\log(p))^{-(1-\gamma)}\tau_{p}$. Last,
with probability at least $(1-o(1/p))$,
$\parallel\hat{M}-M\parallel_{\infty}\leq
C\sqrt{\log(p)}p^{-[\kappa-(1-\vartheta)]/2}$. Inserting these into (3.61)
gives that with probability at least $(1-o(1/p))$,
$w^{\prime}M^{-1}w\geq\frac{1}{4}\biggl{[}\bigl{(}\sqrt{\delta
r}+\frac{(|S_{1}|-|S_{0}|)\vartheta}{\sqrt{\delta
r}}\bigr{)}_{+}\biggr{]}^{2}(2\log(p))+O((\log(p))^{\gamma})$. Since
$\gamma<1$, $O((\log(p))^{\gamma})$ is negligible. We note that
$w^{\prime}M^{-1}w\sim\chi_{|{\cal I}_{0}|}^{2}(0)$. Inserting this back to
(3.60), the left hand side $\leq\epsilon_{p}^{|S_{0}|}P(\chi_{|{\cal
I}_{0}|}^{2}(0)\geq[(\sqrt{\delta r}+(|S_{1}|-|S_{0}|)\vartheta/\sqrt{\delta
r})_{+}]^{2}(\log(p)/2))+o(1/p)$. Assume $\hbox{\rm
sgn}(\beta_{j})\neq\hbox{\rm sgn}(\hat{\beta}_{j})$, and fix all parameters
except $\delta$, $S_{0}$ and $S_{1}$. By arguments similar to the proof of
Lemma 18, the above quantity cannot achieve its maximum in the cases where
$S_{0}=S_{1}$. Hence we only need to consider the cases where $S_{0}\neq
S_{1}$. We also only need to consider the cases where
$\max(|S_{0}|,|S_{1}|)\leq m_{0}$, since the sum of the probabilities of other
cases is controlled by $p^{1-(m_{0}+1)\vartheta}$. The claim follows by the
definitions of $\rho_{j}^{*}$. $\Box$

##  4 Simulations

We conduct a small-scale simulation study to investigate the numerical
performance of Graphlet Screening and compare it with the lasso and the UPS.
The subset selection is not included for comparison since it is
computationally NP hard. We consider the experiments for both random design
and fixed design, where as before, the parameters $(\epsilon_{p},\tau_{p})$
are tied to $(\vartheta,r)$ by $\epsilon_{p}=p^{-\vartheta}$ and
$\tau_{p}=\sqrt{2r\log(p)}$ (we assume $\sigma=1$ for simplicity in this
section).

In random design settings where $p$ is not very large, we follow the spirit of
the refined UPS in Ji and Jin (2011) and propose the iterative Graphlet
Screening algorithm where we iterate Graphlet Screening for a few times ($\leq
5$). The main purpose for the iteration is to denoise the Gram matrix; see Ji
and Jin (2011, Section 3) for more discussion.

Even with the refinement as in Ji and Jin (2011, Section 3), UPS behaves
poorly for most examples presented below. Over close investigations, we find
out that this is due to the threshold choice in the initial $U$-step is too
low, and increasing the threshold largely increases the performance. Note that
the purpose of this step is to denoise the Gram matrix Ji and Jin (2011,
Section 3), not for signal retainment, and so a larger threshold helps.

In this section, we use this improved version of refined UPS, but for
simplicity, we still call it the refined UPS. With that being said, recall
that UPS is unable to resolve the problem of signal cancellation, so it
usually performs poorer than GS, especially when the effect of signal
cancellation is strong. For this reason, part of the comparison is between GS
and the lasso only.

The experiments with random design contain the following steps.

  1. 1.

Fix $(p,\vartheta,r,\mu,\Omega)$ such that $\mu\in\Theta_{p}(\tau_{p})$.
Generate a vector $b=(b_{1},b_{2},\ldots,b_{p})^{\prime}$ such that
$b_{i}\stackrel{{\scriptstyle iid}}{{\sim}}\mathrm{Bernoulli}(\epsilon_{p})$,
and set $\beta=b\circ\mu$.

  2. 2.

Fix $\kappa$ and let $n=n_{p}=p^{\kappa}$. Generate an $n\times p$ matrix with
$iid$ rows from $N(0,(1/n)\Omega)$.

  3. 3.

Generate $Y\sim N(X\beta,I_{n})$, and apply the iterative Graphlet Screening,
the refined UPS and the lasso.

  4. 4.

Repeat 1-3 independently, and record the average Hamming distances or the
Hamming ratio, the ratio of the Hamming distance and the number of the
signals.

The steps for fixed design experiments are similar, except for that $n_{p}=p$,
$X=\Omega^{1/2}$ and we apply GS and UPS directly.

GS uses tuning parameters $(m_{0},{\cal Q},u^{gs},v^{gs})$. We set $m_{0}=3$
for our experiments, which is usually large enough due to signal sparsity. The
choice of ${\cal Q}$ is not critical, as long as the corresponding parameter
$q$ satisfies (2.41), and we use the maximal ${\cal Q}$ satisfying (2.41) in
most experiments. Numerical studies below (e.g., Experiment $5a$) support this
point. In principle, the optimal choices of $(u^{gs},v^{gs})$ depend on the
unknown parameters $(\epsilon_{p},\tau_{p})$, and how to estimate them in
general settings is a lasting open problem (even for linear models with
orthogonal designs). Fortunately, our studies (e.g., Experiment 5b-5d) show
that mis-specifying parameters $(\epsilon_{p},\tau_{p})$ by a reasonable
amount does not significantly affect the performance of the procedure. For
this reason, in most experiments below, assuming $(\epsilon_{p},\tau_{p})$ are
known, we set $(u^{gs},v^{gs})$ as $(\sqrt{2\log(1/\epsilon_{p})},\tau_{p})$.
For the iterative Graphlet Screening, we use the same tuning parameters in
each iteration.

For the UPS and the refined UPS, we use the tuning parameters
$(u^{ups},v^{ups})=(u^{gs},v^{gs})$. For both the iterative Graphlet Screening
and the refined UPS, we use the following as the initial estimate:
$\hat{\beta}_{i}=\hbox{\rm sgn}(\tilde{Y}_{i})\cdot
1\\{|\tilde{Y}_{i}|\geq\tau_{p}\\}$, $1\leq i\leq p$, where
$\tilde{Y}=X^{\prime}Y$. The main purpose of initial estimate is to denoise
the Gram matrix, not for screening. We use glmnet package (Friedman et al,
2010) to perform lasso. To be fair in comparison, we apply the lasso with all
tuning parameters, and we report the Hamming error associated with the “best”
tuning parameter.

The simulations contain $6$ different experiments which we now describe
separately.

Experiment 1. The goal of this experiment is two-fold. First, we compare GS
with UPS and the lasso in the fixed design setting. Second, we investigate the
minimum signal strength levels $\tau_{p}$ required by these three methods to
yield exact recovery, respectively.

Fixing $p=0.5\times 10^{4}$, we let $\epsilon_{p}=p^{-\vartheta}$ for
$\vartheta\in\\{0.25,0.4,0.55\\}$, and $\tau_{p}\in\\{6,7,8,9,10\\}$. We use a
fixed design model where $\Omega$ is a symmetric diagonal block-wise matrix,
where each block is a $2\times 2$ matrix, with $1$ on the diagonals, and $\pm
0.7$ on the off-diagonals (the signs alternate across different blocks).
Recall the $\beta=b\circ\mu$. For each pair of $(\epsilon_{p},\tau_{p})$, we
generate $b$ as $p$ iid samples from $Bernoulli(\epsilon_{p})$, and we let
$\mu$ be the vector where the signs of $\mu_{i}=\pm 1$ with equal
probabilities, and $|\mu_{i}|\stackrel{{\scriptstyle
iid}}{{\sim}}0.8\nu_{\tau_{p}}+0.2h$, where $\nu_{\tau_{p}}$ is the point mass
at $\tau_{p}$ and $h(x)$ is the density of $\tau_{p}(1+V/6)$ with
$V\sim\chi_{1}^{2}$. The average Hamming errors across $40$ repetitions are
tabulated in Table 3. For all $(\vartheta,\tau_{p})$ in this experiment, GS
behaves more satisfactorily than the UPS, which in turn behaves more
satisfactorily than the lasso.

Suppose we say a method yields ‘exact recovery’ if the average Hamming error
$\leq 3$. Then, when $\vartheta=0.25$, the minimum $\tau_{p}$ for GS to yield
exact recovery is $\tau_{p}\approx 8$, but that for UPS and the lasso are much
larger ($\geq 10$). For larger $\vartheta$, the differences are less
prominent, but the pattern is similar.

The comparison between GS and UPS is particularly interesting. Due to the
block structure of $\Omega$, as $\vartheta$ decreases, the signals become
increasingly less sparse, and the effects of signal cancellation become
increasingly stronger. As a result, the advantage of GS over the UPS becomes
increasingly more prominent.

| $\tau_{p}$ | 6 | 7 | 8 | 9 | 10
---|---|---|---|---|---|---
$\vartheta=0.25$ | Graphic Screening | 24.7750 | 8.6750 | 2.8250 | 0.5250 | 0.1250
UPS | 48.5500 | 34.6250 | 36.3500 | 30.8750 | 33.4000
lasso | 66.4750 | 47.7000 | 43.5250 | 35.2500 | 35.0500
$\vartheta=0.40$ | Graphic Screening | 6.9500 | 2.1500 | 0.4000 | 0.0750 | 0.0500
UPS | 7.7500 | 4.0000 | 2.2000 | 2.7750 | 2.4250
lasso | 12.8750 | 6.8000 | 4.3250 | 3.7500 | 2.6750
$\vartheta=0.55$ | Graphic Screening | 1.8750 | 0.8000 | 0.3250 | 0.2250 | 0.1250
UPS | 1.8750 | 0.8000 | 0.3250 | 0.2250 | 0.1250
lasso | 2.5000 | 1.1000 | 0.7750 | 0.2750 | 0.1250
Table 3: Comparison of average Hamming errors (Experiment 1).

Experiment 2. In this experiment, we compare GS, UPS and the lasso in the
random design setting, and investigate the effect of signal cancellation on
their performances. We fix $(p,\kappa,\vartheta,r)=(0.5\times
10^{4},0.975,0.35,3)$, and assume $\Omega$ is blockwise diagonal. We generate
$\mu$ as in Experiment $1$, but to better illustrate the difference between
UPS and GS in the presence of signal cancellation, we generate the vector $b$
differently and allow it to depend on $\Omega$. The experiment contains $2$
parts, $2a$ and $2b$.

In Experiment 2a, $\Omega$ is the block-wise matrix where each block is $2$ by
$2$ matrix with $1$ on the diagonals and $\pm.5$ on the off diagonals (the
signs alternate on adjacent blocks). According to the blocks in $\Omega$, the
set of indices $\\{1,2,\ldots,p\\}$ are also partitioned into blocks
accordingly. For any fixed $\vartheta$ and
$\eta\in\\{0,0.01,0.02,0.03,0.04,0.05,0.06,0.1,0.2\\}$, we randomly choose
$(1-2p^{-\vartheta})$ fraction of the blocks (of indices) where $b$ is $0$ at
both indices, $2(1-\eta)p^{-\vartheta}$ fraction of the blocks where $b$ is
$0$ at one index and $1$ at the other (two indices are equally likely to be
$0$), $2\eta p^{-\vartheta}$ faction of the blocks where $b$ is $1$ on both
indices.

Experiment 2b has similar settings, where the difference is that (a) we choose
$\Omega$ to be a diagonal block matrix where each block is a $4$ by $4$ matrix
(say, denoted by $A$) satisfying $A(i,j)=1\\{i=j\\}+0.4\cdot
1\\{|i-j|=1\\}\cdot\hbox{\rm sgn}(6-i-j)+0.05\\{|i-j|\geq 2\\}\cdot\hbox{\rm
sgn}(5.5-i-j)$, $1\leq i,j\leq 4$, and (b) $(1-4p^{-\vartheta})$ is the
fraction of blocks where $b$ is nonzero in $k=0$ indices,
$4(1-\eta)p^{-\vartheta}$ is that for $k=1$, and $4\eta p^{-\vartheta}$ is
that for $k\in\\{2,3,4\\}$ in total. In a block where $\beta$ is nonzero at
$k$ indices, all configurations with $k=1$ are equally likely, and all those
with $k\in\\{2,3,4\\}$ are equally likely.

The average Hamming ratio results across $40$ runs for two Experiment 2a and
2b are reported in Figure 3, where UPS and GS consistently outperform the
lasso. Additionally, when $\eta$ is small, the effect of signal cancellation
is negligible, so UPS and GS have similar performances. However, when $\eta$
increases, the effects of signal cancellation grows, and the advantage of GS
over UPS becomes increasingly more prominent.

(a) Blockwise diagonal $\Omega$ in $2\times 2$ blocks

(b) Blockwise diagonal $\Omega$ in $4\times 4$ blocks

Figure 3: Hamming ratio results in Experiment 2

Through Experiment 1-2, the comparison of UPS and GS is more or less
understood. For this reason, we do not include UPS for study in Experiment
3-5, but we include UPS for study in Experiment 6 where we investigate
robustness of all three methods.

Experiment 3. In this experiment, we investigate how different choices of
signal vector $\beta$ affect the comparisons of GS and the lasso. We use a
random design model, and $\Omega$ is a symmetric tri-diagonal correlation
matrix where the vector on each sub-diagonal consists of blocks of
$(.4,.4,-.4)^{\prime}$. Fix $(p,\kappa)=(0.5\times 10^{4},0.975)$ (note
$n=p^{\kappa}\approx 4,000$). We let $\epsilon_{p}=p^{-\vartheta}$ with
$\vartheta\in\\{0.35,0.5\\}$ and let $\tau_{p}\in\\{6,8,10\\}$. For each
combination of $(\epsilon_{p},\tau_{p})$, we consider two choices of $\mu$.
For the first choice, we let $\mu$ be the vector where all coordinates equal
to $\tau_{p}$ (note $\beta$ is still sparse). For the second one, we let $\mu$
be as in Experiment $1$. The average Hamming ratios for both procedures across
$40$ repetitions are tabulated in Table 4.

$\tau_{p}$ | 6 | 8 | 10
---|---|---|---
Signal Strength | Equal | Unequal | Equal | Unequal | Equal | Unequal
$\vartheta=0.35$ | Graphic Screening | 0.0810 | 0.0825 | 0.0018 | 0.0034 | 0 | 0.0003
lasso | 0.2424 | 0.2535 | 0.1445 | 0.1556 | 0.0941 | 0.1109
$\vartheta=0.5$ | Graphic Screening | 0.0315 | 0.0297 | 0.0007 | 0.0007 | 0 | 0
lasso | 0.1107 | 0.1130 | 0.0320 | 0.0254 | 0.0064 | 0.0115
Table 4: Hamming ratio results of Experiment $3$, where “Equal” and “Unequal”
stand for the first and the second choices of $\mu$, respectively.

Experiment 4. In this experiment, we generate $\beta$ the same way as in
Experiment $1$, and investigate how different choices of design matrices
affect the performance of the two methods. Setting
$(p,\vartheta,\kappa)=(0.5\times 10^{4},0.35,0.975)$ and
$\tau_{p}\in\\{6,7,8,9,10,11,12\\}$, we use Gaussian random design model for
the study. The experiment contains $3$ sub-experiments $4a$-$4c$.

In Experiment $4a$, we set $\Omega$ as the symmetric diagonal block-wise
matrix, where each block is a $2\times 2$ matrix, with $1$ on the diagonals,
and $\pm 0.5$ on the off-diagonals (the signs alternate across different
blocks). The average Hamming ratios of $40$ repetitions are reported in Figure
4.

In Experiment $4b$, we set $\Omega$ as a symmetric penta-diagonal correlation
matrix, where the main diagonal are ones, the first sub-diagonal consists of
blocks of $(.4,.4,-.4)^{\prime}$, and the second sub-diagonal consists of
blocks of $(.05,-.05)^{\prime}$. The average Hamming ratios across $40$
repetitions are reported in Figure 4.

In Experiment $4c$, we generate $\Omega$ as follows. First, we generate
$\Omega$ using the function _sprandsym(p,K/p)_ in _matlab_. We then set the
diagonals of $\Omega$ to be zero, and remove some of entries so that $\Omega$
is $K$-sparse for a pre-specified $K$. We then normalize each non-zero entry
by the sum of the absolute values in that row or that column, whichever is
larger, and multiply each entry by a pre-specified positive constant $A$.
Last, we set the diagonal elements to be 1. We choose $K=3$ and $A=0.7$, draw
$5$ different $\Omega$ with this method, and for each of them, we draw
$(X,\beta,z)$ $10$ times independently. The average Hamming ratios are
reported in Figure 4. The results suggest that GS is consistently better than
the lasso.

Figure 4: $x$-axis: $\tau_{p}$. $y$-axis: Hamming ratios. Left to right:
Experiment $4a$, $4b$, and $4c$.

Experiment 5. In this experiment, we investigate how sensitive GS is with
respect to the tuning parameters. The experiment contains $4$ sub-experiments,
$5a$-$5d$. In Experiment $5a$, we investigate how sensitive the procedure is
with respect to the tuning parameter $q$ in ${\cal Q}$ (recall that the main
results hold as long as $q$ fall into the range given in (2.41)), where we
assume $(\epsilon_{p},\tau_{p})$ are known. In Experiment $5b$-$5d$, we mis-
specify $(\epsilon_{p},\tau_{p})$ by a reasonably small amount, and
investigate how the mis-specification affect the performance of the procedure.
For the whole experiment, we choose $\beta$ the same as in Experiment $1$, and
$\Omega$ the same as in Experiment $4b$. We use a fixed design model in
Experiment $5a$-$5c$, and a random design model in Experiment $5d$. For each
sub-experiment, the results are based on $40$ independent repetitions. We now
describe the sub-experiments with details.

In Experiment $5a$, we choose $\vartheta\in\\{0.35,0.6\\}$ and
$r\in\\{1.5,3\\}$. In GS, let $q_{max}=q_{max}(\hat{D},\hat{F})$ be the
maximum value of $q$ satisfying (2.41). For each combination of
$(\vartheta,r)$ and $(\hat{D},\hat{F})$, we choose
$q(\hat{D},\hat{F})=q_{max}(\hat{D},\hat{F})\times\\{0.7,0.8,0.9,1,1.1,1.2\\}$
for our experiment. The results are tabulated in Table 5, which suggest that
different choices of $q$ have little influence over the variable selection
errors. We must note that the larger we set $q(\hat{D},\hat{F})$, the faster
the algorithm runs.

$q(\hat{F},\hat{D})/q_{max}(\hat{F},\hat{D})$ | 0.7 | 0.8 | 0.9 | 1 | 1.1 | 1.2
---|---|---|---|---|---|---
$(\vartheta,r)=(0.35,1.5)$ | 0.0782 | 0.0707 | 0.0661 | 0.0675 | 0.0684 | 0.0702
$(\vartheta,r)=(0.35,3)$ | 0.0066 | 0.0049 | 0.0036 | 0.0034 | 0.0033 | 0.0032
$(\vartheta,r)=(0.6,1.5)$ | 0.1417 | 0.1417 | 0.1417 | 0.1417 | 0.1417 | 0.1417
$(\vartheta,r)=(0.6,3)$ | 0.0089 | 0.0089 | 0.0089 | 0.0089 | 0.0089 | 0.0089
Table 5: Hamming ratio results in Experiment $5a$.

In Experiment $5b$, we use the same settings as in Experiment $5a$, but we
assume $\vartheta$ (and so $\epsilon_{p}$) is unknown (the parameter $r$ is
assumed as known, however), and let $\vartheta^{*}$ is the misspecified value
of $\vartheta$. We take
$\vartheta^{*}\in\vartheta\times\\{0.85,0.925,1,1.075,1.15,1.225\\}$ for the
experiment.

In Experiment $5c$, we use the same settings as in Experiment $5b$, but we
assume $r$ (and so $\tau_{p}$) is unknown (the parameter $\vartheta$ is
assumed as known, however), and let $r^{*}$ is the misspecified value of $r$.
We take $r^{*}=r\times\\{0.8,0.9,1,1.1,1.2,1.3\\}$ for the experiment.

In Experiment $5b$-$5c$, we run GS with tuning parameters set as in Experiment
$1$, except $\vartheta$ or $r$ are replaced by the misspecified counterparts
$\vartheta^{*}$ and $r^{*}$, respectively. The results are reported in Table
6, which suggest that the mis-specifications have little effect as long as
$r^{*}/r$ and $\vartheta^{*}/\vartheta$ are reasonably close to $1$.

$\vartheta^{*}/\vartheta$ | 0.85 | 0.925 | 1 | 1.075 | 1.15 | 1.225
---|---|---|---|---|---|---
$(\vartheta,r)=(0.35,1.5)$ | 0.0799 | 0.0753 | 0.0711 | 0.0710 | 0.0715 | 0.0746
$(\vartheta,r)=(0.35,3)$ | 0.0026 | 0.0023 | 0.0029 | 0.0030 | 0.0031 | 0.0028
$(\vartheta,r)=(0.6,1.5)$ | 0.1468 | 0.1313 | 0.1272 | 0.1280 | 0.1247 | 0.1296
$(\vartheta,r)=(0.6,3)$ | 0.0122 | 0.0122 | 0.0139 | 0.0139 | 0.0130 | 0.0147
$r^{*}/r$ | 0.8 | 0.9 | 1 | 1.1 | 1.2 | 1.3
$(\vartheta,r)=(0.35,1.5)$ | 0.0843 | 0.0731 | 0.0683 | 0.0645 | 0.0656 | 0.0687
$(\vartheta,r)=(0.35,3)$ | 0.0062 | 0.0039 | 0.0029 | 0.0030 | 0.0041 | 0.0054
$(\vartheta,r)=(0.6,1.5)$ | 0.1542 | 0.1365 | 0.1277 | 0.1237 | 0.1229 | 0.1261
$(\vartheta,r)=(0.6,3)$ | 0.0102 | 0.0076 | 0.0085 | 0.0059 | 0.0051 | 0.0076
Table 6: Hamming ratio results in Experiment $5b$ (top) and in Experiment $5c$
(bottom).

In Experiment $5d$, we re-examine the mis-specification issue in the random
design setting. We use the same settings as in Experiment $5b$ and Experiment
$5c$, except for (a) while we use the same $\Omega$ as in Experiment $5b$, the
design matrix $X$ are generated according to the random design model as in
Experiment $4b$, and (b) we only investigate for the case of $r=2$ and
$\vartheta\in\\{0.35,0.6\\}$. The results are summarized in Table 7, which is
consistent with the results in $5b$-$5c$.

$\vartheta^{*}/\vartheta$ | 0.85 | 0.925 | 1 | 1.075 | 1.15 | 1.225
---|---|---|---|---|---|---
$(\vartheta,r)=(0.35,2)$ | 0.1730 | 0.1367 | 0.1145 | 0.1118 | 0.0880 | 0.0983
$(\vartheta,r)=(0.6,2)$ | 0.0583 | 0.0591 | 0.0477 | 0.0487 | 0.0446 | 0.0431
$r^{*}/r$ | 0.8 | 0.9 | 1 | 1.1 | 1.2 | 1.3
$(\vartheta,r)=(0.35,2)$ | 0.1881 | 0.1192 | 0.1275 | 0.1211 | 0.1474 | 0.1920
$(\vartheta,r)=(0.6,2)$ | 0.0813 | 0.0515 | 0.0536 | 0.0397 | 0.0442 | 0.0510
Table 7: Hamming ratio results in Experiment $5d$.

Experiment $6$. In this experiment, we investigate the robustness of all three
methods for the mis-specification of the linear model (1.1). We use the random
design setting as in Experiment $4b$, except that we fix
$(\vartheta,r)=(0.35,3)$. The experiment contains $3$ sub-experiments,
$6a$-$6c$, where we consider three scenarios where the linear model (1.1) is
in question: the presence of nonGaussianity, the presence of missing
predictors, and the presence of non-linearity, correspondingly.

In Experiment $6a$, we assume the noise vector $z$ in Model (1.1) is
nonGaussian, where the coordinates are iid samples from a $t$-distribution
with the same degree of freedom (df) (we assume that $z$ is normalized so each
coordinate has unit variance), where the df range in
$\\{3,4,5,6,7,8,9,10,30,50\\}$. Figure 5(a) shows how the Hamming ratios
(based on $40$ independent repetitions) change when the $df$ decreases. The
results suggest that all three methods are reasonably robust against
nonGaussianity, but GS continues to have the best performance.

In Experiment $6b$, we assume that the true model is $Y=X\beta+z$ where
$(X,\beta,z)$ are generated as in $4b$, but the model that is accessible to us
is a misspecified model where the some of the true predictors are missing. Fix
$\eta\in(0,1)$, and let $S(\beta)$ be the support of $\beta$. For each $i\in
S(\beta)$, we flip a coin that lands on head with probability $\eta$, and we
retain $i$ if and only if the coin lands on tail. Let $S^{*}\subset S(\beta)$
be the set of retained indices, and let $R=S^{*}\cup S^{c}$. The misspecified
model we consider is then $Y=X^{\otimes,R}\beta^{R}+z$.

For the experiment, we let $\eta$ range in
$0.02\times\\{0,1,2,3,4,5,6,7,8,9,10\\}$. The average Hamming ratios (based on
$40$ independent repetitions) are reported in Figure 5(b). The results suggest
that all three results are reasonably robust to missing predictors, with the
lasso being the most robust. However, as long as the proportion of true
predictors that are missing is reasonably small (say, $\eta\leq.1$), GS
continues to outperform UPS and the lasso.

In Experiment $6c$, for $i=1,\ldots,n$, the true model is an additive model in
the form of $Y_{i}=\sum_{j=1}^{p}f_{j}(X_{ij})\beta_{j}+z_{i}$, but what is
accessible to us is the linear model
$Y_{i}=\sum_{j=1}^{p}X_{ij}\beta_{j}+z_{i}$ (and thus misspecified; the true
model is non-linear). For experiment, we let $(X,\beta,z)$ be generated as in
$4b$, and $S(\beta)$ be the support of $\beta$. Fixing $\eta\in(0,1)$, for
each $i\in S(\beta)$, we flip a coin that lands on head with probability
$\eta$, and let $S_{nl}\subset S(\beta)$ be all indices of the heads. We then
randomly split $S_{nl}$ into two sets $S_{1}$ and $S_{2}$ evenly. For
$j=1,\ldots,p$, we define $f_{j}(x)=[\hbox{\rm sgn}(x)x^{2}\cdot 1\\{j\in
S_{1}\\}+(e^{\sqrt{n}x}-a_{j})\cdot 1\\{j\in S_{2}\\}+x\cdot 1\\{j\in
S_{nl}^{c}\\}]/c_{j}$, where $a_{j}$ and $c_{j}$ are constants such that
$\\{f_{j}(X(i,j))\\}_{i=1}^{n}$ has mean $0$ and variance $1/n$.

For the experiment, we let $\eta$ range in $.05\times\\{0,1,2,3,4,5,6,7,8\\}$.
The average Hamming ratios (based on $40$ independent repetitions) are
reported in Figure 5(c). The results suggest that all three methods are
reasonably robust to the presence of nonlinearity, and GS continues to
outperform UPS and the lasso when the degree of nonlinearly is moderate (say,
$\eta<.2$).

(a) nonGaussian noise

(b) Missing predictors

(c) Non-linear predictors

Figure 5: Hamming ratio results in Experiment $6$

##  5 Connection to existing literature and possible extensions

Our idea of utilizing graph sparsity is related to the graphical lasso
(Meinshausen and Buhlmann, 2006; Friedman et al, 2008), which also attempts to
exploit graph structure. However, the setting we consider here is different
from that in Meinshausen and Buhlmann (2006); Friedman et al (2008), and our
emphasis on precise optimality and calibration is also very different. Our
method allows nearly optimal detection of very rare and weak effects, because
they are based on careful analysis that has revealed a number of subtle high-
dimensional effects (e.g., phase transitions) that we properly exploit.
Existing methodologies are not able to exploit or capture these phenomena, and
can be shown to fail at the levels of rare and weak effects where we are
successful.

The paper is closely related to the recent work by Ji and Jin (Ji and Jin,
2011) (see also Fan and Lv (2008); Genovese et al (2012)), and the two papers
use a similar rare and weak signal framework and a similar random design
model. However, they are different in important ways, since the technical
devise developed in Ji and Jin (2011) can not be extended to the current
study. For example, the lower bound derived in this paper is different and
sharper than that in Ji and Jin (2011). Also, the procedure in Ji and Jin
(2011) relies on marginal regression for screening. The limitation of marginal
regression is that it neglects the graph structure of GOSD for the regularized
Gram matrix (1.5), so that it is incapable of picking variables that have weak
marginal correlation but significant joint correlation to $Y$. Correct
selection of such hidden significant variables, termed as the challenge of
signal cancellation (Wasserman and Roeder, 2009), is the difficulty at the
heart of the variable selection problem. One of the main innovation of GS is
that it uses the graph structure to guide the screening, so that it is able to
successfully overcome the challenge of signal cancellation.

Additionally, two papers have very different objectives, and consequently the
underlying analysis are very different. The main results of each of these two
papers can not be deduced from the other. For example, to assess optimality,
Ji and Jin (2011) uses the criterion of the partition of the phase diagram,
while the current paper uses the minimax Hamming distance. Given the
complexity of the high dimensional variable selection, one type of optimality
does not imply the other, and vice versa. Also, the main result in Ji and Jin
(2011) focuses on conditions under which the optimal rate of convergence is
$L_{p}p^{1-(\vartheta+r)^{2}/(4r)}$ for the whole phase space. While this
overlaps with our Corollaries 10 and 11, we must note that Ji and Jin (2011)
deals with the much more difficult cases where $r/\vartheta$ can get arbitrary
large; and to ensure the success in that case, they assume very strong
conditions on the design matrix and the range of the signal strength. On the
other hand, the main focus of the current paper is on optimal variable
selection under conditions (of the Gram matrix $G$ as well as the signal
vector $\beta$) that are as general as possible.

While the study in this paper has been focused on the Random Design model
RD$(\vartheta,\kappa,\Omega)$, extensions to deterministic design models are
straightforward (in fact, in Corollary 7, we have already stated some results
on deterministic design models), and the omission of discussion on the latter
is largely for technical simplicity and the sake of space. In fact, for models
with deterministic designs, since the likelihood ratio test in the derivation
of the lower bound matches the penalized MLE in the cleaning step of GS, the
optimality of GS follows from the Sure Screening and Separable After Screening
properties of GS. The proof of these properties, and therefore the optimality
of GS, follows the same line as those for random design as long as
$\max_{j}|\sum_{i}\beta_{i}G(i,j)I\\{\Omega^{*,\delta}(i,j)=0\\}|/\tau_{p}$ is
small. This last condition on $G$ holds when
$p^{1-\vartheta}\|G-\Omega\|_{\infty}=o(1)$ with a certain $\Omega\in{\cal
M}^{*}_{p}(\gamma,c_{0},g,A)$. Alternatively, this condition holds when
$p^{1-\vartheta}\|G-\Omega\|_{\infty}^{2}\log p=o(1)$ with $\Omega\in{\cal
M}^{*}_{p}(\gamma,c_{0},g,A)$, provided that $\hbox{\rm sgn}(\beta_{j})$ are
iid symmetric random variables as in Candes and Plan (2009).

In this paper, we assume the signal vector $\beta$ is independent of the
design matrix $X$, and that $\beta$ is modeled by a Bernoulli model through
$\beta=b\circ\mu$. Both assumptions can be relaxed. In fact, in order for GS
to work, what we really need is some decomposability condition similar to that
in Lemma 1, where except for negligible probabilities, the maximum size of the
graphlets $m_{0}^{*}=m_{0}^{*}(S(\beta),G,\delta)$ is small. In many
situations, we can show that $m_{0}^{*}$ does not exceed a fixed integer. One
of such examples is as follows. Suppose for any fixed integer $m\geq 1$ and
size-$m$ subset $S$ of $\\{1,2,\ldots,p\\}$, there are constants $C>0$ and
$d>0$ such that the conditional probability $P(\beta_{j}\neq 0,\forall j\in
S|X)\leq Cp^{-dm}$. In fact, when such a condition holds, the claim follows
since ${\cal G}^{*,\delta}$ has no more than $C(eK)^{m}$ size-$m$ connected
subgraphs if it is $K$-sparse. See the proof of Lemma 1 for details. Note that
when $\epsilon_{p}=p^{-\vartheta}$ as in the ARW, then the condition holds for
the Bernoulli model in Lemma 1, with $d=\vartheta$. Note also that the
Bernoulli model can be replaced by some Ising models.

Another interesting direction of future research is the extension of GS to
more general models such as logistic regression. The extension of the lower
bound in Theorem 6 is relatively simple since the degree of GOLF can be
bounded using the true $\beta$. This indicates the optimality of GS in
logistic and other generalized linear models as long as proper generalized
likelihood ratio or Bayes tests are used in both the $GS$\- and $GC$-steps.

##  6 Proofs

In this section, we provide all technical proofs. We assume $\sigma=1$ for
simplicity.

###  6.1 Proof of Lemma 1

When ${\cal G}_{S}^{*,\delta}$ contains a connected subgraph of size $\geq
m_{0}+1$, it must contain a connected subgraph with size $m_{0}+1$. By Frieze
and Molloy (1999), there are $\leq p(eK)^{m_{0}+1}$ connected subgraph of size
$m_{0}+1$. Therefore, the probability that ${\cal G}_{S}^{*,\delta}$ has a
connected subgraph of size $(m_{0}+1)$ $\leq
p(eK)^{m_{0}+1}\epsilon_{p}^{m_{0}+1}$. Combining these gives the claim.
$\Box$

###  6.2 Proof of Theorem 6

Write for short $\rho_{j}^{*}=\rho_{j}^{*}(\vartheta,r,a,\Omega)$. Without
loss of generality, assume
$\rho_{1}^{*}\leq\rho_{2}^{*}\leq\ldots\leq\rho_{p}^{*}$. We construct indices
$i_{1}<i_{2}<\ldots<i_{m}$ as follows. (a) start with $B=\\{1,2,\ldots,p\\}$
and let $i_{1}=1$, (b) updating $B$ by removing $i_{1}$ and all nodes $j$ that
are neighbors of $i_{1}$ in GOLF, let $i_{2}$ be the smallest index, (c)
defining $i_{3},i_{4},\ldots,i_{m}$ by repeating (b), and terminates the
process when no indices is left in $B$. Since each time we remove at most
$d_{p}({\cal G}^{\diamond})$ nodes, it follows that

$\sum_{j=1}^{p}p^{-\rho_{j}^{*}}\leq d_{p}({\cal
G}^{\diamond})\sum_{k=1}^{m}p^{-\rho_{i_{k}}^{*}}.$ (6.62)

For each $1\leq j\leq p$, as before, let $(V_{0j}^{*},V_{1j}^{*})$ be the
least favorable configuration, and let
$(\theta_{*j}^{(0)},\theta_{*j}^{(1)})=\mathrm{argmin}_{\\{\theta^{(0)}\in
B_{V_{0j}^{*}},\theta^{(1)}\in B_{V_{1j}^{*}},\hbox{\rm
sgn}(\theta^{(0)})\neq\hbox{\rm
sgn}(\theta^{(1)})\\}}\alpha(\theta^{(0)},\theta^{(1)};\Omega)$. By our
notations, it is seen that

$\rho_{j}^{*}=\eta(V_{0j}^{*},V_{1j}^{*};\Omega),\qquad\alpha^{*}(V_{0j}^{*},V_{1j}^{*};\Omega)=\alpha(\theta_{*j}^{(0)},\theta_{*j}^{(1)};\Omega).$
(6.63)

In case $(\theta_{*j}^{(0)},\theta_{*j}^{(1)})$ is not unique, pick one
arbitrarily. We construct a $p\times 1$ vector $\mu^{*}$ as follows. Fix
$j\in\\{i_{1},\cdots,i_{m}\\}$. For all indices in $V_{0j}^{*}$, set the
constraint of $\mu^{*}$ on these indices to be $\theta_{*j}^{(0)}$. For any
index $i\notin\cup_{k=1}^{m}V_{0i_{k}}^{*}$, set $\mu^{*}_{i}=\tau_{p}$. Since

$\mathrm{Hamm}_{p}^{*}(\vartheta,\kappa,r,a,\Omega)\geq\inf_{\hat{\beta}}H_{p}(\hat{\beta};\epsilon_{p},n_{p},\mu^{*},\Omega)=\inf_{\hat{\beta}}\sum_{i=1}^{p}P(\hbox{\rm
sgn}(\hat{\beta}_{j})\neq\hbox{\rm sgn}(\beta_{j})),$ (6.64)

it follows that

$\mathrm{Hamm}_{p}^{*}(\vartheta,\kappa,r,a,\Omega)\geq\sum_{k=1}^{m}\sum_{j\in
V_{0i_{k}}\cup V_{1i_{k}}}P(\hbox{\rm sgn}(\hat{\beta}_{j})\neq\hbox{\rm
sgn}(\beta_{j})),$ (6.65)