Signed Support Recovery for Single Index Models in
High-Dimensions
arXiv:1511.02270v2 [] 23 Jun 2016

Matey Neykov∗

Qian Lin†

Jun S. Liu‡

Abstract
In this paper we study the support recovery problem for single index models Y = f (X | β, ε),
where f is an unknown link function, X ∼ Np (0, Ip ) and β is an s-sparse unit vector such that
β i ∈ {± √1s , 0}. In particular, we look into the performance of two computationally inexpensive
algorithms: (a) the diagonal thresholding sliced inverse regression (DT-SIR) introduced by
Lin et al. (2015); and (b) a semi-definite programming (SDP) approach inspired by Amini &
Wainwright (2008). When s = O(p1−δ ) for some δ > 0, we demonstrate that both procedures
n
can succeed in recovering the support of β as long as the rescaled sample size Γ = s log(p−s)
is
larger than a certain critical threshold. On the other hand, when Γ is smaller than a critical
value, any algorithm fails to recover the support with probability at least 12 asymptotically. In
other words, we demonstrate that both DT-SIR and the SDP approach are optimal (up to a
scalar) for recovering the support of β in terms of sample size. We provide extensive simulations,
as well as a real dataset application to help verify our theoretical observations.

Keywords: Single index models, Sliced inverse regression, Sparsity, Support recovery, Highdimensional statistics, Semidefinite programming

1

Introduction

Due to the recent advances in technology, collecting data becomes a routine. The ‘small n, large p’
characteristic of modern data brings new challenges in interpreting and processing it. Dimension
reduction and variable selection procedures become an indispensable step in data exploration.
Regrettably, the majority of the classical algorithms were developed to work in the regime p  n,
and hence one needs to exercise caution when applying existing methods in a high dimensional
setting. A firm understanding of the limitations of classical dimension reduction procedures in the
modern p  n regime, will help facilitate their appropriate application.
For example, the archetypical unsupervised dimension reduction procedure — principal component analysis (PCA), has been widely and successfully applied in a range of scientific problems
(Price et al. 2010, Brenner et al. 2000, e.g.). The behavior of PCA in high dimensions has been
well studied in the recent years. In Johnstone & Lu (2004, 2009), Paul (2007), it was shown that
in the spiked covariance model, PCA succeeds if and only if lim np 6= 0. This result stimulated the
statistical community to discuss the minimax rate of estimating the principal space under sparsity
∗

Department of Operations Research and Financial Engineering, Princeton University, Princeton, NJ 08544
Center of Mathematical Sciences and Applications, Harvard University, Cambridge, MA 02138
‡
Department of Statistics, Harvard University, Cambridge, MA 02138
†

1

constraints (see e.g., Cai et al. (2013), Vu & Lei (2012) and references therein) and the tradeoff
between the statistical and computational efficiency (see e.g., Berthet & Rigollet (2013)).
Another line of dimension reduction research, studies the so-called sufficient dimension reduction
(SDR). The goal of SDR is to find the minimal subspace S such that Y ⊥⊥ X PS X, where X ∈ Rp
is the predictor vector and Y ∈ R is the response (Cook et al. 2004, Li 1991, Cook & Ni 2005, e.g.).
Unlike its unsupervised counterpart, much less attention has been paid to how SDR algorithms
behave in a high dimensional setting. The optimal estimation rates of SDR algorithms in terms of
sparsity (s), dimensionality (p), and sample size (n) are unclear.
Sliced inverse regression, proposed by Li (1991), is one of the most popular SDR methods for
estimating the space S. When the dimensionality p is larger than or comparable to the sample size
n, sparsity assumptions are often imposed on the loading vector β (Li & Nachtsheim 2006, e.g.).
b β)] > 0 if ρ = lim p 6= 0 and sin ∠(β,
b β) = 0 when
Lin et al. (2015) proved that in fact E[∠(β,
n
b is the SIR estimator of β. In other words, the SIR estimator β
b is consistent (up to
ρ = 0 where β
p
a sign) if and only if ρ = lim n = 0. One implication of this result is that in order to estimate β,
structural assumptions such as sparsity are necessary in the high dimensional setting.
In the present paper, inspired by Amini & Wainwright (2008), we investigate the support
recovery problem for single index models (SIM) (1), under the sparsity assumption kβk0 = s in the
regime s = O(p1−δ ) for some δ > 0. More formally we study the models:
Y = f (X | β, ε) with X ∼ Np (0, Ip ),
where the noise ε is independent of X and f, ε, β belong to the class:
n
FA = (f, ε, β) : Var(E[Z|f (Z, ε)]) ≥ A where Z ∼ N (0, 1),
n
o
1 o
β i ∈ ± √ , 0 and (f, ε) is sliced stable (see Definition 1) .
s

(1)

(2)

Notice that no generality is lost in assuming that the vector β is a unit vector since otherwise
model (1) is not identifiable. Model class (2), further assumes the idealized setting where all nonzero coordinates of β have a signal strength of exactly the same magnitude, which simplifies our
presentation while preserving the inherent complexity of the support recovery problem. We believe
that studying support recovery in SIM, would bring us insightful understanding of SIR and other
SDR algorithms.
In this paper, we study two procedures for signed support recover of SIM (1): the DT-SIR
introduced by Lin et al. (2015) and the SDP approach inspired by Amini & Wainwright (2008). We
n
let Γ = s log(p−s)
be the rescaled sample size. Our main contribution is to establish the existence of
constants ω > 0 and Ω > 0 such that when Γ > Ω both DT-SIR and SDP approaches recover the
signed support of β correctly, with probability converging to 1 asymptotically. Conversely, we show
that when Γ < ω any algorithm fails to recover the support of β with probability at least 1/2. In
other words, we show that the optimal sample size of the support recovery problem of model (1) is
of the order s log(p − s). To the best of our knowledge, this optimality result, regarding the sample
size of SIM (1), has not been previously discussed in the literature. Our second contribution is,
to establish a sliced stability conjecture formulated by Lin et al. (2015), under the SIM case. We
demonstrate that classical conditions of Hsing & Carroll (1992) imply sliced stability in Section
2.1. This technical result might be of independent interest, especially when one wants to discuss
the optimal rate problem for other SDR algorithms such as SAVE. We further develop a novel tool
2

along the way of our analysis, which may represent further interest — a concentration inequality
stated under Lemma 7.

1.1

Related Work

In the fixed p setting, the first asymptotic results on SIR appeared in the seminal papers (Duan &
Li 1991, Hsing & Carroll 1992). Later on Zhu et al. (2006) allowed p to diverge slowly with n and
established asymptotics in the regime p = o(n1/2 ). In the super high-dimensional setting where
p  n, several algorithms, hinging on regularization such as LASSO (Tibshirani et al. 1997) and
Dantzig Selector (Candes & Tao 2007) were proposed by Li & Nachtsheim (2006), Yu et al. (2013),
but these algorithms are not concerned with support recovery. Moreover, the algorithm suggested
by Li & Nachtsheim (2006) did not come with theoretical guarantees, and in Yu et al. (2013) it is
not allowed for s to scale with p and n. A generic variable selection procedure was suggested in
Zhong et al. (2012), with guarantees of support recovery, in a more general setting than our present
paper, but with a much more restrictive relationship (p = o(n1/2 )) than the one we consider.
In the parallel line of research on sparse PCA there have been more developments. In Johnstone
& Lu (2004) the algorithm Diagonal Thresholding (DT) was suggested, to deal with the spikedcovariance model. It was later analyzed by Amini & Wainwright (2008), who showed that support
recovery is achieved by DT in the sparse spiked covariance model, provided that n & s2 log(p).
Amini & Wainwright (2008) further showed an information theoretic obstruction, in that no algorithm can recover the support of the principal eigenvector if n . s log(p). A computationally
inefficient algorithm that succeeds in support
recovery with high probability as long as n & s log(p)

p
is exhaustively scanning through all s subsets of the coordinates of the principal eigenvector. In
order to find feasible procedures, Amini & Wainwright (2008) studied a semidefinite programming
(SDP) estimator originally suggested in d’Aspremont et al. (2008) — and showed that if n & s log(p)
and the SDP has a rank 1 solution, this solution can recover the signed support with high probability. Surprisingly however, Krauthgamer et al. (2013) showed that the rank 1 condition, does not
hold if s2 log(p) & n & s log(p). In contrast to the PCA case, our current paper argues that if one
is concerned with the support recovery for SIM in the class FA , no computational and statistical
tradeoff exists in the regime s = O(p1−δ ), since the computationally tractable algorithms DT-SIR
and SDP algorithms solve the support recovery problem of SIM with optimal sample size.

1.2

Preliminaries and Notation

We first briefly recall the SIR procedure for SIM. Suppose we observe n = Hm independent and
identically distributed (i.i.d.) samples (Yi , X i ) from model (1). SIR proceeds to sort and divide
the data into H slices of equal size, according to the order statistics Y(i) of Yi . Let the concomitant
of Y(i) be X(i) . We will also use the double subscript notation Xh,i for X((h−1)m+i) . Let Sh =
(Y((h−1)m) , Y(hm) ], 1 ≤ h < H, and SH = (Y((H−1)m) , +∞) denote the random intervals partitioned
by the points Y(jm) (with Y(0) = −∞), j ≤ H − 1. Furthermore let mj (Y ) = E[X j |Y ] denote the j th
coordinate of the centered inverse regression curve E[X|Y ], and µjh denotes E[X j |Y ∈ Sh ]. Note
that conditionally on the values Y((h−1)m) and Y(hm) the quantities Sh and µjh become deterministic.
j
j
1 P
For any integer k ∈ N put [k] = {1, . . . , k} for brevity. Let X h,S = |S|
i∈S X h,i , where S ⊂ [m],
j
j
j
j
1 Pm
1 PH
j ∈ [p]. If S = [m] we omit it from the notation, i.e., X h = m
i=1 X h,i . Put X = H
h=1 X h
for the global average. In terms of this notation, the SIR estimator V of the conditional covariance
3

matrix is given by
V

jk

H
1 X j k
X hX h,
=
H

(3)

h=1

b of β is defined as the principal eigenvector of V. We will denote
entrywise. The SIR estimator β
the support of the vector β by Sβ , i.e. we set Sβ := supp(β) = {j : β j 6= 0}.
Finally, we need some standard notations. For a vector v, let kvkp denote the usual `p
norm for 1 ≤ p ≤ ∞ and kvk0 = | supp(v)|. For a real random variable X, let kXkψ2 =
supp≥1 p−1/2 (E|X|p )1/p , kXkψ1 = supp≥1 p−1 (E|X|p )1/p . Recall that a random variable is called
sub-Gaussian if kXkψ2 < ∞ and sub-exponential if kXkψ1 < ∞. For a matrix M ∈ Rd×d and
two sets S1 , S2 ⊂ [d], by double indexing MS1 ,S2 we denote the sub-matrix of M with entries
Mij for i ∈ S1 , j ∈ S2 . Furthermore, for a d × d matrix Md×d , let kMkmax = maxjk |Mjk | and
kMkp,q = supkvkp =1 kMvkq . In particular, we have kMk2,2 = maxi∈[d] {σi (M)}, where σi (M) is
P
P
the ith singular value of M, and kMk∞,∞ = maxi∈[d] dj=1 |Mij |. Let Fn (x) = n1 ni=1 I(Yi ≤ x)
denote the empirical distribution of the Y sample, and Φ (resp. φ) denote the cdf (resp. pdf) of a
standard normal random variable. We will also occasionally use the abbreviation WLOG to stand
for “without loss of generality”.

1.3

Organization

The paper is organized as follows. The main results, and confirmatory numerical studies and a real
data analysis are presented in Section 2 and Section 3 respectively. Proofs of our main results are
included in Sections 4 and 5 while technical lemmas are deferred to Appendix A. A brief discussion
on the potential directions is included in Section 6. In Appendix B we provide several concrete
examples of SIM.

2

Main Results

In this section, we first prove a conjecture regarding the coordinate-wise sliced stability conditions
introduced by Lin et al. (2015). Second, we establish the optimal rate of support recovery in terms
of the sample size. Throughout the remainder of the paper, we assume that Y is a continuously
distributed random variable.

2.1

Sliced Stability

Our proofs of signed support recovery rely on the following property of the inverse regression curve
of a SIM.
Definition 1 (Sliced Stability). We call the pair (f, ε) sliced stable iff there exist constants 0 <
l < 1, 1 < K, 0 < M , such that for any H ∈ N, H > M , and all partitions of R = {a1 =
−∞, . . . , aH+1 = +∞} with Hl ≤ P(ah < Y ≤ ah+1 ) ≤ K
H there exist two constants 0 ≤

4

κ(l, K, M ) < 1, C(l, K, M ) > 0 such that for all j ∈ Sβ 1 :
H
X

Var[mj (Y )|ah < Y ≤ ah+1 ] ≤ C(l, K, M )H κ(l,K,M ) Var[mj (Y )].

(4)

h=1

The sliced stability assumption, is an implicit assumption on the function f and the error
distribution ε. If κ = 0, the condition means that the cumulative relative variability of the inverse
regression curve is bounded for all slicing schemes with sufficiently small slices. If κ > 0 the
cumulative relative variability of the inverse regression curve is allowed to scale sub-linearly with the
number of slices (see also Figure 1). Intuitively, sliced stability allows us to ensure that the estimates
Vjj of Var[mj (Y )] become increasingly accurate the more slices we introduce (see Lemma 3 for a
more rigorous treatment). Relying on the subsequent developments of this section, in Example
2 and Remark 3 of Appendix B we demonstrate that models of the form Y = G(h(X | β) + ε),
where G, h are continuous and monotone and ε is a log-concave random variable, satisfy the sliced
stability assumption.
Figure 1: Below, with a solid line we plot the standardized inverse regression curve m(y) :=
sign(βj )mj (y)
√
, j ∈ Sβ , for the model Y = 2 atan(X | β) + N (0, 1) (see also (14)). The different colored
Var(mj (Y ))

parts of the dashed curve represent the conditional densities of Y |Y ∈ Sh , where Sh for h ∈ [H]
(with H = 7) are the slices illustrated by punctured vertical black lines. Finally the seven points’
vertical-axis values represent the variances Var[m(Y )|Y ∈ Sh ]. Sliced stability ensures that the
average value of the variances Var[m(Y )|Y ∈ Sh ] has a decay rate of the order H κ−1 .

Y
1

Observe that (4) automatically holds for j 6∈ Sβ in our case, since both the LHS and the RHS of (4) are 0 in this
case.

5

Definition 1 is the sliced stability definition from Lin et al. (2015) restated in terms of the
SIM. Lin et al. (2015) conjectured that the sliced stability condition could be implied from the
well accepted conditions proposed by Hsing & Carroll (1992), which we state below with a slight
modification. The first contribution of this paper is proving this conjecture in the case of SIM.
Let AH (l, K), with 1 < K, 0 < l < 1, denote all partitions of R of the sort {−∞ = a1 ≤
a2 ≤ . . . ≤ aH+1 = +∞}, such that Hl ≤ P(ah ≤ Y ≤ ah+1 ) ≤ K
H . Moreover, for any fixed
B ∈ R, let Πr (B) denote all possible partitions of the closed interval [−B, B] into r points −B ≤
b1 ≤ b2 ≤ . . . ≤ br ≤ B. Define the normalized version of the centered inverse regression curve
sign(β )m (y)
m(y) := √ j j , j ∈ [p]2 , and let m satisfy the following smoothness condition:
Var(mj (Y ))

lim

sup r

−1/(2+ξ)

r
X

r→∞ b∈Π (B)
r

|m(bi ) − m(bi−1 )| = 0,

(5)

i=2

for any B > 0 for some fixed ξ > 0. Note that as mentioned in Hsing & Carroll (1992), assumption
(5) is weaker than assuming that m is of bounded variation, and furthermore the bigger the ξ
the more stringent this assumption becomes. In addition, assume that there exists B0 > 0 and a
e : (B0 , ∞) 7→ R, such that:
non-decreasing function m
e
e
|m(x) − m(y)| ≤ |m(|x|)
− m(|y|)|,
for x, y ∈ (−∞, −B0 ) or (B0 , +∞),

(6)

e
e
and moreover, E[|m(|Y
|)|(2+ξ) ] < ∞ (where in the expectation we set m(y)
= 0 for |y| ≤ B0 ). We
are now in a position to formulate the following:
Proposition 1. Assume that the standardized centered inverse regression curve satisfies properties
(5) and (6) for some ξ > 0. Then we have that for any fixed 0 < l < 1 < K:
lim

sup

H→∞ a∈A (l,K)
H

1
H 2/(2+ξ)

H
X

Var[m(Y )|ah < Y ≤ ah+1 ] → 0.

(7)

h=1

We defer the proof of Proposition 1 to Appendix A. It is clear however that (7) implies the
2
< 1.
existence of constants M , C(l, K, M ) such that (4) holds, with κ = 2+ξ

2.2

Optimal sample size for support recovery

n
Recall that Γ = s log(p−s)
is the rescaled sample size. First we establish an information theoretic
barrier of the support recovery problem, i.e. we show that there exists a positive constant ω such
that, when Γ < ω every algorithm fails with probability at least 1/2. Then we prove that two
algorithms — DT-SIR and an SDP based procedure achieve this bound, i.e. we demonstrate that
there exists a positive constant Ω, such that, when Γ > Ω, these two algorithms successfully recover
the signed support with probability 1 as n → ∞.

2.2.1

Information theoretic barrier

Intuitively, when the sample size is small, no algorithm is expected to be able to recover the support
successfully. In fact, let us consider the simple linear regression model:
Y = X | β + ε where X ∼ Np (0, Ip ), ε ∼ N (0, σ 2 ),
2

By symmetry sign(βj )mj (y) = sign(βi )mi (y) for i, j ∈ Sβ

6

(8)

where β is a unit vector such that β i ∈ {± √1s , 0}. This model belongs to the class F(σ2 +1)−1 (to
verify sliced stability refer to Example 2 and Remark 3 of Appendix B). The following result, which
is obtained by a slight modification of the arguments in (Wainwright 2009), shows that support
recovery is infeasible if Γ is small, by effectively arguing that no algorithm works even for model
(8).
Proposition 2. Suppose there are n observations from the model (8). Then there exists a positive
constant ω, such that if
n
Γ=
< ω,
s log(p − s)
any algorithm for support support recovery will have errors with probability at least
2.2.2

1
2

asymptotically.

Optimality of DT-SIR

Recall that we are working with the class of models FA (2). One key observation is the following
signal strength result.
Lemma 1. For any j ∈ Sβ = supp(β), one has:
A
1
≤ Var(mj (Y )) ≤ ,
s
s

(9)

and Var(mj (Y )) = 0 for j 6∈ Sβ .
By Lemma 1 WLOG we can assume that there exists a CV > 0 such that Var[mj (Y )] = CsV for
any j ∈ Sβ . Now we are ready to discuss the properties of the DT algorithm which we formulate
below. In this paper, we assume that the sparsity s is known.
Algorithm 1: DT algorithm
input: (Yi , X i )ni=1 : data, H: number of slices, s: the sparsity of β
1. Calcluate Vjj , j ∈ [p] – according to formula (3);
b
2. Collect the s highest Vjj into the set S;
b
3. Output the set {j : Vjj ∈ S}.
Recall the definitions (see (3)) of V and Vjj the estimates of Cov[E[X|Y ]] and Var[mj (Y )]
respectively. To obtain the result on sample size optimality of Algorithm 1, a key point is to
establish a concentration inequality for Vjj . When j 6∈ Sβ , a standard deviation inequality for the
χ2 distribution is applicable. For j ∈ Sβ , we need to pay extra effort to obtain the appropriate
deviation inequality, which is also the main technical contribution of this paper. Once we have
established these deviation inequalities, we can show the following theorem.
Theorem 1. Suppose s = O(p1−δ ) for some δ > 0. There exists a positive constant Ω such that,
for any
n
≥ Ω,
(10)
Γ=
s log(p − s)
the support S is recovered by the DT algorithm (i.e., Sb = S) with probability converging to 1 as n
increases. Additionally, the number of slices can be held large enough but fixed (depending solely
on C, l, K, M, κ, CV ).
7

Remark 1 (Choice of H). An interesting by-product of our analysis is that we can choose the
number of slices H large enough but finite. When dimension p is fixed, this has already been
observed in the literature (Li 1991), however in high dimensional setting to the best of our knowledge
this property has not been discussed.
Clearly the DT algorithm does not recover the signed support of β standalone. One can apply
DT first to select the variables and then apply the SIR procedure to obtain an estimate of the
principal eigenvector. To obtain the signed support of β, take the sign of the principal eigenvector.
We summarize this algorithm in the following:
Algorithm 2: DT-SIR
input: (Yi , X i )ni=1 : data, H: number of slices, s: the sparsity of β
b
1. Perform Algorithm 1 to obtain the set {j : Vjj ∈ S}.
2. Evaluate v — the principal eigenvector of the matrix VS,
bS
b.
3. Output sign(v).
We note that this algorithm recovers the signed support, up to multiplication by ±1. The
following result is a direct corollary of Theorem 1, whose proof is omitted.
Proposition 3. Under the assumptions of Theorem 1, with a potentially bigger value of Ω, applying
Algorithm 2 restores the signed support (up to a sign) with asymptotic probability converging to 1.
b where the jk th entry of
In fact Algorithm 2 works when X ∼ Np (µ, Ip ) with V substituted with V
P
b is defined as V
b jk = 1 H (X j − X j )(X k − X k ).
V
h
h
h=1
H
2.2.3

Optimality of SDP

Algorithm 2 is a two-step procedure, selecting variables by applying DT and then obtaining the
SIR estimate of β. Recent advances of optimization theory provide us with a more sophisticated
one-step approach to obtain a sparse principal eigenvector. It is well known that the principal
eigenvector of a symmetric positive definite matrix A is given by:
b
z = argmax z| Az
z∈R:kzk2 =1

and the principal eigenvalue is b
z| Az. When z is sparse, i.e., kzk0 ≤ s , the above optimization problem with an additional sparsity constraint is computationally expensive. To remedy this difficulty,
d’Aspremont et al. (2007) proposed an SDP approach to solve the sparse optimization problem.
They suggested to solve the following convex program:
b=
Z

argmax
tr(Z)=1,Z∈Sp+

tr(AZ) − λn

p
X

|Zij |,

(11)

i,j=1

where Sp+ is the set of all the p × p positive semi-definite matrices. If the solution happens to be a
b =b
zb
z| , and hence we can easily obtain an estimate of the
rank 1 solution, then it is of the form Z
principal eigenvector.

8

In the following algorithm we summarize the SDP approach, tailored for the signed support
recovery of SIM. We remind the reader, that this algorithm recovers the signed support, up to
multiplication by a global constant equal to ±1.
Algorithm 3: SDP algorithm for SIR
input: (Yi , X i )ni=1 : data, H: number of slices, s: the sparsity of β
1. Calcluate the matrix V — as given in (3);
b by solving (11), with A = V;
2. Obtain the matrix Z
b
3. Find the principal eigenvector b
z of Z;
4. Output sign(b
z).
The performance of this algorithm is guaranteed by the following theorem.
Theorem 2. Suppose s = O(p1−δ ) for some δ > 0, then there is a positive constant Ω such that,
for any
Γ=

n
≥ Ω,
s log(p − s)

(12)

with a properly chosen the tuning parameter λn , Algorithm 3 recovers the signed support with
probability converging to 1, i.e. P(sign(b
z) = sign(β)) → 13 .
b where the
Remark 2. In fact Algorithm 3 works when X ∼ Np (µ, Ip ) with V substituted with V
P
j
j
k
k
H
b is defined as V
b jk = 1
jk th entry of V
h=1 (X h − X )(X h − X ). The proof of this fact is trivial
H
and is omitted.

3

Numerical Experiments and Data Analysis

We open this section with extensive numerical studies and in Section 3.2 we apply our algorithms
to a real dataset.

3.1

Simulations

In this section we compare Algorithms 2 and 3 in terms of signed support recovery. We consider
the following scenarios:
Y = X | β + sin(X | β) + N (0, 1),
|

Y = 2 atan(X β) + N (0, 1),
|

3

Y = (X β) + N (0, 1),
|

Y = sinh(X β) + N (0, 1).

(13)
(14)
(15)
(16)

√
In each simulation we create a sample of size n with dimensionality p and sparsity levels s = p
(s = log(p) resp.) and we vary the rescaled sample size Γ ∈ [0, 30] for the DT-SIR and Γ ∈ [0, 40]
3

We understand that b
z is selected so that b
z| β ≥ 0 by convention.

9

for the SDP approach. The vector β is selected in the following manner:
1
1
βj = √ , j ∈ [s − 1], βs = − √ , and βj = 0, for j > s.
s
s
Each simulation is repeated 500 times, and we report the proportion of correctly recovered signed
supports of the vector β. We remark that in general, it should not be expected that the phase
transition described in Theorems 1 and 2 to occur at the same places for these 4 models.
We first explore the predictions of Proposition 3 and Algorithm 2. Even though we provide
theoretical values of the constants H and m, we ran all simulations with H = 10 slices. We believe
this scenario, is still reflective of the true nature of the DT-SIR algorithm, as the theoretical value
of H we provide is not optimized in any fashion. In Figure 2, we present DT-SIR results from plots
√
for different p values in the regime s = p.

●

●

●

●

●

●

●

√

n
s log(p−s)

p, Γ =

∈ [0, 30]

1.0

1.0

Figure 2: Efficiency Curves for DT-SIR, s =

●

●

●

●

●

●

●

●

●

●

●

p = 100
p = 200
p = 300
p = 600
p = 1200

● ●

●

●

0

5

10

15
n

20

25

0.8
0.4

0.6

●

●

●

30

● ●

●

●

0

5

10

15
n

slog(p − s)

●

●

●

●

●

1.0

1.0

●

●

●

●

●

●

●

●

●

●

●

15
n

20

25

0.8
0.6

●

0.4

^
P(supp(sign(β)) = supp(sign(β)))

●
●

0.0

●
●

p = 100
p = 200
p = 300
p = 600
p = 1200

●

0.2

0.8
0.6
0.4
0.2

●
●

0.0

^
P(supp(sign(β)) = supp(sign(β)))

●

10

30

●

●

5

25

(b) Model (14)

●

0

20

slog(p − s)

(a) Model (13)

● ●

p = 100
p = 200
p = 300
p = 600
p = 1200

●

0.0

●

●

0.2

^
P(supp(sign(β)) = supp(sign(β)))

0.8
0.2

0.4

0.6

●

0.0

^
P(supp(sign(β)) = supp(sign(β)))

●

30

● ●

0

slog(p − s)

●

●

5

10

15
n

20

slog(p − s)

(c) Model (15)

(d) Model (16)

10

25

p = 100
p = 200
p = 300
p = 600
p = 1200
30

n
On the X-axis we have plotted the rescaled sample size s log(p−s)
and on the Y -axis is the
estimated probability of successful signed support recovery. We would refer to these curves as
√
efficiency curves (EC). The EC plots in the case s = p, are very similar to the corresponding
plots in the regime s = log(p), which can be seen in Figure 4 in Appendix C. Both EC plots are
in concordance with the predictions from our theoretical results. We can distinctly see the phase
transition occurring in approximately the same place regardless of the values of the dimension p.
Next, we present the corresponding ECs for Algorithm 3. We used the code from an efficient
implementation of program (11), as suggested in Zhang & Ghaoui (2011). The code was kindly
provided to us by the authors of Zhang & Ghaoui (2011). In Figure 3 we plot the ECs for the four
√
models in the case when s = p. Due to running time limitations we only show scenarios where
p ∈ {100, 200, 300}.

●
●

●

●

●

●

●

√

p, Γ =

n
s log(p−s)

∈ [0, 40]

1.0

1.0

Figure 3: Efficiency Curves for SDP, s =

●

●

●

●

●

●

●

●

●

●

●

●
●

●

●●

●

●

●

0

p = 100
p = 200
p = 300

●

●

10

20
n

30

0.8
0.6

●

●

40

●●

●

●

●

0

10

20
n

slog(p − s)

(b) Model (14)

●
●

●

●

●

●

1.0

1.0

●

●

●

●

●

●

●

●

●

●

●

20
n

30

p = 100
p = 200
p = 300

0.8

●

0.6

●

●

●

●

0.0

●
●
●

●

0.4

^
P(supp(sign(β)) = supp(sign(β)))

●

0.2

●

●

0.2

0.8
0.4

0.6

●

0.0

^
P(supp(sign(β)) = supp(sign(β)))

●

10

●
●

●

0

40

●

●

●

30

slog(p − s)

(a) Model (13)

●●

p = 100
p = 200
p = 300

●
●

0.0

0.0

●

●

●

0.4

0.4

●

●

0.2

^
P(supp(sign(β)) = supp(sign(β)))

0.8
0.6

●

0.2

^
P(supp(sign(β)) = supp(sign(β)))

●

40

●●

0

slog(p − s)

●

●

●

10

20
n
slog(p − s)

(c) Model (15)

(d) Model (16)

11

30

p = 100
p = 200
p = 300
40

Here we have again used H = 10 in all scenarios, for simplicity. We observe that phase transitions are occurring in all of the models, and the signed support is being correctly recovered for
n
large enough values of Γ = s log(p−s)
. Plots for the setting s = log(p) are provided in Figure 5.
An important empirical observation is that the constant Ω needs to be higher when running the
SDP algorithm compared to when running the DT-SIR. This fact coupled with the much slower
run-times of the SDP algorithm speak in favor of using the DT-SIR algorithm in practice.
In addition to the above simulations, we also performed numerical studies over the same set of
models (13), (14), (15) and (16), where the coefficients of the vector β were not of equal magnitude.
e
e were drawn in each
In our second set of simulations we chose β = β , where the entries of β
e 2
kβk

simulation as
β̃j = Uj for j ≤ bs/2c, β̃j = −Uj for s/2 < j ≤ s, and β̃j = 0 otherwise,

(17)

where Uj for j ∈ [s] are i.i.d. uniform U ( 12 , 1) random variables. The results are similar to the
ones reported above, although the efficiency curves required somewhat higher values of the rescaled
sample size Γ, and extra variability can be seen in the curves due to the random choice of β. This is
to be expected, since under such a generation scheme, the lowest signal will be of lower magnitude
√
compared to the uniform signal case. We show the results for the case s = p in Appendix C under
Figures 6 and 7. Due to space considerations, we do not show results for the case s = log p, as the
results are comparable to the results plotted in Figures 4 and 5.

3.2

Real Data Example

In this sub-section we will illustrate the practical benefits of our algorithm on a real dataset. Our
dataset is the mouse embryonic stem cells (ESCs) dataset, and has been previously analyzed by
Zhong et al. (2012), Jiang et al. (2014).
The outcome variable is the expression levels of 12, 408 genes, and is derived using RNA-seq
technology in mouse ESCs (Cloonan et al. 2008). The predictors correspond to 312 transcription
factors (TFs). From the predictors, 12 TFs are known to be associated with different roles in
the ES-cell biology as constituents of important signaling pathways, self renewal regulators, and
key reprogramming factors, and the remaining 300 are putative mouse TFs compiled from the
TRANSFAC database. For each gene and TF, the predictor matrix contains a score which is the
transcription factor association strength score (Chen et al. 2008, Ouyang et al. 2009) for the 12
known TFs, and motif matching scores for the remaining 300 genes which were introduced by Zhong
et al. (2012). Hence, the design matrix X = (X |i )|i∈[n] is a n × p matrix, where n = 12, 408 and
p = 312, each entry of which represents the score of a TF for the corresponding gene. The EFCs
dataset is not truly a high-dimensional dataset in the sense n < p, but serves to illustrate that our
algorithms can successfully perform variable selection. Since X is coming from a real dataset with
non-homogeneous scores across the predictors, we wouldn’t expect it to be centered with identity
e jk in place of Vjk where the matrix
covariance. To deal with this problem we use the statistics V
e
V is defined as:
H
X
e := Σ
b −1/2 1
b −1/2 ,
V
(X h − X)(X h − X)| Σ
H
h=1

|
b −1/2 := n−1 P
and Σ
i∈[n] (X i − X)(X i − X)

h

i−1/2

12

is the symmetric square root of the covariance

b −1/2 can be estimated without the application of sparse
matrix of X. Notice that in this dataset Σ
inverse procedures since n > p. As in the simulation section we used H = 10 slices.
We ran the DT and SDP procedures (i.e. Algorithms 1 and 3) on the ESCs dataset, and the DT
procedure selected 28 TFs which it found associated with the outcome, while the SDP procedure
selected 36. In Table 1 we compare the rankings of the 12 known ES-cell TFs within the TFs
selected by DT and SDP algorithms to the same rankings produced by the procedures SIRI-AE
(Jiang et al. 2014) and COP (Zhong et al. 2012).
Table 1: Rankings of the 12 known TFs among the selected TFs by different algorithms
TFs names
E2f1
Zfx
Mycn
Klf4
Myc
Esrrb
Oct4
Tcfcp2l1
Nanog
Stat3
Sox2
Smad1

DT
1
2
3
7
4
14
20
8
11
–
24
25

SDP
1
3
2
4
7
10
20
8
17
23
–
25

SIRI-AE
1
3
4
5
6
8
9
10
14
17
18
32

COP
1
3
10
19
–
–
11
36
–
20
–
13

The top 10 highly ranked genes of SIRI-AE, DT, SDP and COP contain 8, 6, 7, 3 known ES-cell
TFs respectively. The DT algorithm is the most parsimonious out of all procedures (SIRI-AE and
COP have reportedly selected 34 and 42 TFs correspondingly), but nevertheless includes all of the
known TFs with the exception of the Stat3. The DT algorithm is additionally able to capture both
Nanog and Sox2 which are generally believed to be the master ESC regulators. These TFs are
completely missed by COP and the Sox2 TF was omitted by SDP. The SIRI-AE algorithm is the
only algorithm which is able to recognize all of the 12 known TFs. We would like to mention that
compared to our work the SIRI-AE algorithm is designed to work in a lower dimensional setup,
and is able to capture interactions beyond the single-index models, which could help explain its
better performance on the ESCs dataset.
All in all, we believe the results of our algorithms are satisfactory, and illustrate that both DT
and SDP can be successfully applied in practice. As a remark, in a truly high-dimensional setting
b −1/2 , in practice one might resort to sparse estimators
when n < p, instead of the estimate Σ
e in
produced by procedures such as CLIME (Cai et al. 2011), and use the corrected statistics V
place of V.

4

Proof of Theorem 1

We will prove this theorem under slightly more general conditions:
i. Let X j , j ∈ [p] be centered, sub-Gaussian random variables with maxj∈[p] kX j kψ2 ≤ K.

13

ii. For j ∈ Sβ , we assume that Var[mj (Y )] ≥

CV
s

.

iii. For j ∈ Sβc , we assume that mj (Y ) = E[X j |Y ] = 0 a.s.
It is easy to see for that for SIM in FA , the above conditions are satisfied in the case when
X ∼ Np (0, Ip ). We start with a high level outline of the proof. As we pointed out, the key
argument is to show a deviation inequality for Vjj , i.e. we will show that |Vjj − Var[mj (Y )]| is
n
small (e.g. < 31 Var[mj (Y )]) with high probability provided that Γ = s log(p−s)
is large enough. Note
that:
Z
H
1 X  j 2
jj
(18)
|V − Var[mj (Y )]| =
X h − m2j (y)pY (y)dy .
H
h=1

We will show below (see (22)), that the above expression is well approximated by:
H
H
1 X  j 2 X j 2
Xh −
(µh ) P(Y ∈ Sh ) ,
H
h=1

(19)

h=1

j

under sliced stability (where Y(0) = −∞). Intuitively, we have 1) P(Y ∈ Sh ) ≈ H1 and 2) X h ≈ µjh ,
which shows that (19) should be close to 0. To rigorously validate this intuition, we use the
Dvoretzky-Kiefer-Wolfowitz (DKW) inequality (Massart 1990) for 1) and concentration inequalities
for 2).
Note that the probability P(Y ∈ Sh ) is a random variable, where the randomness comes from the
two endpoints — Y(m(h−1)) and Y(mh) of the interval Sh . Recall that Fn is the empirical distribution
function of Y , based on the sample Yi . By the DKW inequality, we have that P(supy |Fn (y)−F (y)| >
) ≤ 2 exp(−2n2 ), which in addition to the fact that Y has a continuous distribution, implies that
for all h we have:




1
h−1
h
h−1
h
1
− 2 ≤ P
< Fn (Y ) <
≤P
≤ Fn (Y ) ≤
≤
+ 2,
(20)
H
H
H
H
H
H
on an event with probability at least 1 − 2 exp(−2n2 ). Denote the event where (20) holds with E.
Next, notice that for h ∈ [H − 2] and a random permutation π : [m − 1] 7→ [m − 1], the
collection random variables X jh,πi , i ∈ [m − 1] are conditionally i.i.d. given (Y(m(h−1)) , Y(mh) ), and

for h = H − 1 and a random permutation π : [m] 7→ [m] the random variables X jh,πi , i ∈ [m] are
conditionally i.i.d. given Y(m(H−1)) , as they can be generated via simple rejection sampling. Hence
j

conditionally on (Y(mh) )h∈[H−1] the sample means X h,1:(m−1) have corresponding means µjh when
j

h ∈ [H − 1] and X H has a mean of µjH . For consistency of notation we will discard a point at
random in the H th slice, and WLOG (upon re-labeling the points) we assume that the discarded
point is X H,m . Next, we formulate the following concentration result for the sliced means, which
we show in Appendix A:
Lemma 2. Let G ⊆ [p]. On the event E, for η > 0 we have the following:




η 2 (m − 1)
j
P
max X h,1:(m−1) − µjh > η ≤ 2|G|H exp −
,
C1 K2 q −1 + C2 Kη
j∈G,h∈[H]
where q =

1
H

− 2, for some absolute constants C1 , C2 > 0.
14

(21)

e the event on which we have
Set G = [p] and denote with E
j

max
j∈[p],h∈[H]

X h,1:(m−1) − µjh ≤ η.

By (21), (20) and the union bound we have that:


η 2 (m − 1)
e
− 2 exp(−2n2 ).
P(E) ≥ 1 − 2pH exp −
C1 K2 q −1 + C2 Kη
Next we move on, to show that (18) is close to (19) on the event E, as well as we collect two
straightforward inequalities in the following helpful:
Lemma 3. Assume that the sliced stability condition (4) holds. Then we have the following inequalities holding on the event E, for large enough H, and small enough :
H
H
1 X  j 2 X j 2
Xh −
(µh ) P(Y ∈ Sh )
H
h=1
h=1


CH κ 1
+ 2 ,
+
s
H
|
{z
}

|Vjj − Var[mj (Y )]| ≤

(22)

B1

H
X

(µjh )2 ≤

h=1

CV
s
1
H

|
H
X

+ B1
,
− 2
{z }

(23)

B2

q
|µjh |

CV
s + B1

1
H − 2

≤

h=1

|

{z

B3

.

(24)

}

Note. We refer to the constants from (4) as C and κ, dropping the dependence on K and M for
brevity, and in fact C = C(l, K, M )CV .
e we get:
Note that by an elementary calculation — using (20) and Lemma 3, on the event E
H
1 X  j 2 (m − 1)2 j 2
Xh −
(µh )
|V − Var[mj (Y )]| ≤
H
m2
h=1


1
(m − 1)2
+ 2 +
−
B2 + B1 ,
|{z}
H
Hm2
|
{z
}
I2
jj

I1

where we used (22), the triangle inequality and (23). Consider the following:

15

(25)
(26)

e
e
e ⊂ E
e such that P(E
e\E
e ) ≤ p exp(−c τ 2 n4 ), for some fixed
Lemma 4. There exists a subset E
4K
τ ∈ [0, 2K2 ) on which we have the following bound for all j ∈ [p]:
r
√
H
1 X  j 2 (m − 1)2 j 2
B2
(2K2 + τ ) 2 2K2 + τ
√
Xh −
(µh ) ≤
2η 2 + 2
+
2
H
m
m }
H
m
| {z
h=1
|
{z
}
I3

(27)

I4

B3
.
+ η 2 + 2η
|{z} | {zH}
I5

I6

We defer the proof of this lemma to the appendix. Next, we provide exact constants, such that
− Var[mj (Y )]| ≤ CρsV for some constant ρ > 0 and all j ∈ [p], so that the probability of the
e
e still converges to 1. The remarkable phenomenon here is that the number of slices H, can
event E
be selected so that it is a constant, which might seem counterintuitive. Select the constants in the
following manner:
( 
 1 )
γCK 1−κ
H = max M,
,
(28)
CV


K −1 1−l
l
 = min
,
,
,
(29)
2H
2H 4(γ + 1)H
√
l CV
(30)
η= p
√ ,
2 γ(γ + 1) s


8K4 log(p − s)
8γ(γ + 1 + l)(2K2 + τ )s 4(γ + 1) 2s log(p − s)
m ≥ max
,
,
+ 1,
(31)
lCV
l
Υ
τ 2 cH
|Vjj

where γ ≥ 12 is a positive constant and Υ =

l2 CV

.
√
C2 Kl CV
4(γ(γ+1)) K2 C1 Hl−1 + √
2

Simple algebra shows

γ(γ+1)

that selecting these constants ensures the following inequality:
max(I1 , I2 , I3 , I4 , I5 , I6 ) ≤

CV
.
γs

By combining (25) and (27) we arrive at:
|Vjj − Var[mj (Y )]| ≤

6CV
,
γs

(32)

e
e converges
as promised for ρ = γ6 . Next we proceed to show that the probability of the event E
to 0. To achieve this, note that in (29) we chose  to be a constant, so evidently exp(−n2 ) → 0.
2
e
p
e
Moreover p exp(−c τ mH
4 ) ≤
2 → 0 by (31). Finally to show that P(E ) → 1, we need:
4K

(p−s)

η 2 (m − 1)
− log(p) − log(H) → +∞.
C1 (H −1 − 2)−1 K2 + C2 Kη

16

√

Noting that H −1 − 2 ≥ lH −1 , η ≤ √l
2

CV
γ(γ+1)

and that by (28) H is fixed, the above expression is

implied by:
(m − 1)
− log(p) ≥ 2 log(p − s) − log(p) → +∞,
s
where we used (31) and the fact that 2 log(p − s) ≥ log(p) asymptotically (since s = O(p1−δ ) for
e
e ) → 1.
δ > 0). Hence we have guaranteed that P(E
Now, note that for variables j 6∈ Sβ we have mj (Y ) = E[X j ] = 0, which implies that for all
h ∈ [H], j ∈ Sβc we have µjh = 0. Using (32) and selecting γ = 18, we conclude that:
Υ

Vjj ≥

2 CV
CV
, j ∈ Sβ and Vjj ≤
, j ∈ Sβc ,
3 s
3s

and hence separation of the signals is asymptotically possible, as we claimed.

5

Proof of Theorem 2

In this section we show that the SDP relaxation will have a rank 1 solution with high probability
and moreover this solution will recover the signed support of the vector β. In contrast to Section
4, here we make full usage of the fact that X ∼ Np (0, Ip ). One simple implication of this, is for
example the fact that maxj kX j kψ2 < 1, and hence we can set K = 1. For the analysis of the
algorithm we set the regularization parameter λn = C2sV .
To this end, we restate Lemma 5 from Amini & Wainwright (2008), which provides a sufficient
condition for a global solution of the SDP problem:
Lemma 5. Suppose there exists a matrix U satisfying:
(
sign(b
zi ) sign(b
zj ), if b
zi b
zj 6= 0;
Uij =
∈ [−1, 1],
otherwise.

(33)

Then if b
z is the principal eigenvector of the matrix A − λn U, b
zb
z| is the optimal solution to problem
(11).
Recall that the SIR estimate of the variance-covariance matrix has entries:
Vjk =

H
1 X j k
X hX h.
H
h=1

e = V − λn U, where U is a to be defined sign matrix from Lemma 5. We furtherDenote with V
e into blocks — V
e S ,S , V
e S c ,S , V
e S c ,S c . Here, these three
more consider the decomposition of V
β β
β β
β β
e
matrices are sub-matrices of the matrix V restricted to entries with indexes in the sets Sβ or Sβc
correspondingly. We observe that USβ ,Sβ = sign(β Sβ ) sign(β Sβ )| .
We first focus on the VSβ ,Sβ matrix. We calculate the value of the covariance of two coordinates
j, k ∈ Sβ :
Cov[mj (Y ), mk (Y )] = E[mj (Y )mk (Y )] = sign(β j ) sign(β k )E[m2k (Y )]
= β j β k CV ,
17

(34)

where we used that sign(β j )mj (Y ) = sign(β k )mk (Y ), which follows by noticing that the distribution
of X j |Y is the same as the distribution of X k |Y except the potential difference in the signs of the
coefficients, because of the symmetry in the problem.
Next, let G = Sβ (hence |G| = s) in Lemma 2 to obtain that:
max
j∈Sβ ,h∈[H]

j

X h,1:(m−1) − µjh ≤ η,

(35)



2 (m−1)
eS ⊂ E, be the event
− 2 exp(−2n2 ). Let E
with probability at least 1 − 2sH exp − Cη1 q−1
β
+C2 η
where (35) holds. We proceed with formulating a bound similar to Lemma 4, but for the covariance:
e
e
τ 2n
e
eS \ E
e S with P(E
Lemma 6. There exists an event E
Sβ ) ≤ s exp(−c 4 ) (c > 0 is an absolute
β
β
constant), for some fixed τ ∈ (0, 2], such that for all j, k ∈ Sβ , we have the following inequality:


1
(m − 1)2
B3
jk
V − Cov(mj (Y ), mk (Y )) ≤ 2 +
−
B2 + B1 + 4η
2
H
Hm
H
r
√
4(2 + τ ) 4 2 + τ
B2
+
+ √
+ 4η 2 .
(36)
2η 2 + 2
m
H
m
Let H 0 , 0 , η 0 are constants selected according to (28) — (30), with K = 1, correspondingly. Set
0
H = H 0 ,  = 0 , η = η2 , and take


8 log(p − s)
8γ(γ + 1 + l)(2 + τ )s 8s log(p − s)
,
+ 1,
,
(37)
m ≥ max 32
lCV
Υ
τ 2 cH
where Υ is specified as before with K = 1. Similarly to Section 4, we can show that the following
inequality:
sup Vjk − Cov(mj (Y ), mk (Y )) ≤
j,k∈Sβ

6CV
,
γs

(38)

e
e
e S tending to 1. To get the bound in (38), one can observe
e S , with the probability of E
holds on E
β
β
that with the choices of constants as above all 6 terms in (36) are guaranteed to be smaller than
CV
γs . Here similarly to Section 4, H is large enough but fixed. Having in mind the above inequality
e S ,S :
we consider the matrix V
β

β

e S ,S = CV β S β | + N,
V
β β
β Sβ
2
where N is some symmetric noise matrix. Note that by (34), (38) gives a bound on kNkmax . Next
we make usage of Lemma 8, which is a restatement of Lemma 6 in Amini & Wainwright (2008) and
can be found in the Appendix for the reader’s convenience. We start by verifying that condition
(53) indeed holds for the matrix N. We have that:
kNkmax = kVSβ ,Sβ − CV β Sβ β |Sβ kmax ≤

6CV
,
γs

with the last inequality following from (38). Selecting γ = 240 bounds kNkmax by
in (53). Thus, by Lemma 8 we conclude that:
18

(39)
CV
40s ,

as required

e and the second largest in magnitude eigenvalue of V
e we have γ1 > |γ2 |.
(a) For γ1 := λmax (V)
e —e
(b) The corresponding principal eigenvector of V
zSβ satisfies the following inequality:
e
z Sβ − β Sβ

∞

1
≤ √ .
2 s

Next we show that the rest of the sign matrix U, i.e. USβc ,Sβ and USβc ,Sβc can be selected in such
e S c ,S c are 0. For this purpose we select US c ,S = 1 VS c ,S
e S c ,S and V
a way, so that the blocks V
λn
β β
β β
β β
β β
and USβc ,Sβc = λ1n VSβc ,Sβc . Since it is clear that the vector (e
zSβ , 0Sβc ) is the principal eigenvector
e if U is a sign matrix, Lemma 5 will conclude that — (e
z| , 0| c ) is the optimal
of V,
z| , 0| c )| (e
Sβ

Sβ

Sβ

Sβ

solution to the optimization problem, which will in turn conclude our claim.
It remains to show that the specified U is indeed a sign matrix. Note that by Cauchy-Schwartz
for k ∈ Sβc and any j, we have:
√
√
Vjk ≤ Vjj Vkk .
(40)
From (38) if j ∈ Sβ , we have that high probability:
Vjj ≤

6CV
(γ + 6)CV
CV
+
=
.
s
γs
γs

Hence, it is sufficient to select m, H large enough so that:
Vkk ≤

γCV
.
4(γ + 6)s

To achieve the above bound, we make usage of the following tail inequality, for χ2 random variables
which we take from Laurent & Massart (2000) (see Lemma 1):
r
 2

χH
x
2x
P
≥1+2
+
≤ exp(−x).
H
H
H
Note that, Vkk ∼

1
2
mH χH

for k ∈ Sβc . Thus applying the bound above we have
r
1 2
1
2
x
2x
2
3x
χH ≤
+
+
≤
+
,
mH
m m H
mH
m mH

(41)

with probability at least exp(−x). Applying (41) it can be easily seen that by selecting:
x=
we can ensure (after using (37)) that Vkk ≤

nγCV
,
24(γ + 6)s
γCV
4(γ+6)s

for all k ∈ Sβc . Requiring:

nγCV
≥ 2 log(p − s),
24(γ + 6)s

(42)

ensures that the probability of the event is asymptotically 1 from the union bound. This combined
with (40) shows that the so defined matrix U is indeed a sign matrix, which concludes the proof.
19

6

Discussion

Sliced inverse regression has been widely applied in various scientific problems. Though it is a
successful tool in terms of data visualization, and provably works when the dimension p is not
large, its behavior in the high-dimensional regime p  n is much less well understood.
In this paper, we studied the support recovery of SIM within the class of models FA . We
demonstrated that the optimal sample size of this problem is of the order s log(p − s). Two
unforeseeable results of our analysis might be of particular interest for future investigations.
Recall that a central subspace of a pair of random variables (Y, X) is the minimal subspace S
such that Y ⊥
⊥ X PS X. The first implication of our results, as we hinted in Remark 1, is that
if we focus on estimating the support of the central subspace rather than consistently estimating
the intermediate matrix, better convergence rate might be expected. For instance, the results in
Lin et al. (2015) show that,
q under mild conditions, the convergence rate of the SIR estimate V of

Var(E[X|Y ]) is OP ( H1 + np ). In order to get a consistent estimation of the intermediate matrix
Var(E[X|Y ]), the slices number H has to be diverging. The result regarding the choice of H in
this paper suggests that it would be possible
q to get an improved convergence rate of estimating the
central subspace which is expected to be np .
The second consequence is that estimating the central subspace might be easier than that of
sparse PCA in terms of computational cost. At first glance, the unknown nonlinear link function f
could bring in difficulties in determining the optimal rate, and it might be reasonable to expect that
this difficulty will increase the computational cost in general. The results in this paper however,
provide the opposite evidence — as long as s = O(p1−δ ) for some δ > 0, the computationally efficient
algorithms DT-SIR and SDP approach can solve the support recovery problem. In other words,
the tradeoff between computational and statistical efficiency for estimating the central subspace
might occur in a more subtle regime, where s ∝ p and n ∝ p. Finally combining our results with
the results of Lin et al. (2015), calculating the minimax rate of an estimator for the SDR space in
model (1) appears to be a plausible exercise, and is left for future work.

Acknowledgments
This research was conveyed while the first author was a graduate student in the Department of
Biostatistics at Harvard University. Lin’s research is supported by the Center of Mathematical
Sciences and Applications at Harvard University. Liu’s research is supported by the NSF Grant
DMS-1120368 and NIH Grant R01 GM113242-01. We thank Professor Tianxi Cai, Harvard University and Professor Noureddine El Karoui, University of California, Berkeley for valuable discussions
which led to improvements of the present manuscript. In addition, the authors are grateful to Professor Laurent El Ghaoui, University of California, Berkeley and Dr Youwei Zhang for providing the
Matlab code used for solving the SDP problem. The authors also express gratitude to the referee
and anonymous reviewers for raising important points, which further bettered the manuscript.

A

Technical Proofs

Proof of Lemma 1. First take any vector b ∈ Rp such that b ⊥ β. We have:
E[X|Y ]| b = E[X | b] = 0,
20

since X | b is independent of X | β and ε. This implies that
E[X|Y ] = c(Y )β,

(43)

for some real valued function c. Since β is a unit vector it follows that c(Y ) = E[X | β|Y ]. Take
j ∈ Sβ and apply (43) to get:
Var[mj (Y )] =

Var[E[X | β|Y ]]
.
s

(44)

Combining the observation above with the following two inequalities:
A ≤ Var(E[X | β|Y ]) ≤ Var(X | β) = 1,
gives the desired result.
Proof of Proposition 1. We first note that without loss of generality we can consider the function
e to be non-negative, at the price of potentially shrinking the interval (B0 , ∞) to (B0 + η, ∞) by
m
e 0 (x) = m(x)
e
e 0 + η) for x ∈ (B0 + η, ∞).
any η > 0. To see this fix an  > 0, and define m
− m(B
Then since (6) holds on (−∞, −B0 ) ∪ (B0 , ∞), clearly:
e 0 (|x|) − m
e 0 (|y|)|, for x, y ∈ (−∞, −B0 − η) or (B0 + η, +∞).
|m(x) − m(y)| ≤ |m
2+ξ + |m(B
e 0 (x)2+ξ ≤ 21+ξ (|m(x)|
e
e 0 + η)|2+ξ )
By the convexity of the map x 7→ x2+ξ we have m
0
(2+ξ)
0
e (|Y |)|
e is non-negative and non-decreasing on
and hence E[|m
] < ∞. Finally by definition m
(B0 + η, ∞).
Next note that if Y has a bounded support, this proposition clearly follows from assumption
(5) alone. Thus, without loss of generality we assume that Y has unbounded support (from both
sides, as if one of them is bounded we can handle it in much the same way as the proof below).
e0 = B0 + η, for some small fixed η > 0. Fix any partition a ∈ AH (l, K). Let S0 = {h :
Let B
e0 , B
e0 ]}, and let hm = min S0 , hM = max S0 . Note that the following simple inequality
ah ∈ [−B
holds for any 2 ≤ h ≤ hm − 2 or hM + 1 ≤ h ≤ H − 1:

Var[m(Y )|ah < Y ≤ ah+1 ] ≤

inf
t∈(ah ,ah+1 ]

≤

E[(m(Y ) − m(t))2 |ah < Y ≤ ah+1 ]

sup
y,t∈(ah ,ah+1 ]

e h |) − m(|a
e h+1 |))2 .
(m(y) − m(t))2 ≤ (m(|a

This gives us the following inequality:
hX
m −2

Var[m(Y )|ah < Y ≤ ah+1 ] ≤

h=2

hX
m −2

e h |) − m(|a
e h+1 |))2
(m(|a

(45)

h=2

e 2 |) − m(|a
e hm −1 |))2 ,
≤ (m(|a
e is non-decreasing. Similar inequality holds for the other tail
where the last inequality holds since m
as well.

21

e0 , B
e0 ]:
Using a similar technique we obtain the following bound on the interval: [−B
hX
M −1

Var[m(Y )|ah < Y ≤ ah+1 ] ≤

h=hm

hX
M −1

E[(m(Y ) − m(ah ))2 |ah < Y ≤ ah+1 ]

h=hm

≤

hX
M −1

(m(y) − m(ah ))2 .

sup

h=hm y∈(ah ,ah+1 ]

Notice further that:
Var[m(Y )|ahm −1 < Y ≤ ahm ] ≤

e0 ))2
(m(y) − m(−B

sup
y∈(ahm −1 ,ahm ]

≤

e0 ))2 +
(m(y) − m(−B

sup

e0 ))2 .
(m(y) − m(−B

sup
e0 ,ah ]
y∈[−B
m

e0 ]
y∈(ahm −1 ,−B

And a similar inequality holds for Var[m(Y )|ahM < Y ≤ ahM +1 ]. Thus:
hM
X

Var[m(Y )|ah < Y ≤ ah+1 ] ≤

hX
M −1

h=hm y∈(ah ,ah+1 ]

h=hm −1

|
+

{z

|

}

I2

e0 ))2
(m(y) − m(−B

sup
|

e0 ))2 +
(m(y) − m(B

sup

}

I1

e0 ,ah ]
y∈[−B
m

e0 ]
y∈(ahm −1 ,−B

e0 ,ah +1 ]
y∈[B
M

|

{z

e0 ))2 +
(m(y) − m(−B

sup

+

(m(y) − m(ah ))2

sup

{z

}

I3

e0 ))2 .
(m(y) − m(B

sup
e0 ]
y∈(ahM ,B

{z

}

I4

|

{z
I5

}

We have:
2|S0 |+3

I1 + I3 + I5 ≤

sup
e0 )
b∈Π2|S0 |+3 (B

X

(m(bi ) − m(bi−1 ))2

i=2



2|S0 |+3

≤

sup

(46)

X


e0 )
b∈Π2|S0 |+3 (B

2
|m(bi ) − m(bi−1 )| .

(47)

i=2

e0 , b3 = ah , . . . , b2|S |+1 = ah , b2|S |+3 =
To see this, consider a partition containing the points b1 = −B
m
M
0
0
2
e
B0 , and b2k = argmaxy∈(b2k−1 ,b2k+1 ] (m(y)−m(b2k−1 )) , k ∈ [|S0 |] and b2|S0 |+2 = argmaxy∈(b
e0 ] (m(y)−
,B
2|S0 |+1

e0 ))2 (note that if the max doesn’t exist we can take a limit of partitions converging to it).
m(B
Next, we control I2 :
I2 =

sup

e0 ))2 ≤ (m(|a
e0 ))2 .
e hm −1 |) − m(
e B
(m(y) − m(−B

e0 ]
y∈(ahm −1 ,−B

with the last inequality following from (6). Combining this bound with (45) we get:
e0 ))2 .
e 2 |) − m(|a
e hm −1 |))2 + I2 ≤ (m(|a
e 2 |) − m(
e B
(m(|a
22

(48)

Similarly, for I4 and the other bound in (45) we have:
e0 ))2 .
e H |) − m(|a
e hM +1 |))2 + I4 ≤ (m(|a
e H |) − m(
e B
(m(|a

(49)

Finally, we deal with the tail part:
Var[m(Y )|Y ≤ a2 ] ≤ E[(m(Y ) − m(a2 ))2 |Y ≤ a2 ]

(50)

e
e 2 |))2 |Y ≤ a2 ]
≤ E[(m(|Y
|) − m(|a
e
e
≤ 4E[(m(|Y
|))2 |Y ≤ a2 ] ≤ 4(E[|m(|Y
|)|2+ξ |Y ≤ a2 ])2/(2+ξ)
2/(2+ξ)
Z a2
2+ξ
−1
e
=4
|m(|y|)|
dP(Y ≤ y)P(Y ≤ a2 )
−∞

= o(1)P(Y ≤ a2 )−2/(2+ξ) .
e
where we used the fact that E[|m(|Y
|)|2+ξ ] is bounded by assumption, and the o(1) is in the sense
of |a2 | → ∞ . We can show a similar inequality for the other tail — Var[m(Y )|Y ≥ aH ].
Combining (45), (48), (49), (46) and (50) we have:

2
2|S0 |+3
H
X
X

Var[m(Y )|ah < Y ≤ ah+1 ] ≤
sup
|m(bi ) − m(bi−1 )|
e0 )
b∈Π2|S0 |+3 (B

h=1

i=2

+ o(1)P(Y ≥ aH )−2/(2+ξ) + o(1)P(Y ≤ a2 )−2/(2+ξ)
e0 ))2 + (m(|a
e0 ))2 .
e 2 |) − m(
e B
e H |) − m(
e B
+ (m(|a
e0 ))2 ≤ 4(m(|a
e 2 |)2+ξ P(Y ≤ a2 ) → 0,
e 2 |)− m(
e B
e 2 |))2 , and we know that m(|a
e 2 |)2+ξ Hl ≤ m(|a
Since (m(|a
1
1
e 2 |)2 H 2/(2+ξ)
this means that m(|a
→ 0. Furthremore o(1)P(Y ≤ a2 )−2/(2+ξ) H 2/(2+ξ)
= o(1). Finally
we recall that by (5) we have:

2
2|S0 |+3


X

sup
|m(bi ) − m(bi−1 )| = o |S0 |2/(2+ξ) .
e0 )
b∈Π2|S0 |+3 (B

i=2

e0 ≤ Y ≤ B
e0 )H/l + 1 and thus:
However |S0 | ≤ P(−B

2
2|S0 |+3


X

sup
|m(bi ) − m(bi−1 )| = o H 2/(2+ξ) ,
e0 )
b∈Π2|S0 |+3 (B

i=2

which finishes the proof.
Proof of Lemma 2. Before we go to the main proof of the lemma we first formulate a simple but
useful concentration inequality.
e be a sub-Gaussian random variable with kXk
e ψ ≤ K. Let A(X,
e ν) ∈ {0, 1} be
Lemma 7. Let X
2
any (randomized) acceptance rule such that P(A = 1) ≥ q, with ν being any random variable. Let
e
e ν) = 1. Denote with µ = E[Xi ]. Then
X1 , . . . , Xr be an i.i.d. samples of the distribution X|A(
X,
there exist some absolute constants C1 , C2 > 0 such that:
!
2 r
p
.
P(|X − µ| > ) ≤ 2 exp −
C1 K2 exp( 1 − log(q)) + C2 K
23

p
As a Corollary to Lemma 7, observe that since supq∈[0,1] exp( 1 − log(q))q = e, we then have
the following:


2 r
P(|X − µ| > ) ≤ 2 exp −
,
C1 K2 q −1 + C2 K
for some absolute constants C1 , C2 > 0.
By (20) we know that P(Y ∈ Sh ) ≥ H1 − 2 on E, thus setting q =
conditionally on {Y(mh) : h ∈ [H − 1]} we have for all j ∈ G and all h:

P

j
X h,1:(m−1)

−

µjh


> η ≤ 2 exp −


η 2 (m − 1)
C1 K2 q −1 + C2 Kη

1
H

− 2, by Lemma 7


.

Note that Lemma 7 is applicable in this case, since the statistics X jh,πi are i.i.d. conditionally on
Y(m(h−1)) and Y(mh) , where π is a random permutation, as we noticed in the main text. Therefore
X jh,πi = X j |A(X j , ν) for the acceptance rule A(X j , ν) := A(X, ε) = 1(f (X | β, ε) ∈ Sh ). Furthermore, notice that the above inequality holds regardless of the values of {Y(mh) : h ∈ [H − 1]}, on
the event E.
Finally, using union bound across the slices and the indexes j ∈ G, we have that this holds for
all slices or in other words




η 2 (m − 1)
j
j
P
max X h,1:(m−1) − µh > η ≤ 2|G|H exp −
,
C1 K2 q −1 + C2 Kη
j∈G,h∈[H]
d

on the event E. This is precisely what we wanted to show.
d e
e ν) satisfies the following inProof of Lemma 7. Observe that the random variable X = X|A(
X,
equality:

e
e
e ≥ t|A(X,
e ν) = 1) = P(|X| ≥ t, A(X, ν) = 1)
P(|X| ≥ t) = P(|X|
e ν) = 1)
P(A(X,
e ≥ t) ≤ q −1 e exp(−ct2 /K2 ),
≤ q −1 P(|X|
e is assumed to be
where c is an absolute constant, and the last inequality follows by the fact that X
sub-Gaussian.
Clearly the above bound
can be substituted with the trivial bound 1, for values of
q
q

q
2
2
t ≤ −K /c log e . Let f (q) := −K /c log qe . Next, for any positive integer j ∈ N, we have:
j

Z

∞

P(|X| ≥ t)jt

E[|X| ] =

j−1

0

Z
dt ≤

f (q)

jt

j−1

Z

∞

dt +

0

q −1 e exp(−ct2 /K2 )jtj−1 dt.

f (q)

Multiplying by λj /j! and summing over j = 1, 2, . . ., we obtain:
Z
E[exp(λ|X|)] ≤ 1 + λ

f (q)

Z

∞

exp(λt)dt + λ
0

f (q)

q −1 e exp(−ct2 /K2 ) exp(λt)dt




√
λe πK
f (q) − K2 λ/(2c)
2 2
√
√ exp(K λ /4c) 1 − Φ
= exp(λf (q)) +
.
q c
K/ 2c

24

Assuming that λ is small enough so that f (q) − K2 λ/(2c) > 0, we can use the well known tail
bound 1 − Φ(x) ≤ φ(x)/x, for x > 0 to get:


√
λeK2 π exp(K2 λ2 /4c)
f (q) − K2 λ/(2c)
√ φ
√
E[exp(λ|X|)] ≤ exp(λf (q)) +
q(f (q) − K2 λ/(2c)) 2c
K/ 2c


λK2
= exp(λf (q)) 1 +
.
2c(f (q) − K2 λ/(2c))
Next, by the triangle inequality and the convexity of the exponent and the absolute value we have:
E[exp(λ|X − µ|)] ≤ E[exp(λ|X|)] exp[E[λ|X|]] ≤ E[exp(2λ|X|)].
Set Z := X − µ. We have showed the following:

E[exp[λ|Z|] − 1 − λ|Z|] ≤ E[exp[λ|Z|]] ≤ exp(2λf (q)) 1 +

λK2
c(f (q) − K2 λ/c)


,

(51)

p
for a λ such that f (q) − K2 λ/c > 0. By selecting λ := 21 Kc2 , one can easily verify that (f (q) −
λK2
K2 λ/c) ≥ λ−1 /4, and hence c(f (q)−K
2 λ/c) ≤ 1, for any q ≤ 1. With this choice of λ (51) becomes:
−2

E[exp[λ|Z|] − 1 − λ|Z|]λ

8K2
exp
≤
c

r


r

8K2
c
q
f (q) =
exp
.
− log
K2
c
e

Recall that a version Bernstein’s inequality (see Lemma 2.2.11 Van Der Vaart & Wellner (1996))
states that if the following moment condition E|Z|j ≤ j!λ−(j−2) v/2 is met for all j ≥ 2 and
Zi ∼ Z, i ∈ [r] are mean 0, i.i.d. then:
!


r
X
2
1
P
.
Zi ≥  ≤ 2 exp −
2 rv + λ−1 
i=1

Observe that by a Taylor expansion the condition E[exp[λ|Z|] − 1 − λ|Z|]λ−2 ≤ 12 v implies the
moment condition from above. Hence we conclude:
!
!
r
1
2 r
1X
p
Xi − µ ≥  ≤ 2 exp −
,
P
√
r
2 16K2 exp( 1 − log(q))/c + 2K/ c
i=1

which completes the proof.
Proof of Lemma 3. Using the sliced stability condition, for large H we get:
Var[mj (Y )] −

H
X

(µjh )2 P(Y

∈ Sh ) =

h=1

H
X

Var[mj (Y )|Y ∈ Sh ]P(Y ∈ Sh )

h=1

CH κ
≤
s




1
+ 2 .
H

This shows (22). Consequently we have:


X


H
H
X
1
CV
CH κ 1
− 2
(µjh )2 ≤
(µjh )2 P(Y ∈ Sh ) ≤
+
+ 2 .
H
s
s
H
h=1

h=1

25

This yields (23). To get (24) we proceed as follows:



v
v
uH
uH
X
H
H
uX
u
X
X
1
j
j
− 2
|µh | ≤
P(Y ∈ Sh )t (µjh )2 P(Y ∈ Sh )
|µh |P(Y ∈ Sh ) ≤ t
H
h=1
h=1
h=1
h=1
s


CV
CH κ 1
≤
+
+ 2 ,
s
s
H

and we are done.
e we have the following chain of inequalities:
Proof of Lemma 4. Note that on the event E
2
H 
1 X
1 j
m−1 j
(m − 1)2 j 2
X h,m +
X h,1:(m−1) −
(µh )
H
m
m
m2

(52)

h=1

H
H
H
1 X j 2 2(m − 1) X j
1 (m − 1)2 X
j
≤
(X h,m ) +
|X h,m |(η + |µh |) +
η(2|µjh | + η)
Hm2
Hm2
H
m2
h=1
h=1
h=1
v
v
u
u
H
n
n
u X
X
η
1 1X j 2
2 u
t (X j )2 t2
(η 2 + (µjh )2 ) + η 2 + 2 B3 .
≤
(X r ) +
r
mn
Hm
H
r=1

r=1

h=1

e in the first inequality, and (24), Cauchy-Schwartz and
where we used that we are on P
the event E
P
j
H
m−1
1
the trivial bounds m < 1, n h=1 (X h,m )2 ≤ n1 nr=1 (X jr )2 in the second one.
To this end, observe that (X jr )2 , r ∈ [n] are i.i.d. random variables with sub-exponential
distributions. This can be seen from the standard inequality:
k(X jr )2 kψ1 ≤ 2kX jr k2ψ2 ≤ 2K2 .
Clearly we also have E[(X jr )2 ] ≤ 2K2 . Denote the mean E[(X jr )2 ] = ν. Now we are in a position to
use a Bernstein type of deviation inequality (see Proposition 5.16 in Vershynin (2010)). We obtain:
!

 2

n
1X j 2
τ n τn
(X r ) > ν + τ ≤ exp −c min
,
,
P
n
4K4 2K2
r=1

e
e ⊂E
e
for some absolute constant c > 0. Hence we infer that, when
τ ≤ 2K2 , there exists a set E
j 2
τ 2n
1 Pn
2
failing with probability at most p exp(−c 4K4 ), such that n r=1 (X r ) ≤ ν + τ ≤ 2K + τ for all
e
e , we get:
j ∈ [p]. Therefore continuing the bound on the event E
r
√
(2K2 + τ ) 2 2K2 + τ
B2
B3
√
(52) ≤
+
2η 2 + 2
+ η 2 + 2η ,
m
H
H
m
where we used (23). This finishes the proof.

26

Proof of Lemma 6. Fix any τ ∈ (0, 2]. Define the event:
e
eS ∩ { sup n−1
eS = E
E
β
β
j∈Sβ

n
X

(X jr )2 ≤ 2 + τ }

i=1

We first note that the following inequality holds:
V

jk

H
CV
1 X j k
CV
− sign(β j ) sign(β k )
=
X h X h − sign(β j ) sign(β k )
s
H
s
h=1

1
H

≤

H
X

j

(X h )2 −

h=1

CV
1
+
s
H

H
X

j

j

k

X h sign(β j )X h − sign(β k )X h ,

h=1

where the LHS equals, the LHS of (36) after using (34). Fortunately (25) and (27) already give
e
e S , which can be seen with identical arguments to Lemmas
bounds on the first term on the event E
β
3 and 4. We now show that the second term is small on the same event. Note that the following
identity holds:
H
1 X
j
j
k
X h sign(β j )X h − sign(β k )X h
H
h=1

H
1 X m−1 j
1
m−1
j
=
X h,1:(m−1) + X jh,m sign(β j )X h −
sign(β k )µjh
H
m
m
m
h=1

+

m−1
k
sign(β k )µkh − sign(β k )X h .
m

eS :
Thus on the event E
β
H
1 X
j
j
k
X h sign(β j )X h − sign(β k )X h
H
h=1


H 
1
1 X
1
1
j
j
k
≤
µh + η +
X h,m
2η +
X h,m +
X h,m ,
H
m
m
m
h=1

where µh = |µjh | = |µkh |, and we used that

m−1
m

eS we have |X j
< 1, and the fact that on E
β
h,1:(m−1) −

µjh | ≤ η and similarly |X kh,1:(m−1) − µkh | ≤ η. Next we have:


H 
1 X
1
1
1
j
j
k
µh + η +
X h,m
2η +
X h,m +
X h,m
H
m
m
m
h=1

H
H
η X
1 X
≤2
µh +
(µh + η)(|X jh,m | + |X kh,m |) + 2η 2
H
Hm

+

2η
mH

h=1
H
X
h=1

h=1

|X jh,m | +

H
H
H
1 X j 2
1 X j 2
1 X k 2
(X
)
+
(X
)
+
(X h,m ) ,
h,m
h,m
m2 H
2m2 H
2m2 H
h=1

h=1

27

h=1

where we used the simple inequality ab ≤ (a2 + b2 )/2. Luckily we have already controlled all of the
above quantities. Using Lemma 4 and (25) we get:
H
1 X
j
j
k
X h sign(β j )X h − sign(β k )X h
H
h=1
r
√
√
2 2+τ
B2
η 2+τ
2+τ
η
2
2
2η + 2
+ 2η + 2 √
,
+2
≤ 2 B3 + √
H
H
m
m
m

e
e S we have 1 Pn (X j )2 ≤ 2 + τ , the rest
where we heavily relied on the fact that on E
r
β
r=1
mH
j
1 PH
of the bounds can be seen in the proof of Lemma 4. (For the term note mH
h=1 |X h,m | ≤
qP
qP
√
j
H
n
η √2+τ
j 2
1
2
√1
≤ η 2 + 2+τ
h=1 (X h,m ) /H ≤ m
r=1 (X r ) /mH). Finally noting that 2
m
m gives the
m
desired result.
Lemma 8. Let N be a s × s symmetric “noise” matrix satisfying:
CV
.
40s

kNkmax ≤

(53)

Then the following occurs:
e S ,S =
(a) Let γ1 be the maximum eigenvalue of V
β β

CV
2

e S ,S ).
β Sβ β |Sβ + N, i.e. γ1 := λmax (V
β β

Then we have |γ1 − C2V | ≤ C40V ; Furthermore if γ2 is the second largest in magnitude eigenvalue
e S ,S , we have |γ2 | ≤ CV , and hence γ1 > |γ2 |.
of V
β β
40
e S ,S — e
(b) The corresponding principal eigenvector of V
zSβ 4 satisfies the following inequality:
β β
e
z Sβ − β Sβ

∞

1
≤ √ .
2 s

Proof. This Lemma is a re-statement of Lemma 6 from Amini & Wainwright (2008), the difference
being that we require precise bounds, rather than considering simply the asymptotics. To see part
(a), observe that (53) implies:
kNk2,2 = sup |v| Nv| ≤ kvk21 kNkmax ≤
kvk2 =1

CV
.
40

(54)

Hence by Weyl’s inequality we have that:
|γ1 −

CV
CV
|≤
,
2
40

|γ2 − 0| ≤

CV
,
40

which is exactly part (a). For the second part observe that
X
CV
CV
kNk∞,∞ = max
|Nij | ≤ s
=
.
i
40s
40

(55)

j

The proof of part (b) can then be seen as in Lemma 6 of Amini & Wainwright (2008), by carefully
exploiting (54) and (55).
4

Here we mean the principal eigenvector oriented so that e
z|Sβ β Sβ ≥ 0.

28

Proof of Proposition 2. The proof is based on an application of Fano’s inequality, which in turn
is a standard approach for showing minimax lower bounds (e.g. see Cover & Thomas (2012),
Wainwright (2009), Yang & Barron (1999), Yu (1997) among others). In particular we base our
proof on ideas from Amini & Wainwright (2008), Wainwright (2009).
For two probability measures P and Q, which are absolutely Rcontinuous with respect to a third
dQ
probability measure µ define their KL divergence by D(P kQ) = p log pq dµ, where p = dP
dµ , q = dµ .
We proceed with the following lemma:
Lemma 9. Let us have n observations from a SIM Y = f (X | β, ε), X ∼ Np (0, Ip ). Assume that
for any fixed u, v ∈ R the following regularity condition holds for the random variables f (u, ε) and
f (v, ε):
D(p(f (u, ε))kp(f (v, ε))) ≤ exp(Ξ(u − v)2 ) − 1,

(56)

where Ξ is a positive constant. Then if
n<

s log(p − s + 1)
,
8Ξ

and s ≥ 8Ξ, any algorithm recovering the support in our model will have errors with probability at
least 12 asymptotically.
To finish the proof simply apply Example 3 with P (x) = x/(2σ 2 ) and G = h = Id.
Proof of Lemma 9. For simplicity of the exposition we will assume that the vector β has only nonnegative entries (i.e. all non-zero entries are √1s ). The proof extends in exactly the same way in
the case when the entries of β are not restricted to be positive.

Let S ⊂ 2[p] , the set of all subsets of [p] with s elements. Clearly, |S| = ps . Let Sb : (Rp+1 )n → S
be any potentially random function, which is used to recover the support of β, based on the sample
{(Yi , X i )}ni=1 . Under the 0-1 loss the risk equals the probability of error:
1 X

PSβ (Sb 6= Sβ ),
p
s

(57)

Sβ ∈S

where by PSβ we are measuring the probability under a dataset generated with supp(β) equal to
the index of the measure PSβ .
Instead of directly dealing with the sum above, we first consider the p − s + 1 element set
e
S = {S ∈ S : [s − 1] ⊂ S, |S| = s}, and we bound the probability of error, on any function Sb (even
if given the knowledge that the true support is drawn from e
S). Let J be a uniformly distributed in
e
S. By Fano’s inequality that:
P(error) ≥ 1 −

I(J; (Y, X)n ) + log(2)
,
log |e
S|

(58)

where I(J; (Y, X)n ) is the mutual information between the sample J and the sample (Y, X)n . Note
now that for the mutual information we have
I(J; (Y, X)n ) = I(J; (f (X | β, ε), X)n )
≤ nH((f (X | β, ε), X)) − nH((f (X | β, ε), X)|J), 5
29

where the last inequality follows from the chain inequality of entropy. We therefore need an upper
bound on the last expression.
1(j=i)
√
Let i ≥ s and set β i := (β1i , β2i , . . . , βpi )| the vector such that that βji = 1(j∈[s−1])+
. When
s
J is unknown the distribution is a mixture, and hence due to the convexity of − log we have the
following inequality:
H((f (X | β, ε), X)) − H((f (X | β, ε), X))|J)
X
p(f (X | β i , ε), X)
1
| i
p((f
(X
β
,
ε),
X))
log
≤
(p − s + 1)2
p(f (X | β j , ε), X)
i,j≥s
=

p−s
D(p((f (X | β s , ε), X))kp((f (X | β s+1 , ε), X))),
p−s+1

where the last equality follows by a symmetric argument. Since the KL divergence is invariant under
changing variables, setting U = X | β s , V = X | β s+1 and W = P{βs ,βs+1 }⊥ X, where P{βs ,βs+1 }⊥
denotes the orthogonal projection on the space span{β s , β s+1 }⊥ . We get:
H((f (X | β, ε), X)) − H((f (X | β, ε), X))|J)
p−s
D(p((f (U, ε), U, V ))kp((f (V, ε), U, V ))),
≤
p−s+1
where we used the fact that W is independent of U, V, ε. Applying assumption (56) we get:
D(p((f (U, ε), U, V ))kp((f (V, ε), U, V ))) ≤ E exp(Ξ(U − V )2 ) − 1
r
s
4Ξ
−1≤
.
=
s − 4Ξ
s
where the first inequality can be obtained by conditioning on U, V , in the equality we used the fact
that U − V ∼ N (0, 2s ), and we assume that the value of s is large enough so that s ≥ 8Ξ. The
inequality in the preceding display helps us to conclude that for large values of s:
I(J; (Y, X)n ) ≤

4Ξn
4Ξ(p − s)n
<
.
s(p − s + 1)
s

1
Consequently by (58) if n < s log(p−s+1)
we will have errors with probability at least
8Ξ
2 , asymp
p
terms, by the
totically. To finish the conclusion, note that the sum (57), can be split into s−1
following operation:

1. Repeat each set in S — s times, and denote this superset by s × S

p
2. For each S of the s−1
, subsets of [p] with s − 1 elements, collect p − s + 1 distinct elements
of s × S containing S
3. Apply the

1
2

error bound obtained from above to this local sum.
p

In the end we get that the probability of error by selecting S ⊂ S uniformly is at least:

s

s + 1) 12 = 12 .
5

1 (s−1)
s (p) (p−

Here we use H to denote the entropy, not to be confused with the number of slices H.

30

B

SIM Examples

In this section we look into models of the type Y = G(h(X | β) + ) and show that under certain
sufficient conditions they belong to a class FA for some A. In addition, we provide a sub-class of
these models, in which support recovery is impossible with probability at least 21 when Γ is small.
Example 1. For SIM Y = G(h(X | β) + ε) with strictly monotone h and G there exists A > 0
(depending on h, G) such that Var(E[X | β|Y ]) ≥ A > 0.
Proof. We will show more generally that if Y = f (X | β, ε), where f, ε satisfy the condition that
the function g(z) = E[f (Z, ε)|Z = z] is strictly monotone, we have: Var(E[Z|f (Z, ε)]) > 0.
Note that the condition ∃A > 0 : Var(E[Z|f (Z, ε)]) ≥ A is equivalent E[Z|f (Z, ε)] to not
being a constant. We argue that the latter is clearly implied if, for example, E[Zf (Z, ε)] 6= 0.
To see this assume the contrary, i.e. E[Z|f (Z, ε)] = 0 a.s., but E[Zf (Z, ε)] 6= 0. Then we have
E[Zf (Z, ε)] = E[E[Z|f (Z, ε)]f (Z, ε)] = 0, which is a contradiction.
Next we show that our condition implies E[Zf (Z, ε)] 6= 0. WLOG assume that g is strictly
increasing. Observe that since Z ∼ N (0, 1) is a continuous random variable, by Chebyshev’s
association inequality (Boucheron et al. 2013) we have:
E[Zf (Z, ε)] = E[Zg(Z)] > E[Z]E[g(Z)] = 0,
and hence as we argued earlier it follows that Var(E[Z|f (Z, ε)]) > 0. Due to the independence
of X and ε models with a coordinate-wise monotone f function (strictly monotone in the first
coordinate) belong to this class.
Example 2. Let Y = G(h(X | β) + ε), where G, h are strictly monotone continuous functions.
Furthermore let ε be such that the function m(y) is continuous and monotone in y. We argue that
such models are sliced stable in the sense of Definition 1.
Proof. Since m is a continuous and monotone function in y, it follows that it is of bounded variation
on any
p closed interval, which in turn implies (5). Furthermore, by Example 1, we are guaranteed
that Var(mj (Y ))s = CV for some CV > 0. Notice also that in the case of a continuous and
e
monotone m one can take m(y)
:= |m(y) − m(−y)| < |m(y)| + |m(−y)|, for y > 0. Finally we need
e
to show that E[|m(|Y
|)|(2+ξ) ] is finite. Due to the last inequality it suffices to control:
e )|2+ξ ≤ s(2+ξ)/2 E|E[X j |Y ]|2+ξ = s(2+ξ)/2 E|E[s−1/2 X | β|Y ]|2+ξ
CV2+ξ E|m(Y
≤ s(2+ξ)/2 s−(2+ξ)/2 EE[|X | β|2+ξ |Y ] = E|Z|2+ξ < ∞,
where Z ∼ N (0, 1) and in the first equality we used the fact that X j − s−1/2 X | β ⊥⊥ X | β and
hence E[X j − s−1/2 X | β|Y ] = E[X j − s−1/2 X | β] = 0. This completes the proof.
Remark 3. To see that the monotonicity condition on m(y) is not vacuous, we give concrete
examples of SIM satisfying it, which include the simple linear regression model (8) as a special
case. Consider the models from Example 2 and let ε have a density function satisfying pε (x) ∝
exp(−P (x2 )), where P is any nonzero polynomial with non-negative coefficients such that P (0) = 0.
To see that m is monotonic and continuous, we start by obtaining a precise expression of m. It is
simple to see that:
E[Z|h(Z) + ε = G−1 (y)]
p
m(y) =
,
Var(mj (Y ))s
31

p
for almost every y ∈ supp(Y ). By Example 1, we have that Var(mj (Y ))s = CV for some
CV > 0. Next we argue that y 7→ m(y) is monotone. This is equivalent to showing that m(G(y))
is monotone. To this end we apply Lemma A.2 of Duan & Li (1991), which implies that it suffices
for the random variable h(Z) + ε|Z to have a monotone likelihood ratio in order for m(G(y)) to
be monotone. Since the family of random variables h(z) + ε is a location family, the normalizing
constants of their densities do not change with z. This in conjuction with the fact that h is
increasing, implies that the monotonicity of the likelihood ratio will be implied if we show that the
function x 7→ P (x2 ) − P ((x − c)2 ) is increasing in x for any fixed c > 0. Notice that since P (x2 ) is a
(y 2 )
differentiable convex function by construction, we have that P (x2 )−P ((x−c)2 ) ≥ c dPdy
|y=x−c > 0.
It is worth noting that the same argument applies more generally to the case where ε is a logconcave random variable (i.e. pε (x) = exp(−ϕ(x)) where ϕ is a convex function). The fact that m
is continuous follows by the continuity of G and h.
Finally, with the help of Lemma 9 we demonstrate that some models discussed in Remark 3
meet the information theoretic barrier described in Proposition 2, and hence their support cannot
be recovered by any algorithm unless Γ is large enough.
Example 3. Let Y = G(h(X | β) + ε), where G, h are strictly monotone continuous functions and
in addition h is an L-Lipschitz function. Furthermore let ε has a density as specified in Remark
3, i.e. pε (x) ∝ exp(−P (x2 )), where P is any nonzero polynomial with non-negative coefficients
such that P (0) = 0. Then if n < s log(p−s+1)
for some constant C > 0 (depending on P, G, h)
C
and s sufficiently large, any algorithm recovering the support in our model will have errors with
probability at least 12 asymptotically.
Proof. Note that all moments of the random variable ε exist. Next we verify that condition (56)
holds in this setup. Since G is 1-1 and KL divergence is invariant under changes of variables WLOG
we can assume our model is simply Y = h(X | β) + ε or in other words f (u, ε) = h(u) + ε. This is
a location family for u ∈ R and thus the normalizing constant of the densities will stay the same
regardless of the value of u. Direct calculation yields:
D(p(f (u, ε))kp(f (v, ε))) = E[P ((ξ + h(u) − h(v))2 ) − P (ξ 2 )]
= Pe((h(u) − h(v))2 ),
where ξ has a density pξ (x) ∝ exp(−P (x2 )), and Pe is another nonzero polynomial with nonnegative
coefficients, with Pe(0) = 0 of the same degree as P . The last equality follows from the fact that
all odd moments of ξ are 0, since ξ is a symmetric about 0 distribution. Since h is L-Lipschitz we
conclude that:
D(p(f (u, ε))kp(f (v, ε))) ≤ Pe(L2 (u − v)2 ).
The last can be clearly dominated by exp(C(u − v)2 ) − 1 for a large enough constant C.

References
Amini, A. A. & Wainwright, M. J. (2008), High-dimensional analysis of semidefinite relaxations
for sparse principal components, in ‘Information Theory, 2008. ISIT 2008. IEEE International
Symposium on’, IEEE, pp. 2454–2458.
32

Berthet, Q. & Rigollet, P. (2013), Complexity theoretic lower bounds for sparse principal component
detection, in ‘Conference on Learning Theory’, pp. 1046–1066.
Boucheron, S., Lugosi, G. & Massart, P. (2013), Concentration inequalities: A nonasymptotic theory
of independence, OUP Oxford.
Brenner, N., Bialek, W. & Van Steveninck, R. d. R. (2000), ‘Adaptive rescaling maximizes information transmission’, Neuron 26(3), 695–702.
Cai, T., Liu, W. & Luo, X. (2011), ‘A constrained l1 minimization approach to sparse precision
matrix estimation 1 minimization approach to sparse precision matrix estimation’, Journal of
the American Statistical Association 106(494), 594–607.
Cai, T. T., Ma, Z., Wu, Y. et al. (2013), ‘Sparse pca: Optimal rates and adaptive estimation’, The
Annals of Statistics 41(6), 3074–3110.
Candes, E. & Tao, T. (2007), ‘The dantzig selector: Statistical estimation when p is much larger
than n’, The Annals of Statistics pp. 2313–2351.
Chen, X., Xu, H., Yuan, P., Fang, F., Huss, M., Vega, V. B., Wong, E., Orlov, Y. L., Zhang, W.,
Jiang, J. et al. (2008), ‘Integration of external signaling pathways with the core transcriptional
network in embryonic stem cells’, Cell 133(6), 1106–1117.
Cloonan, N., Forrest, A. R., Kolle, G., Gardiner, B. B., Faulkner, G. J., Brown, M. K., Taylor,
D. F., Steptoe, A. L., Wani, S., Bethel, G. et al. (2008), ‘Stem cell transcriptome profiling via
massive-scale mrna sequencing’, Nature methods 5(7), 613–619.
Cook, R. D. & Ni, L. (2005), ‘Sufficient dimension reduction via inverse regression’, Journal of the
American Statistical Association 100(470).
Cook, R. D. et al. (2004), ‘Testing predictor contributions in sufficient dimension reduction’, The
Annals of Statistics 32(3), 1062–1092.
Cover, T. M. & Thomas, J. A. (2012), Elements of information theory, John Wiley &amp; Sons.
d’Aspremont, A., Bach, F. & Ghaoui, L. E. (2008), ‘Optimal solutions for sparse principal component analysis’, The Journal of Machine Learning Research 9, 1269–1294.
d’Aspremont, A., El Ghaoui, L., Jordan, M. I. & Lanckriet, G. R. (2007), ‘A direct formulation for
sparse pca using semidefinite programming’, SIAM review 49(3), 434–448.
Duan, N. & Li, K.-C. (1991), ‘Slicing regression: a link-free regression method’, The Annals of
Statistics pp. 505–530.
Hsing, T. & Carroll, R. J. (1992), ‘An asymptotic theory for sliced inverse regression’, The Annals
of Statistics pp. 1040–1061.
Jiang, B., Liu, J. S. et al. (2014), ‘Variable selection for general index models via sliced inverse
regression’, The Annals of Statistics 42(5), 1751–1786.
Johnstone, I. M. & Lu, A. Y. (2004), ‘Sparse principal components analysis’, Unpublished
manuscript .
33

Johnstone, I. M. & Lu, A. Y. (2009), ‘On consistency and sparsity for principal components analysis
in high dimensions’, Journal of the American Statistical Association 104(486).
Krauthgamer, R., Nadler, B. & Vilenchik, D. (2013), ‘Do semidefinite relaxations really solve sparse
pca?’, arXiv preprint arXiv:1306.3690 .
Laurent, B. & Massart, P. (2000), ‘Adaptive estimation of a quadratic functional by model selection’, Annals of Statistics pp. 1302–1338.
Li, K.-C. (1991), ‘Sliced inverse regression for dimension reduction’, Journal of the American Statistical Association 86(414), 316–327.
Li, L. & Nachtsheim, C. J. (2006), ‘Sparse sliced inverse regression’, Technometrics 48(4).
Lin, Q., Zhao, Z. & Liu, S. J. (2015), ‘On consistency and sparsity of sliced inverse regression in
high dimensions’, arXiv preprint arXiv:1507.03895 .
Massart, P. (1990), ‘The tight constant in the dvoretzky-kiefer-wolfowitz inequality’, The Annals
of Probability pp. 1269–1283.
Ouyang, Z., Zhou, Q. & Wong, W. H. (2009), ‘Chip-seq of transcription factors predicts absolute
and differential gene expression in embryonic stem cells’, Proceedings of the National Academy
of Sciences 106(51), 21521–21526.
Paul, D. (2007), ‘Asymptotics of sample eigenstructure for a large dimensional spiked covariance
model’, Statistica Sinica 17(4), 1617.
Price, A. L., Zaitlen, N. A., Reich, D. & Patterson, N. (2010), ‘New approaches to population
stratification in genome-wide association studies’, Nature Reviews Genetics 11(7), 459–463.
Tibshirani, R. et al. (1997), ‘The lasso method for variable selection in the cox model’, Statistics
in medicine 16(4), 385–395.
Van Der Vaart, A. W. & Wellner, J. A. (1996), Weak Convergence, Springer.
Vershynin, R. (2010), ‘Introduction to the non-asymptotic analysis of random matrices’, arXiv
preprint arXiv:1011.3027 .
Vu, V. Q. & Lei, J. (2012), ‘Minimax rates of estimation for sparse pca in high dimensions’, arXiv
preprint arXiv:1202.0786 .
Wainwright, M. J. (2009), ‘Information-theoretic limits on sparsity recovery in the high-dimensional
and noisy setting’, Information Theory, IEEE Transactions on 55(12), 5728–5741.
Yang, Y. & Barron, A. (1999), ‘Information-theoretic determination of minimax rates of convergence’, Annals of Statistics pp. 1564–1599.
Yu, B. (1997), Assouad, fano, and le cam, in ‘Festschrift for Lucien Le Cam’, Springer, pp. 423–435.
Yu, Z., Zhu, L., Peng, H. & Zhu, L. (2013), ‘Dimension reduction and predictor selection in
semiparametric models’, Biometrika p. ast005.

34

Zhang, Y. & Ghaoui, L. E. (2011), Large-scale sparse principal component analysis with application
to text data, in ‘Advances in Neural Information Processing Systems’, pp. 532–539.
Zhong, W., Zhang, T., Zhu, Y. & Liu, J. S. (2012), ‘Correlation pursuit: forward stepwise variable selection for index models’, Journal of the Royal Statistical Society: Series B (Statistical
Methodology) 74(5), 849–870.
Zhu, L., Miao, B. & Peng, H. (2006), ‘On sliced inverse regression with high-dimensional covariates’,
Journal of the American Statistical Association 101(474).

Extra Numerical Studies Figures

●

●

●

●

●

●

●

n
s log(p−s)

1.0

1.0

Figure 4: Efficiency Curves for DT-SIR, s = log p, Γ =

●

●

0.0

●

10

●

●

●

●

15
n

20

25

0.8
0.6

●

●

30

● ●

●

●

0

5

10

15
n

20

25

30

slog(p − s)

●

●

(b) Model (14)

●

●

●

●

1.0

(a) Model (13)

●

p = 100
p = 200
p = 300
p = 600
p = 1200

●

slog(p − s)

1.0

●

●

0.0

p = 100
p = 200
p = 300
p = 600
p = 1200

●
●

5

●

0.4

^
P(supp(sign(β)) = supp(sign(β)))

0.6
0.4

●

0

●

0.2

0.8

●

●

● ●

∈ [0, 30]

●

0.2

^
P(supp(sign(β)) = supp(sign(β)))

●

●

●

●

●

●

●

●

●

●

●

●

●
● ●

0

●

5

10

15
n

20

25

p = 100
p = 200
p = 300
p = 600
p = 1200

0.8
0.6

●
●

●

0.0

●

●

0.4

0.4

●

●

0.2

^
P(supp(sign(β)) = supp(sign(β)))

0.6

0.8

●

0.2

^
P(supp(sign(β)) = supp(sign(β)))

●

0.0

C

30

● ●

0

slog(p − s)

●

5

10

15
n

20

slog(p − s)

(c) Model (15)

(d) Model (16)

35

25

p = 100
p = 200
p = 300
p = 600
p = 1200
30

n
s log(p−s)

●
●

●

●

●

∈ [0, 40]

1.0

1.0

Figure 5: Efficiency Curves for SDP, s = log p, Γ =

●

●

●

●

●

●

●

●

●

●

0.4

●

0.8
0.6

●
●

●

●

●

0.4

0.6

●

●

●

0.2

0.8

^
P(supp(sign(β)) = supp(sign(β)))

●

●

0.2

^
P(supp(sign(β)) = supp(sign(β)))

●

●

p = 100
p = 200
p = 300

0.0

●●

●

●

●

0

●

10

20
n

30

●

0.0

●
●

40

●●

●

●

●

0

10

20
n

slog(p − s)

30

40

slog(p − s)

●

●

●

●

●

●

●

1.0

(b) Model (14)

1.0

(a) Model (13)

●

●

●

●

●

●

●

●

●

●

●
●

0.4

●

0.2

●

0.8

●

0.6

0.6

●

●

●

●

0.4

●

0.2

^
P(supp(sign(β)) = supp(sign(β)))

0.8

●

●

●●

0

●

●

●

10

20
n

30

●

p = 100
p = 200
p = 300

●
●

0.0

●
●

0.0

^
P(supp(sign(β)) = supp(sign(β)))

p = 100
p = 200
p = 300

●

●

40

●●

0

slog(p − s)

●

●

●

10

20
n
slog(p − s)

(c) Model (15)

(d) Model (16)

36

30

p = 100
p = 200
p = 300
40

●
●

●

√

p, Γ =

n
s log(p−s)

e βk
e 2 according to
∈ [0, 30], β = β/k

1.0

1.0

Figure 6: Efficiency Curves for DT-SIR, s =
(17)

●

●
●

●

●
●

●

p = 100
p = 200
p = 300
p = 600
p = 1200

0.0

●
● ●

●

0

5

10

15
n

20

25

0.8
0.6

●

●

p = 100
p = 200
p = 300
p = 600
p = 1200

●

●
●

0.0

●

●

●

●

0.4

^
P(supp(sign(β)) = supp(sign(β)))

0.4

●

●

●

0.2

0.8
0.6

●

0.2

^
P(supp(sign(β)) = supp(sign(β)))

●
●
●

30

● ●

●

0

●

5

10

15
n

slog(p − s)

20

25

30

slog(p − s)

(b) Model (14)

●

●

1.0

1.0

(a) Model (13)

●

●

●

●
●

●

●

●

●
●
●
● ●

0

●

●

5

10

15
n

20

25

p = 100
p = 200
p = 300
p = 600
p = 1200

0.8
0.6

●

●

●

0.0

0.2

0.4

●

●
●
●

0.4

●

●

0.2

^
P(supp(sign(β)) = supp(sign(β)))

0.8
0.6

●

0.0

^
P(supp(sign(β)) = supp(sign(β)))

●

●

30

● ●

0

slog(p − s)

●

●

5

●

10

15
n

20

slog(p − s)

(c) Model (15)

(d) Model (16)

37

25

p = 100
p = 200
p = 300
p = 600
p = 1200
30

●
●

●

●

√

p, Γ =

n
s log(p−s)

●

●

●

●

●

●

●

●
●

●

●

●

●

●

p = 100
p = 200
p = 300

●

●

●

0

●

●

10

20
n

30

0.8

●

●

●

p = 100
p = 200
p = 300

●
●

0.0

0.0

●

●●

0.6

●

0.4

●

●

●

0.4

0.6

●

●

0.2

^
P(supp(sign(β)) = supp(sign(β)))

0.8

●

0.2

^
P(supp(sign(β)) = supp(sign(β)))

e βk
e 2 according to (17)
∈ [0, 40], β = β/k

1.0

1.0

Figure 7: Efficiency Curves for SDP, s =

40

●●

●

●

●

0

●

10

20
n

slog(p − s)

30

40

slog(p − s)

●

(b) Model (14)

●

●

●

●

●

●

1.0

1.0

(a) Model (13)

●

●

●

●

●

●

●

●

●

0.4

●

0.2

●

0

●

●

●

10

20
n

30

p = 100
p = 200
p = 300

0.8
0.6

●
●

0.0

●

●

●

●

●

●●

●
●

0.4

0.6

●

●

●

0.2

^
P(supp(sign(β)) = supp(sign(β)))

0.8

●

●

0.0

^
P(supp(sign(β)) = supp(sign(β)))

●
●

40

●●

0

slog(p − s)

●

●

●

10

20
n
slog(p − s)

(c) Model (15)

(d) Model (16)

38

30

p = 100
p = 200
p = 300
40

