Partial Consistency with Sparse Incidental Parameters
Jianqing Fan1 , Runlong Tang2 , and Xiaofeng Shi1

arXiv:1210.6950v3 [] 2 Aug 2017

1
2

Princeton University

Johns Hopkins University
March 14, 2018

Abstract
Penalized estimation principle is fundamental to high-dimensional problems. In the literature, it has been extensively and successfully applied to various models with only structural
parameters. As a contrast, in this paper, we apply this penalization principle to a linear regression model with a finite-dimensional vector of structural parameters and a high-dimensional
vector of sparse incidental parameters. For the estimators of the structural parameters, we derive their consistency and asymptotic normality, which reveals an oracle property. However, the
penalized estimators for the incidental parameters possess only partial selection consistency but
not consistency. This is an interesting partial consistency phenomenon: the structural parameters are consistently estimated while the incidental ones cannot. For the structural parameters,
also considered is an alternative two-step penalized estimator, which has fewer possible asymptotic distributions and thus is more suitable for statistical inferences. We further extend the
methods and results to the case where the dimension of the structural parameter vector diverges
with but slower than the sample size. A data-driven approach for selecting a penalty regularization parameter is provided. The finite-sample performance of the penalized estimators for the
structural parameters is evaluated by simulations and a real data set is analyzed.
Keywords: Structural Parameters, Sparse Incidental Parameters, Penalized Estimation, Partial Consistency,
Oracle Property, Two-Step Estimation, Confidence Intervals

1

Introduction

Since the pioneering papers by Tibshirani (1996) and Fan and Li (2001), the penalized estimation
methodology for exploiting sparsity has been studied extensively. For example, Zhao and Yu (2006)
1

provides an almost necessary and sufficient condition, namely Irrepresentable Condition, for the
LASSO estimator to be strong sign consistent. Fan and Lv (2011) shows that an oracle property
holds for the folded concave penalized estimator with ultrahigh dimensionality. For an overview on
this topic, see Fan and Lv (2010).
All the aforementioned papers consider models only with the so-called structural parameters,
which are related to every data point. By contrast, we consider in this paper another type of models
where there are not only the structural parameters but also the so-called incidental parameters,
each of which is related to only one data point. Specifically, suppose data {X i , Yi }ni=1 are from the
following linear model:
Yi = µ?i + X Ti β ? + i ,

(1.1)

where the vector of incidental parameters µ? = (µ?1 , · · · , µ?n )T is sparse, the vector of structural
parameters β ? = (β1? , · · · , βd? )T is of main interest, and for each i, X i is a d-dimensional covariate
vector, and i is a random error. Let ν = (µ?T , β ?T )T . Then, in model (1.1), a different data point
(X i , Yi ) depends on a different subset of ν, that is, µ?i and β ? .
Model (1.1) arises as a working model for estimation from Fan et al. (2012b), which considers a
large-scale hypothesis testing problem under arbitrary dependence of test statistics. By Principal
Factor Approximation, a method proposed by Fan et al. (2012b), the dependent test statistics
Z = (Z1 , · · · , Zp )T ∼ N (µ, Σ) can be decomposed as Zi = µi +bTi W +Ki , where µ = (µ1 , · · · , µp )T
and bi is the ith row of the first k unstandardized principal components, denoted by B, of Σ
and K = (K1 , · · · , Kp )T ∼ N (0, A) with A = Σ − BB T . The common factor W drives the
dependence among the test statistics. This realized but unobserved factor is critical for False
Discovery Proportion (FDP) estimation and power improvements by removing the common factor
{bTi W } from the test statistics. Hence, an important goal is to estimate W with given {bi }ni=1 .
In many applications on large-scale hypothesis testing, the parameters {µi }pi=1 are sparse. For
example, genome-wide association studies show that the expression level of gene CCT8 is highly
related to the phenotype of Down Syndrome. It is of interest to test the association between each
of millions of SNP’s and the CCT8 gene expression level. In the framework of Fan et al. (2012b),
each µi stands for such an association. That is, if µi = 0, the ith SNP is not associated with the
CCT8 gene expression level; otherwise, it is associated. Since most of the SNP’s are not associated
the CCT8 gene expression level, it is reasonable to assume {µi }pi=1 are sparse. Replacing Zi , µi , bi ,
W , k, p, and Ki with Yi , µ?i , X i , β ? , d, n, and i respectively, we obtain model (1.1) formally. It
is of interest to study this model independently with simplifications.
2

Although model (1.1) emerges from a critical component of estimating FDP in Fan et al.
(2012b), it stands with its own interest. For example, in some applications, there are only few
signals (nonzero µ?i ’s) and what is interesting is to learn about β ? , which reflects the relationship
between the covariates and response. For another example, those few nonzero µ?i ’s might be some
measurement or recording errors of the responses {Yi }. In these cases, model (1.1) is suitable for
modeling data with contaminated responses and a method producing a reliable estimator for β ? is
essentially a robust replacement for ordinary least squares estimate, which is sensitive to outliers.
Several models with structural and incidental parameters have first been studied in a seminal
paper by Neyman and Scott (1948), which points out the inconsistency of the maximum likelihood
estimators (MLE) of structural parameters in the presence of a large number of incidental parameters and provides a modified MLE. However, their method does not work for model (1.1) due to
no exploration of the sparsity of incidental parameters. Kiefer and Wolfowitz (1956) shows the
consistency of the MLE of the structural parameters when the incidental parameters are assumed
to be from a common distribution. That is, they eliminate the essential high-dimensional issue
of the incidental parameters by randomizing them. In contrast, this paper considers deterministic incidental parameters and handles the high-dimensional issue by penalization with a sparsity
assumption. Basu (1977) considers the elimination of nuisance parameters via marginalizing and
conditioning methods and Moreira (2009) solves the incidental parameter problem with an invariance principle. For a review of the incidental parameter problems in statistics and economics, see
Lancaster (2000).
Without loss of generality, suppose the first s incidental parameters {µ?i }si=1 are nonvanishing
and the remaining are zero. Then, model (1.1) can be written in a matrix form as Y = Xν + ,
where


X=

Is

X T1,s

0

0

X Ts+1,n

I n−s


,

X Ti,j = (X i , X i+1 , · · · , X j )T , I k is a k × k identity matrix, 0 is a generic block of zeros and
ν = (µ?1 , · · · , µ?s , β T , µ?s+1 , · · · , µ?n )T . Although this is a sparse high-dimensional problem, the
matrix X does not satisfy the sufficient conditions of the theoretical results in Zhao and Yu (2006)
and Fan and Lv (2011) due to the inconsistency of the estimation of the incidental parameters in
ν and the penalty should not be simply placed on all parameters. For details, see Supplement C.

3

In this paper, we investigate a penalized estimator of (µ? , β ? ) defined by
(µ̂, β̂) = argmin

n
X

(Yi − µi −

(µ,β)∈Rn+d i=1

X Ti β)2

+

n
X

pλ (|µi |),

(1.2)

i=1

where pλ is a penalty function with a regularization parameter λ. Since only the incidental parameters are sparse, the penalty is imposed on them. An iterative algorithm is proposed to numerically
compute the estimators. The estimator β̂ possesses consistency, asymptotic normality and an oracle property. On the other hand, the nonvanishing elements of µ? cannot be consistently estimated
even if β ? were known. So, there is a partial consistency phenomenon.
Penalized estimation (1.2) is a one-step method. For the estimation of β ? , We also propose a twostep method whose first step is designed to eliminate the influence of the data with large incidental
parameters. The estimator β̃ from the two-step method has fewer possible asymptotic distributions
than β̂ and thus is more suitable for constructing confidence regions for β ? . It is asymptotically
equivalent to the one-step estimator β̂ when the sizes of the nonzero incidental parameters are
small enough, that is, when the incidental parameters are really sparse. Also, the two-step method
improves the convergence rate and efficiency over the one-step method for challenging situations
where large nonzero incidental parameters increase the asymptotic covariance or even reduce the
convergence rate for the one-step method.
The rest of the paper is organized as follows. In Section 2, the model and penalized estimation
method are formally introduced and the corresponding penalized estimators are characterized. In
Section 3, asymptotic properties of the penalized estimators are derived; a penalized two-step
estimator is proposed and its theoretical properties are obtained; we also provide a data-driven
approach for selecting the regularization parameter. In Section 4, we consider the case where
the number of covariates grows with but slower than the sample size. In Section 5, we present
simulation results and analyze a read data set. Section 6 concludes this paper with a discussion
and all the proofs and some theoretical results are relegated to the appendix and supplements.

2

Model and Method

The matrix form of model (1.1) is given by
Y = µ? + Xβ ? + ,

(2.1)

where Y = (Y1 , Y2 , · · · , Yn )T , X = (X 1 , X 2 , · · · , X n )T , and  = (1 , 2 , · · · , n )T . The covariates
{X i }ni=1 are independent and identically distributed (i.i.d.) copies of X 0 ∈ Rd , which is a random
4

vector with mean zero and a covariance matrix ΣX > 0. They are independent of the random
errors {i }, which are i.i.d. copies of 0 , which is a random variable with mean zero and variance
σ 2 > 0. Denote an  bn and an  bn if an = o(bn ) and bn = o(an ), respectively. There is an
assumption on the covariates and random errors.
Assumption (A): There exist positive sequences κn 

√

n, γn 

√

n such that

P ( max kX i k2 > κn ) → 0 and P ( max |i | > γn ) → 0, as n → ∞,
1≤i≤n

1≤i≤n

(2.2)

where k · k2 stands for the l2 norm of Rd .
Suppose there are three types of incidental parameters in model (1.1) or (2.1): for simplicity on
1
the indexes, the first s1 incidental parameters {µ?i }si=1
are large in the sense that |µ?i |  max{κn , γn }

for 1 ≤ i ≤ s1 ; the next s2 ones {µ?i }si=s1 +1 are nonzero and bounded by γn with s = s1 + s2 ; the
last n − s ones {µ?i }ni=s+1 are zero. Note that it is unknown to us which µ?i ’s are large, bounded or
zero. The sparsity of µ? is understood by s1 + s2  n, i.e. s1 + s2 = o(n). Denote the vectors of
the three types of incidental parameters µ?1 , µ?2 , and µ?3 , respectively.
The penalized estimation (1.2) can be written as
(µ̂, β̂) = argmin L(µ, β),

L(µ, β) = kY − µ −

Xβk22

+

(µ,β)

n
X

pλ (|µi |).

(2.3)

i=1

The penalty function pλ can be the soft (i.e. L1 or LASSO), hard, SCAD or a general folded
concave penalty function (Fan and Li, 2001). For simplicity, we next consider only the soft penalty
function, that is, pλ (|µi |) = 2λ|µi |. The cases with the hard and SCAD penalties can be considered
in a similar way.
By Lemma D.1 in Supplement D, A necessary and sufficient condition for (µ̂, β̂) to be a minimizer of L(µ, β) is that β̂ = (X T X)−1 X T (Y − µ̂), Yi − µ̂i − X Ti β̂ = λSign(µ̂i ) for i ∈ Iˆ0c and
|Yi − X Ti β̂| ≤ λ for i ∈ Iˆ0 , where Sign(·) is a sign function and Iˆ0 = {1 ≤ i ≤ n : µ̂i = 0}.
Numerically, the special structure of L(µ, β) suggests a marginal decent algorithm for the minimization problem in (2.3), which iteratively computes µ(k) = argmin L(µ, β (k−1) ) and β (k) =
µ∈Rn

argmin

L(µ(k) , β)

until convergence. The advantage of this algorithm is that there exist analytic

β∈Rd

solutions to the above two minimization problems. They are respectively the soft-threshold estimators with residuals {Yi − X Ti β (k−1) } and ordinary least-squares estimator with responses Y − µ(k) .
In this section and the next, the number of the covariates d is assumed to be a fixed integer. A
case where d diverging to infinity will be considered in Section 4. For the case where d is finite, we
make the following assumption on λ.
5

Assumption (B): The regularization parameter λ satisfies
κn  λ, αγn ≤ λ, and λ  min{µ? ,

√

n},

(2.4)

where κn and γn are defined in (2.2), α is a constant greater than 2, and µ? = min1≤i≤s1 |µ?i |.
For simplicity, abbreviate “with probability going to one” to “wpg1”. A stopping rule for the above
algorithm is based on the successive difference kβ (k+1) − β (k) k2 . By Proposition D.2 in Supplement
D, wpg1, the iterative algorithm stops at the the second iteration, given the initial estimator is
bounded wpg1.
Suppose {β (k) } has a theoretical limit β (∞) , corresponding to which, there is a limit estimator
µ(∞) . Then, (µ(∞) , β (∞) ) is a solution of the following system of nonlinear equations
β = (X T X)−1 X T (Y − µ),

(2.5)

and, with soft-threshold applied to each component, it follows
µ = (|Y − X T β| − λ)+ Sign(Y − X T β),

(2.6)

where (·)+ returns the maximum value of the input and zero. By Lemma D.3 in Supplement D, a
necessary and sufficient condition for (µ̂, β̂) to be a minimizer of L(µ, β) is that it is a solution to
equations (2.5) and (2.6). Hence, (µ(∞) , β (∞) ) is a minimizer of L(µ, β) and can also be denoted
as (µ̂, β̂).
Note that β̂ is also the minimizer of the profiled loss function L̃(β) = L(µ(β), β), where
µ(β) as a function of β is given by (2.6) . Interestingly, this profiled loss function is a criterion
function equipped with the famous Huber loss function (see Huber (1964) and Huber (1973)).
Pn
T
Specifically, the profiled loss function can be expressed as L̃(β) =
i=1 ρ(Yi − X i β), where
ρ(x) = x2 I(|x| ≤ λ) + (2λx − λ2 )I(|x| > λ) is exactly the Huber loss function, which is optimal in a
minimax sense. This equivalence between the penalized estimation and Huber’s robust estimation
indicates that the penalization principle is versatile and can naturally produce an important loss
function in robust statistics. This equivalence also provides a formal endorsement of the least
absolute deviation robust regression (LAD) in Fan et al. (2012b) and indicates that it is better
to use all data with LAD regression rather than 90% of them. It is worthwhile to note that the
penalized estimation is only formally equal to the Huber’s. Our model (2.1) considers deterministic
sparse incidental parameters µ?i ’s, while the model in Huber’s works assumes random contamination
as in Kiefer and Wolfowitz (1956). Recently, there appear a few papers on robust regression in
6

high-dimensional settings, see, for example, Chen et al. (2010), Lambert-Lacroix and Zwald (2011),
? and Bean et al. (2012). Portnoy and He (2000) provide a high level review of literature on robust
statistics.
From the equations (2.5) and (2.6), β̂ is a solution to
ϕn (β) = 0, where ϕn (β) = β − (X T X)−1 X T (Y − µ(β)).

(2.7)

In general, this is a Z-estimation problem. The following theoretical analysis is based on this
characterization of β̂.
At the end of this section, we provide for further analysis some notations and an expansion of
P
P
P
P
P
ϕn (β). Let S = ni=1 X i X Ti , SS = i∈S X i X Ti , SµS = i∈S X i µ?i , SS = i∈S X i i , S = ni=1 X i
P
and SS = i∈S X i , where S is a subset of {1, 2, · · · , n}. It is straightforward to show
ϕn (β) = (SS10 + SS11 + SS12 )(β − β ? ) − (SµS11 + SµS12 )
− (SS10 + SS11 + SS12 ) − λ(SS20 + SS21 + SS22 − SS30 − SS31 − SS32 ),

(2.8)

where the index sets S10 = {s + 1 ≤ i ≤ n : |X Ti (β ? − β) + i | ≤ λ}, S11 = {1 ≤ i ≤ s1 :
|µ?i + X Ti (β ? − β) + i | ≤ λ} and S12 = {s1 + 1 ≤ i ≤ s : |µ?i + X Ti (β ? − β) + i | ≤ λ}; S20 , S21
and S22 are defined similarly except that the absolute operation is omitted and “≤” is replaced by
“>”; S30 , S31 and S32 , are defined similarly with S20 , S21 and S22 except that “> λ” is replaced
by “< −λ”. Note that all these index sets depend on β.

3

Asymptotic Properties

In this section, we consider the asymptotic properties of the penalized estimators β̂ and µ̂. Assumption (A), together with Assumption (B), enables the penalized estimation method to distinguish
the large incidental parameters from others, and thus simplifies the asymptotic properties of the
index sets Sij ’s in (2.8) in the sense that they become independent of β wpg1. Denote a hypercube
of β ? by BC (β ? ) = {β ∈ Rd : |βj − βj? | ≤ C, 1 ≤ j ≤ d} with a constant C > 0.
Lemma 3.1 (On Index Sets Sij ’s). Under Assumptions (A) and (B), for every C > 0 and every
? , S
?
?
?
β ∈ BC (β ? ), wpg1, S10 = S10
11 = ∅, S12 = S12 , S20 = ∅, S21 = S21 , S22 = ∅, S30 = ∅, S31 = S31
? = {s + 1, s + 2, · · · , n}, S ? = {s + 1, s + 2, · · · , s},
and S32 = ∅, where the limit index sets S10
1
12
? = {1 ≤ i ≤ s : µ? > 0} and S ? = {1 ≤ i ≤ s : µ? < 0}.
S21
1
1
31
i
i

7

By Lemma 3.1, wpg1, the solution β̂ to (2.7) has an analytic expression:
−1 µ

? + SS ? )
β̂ = β ? + (SS10
[SS ? + (SS10
? + SS ? ) + λ(SS ? − SS ? )],
12
21
31
12
12

(3.1)

from which, we derive asymptotic properties of β̂. Some analysis needs the following assumption.
Assumption (C): There exists some constant δ > 0 such that EkX 0 k2+δ
< ∞ and kµ?2 k2 /kµ?2 k2+δ
2

Ps
? 2+δ 1/(2+δ) .
diverges to infinity, where kµ?2 k2+δ =
i=s1 +1 |µi |
The following result shows the existence of a unique consistent estimator of β ? .
Theorem 3.2 (Existence and Consistency of β̂). Under Assumptions (A) and (B), if either s2 =
o(n/(κn γn )) or Assumption (C) holds, then, for every fixed C > 0, wpg1, there exists a unique
P

estimator β̂ n ∈ BC (β ? ) such that ψn (β̂ n ) = 0 and β̂ n −→ β ? .
In Theorem 3.2, there are two different kinds of sufficient conditions: on is on s2 , which is the
size of bounded incidental parameters µ?2 , and the other is Assumption (C), which is about the
norms of µ?2 . They come from different analysis approaches on the term SµS ? in (3.1). One does
12

not imply the other. For details, see Supplement E. Specially, if s2 =

O(nα2 )

for some α2 ∈ (0, 1)

and κn γn  n(1−α2 ) , then β̂ is consistent by Theorem 3.2.
Next, we consider the asymptotic distributions of the consistent estimator β̂ n obtained in The? = {1 ≤ i ≤ s : µ? > 0}
orem 3.2. Without loss of generality, we assume the sizes of index sets S21
1
i
? = {1 ≤ i ≤ s : µ? < 0} are asymptotically equivalent to as and (1 − a)s with a constant
and S31
1
1
1
i

a ∈ (0, 1). Similar to Theorem 3.2, there are two different sets of conditions on µ?2 corresponding
to two different analysis approaches. Denote ∼ as the asymptotic equivalence and Dn = kµ?2 k2 .
Theorem 3.3 (Asymptotic Distributions on β̂ n ). Under Assumptions (A) and (B), suppose s2 
√
n/(κn γn ) holds or Assumption (C) and Dn2 /n = o(1) hold.
√

d

n(β̂ n − β ? ) −→ N (0, σ 2 Σ−1
X ); [main case]
√
d
+
(2) If s1 ∼ bn/λ2 , then n(β̂ n − β ? ) −→ N (0, (b + σ 2 )Σ−1
X ), for every constant b ∈ R ;
√
d
(3) If s1  n/λ2 , then rn (β̂ n − β ? ) −→ N (0, Σ−1
X ), where rn ∼ n/(λ s1 ).
(1) If s1  n/λ2 , then

When the incidental parameters are really sparse, the size s1 of large incidental parameters is
small and the size s2 or the magnitude Dn of bounded incidental parameters is also small so that the
conditions of case (1) tends to hold. This case is of most interest and we denote it as the main case.
The other cases are presented to provide a relatively complete picture of the asymptotic distributions
8

of β̂. In fact, Theorem E.1 in Supplement E shows more possible asymptotic distributions. Note
that the constant a does not appear in the limit distributions of Theorem 3.3 due to cancelation and
√
that the sub- n consistency emerges in case (3) when s1 is large, because for this case the impact
of the large incidental parameters is too big to be handled efficiently by the penalized estimation.
For case (2), in one direction, as b → 0, its condition and limit distribution become those of case
(1); in the other direction, as b increases, it approaches case (3). This boundary phenomenon was
in spirit similar to that in Tang et al. (2012). Specially, if λ  nα1 , κn γn  nα2 , s1  n1−α1 and
√
d
s2  n1/2−α2 , for some α1 ∈ (0, 1) and α2 ∈ (0, 1/2), then n(β̂ n − β ? ) −→ N (0, σ 2 Σ−1
X ) by the
main case of Theorem 3.3.
Remark 1 (An Oracle Property). Suppose an oracle tells the true µ? . Then, with the adjusted
(O)

responses Y −µ? , the oracle estimator of β ? is given by β̂
= (XX T )−1 X T (Y −µ? ). The limiting
√
(O)
distribution of n(β̂ n − β ? ) is N (0, σ 2 Σ−1
X ). Comparing this with the main case of Theorems 3.3,
it follows that the penalized estimator β̂ n enjoys an oracle property.
Although mainly interested in the estimation of β ? , we also obtain the soft-threshold estimator
µ̂ of µ? : for each i,
µ̂i = µi (β̂) = (|Yi − X Ti β̂| − λ)+ sgn(Yi − X Ti β̂).

(3.2)

Denote E = {µ̂i 6= 0, for i = 1, 2, · · · , s1 ; and µ̂i = 0, for i = s1 + 1, s1 + 2, · · · , n}.
P

Theorem 3.4 (Partial Selection Consistency on µ̂). Under Assumptions (A) and (B), if β̂ −→ β ? ,
then P (E) → 1.
Theorem 3.4 shows that, wpg1, the indexes of µ?1 and µ?3 are estimated correctly, but those of
µ?2 wrongly. We call this a partial selection consistency phenomenon.

3.1

Two-Step Estimation

Theorems 3.3 shows that the penalized estimator β̂ n has multiple different limit distributions,
which complicates the application of these theorems in practice. In addition, the convergence rate
√
of β̂ n is less than the optimal rate n in the challenging cases where the impact of large incidental
parameters is substantial. To address these issues, we propose the following two-step estimation
method: firstly, we apply the penalized estimation (2.3) and let Iˆ0 = {1 ≤ i ≤ n : µ̂i = 0}; secondly,
we define the two-step estimator as
β̃ = (X TIˆ X Iˆ0 )−1 X TIˆ Y Iˆ0 ,
0

0

9

(3.3)

where X Iˆ0 consists of X i ’s whose indexes are in Iˆ0 and Y Iˆ0 consists of the corresponding Yi ’s. The
following theorem shows that β̃ is consistent and Asymptotic Gaussian.
Theorem 3.5 (Consistency and Asymptotic Normality on β̃). Suppose Assumptions (A) and
P

(B) hold. If either s2 = o(n/(κn γn )) or Assumption (C) holds, then β̃ −→ β ? . If either s2 =
√
√
d
o( n/(κn γn )) holds or Assumptions (C) and Dn2 /n = o(1) hold, then n(β̃−β ? ) −→ N (0, σ 2 Σ−1
X ).
Comparing Theorem 3.5 with Theorem 3.3, we see that β̂ has three possible asymptotic distributions but β̃ has only one since for β̃ the conditions on s1 disappear. It is because the two-step
method identifies and removes large incidental parameters by exploiting the partial selection consistency property of µ̂ in Theorem 3.4. Further, the two-step estimator improves the convergence rate
to the optimal one over the one-step estimator for the challenging case with s1  n/λ2 . Because
of these advantages, we suggest to use the two-step method to make statistical inferences.
When the incidental parameters are sparse in the sense that the size or the magnitude of the
√
bounded incidental parameters are small, i.e. s2 = o( n/(κn γn )) or Dn2 /n = o(1), it follows,
√
d
by Theorem 3.5, n(β̃ − β ? ) −→ N (0, σ 2 Σ−1
X ), from which a confidence region with asymptotic
confidence level 1 − α is given by
√
1/2
{β ∈ Rd : σ −1 nkΣX (β̃ − β)k2 ≤ qα (χd )},

(3.4)

where qα (χd ) is the upper α-quantile of χd , the square root of the chi-squared distribution with
degrees of freedom d. For each component βj? of β ? , an asymptotic 1 − α confidence interval is
given by
−1/2

[β̃j ± n−1/2 σΣX
−1/2

where ΣX

(j, j)zα/2 ],

(3.5)

(j, j) is the square root of the (j, j) entry of Σ−1
X and zα/2 is the upper α/2-quantile

of N (0, 1). The confidence region (3.4) and interval (3.5) involve unknown parameters ΣX and σ.
They can be estimated by Σ̂X = (1/n)X T X and
σ̂ = #(Iˆ0 )−1/2 kY Iˆ0 − X TIˆ β̃k2 .
0

(3.6)

By Law of Large Numbers, Σ̂X is consistent. By lemma E.3 in Supplement, σ̂ is also consistent.
Hence, after replacing ΣX and σ in the confidence region (3.4) and interval (3.5) with Σ̂X and σ̂,
the resulting confidence region and interval have the asymptotic confidence level 1 − α.

10

3.2

Theoretical and Data-Driven Regularization Parameters

By Assumption (B), the theoretical regularization parameter λ depends on κn and γn , which are also
crucial to the boundary conditions of the asymptotic properties of the penalized estimators β̂ and
β̃. By Assumption (A), κn and γn is determined by the distributions of X 0 and 0 , respectively. It
is of interest to explicitly derive κn and γn for some typical cases on the covariates and errors. When
√
the covariates are bounded with CX > 0 and the random errors follow N (0, σ 2 ), let κn = dCX
p
and γn = 2σ 2 log(n). They satisfy (2.2) in Assumption (A) and the specification of λ (2.4) in
p
√
Assumption (B) becomes α 2σ 2 log(n) ≤ λ  min{µ? , n}. When X 0 and 0 follow N (0, ΣX )
2 the maximum of diagonal elements of Σ . We can take
and N (0, σ 2 ), respectively. Denote by σX
X
q
p
p
√
2 log(n) and γ =
2σ 2 log(n). Then, (2.4) becomes log(n)  λ  min{µ? , n}. A
κn = 2dσX
n

case on exponentially tailed random variables has all been considered in Supplement E.2.
Although the theoretical specification of λ guaranties desired asymptotic properties, a datadriven specification is of interest in practice. A popular way to specify λ is to use multi-fold
cross-validation. The validation set, however, needs to be made as little contaminated as possible.
We propose the following procedure to identify a data-driven regularization parameter:
Step 1: On the training and testing data sets.
(OLS)

1. Apply ordinary least squares (OLS) to all the data and obtain residuals ˆi
X Ti β̂

(OLS)

= Yi −

for each i.
(OLS)

2. Identify the set of “pure” data corresponding to the npure smallest values in {|ˆ
i
3. Compute the updated OLS estimator β̂
residuals

(OLS,2)
{ˆ
i
}

(OLS,2)

|}.

with the “pure” data and obtain updated

for each i.
(OLS,2)

4. Identify the updated “pure” data set corresponding to the npure smallest {|ˆ
i

|} and

label the remaining as the “contaminated” data set.
5. Randomly select a subset from the updated “pure” data set as a testing set and merge
the remaining “pure” data set and the “contaminated” one into a training set.
Step 2: On the range [λL , λU ] of the regularization parameter.
1. Compute the standard deviation σ̂pure of the residuals of the “pure” data set.
2. Set λL = αl σ̂pure and λU = αu σ̂pure , where αl < αu are positive constants.
Step 3: On the data-driven regularization parameter.
11

1. For each grid point of λ in the interval [λL , λU ], apply a penalized method to the train2
ing set and obtain the estimator β̂ λ,train and the corresponding test error σ̂λ,test
=
P
T
2
testing set (Yi − X i β̂ λ,train ) .
2
2. Identify the data-driven regularization parameter λopt , which minimizes σ̂λ,test
among

the grid points.
This simple data-driven procedure can certainly be further improved. For example, In Step 1, the
sub-steps 3 and 4 can be repeated more times to obtain a better “pure” data set. In Step two, the
(OLS,2)

range for λ can also be obtained from quantiles of {|ˆ
i
(OLS,2)

on σ̂pure and quantiles of {|ˆ
i

|}. We can also hybrid quantities based

|} to determine [λL , λU ].

The good performance of this data-driven regularization parameter will be demonstrated in
Subsection 5.2.

4

Diverging number of structural parameters

In Sections 2 and 3, we have considered model (2.1) under the assumption that the number of
covariates d is a fixed integer. However, when there are a moderate or large number of covariates,
it is appropriate to assume that d diverges to infinity with the sample size. In this section, we
consider model (2.1) with the assumption that d → ∞ and d  n.
Since the number of covariates grows orderly slower than the sample size, we chose to continue
use the penalized estimation (2.3) for (µ? , β ? ) and the penalized two-step estimation (3.3) for β ? .
The corresponding estimators are still denoted as (µ̂, β̂) and β̃, but we should keep it in mind that
their dimensions diverge to infinity with n. The characterizations of β̂ in Lemmas D.1 and D.3
are still valid since they are finite-sample results. The iteration algorithm also wpg1 stops at the
second iteration, which is shown by Proposition F.3 in Supplement F.
As before, it is critical to properly specify the regularization parameter λ for the case with a
diverging number of covariates.
Assumption (B’) The regularization parameter λ satisfies
√

dκn  λ, αγn ≤ λ, and λ  µ? ,

where κn and γn are defined in (2.2) and α > 2.

12

(4.1)

Comparing Assumption (B’) with Assumption (B), the main difference in formation is that κn is
√
changed to dκn . In fact, κn in (4.1) also depends on d, which will be shown in Supplement E.2.
This difference is caused by the assumption that d diverges to ∞.
Lemma 4.1 (On Index Sets Sij ’s). Under Assumptions (A) and (B’), the conclusion of Lemma
3.1 holds.
Thus, wpg1, still valid is the crucial analytic expression of β̂ (3.1), from which we derive its
theoretical properties. They are similar to those in the previous section, with additional technical
complexity caused by the diverging dimension d.
Denote k·kF,d = d−1/2 k·kF , where k·kF is the Frobenius norm, and κX = d−1

Pd

4 1/2 ,
j=1 (E[X0j ])

the average of the square root of the fourth marginal moments of X 0 . We make the following
assumptions on ΣX and κX .
Assumption (D): kΣ−1
X kF,d is bounded.
Assumption (E): κX is bounded.
Theorem 4.2 (Existence and Consistency on β̂). Suppose Assumptions (A), (B’), (D) and (E)
hold. If there exists rd , a sequence of positive numbers depending on d, such that d3 /n → 0,
√
√
(rd d)2 /n → 0, s1 = o(n/(rd dκn λ)) and s2 = o(n/(rd dκn γn )), then, for every fixed C > 0, wpg1,
P

there exists a unique estimator β̂ ∈ BC (β ? ) such that ψn (β̂) = 0 and rd kβ̂ − β ? k2 −→ 0.
Next, we consider the asymptotic distribution on β̂. Since the dimension of β̂ diverges to infinity,
following Fan and Lv (2011), it is more appropriate to study its linear maps. Let An be a q × d
matrix, where q is a fixed integer, Gn = An ATn with the largest eigenvalue λmax (Gn ), and GX,n =
T
2
An Σ−1
X An . Denote by λmin (ΣX ) the smallest eigenvalue of ΣX , σX,max = max1≤j≤d Var[X0j ],
2
σX,min
= min1≤j≤d Var[X0j ] and γX,max = max1≤j≤d E|X0j |3 . Abbreviate “with respect to” by

“wrt”. We assume further
Assumption (D’): λmin (ΣX ) is bounded away from zero, which implies Assumption (D).
Assumption (D”): kΣX kF,d is bounded.
Assumption (F): kAn kF and λmax (Gn ) are bounded and GX,n converges to a q × q symmetric
matrix GX wrt k·kF .
Assumption (G): σX,max > 0; σX,max and γX,min are bounded from above and σX,min is bounded
away from zero.
13

Similar to the main case of Theorem 3.3, a properly scaled β̂ n is asymptotically Gaussian.
Theorem 4.3 (Asymptotic Distribution on β̂). Suppose Assumptions (A), (B’), (D’), (D”), (E),
√
√
√ √
(F) and (G) hold. If d5 log d = o(n), s1 = o( n/(λ dκn )) and s2 = o( n/( dκn γn )), then
√
d
nAn (β̂ − β ? ) −→ N (0, σ 2 GX ).
The penalized estimator µ̂ obtained by (3.2) is partially consistent.
Theorem 4.4 (Partial Selection Consistency on µ̂). Suppose Assumptions (A) and (B’)hold and
√
β̂ is a consistent estimator of β ? wrt rd k·k2 . If rd ≥ 1/ d, then P (E) → 1.
We can construct the penalized two-step estimator β̃ through (3.3) with µ̂. This two-step
estimator is consistent by Theorem F.4 in Supplement F and its asymptotic distribution, as an
extension of the main case in Theorem 3.5, is given by the following theorem.
Theorem 4.5 (Asymptotic Distribution on β̃). Suppose all the assumptions and conditions of The√
d
orem 4.3 hold except that the condition on s1 is not required. Then nAn (β̃−β ? ) −→ N (0, σ 2 GX ).
From Theorems 4.3 and 4.5, Wald-type asymptotic confidence regions of β ? are availabe. For
example, a confidence region based on β̃ with asymptotic confidence level 1 − α is given by
√
−1/2
{β ∈ Rd : σ −1 nkGX,n An (β̃ − β)k2 ≤ qα (χq )}.

(4.2)

−1

Since GX,n involves the unknown ΣX , we estimate it by ĜX,n = An Σ̂X ATn . On the other hand,
σ is estimated by σ̂ in (3.6) as before. After plugging ĜX,n and σ̂ into (4.2), we obtain
√
−1/2
{β ∈ Rd : σ̂ −1 nkĜX,n An (β̃ − β)k2 ≤ qα (χq )}.

(4.3)

By Lemma A.5 in the appendix, the consistency of σ̂ is assured. Then, Theorem A.6 in the appendix
guarantees the asymptotic validity of the confidence region (4.3).

5

Numerical Evaluations and Real Data Analysis

In this section, we evaluate the finite-sample performance of the penalized estimation through
simulations and use it to analyze a real data set. The model for simulations is given by, for
? +  . For simplicity, the deterministic sparse
i = 1, · · · , n, Yi = µ?i + Xi,1 β1? + · · · + Xi,50 β50
i

incidental parameters {µ?i } are generated as i.i.d. copies of µ: µ is 0, U [−c, c] and W (c + Exp(τ ))
with probabilities p0 , p1 , and p2 respectively, where U [−c, c] is a uniform random variable in [−c, c],
14

W takes values 1 and −1 with probabilities pw and 1 − pw , respectively, and Exp(τ ) follows an
exponential distribution with mean 1/τ > 0. Note that c can be viewed as a contamination
parameter. The larger c is, the more contaminated the data. On the other hand, pw determines
? = 1;
the asymmetry of the incidental parameters. The regression coefficients β1? = · · · = β50
i.i.d.

{(Xi,1 , · · · , Xi,50 )} ∼ N (0, ΣX ), where ΣX (i, j) = 2 exp(−|i − j|)), which is a Toeplitz matrix
and the constant 2 is used to inflate the covariance a little; the covariates are independent of
i.i.d.

{i } ∼ N (0, 1); and n = 500; p0 = 0.8, p1 = 0.1, p2 = 0.1, c is 0.5, 1, 3 or 5, pw is 0.5 or 0.75 and
τ = 1.

5.1

Performance of Penalized Methods

The following methods for estimating β ? are evaluated. (i) Oracle method (O): an oracle knows the
index set S of zero µ?i ’s. Its performance is used as a benchmark. (ii) Ordinary least squares method
(OLS): all µ?i ’s are thought as zeros. (iii) Four penalized least squares (PLS) methods, namely, PLS
with soft penalty (PLS.Soft or S), PLS with hard penalty (PLS.Hard or H), two-step PLS with soft
penalty (PLS.Soft.TwoStep or S.TS) and two-step PLS with hard penalty (PLS.Hard.TwoStep or
P
P
(O)
H.TS). More specifically, the oracle estimator of β ? is given by β̂
= ( i∈S X i X Ti )−1 i∈S X i Yi .
The hard penalty function is pλ (|t|) = λ2 −(|t|−λ)2 {|t| < λ} (see Fan and Li (2001)). Each method is
evaluated by the square root of the empirical mean squared error (RMSE). Every penalized method
is evaluated with a grid of values for the regularization parameter λ, ranging from .5 to 5 by .25.
The sequence plot of Figure 1 shows 500 realized incidental parameters µ?i ’s with c = 3 and
pw = 0.75. They are used in data generation of simulations. The scatter plot of Figure 1 shows
the responses Yi ’s against the first covariate Xi1 ’s of a generated data set and the red stars stand
for the contaminated sample points, the ones with nonzero µ?i ’s. With fifty covariates, usually it is
difficult to graphically identify the contaminated data points.
With the above incidental parameters, those six methods are evaluated by simulations with
iteration number 1000. Since ΣX is a Toeplitz matrix with equal diagonal elements, the asymptotic
variances of the estimators of β1? and β2? are different and representative for estimators of other
βi? ’s. So, we only report simulation results on the estimation of β1? and β2? .
Figure 2 shows RMSE’s of six estimators for β1? . RMSE’s for β2? are similar. As expected,
the oracle method has the smallest RMSE and OLS the largest. RMSE of each PLS method as a
function of λ forms a convex curve, which achieves a minimal RMSE significantly below the green
line of OLS and close to the cyan line of O. More specifically, RMSE of PLS.Hard achieves the
15

Marginal Data Plot
40

Sequence Plot of Incidental Parameters
●

●

●
●

●
●

●
● ●

●

●

●
● ● ●

●

●

●
●
●
●●

●

●

●
●

●

●

●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●

●

●

●
●
●

●

●

●

●

●

●

●●

●

●

●

●

●
●

●

●

●
●

●
● ●
●

●

●

0

●

●

−40

−5

●
●

●

●

●

●
●

●

0

●

−20

0

mu

●

●

●

100

200

300

●

●

● ● ●

●

●

●
●
●
●
●
●●
●
●
●●●●● ● ●
●
●
●
●
●
●
●●
●●
●
●
●
●
● ●●
●
●
●●
●
●
●
●●●
● ●
●
●
● ●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●●
● ● ●
●
●
● ●
●
●
●
● ●●
● ●●
● ●
● ●●● ●●
●●
●
●
●
● ●
●
●
●
●
●
●
●
● ●
●
● ●
●●
● ●●●
●
● ●● ● ● ● ●● ●●
●● ●●
●● ●
●
●
●● ● ● ● ●●●
●
●
● ● ●
●
●
●
●
●
●
●
●
●
●●
●
●
●● ●
●●
●
● ●
●
●
● ●
●
●
●
●●
● ●
● ●● ● ●
●
●
● ●
●
● ●
●
● ●
● ●● ●●
● ●
● ● ●●●
●●
●
●
●● ●●●● ●● ●● ●
● ● ●● ● ● ● ●
● ●●
●●
●● ●
●
●●
●
●
●●
●● ●
● ●●●●●
●●
●
● ●
●
●●
●
●
●●●
●
●
●
● ● ●●
●
● ● ● ● ●●●
●
●
●●
●
● ● ●
●
●
●
●
● ●● ●
● ●
●
●
●
●
●● ● ●
●
●
●
●
●
● ●
●
●
●
● ●
●
● ● ●●●
●
●●
●
●
●
● ●●
●
●
● ●
●
●
●●● ●
●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●

●●

●

●

Y

5

●

●

●

20

●

●
●
●●

●
●

●
●
●

●

●

●

●

●

400

500

−4

−2

Index

0

2

4

X1

Figure 1: The sequence plot on the left shows 500 incidental parameters µ?i ’s with c = 3 and pw = 0.75. The
non-zero µ?i ’s are in red. The scatter plot on the right shows the responses Yi ’s against the first covariate
Xi1 ’s of a data set generated with those 500 incidental parameters. The red stars stand for the contaminated
sample points, the ones with nonzero µ?i ’s.

minimal RMSE when λ is around 2.75. On the other hand, RMSE of PLS.Soft decreases a little
till λ is around .75, then increases and stays above RMSE of PLS.Hard. This reflects the fact
that a large λ in a soft-threshold method usually causes bias. PLS.Hard.TwoStep has very similar
performance with PLS.Hard for all λ. PLS.Soft.TwoStep has similar performance with PLS.Soft
when λ is small. However, as λ becomes large, PLS.Soft.TwoStep moves closer to PLS.Hard than
PLS.Soft. It is because PLS.Soft.TwoStep and PLS.Hard.TwoStep have similar estimation when λ
is large. The minimal RMSE of PLS.Soft is slightly larger than those of other PLS Methods.
Table 1 depicts the minimal RMSE’s of the estimators for β1? and β2? with the corresponding optimal λ’s and biases. The biases are ignorable comparing withe the RMSE’s. The optimal λ for PLS.Soft and other PLS methods are around .75 and 2.5, respectively. This indicates
the simple soft threshold method tends to work best with a small λ due to the bias issue. Denote the empirical relative efficiency (ERE) of an estimator A with respect to another estimator
B as RMSE(B)/RMSE(A). Then, for the estimation of β1? , the ERE’s of PLS.Soft, PLS.Hard,
PLS.Soft.TwoStep and PLS.Hard.TwoStep with respect to O are around 78%, 87%, 87% and 87%,
respectively; the ERE’s of the PLS methods with respect to OLS are around 176%, 196%, 196%
and 196%, respectively. The ERE’s for β2? are similar. Thus, in terms of ERE (and RMSE), the
PLS methods perform closely to O and significantly better than OLS.
16

RMSE's for beta1 with Deterministic Sparse Inicidental Parameters
OLS

●

S

H

S.TS

H.TS

0.08

0.09

O

●
●
●

0.07

●

●

●

●
●

●

●

●

●

●

●

●

●

0.04

0.05

●

●

0.06

RMSE

●

1

2

3

4

5

lambda

Figure 2:

RMSE’s of the O(Oracle), OLS, S(PLS.Soft), H(PLS.Hard), S.TH(PLS.Soft.TwoStep) and

H.TH(PLS.Hard.TwoStep) estimators of β1? with the incidental parameters shown in Figure 1. The top and
bottom solid horizontal lines show the RMSE’s for OLS and O, respectively. Other four horizontal lines
indicate the minimal RMSE’s for those four PLS methods and the corresponding four best λ’s are shown by
the vertical lines.

From Table 1, we can also see that the RMSE’s of the estimators for β1? are always smaller than
those for β1? . This is because the first covariate is less correlated with others covariates than the
second one.
In order to examine the performance of the methods with general incidental parameters but
not just those in Figure 1, we generate µ? randomly for each iteration. The iteration number for
each simulation is also 1000.
Figure 3 shows the RMSE’s of six estimators of β1? under two settings: pw = 0.5 and c = 1 or
5. Each plot in Figure 3 presents a similar pattern with Figure 2. When pw is fixed at 0.5, the
RMSE’s of each non-oracle estimator of β1? increases as the contamination parameter c increases
from 1 to 5. This indicates that each non-oracle estimator performs worse as the data becomes more
contaminated. However, the PLS estimators are more robust than OLS, which is very sensitive to
the change of c. We have also done simulations with pw = 0.75 and the RMSE’s of the estimators
of β1? are similar to those with pw = 0.75 so that the corresponding plots are similar to those in
Figure 3. In other words, the RMSE’s of all estimators are stable with respect to pw , which means

17

O

OLS

S

H

S.TS

H.TS

S.P

H.P

LAD

Bias(β̂1 )×104

.8

17.8

1.3

.6

5.9

20.4

9.5

11.5

10.4

RMSE(β̂1 )×102

4.0

8.4

5.1

4.6

4.6

4.6

5.4

5.0

5.4

.75

2.75

2.25

2.5

2.45

2.47

λ
Bias(β̂2 )×104

-2.3

-46.7

24.6

43.5

-21.8

10.0

32.3

13.9

-7.3

RMSE(β̂2 )×102

4.4

9.0

5.3

4.9

5.0

4.9

5.8

5.5

5.9

.75

2.75

2.75

2.5

2.45

2.47

λ

Table 1: RMSE’s of the Oracle, OLS and LAD estimators and The minimal RMSE’s of the penalized
estimators of β1? and β2? with the corresponding optimal or data-driven λ’s and biases when the incidental
parameters shown in Figure 1 are used. For a data-driven method, the data-driven λ is different in each
iteration so that the reported λ’s are averages. The lines over the numbers emphasize the numbers are
averages. The standard deviations for the data-driven λ’s of S.P and H.P are 0.37 and 0.34, respectively.

RMSE's for beta1 with Random Sparse Incidental Parameters

c = 1 , p.w = 0.5

●

1

●

●

●

2

●

●

●

●

●

3

●

●

●

4

●

●

●

0.08
0.06

●
●

●
●

●

●

●
●

●

●

●

●

●

●

●

●

●

0.04

0.06

RMSE

0.04

●

●

●
●
●

RMSE

0.08

O
OLS
S
H
S.TS
H.TS

●

●

c = 5 , p.w = 0.5

5

1

lambda

2

3

4

5

lambda

Figure 3: Similar to Figure 2, these plots show the RMSE’s of the O(Oracle), OLS, S(PLS.Soft),
H(PLS.Hard), S.TH(PLS.Soft.TwoStep) and H.TH(PLS.Hard.TwoStep) estimators of β1? with randomly
generated µ? under two settings with pw = 0.5 and c = 1 or 5.

18

RMSE(βˆ1 )×102

O

OLS

S

H

S.TS

H.TS

S.P

H.P

LAD

(pw , c) = (.5, .5)

4.06

4.06

3.92

3.95

3.90

3.96

4.01

4.38

4.85

1.25

3.5

2.75

2.5

2.08

2.15

4.11

4.27

4.21

4.27

4.01

4.59

2

3.75

4.25

3.5

2.20

2.24

4.81

4.99

4.81

4.80

5.64

5.19

1

2.5

2

2.25

2.42

2.45

4.96

4.63

4.66

4.64

6.36

4.87

1

2.25

2.75

3

2.52

2.58

3.91

3.88

3.98

3.91

4.08

4.41

1.5

5

3.25

4.5

2.08

2.14

4.19

4.15

4.15

4.20

4.16

4.59

3.25

3

4

3.5

2.15

2.23

4.91

4.93

4.80

5.02

5.35

5.06

1

2.25

2

2.5

2.44

2.47

5.01

4.66

4.75

4.66

6.22

4.90

1.5

2.25

2.5

3

2.55

2.60

λ
(pw , c) = (.5, 1)

4.06

4.58

λ
(pw , c) = (.5, 3)

4.12

6.47

λ
(pw , c) = (.5, 5)

4.07

8.50

λ
(pw , c) = (.75, .5)

4.13

4.02

λ
(pw , c) = (.75, 1)

4.17

4.41

λ
(pw , c) = (.75, 3)

4.15

5.99

λ
(pw , c) = (.75, 5)
λ

3.97

8.41

4.75

5.49

5.52

4.79

4.98

5.52

5.76

Table 2: Similar to Table 1, this one shows the RMSE’s and minimal RMSE’s of nine estimators of β1? under
eight settings on randomly generated µ? with pw = 0.5, 0.75 and c = 0.5, 1, 3, 5. The standard deviations of
the data-driven λ’s of S.P and H.P for different settings are between 0.2 and 0.45.

the magnitudes of the nonzero incidental parameters matter most but not their signs. We also note
that some penalized methods perform closely to or even outperform the oracle one when c is small
as showed in the plot with c = 1. This happens because that O ignores all the contaminated data
points, even those with very light contamination, but the penalized methods exploit information
in such points.
Table 2 contains the RMSE’s of the estimators of β1? under eight settings with pw = 0.5 or
0.75 and c = 0.5, 1, 3 or 5. For each pw , as c increases from 0.5 to 5, the RMSE’s of O with the
multiplication factor 102 is almost constantly around 4, those of OLS increases from about 4 to 8.5
and those of PLS ones grow from about 4 to 5, which confirms the robustness of the PLS estimators.
When c ≤ 1 is small with respect to the variance of random error σ = 1, the data points are only
slightly contaminated. OLS and PLS methods perform similar to O. However, when c ≥ 3 is large,
19

which means the data are more contaminated, the RMSE’s of OLS become significantly larger, but
PLS methods perform still closely to O.

5.2

Performance of Data-Driven Penalized Methods

Previous simulations have shown the PLS methods with optimal λ’s have good RMSE’s comparing
with those of the Oracle and OLS ones. In practice, however, these optimal λ’s are unknown. One
approach to obtain a data driven λ has been introduced in Subsection 3.2. Since, as shown in the
previous simulation results, the two-step PLS methods perform similarly with the one-step PLS
methods, i.e. PLS.Soft and PLS.Hard, only the latter are studied by simulations with data-driven
λ’s and denoted as PLS.Soft.Prac (S.P) and PLS.Hard.Prac (H.P), respectively. For estimating
data-driven λ’s, αl = 2 and αu = 7. The size of the pure data set npure is n/2 and that of the
testing data set is npure /2.
Simulations are first run with the deterministic sparse incidental parameters as showed in
Figure 1. We can see in Table 1 that the RMSE’s of estimators of β1? from PLS.Soft.Prac and
PLS.Hard.Prac are around 5.4 and 5.0, slightly larger than the optimal values 5.1 and 4.6.,respectively. However, they are still significantly smaller than RMSE of OLS, which is 8.4. The
observations of the estimators of β2? are similar. As before, we also evaluate the performance of
the data-driven PLS methods with random sparse incidental parameters. Table 2 shows that, for
a given pw , when c is small such as 0.5 and 1, the RMSE’s of PLS.Soft.Prac and PLS.Hard.Prac
are close to those of PLS.Soft and PLS.Hard with the optimal λ’s. In these cases, PLS.Soft.Prac
performs slightly better than PLS.Hard.Prac, and even better than PLS.Soft with the optimal λ
and the Oracle method. On the other hand, for a given pw , when c is large such as 3 and 5, the
RMSE’s of PLS.Soft.Prac and PLS.Hard.Prac are greater than those of PLS.Soft and PLS.Hard,
respectively, but still less than those of OLS. In these cases, the RMSE’s of PLS.Soft.Prac are larger
than those of PLS.Hard.Prac, which indicates the bias issue of the soft threshold method. Thus,
the data-driven regularization parameter works well with penalized estimation. When the data is
slightly contaminated, the soft penalty is preferred; otherwise, the hard penalty is recommended.
Tables 1 and 2 also contain RMSE of the least absolute deviation regression method (LAD)
used in Fan et al. (2012b) with all but not part of the sample points with small residuals. Generally speaking, in both deterministic and random incidental parameter cases, LAD performs
similarly with the PLS methods with data-driven λ’s. More specifically, when c ≤ 1 is small,
PLS.Soft.Prac outperform LAD; otherwise, LAD performs better. For all the cases, LAD is domi20

nated by PLS.Hard.Prac. These observations confirm that LAD is an effective robustness method
and the penalized methods make improvement.

5.3

Data-Driven Confidence Intervals

We next turn to investigate the finite-sample performance of the asymptotic confidence interval (CI)
(3.5) for βj? with j = 1, 2 based on PLS two-step methods. Since CI (3.5) is based on the properties of
the penalized two-step estimator with the soft penalty, we focus on PLS.TS.Soft with a data-driven
regularization parameter λ. The choice of λ in Subsection 3.2 for minimizing RMSE is usually no
longer suitable for constructing confidence intervals, since it is designed to achieve minimal RMSE.
We propose to first obtain σ̂pure as in the data-driven procedure in Subsection 3.2 and then simply
set the data-driven λ be five times of σ̂pure . Since σ̂pure tends to underestimate σ, this data-driven
λ is usually not large with respect to σ. Denote this method as PLS.TwoStage.Soft.Prac or S.TS.P.
−1

After plugging in σ̂ and σ̂j−1 , the square root of the (j, j)th element of Σ̂X , and replace n by
m = #(Iˆ0 ) in the theoretical CI (3.5), we obtain a data-driven CI [β̃j ± m−1/2 σ̂σ̂j−1 zα/2 ], where β̃j
is the PLS.TwoStage.Soft.Prac estimator of βj? for each j.
This data-driven CI is compared with CI’s based on Oracle and OLS methods. More specif(O)

ically, denote the Oracle and OLS estimators of βj? as β̂j
(O)

corresponding CI’s are given by [β̂j

−1/2 (O) −1
σ̂ σ̂j zα/2 ]

± mo

(OLS)

and β̂j

, respectively. Then, the

(OLS)

and [β̂j

± n−1/2 σ̂ (OLS) σ̂j−1 zα/2 ],

where mo is the number of zero incidental parameters and σ̂ (O) and σ̂ (OLS) are the estimators of σ
from O and OLS methods, respectively.
The simulation settings are the same to the previous ones with deterministic sparse incidental
parameters except the following changes. (a) The number of covariates d is reduced to 5 from 50.
This is because when d = 50 and the nominal level is 95%, even the empirical coverage rate (CR)
of the oracle confidence interval for β1? becomes 93.5%, not very close to 95%. (b) The iteration
number is increased from 1000 to 10000 to improve the accuracy of CR’s. (c) The probabilities of
nonzero incidental parameters (p1 , p2 ) are set to be (0.01, 0.01), (0.03, 0.03) and (0.05, 0.05); the
contamination parameter c is increased to 10. In order to achieve good second order asymptotic
approximation, we can either increase the sample size or enlarge the signal noise ratio. Here we
adopt the latter.
Table 3 reports the empirical coverage rates (CR) and average lengths (AL) of the CI’s of β1? and
β2? from O, OLS and PLS.TS.Soft.Prac methods under three different settings on the incidental
parameters. For the oracle method, these three settings are the same and thus only one set of
21

(p1 , p2 )

β1?

O

OLS

S.TS.P

.950

.944

.951

.948

(.05,.05)

.946

(.01,.01)

(.01,.01)
(.03,,03)

(.03,,03)

CR

AL

(.05,.05)

.955

.133

β2?

OLS

S.TS.P

.948

.945

.946

.947

.945

.949

.946

.200

.135

.213

.144

.279

.137

.297

.146

.382

.139

.407

.149

CR

AL

O

.948

.142

Table 3: Coverage rates (CR) and average length (AL) of 95% confidence intervals for β1? and β2? from O,
OLS, PLS.TwoStage.Soft.Prac methods under three settings on deterministic sparse incidental parameters.

simulation results are presented. Table 3 shows that the CR’s of all methods under all settings are
close to the nominal level .95. The OLS treats the deterministic incidental parameters as random
ones and achieves excellent CR’s. However, the AL’s of OLS are significantly larger than those
of O and PLS.TS.Soft.Prac, especially when there are more non-zero incidental parameters. On
the other hand, the AL’s of PLS.TS.Soft.Prac are only slightly larger than those of O. This means
PLS.TS.Soft.Prac has excellent efficiency in terms of AL’s given excellent CR’s. Also note that
the AL’s for β1? are less than those for β2? . This is because the asymptotic variance of β̂1 is less
than that of β̂2 when the covariance matrix ΣX is a Toeplitz matrix. Simulations with random
incidental parameters under the same settings have also been done and the results are similar to
those in Table 3 with slightly inflated AL’s for OLS and PLS.TS.Soft.Prac due to the randomness
of the incidental parameters.

5.4

Real Data Analysis

We implement the penalized estimation with the soft penalty in the method of estimating false
discovery proportion of a multiple testing procedure proposed by Fan et al. (2012b) for investigating the association between the expression level of gene CCT8, which is closely related to Down
Syndrome phenotypes, and thousands of SNPs. The data set consists of three populations: 60 Utah
residents (CEU), 45 Japanese and 45 Chinese (JPTCHB) and 60 Yoruba (YRI). More details on
the data set can be found in Fan et al. (2012b).
In the testing procedure by Fan et al. (2012b), a filtered least absolute deviation regression
(LAD) is used to estimate the loading factors with 90% of the cases (SNPs) whose test statistics
are small and thus the resulting estimator is statistically biased. We upgrade this step with S.P
22

10

●

●
●
●

●
●
●
●
●

3.0

3.5

●
●

2.0

2.5

3.0

●

0.8

50

●

0

0

2.5

●

●
●

2.0

●
●
● ●

●

●

●
●

●
●
●
●
●

0.4

● ●
●

●● ●
●

Estimated FDP(t)

●
●

●● ●
●●●
●
●●
●●

LAD
S.P

●

●

0.0

●
● ●

30

R(t)

●

●
●
●
●
●
●
●
●

30

Estimated V(t)

●
●
●

10

50

CEU
●
●
●●
●
●●
●●
●
●●

3.5

2.0

2.5

3.0

3.5

50

●
●● ●
●

●

●
●
●

20

●

●

●
●

10

●

●

●
●
●
●

●●
●
●
●

●

●

●

●

● ●●
●

●
●

0.8

●

●
●
●
●

●
●

●

0.4

●

Estimated FDP(t)

●

Estimated V(t)

30

●

●
●

30

40

●
●
●

40

●
●

●

●

0.0

R(t)

●●
●●

20

●

10

50

JPTCHB
●
●
●

●

4

5

6

7

8

●

0

●

9

4

5

6

7

8

9

4

5

6

7

8

9

●●
●

●

●

● ●
●

●

●●
●

0

5

6
−log(t)

7

●
●●
●●

●●
●●

●

● ●
●

●

8

●●

●● ●

●

●●

4

5

6
−log(t)

7

8

●

●●

0.8

●
●●
● ●
●
●● ●
●
●
●
●

● ● ●

●
●
●●
●
●
●
●
●

0.4

80 100
60

●
●

●
●●
●
●

4

Estimated FDP(t)

●
● ●
●
●●
●●

●
●
●
●
●
●
●
●
●

●●●
●
●●
●
●
●
●● ●
●
●●
●
●
●●
●
●●
●
●
●
●
●
●

●

●

0.0

●
●

●
●
●
●
●
●

40

Estimated V(t)

●
●
●
●●
●
●
●

●
●
●

0

40

60

●
●●
●
●

●●
●
●
●
●●
●

20

●
●

20

R(t)

80 100

YRI
●●
●
●
●
●●
●●
●
●●

4

5

6

7

8

−log(t)

Figure 4: Discovery number R(t), estimated false discovery number V (t) and estimated false discovery
proportion FDP(t) as functions of a threshold on t for populations CEU, JTPCHB and YRI. The x-axis is
− log10 (t).

23

Population

t

R(t)

[
FDP(t)
with LAD

[
FDP(t)
with S.P

CEU

6.12 × 10−4

4

.810

.845

JPTCHB

1.51 × 10−9

5

.153

.373

YRI

2.54 × 10−9

2

.227

.308

ˆ
Table 4: Discovery numbers R(t) and estimated false discover proportions FDP(t)s
from methods with LAD
and S.P for specific values of threshold t.

described in Subsections 3.2 and 5.2 and re-estimate the number of false discoveries V (t) and the
false discovery proportion FDP(t) as functions of − log10 (t), where t is a thresholding value. Figure
[
4 shows the number of total discoveries R(t), V̂ (t) and FDP(t)
from procedures using filtered LAD
[
and S.P. It is clear that V̂ (t) and FDP(t)
with S.P are uniformly larger than but reasonably close
[
to those with filtered LAD. Table 4 contains R(t) and FDP(t)
with filtered LAD and S.P for several
specific thresholds. The estimated FDPs with S.P for CEU and YRI are slightly larger than those
[ for JPTCHB with S.P is more than double of that with filtered LAD. This
with LAD and FDP
suggests that the estimation of FDP with filtered LAD might tend to be optimistic.

6

Conclusion and Discussion

This paper considers the estimation of structural parameters with a finite or diverging number
of covariates in a linear regression model with the presence of high-dimensional sparse incidental
parameters. By exploiting the sparsity, we propose an estimation method penalizing the incidental
parameters. The penalized estimator of the structural parameters is consistent and asymptotically
Gaussian and achieves an oracle property. On the contrary, the penalized estimator of the incidental
parameters possesses only partial selection consistency but not consistency. Thus, the structural
parameters are consistently estimated while the incidental parameters not, which presents a partial
consistency phenomenon. Further, in order to construct better confidence regions for the structural
parameters, we propose a two-step estimator, which has fewer possible asymptotic distributions and
can be asymptotically even more efficient than the one-step penalized estimator when the size and
magnitude of nonzero incidental parameters are substantially large.
Simulation results show that the penalized methods with best regularization parameters achieve
significantly smaller mean square errors than the ordinary least squares method which ignores the
incidental parameters. Also provided is a data-driven regularization parameter, with which the

24

penalized estimators continue to significantly outperform ordinary least squares when the incidental
parameters are too large to be neglected. In terms of average length together with excellent
coverage rates, the advantage of the confidence intervals based on the two-step estimator with an
alternative data-driven regularization parameter is verified by simulations. A data set on genomewide association is analyzed with a multiple testing procedure equipped with a data-driven penalized
method and false discovery proportions are estimated.
In econometrics, a fixed effect panel data model is given by, for 1 ≤ i ≤ n and 1 ≤ t ≤ T ,
Yit = µ?i + X Tit β ? + it ,

(6.1)

where µ?i ’s are unknown fixed effects. When T diverges, the fixed effects can be consistently
estimated. When T is finite and greater than or equal to 2, although the fixed effects can no
longer be consistently estimated, they can be removed by a within-group transformation: for each
i, Yit − Ȳi = (X it − X̄ i )T β ? + it − ¯i , where Ȳi , X̄ i and ¯i are the averages of Yit ’s, X it ’s and it ’s,
respectively. When T is equal to 1, however, the within-group transformation fails. Note that,
with T = 1, Model (6.1) becomes Model (1.1) so that the proposed penalized estimations provide
a solution under the sparsity assumption on the fixed effects.
Although this paper only illustrates the partial consistency phenomenon of a penalized estimation method for a linear regression model, such a phenomenon shall universally exist for a general
parametric model, which contains both a structural parameter and a high-dimensional sparse incidental parameter. For example, consider a panel data logistic regression model: P (Yit = 1|X it ) =
(1 + exp{−(µ?i + X Tit β ? )})−1 . When T is finite, the fixed effects µ?i ’s cannot be removed by the
within-group transformation as in the panel data linear model (6.1). However, the proposed penalized estimations can still provide a solution.
Further, if the structural parameter has a dimension diverging faster than the sample size and
is sparse, it is expected that the partial consistency phenomenon will continue to appear when
sparsity penalty is imposed on both the structural and incidental parameters.

Acknowledgement
The authors thank the Editor, an associate editor, and three referees for their many helpful comments that have resulted in significant improvements in the article. The research was partially
supported by NSF grants DMS-1206464 and DMS-0704337 and NIH Grants R01-GM100474-01
and R01-GM072611-08.
25

A

Appendix

In this appendix, we provide the proofs of the theoretical results in Section 4. The proofs of the
results in Sections 2 and 3 are in Supplements D and E.
Denote Sk,l = S{k,k+1,··· ,l} and Sk,l = S{k,k+1,··· ,l} . Let B = {maxs+1≤i≤n kX i k2 ≤ κn } and
T
D = ni=1 {−γn ≤ i ≤ γn }. Then P (B) → 1 and P (D) → 1 by (2.2).
Proof of Lemma 4.1. We first consider Si0 ’s, then Si1 ’s, and finally Si2 ’s with i = 1, 2, 3. Consider
? }. Note that P (A) ≥ P (A|B)P (B) and P (B) → 1. It
S10 , S20 and S30 . Let A = {S10 = S10
√
suffices to show that P (A|B) → 1. By λ  dκn , it follows P (A|B) ≥ P ({s + 1 ≤ i ≤ n :
√
√
? |B) ≥ P ({s + 1 ≤ i ≤
−λ + maxs+1≤i≤n kX i k2 dC ≤ i ≤ λ − maxs+1≤i≤n kX i k2 dC} ⊃ S10
√
√
?
? ) ≥ P (D) → 1. Thus, wpg1, S
n : −λ + κn dC ≤ i ≤ λ − κn dC} ⊃ S10
10 = S10 . From
? , it follows that, wpg1, S
S10 ∪ S20 ∪ S30 = S10
20 = S30 = ∅. Consider S21 , S31 and S11 . Recall
√
that µ? = min{|µ?i | : 1 ≤ i ≤ s1 } and note that λ − µ? + dCκn < −γn when n is large. Let
? and S
?c
?
S211 = S21 S21
212 = S21 S21 . We will show P (S211 = S21 ) → 1 and P (S212 = ∅) → 1.
? ) → 1. Denote A = {S
?
Then P (S21 = S21
1
211 ⊃ S21 }. On the event B, S211 ⊃ {1 ≤ i ≤ s1 :
√
i > λ − µ? + dCκn and µ?i > 0} ⊃ {1 ≤ i ≤ s1 : i > −γn and µ?i > 0}. Then, P (A1 ) ≥
? )P (B) → 1 · 1 = 1. It follows that,
P (A1 |B)P (B) ≥ P ({1 ≤ i ≤ s1 : i > −γn and µ?i > 0} ⊃ S21
? . Note that S
?
?
wpg1, S211 ⊃ S21
211 ⊂ S21 . Then, wpg1, S211 = S21 . Denote A2 = {S212 = ∅}.
√
On the event B, S212 ⊂ {1 ≤ i ≤ s1 : i > λ + µ? − dCκn and µ?i < 0}, which contains

{1 ≤ i ≤ s1 : i > γn }. Then, P (A2 ) ≥ P (A2 |B)P (B) ≥ P ({1 ≤ i ≤ s1 : i > γn } = ∅)P (B) →= 1.
? ) → 1. Similarly, we can show, wpg1, S
?
Then, wpg1, S212 = ∅. Thus, P (S21 = S21
31 = S31 .
? ∪ S ? . Then, wpg1, S
Note that S11 , S21 and S31 are disjoint and their union is S21
11 = ∅.
31
√
? }. Note that −λ − µ? + dCκ < −γ and
Consider S12 , S22 and S32 . Denote A = {S12 = S12
n
n
i
√
λ − µ?i − dCκn > γn when n is large for s1 + 1 ≤ i ≤ s. On the event B, S12 ⊃ {s1 + 1 ≤ i ≤
√
√
s : −λ − µ?i + dCκn ≤ i ≤ λ − µ?i − dCκn }, which contains {s1 + 1 ≤ i ≤ s : −γn ≤ i ≤ γn }.
? )P (B) → 1. Thus,
Then, P (A) ≥ P (A|B)P (B) ≥ P ({s1 + 1 ≤ i ≤ s : −γn ≤ i ≤ γn } = S12
? . Note that S , S
?
wpg1, S12 = S12
12
22 and S32 are disjoint and their union is S12 . Then, wpg1,

S22 = S32 = ∅.
2 = (1/d)
Before proceeding to the proofs of Theorems 4.2 to A.6, we denote σ̄X
P
P
2
and σ̄XX
= (1/d2 ) dk=1 dl=1 Var[X0k X0l ] and make the following assumptions.
2 is bounded.
Assumption (E1): σ̄X

26

Pd

j=1 Var[X0j ]

2
Assumption (E2): σ̄XX
is bounded.

Assumption (E) in Section 4 implies Assumptions (E1) and (E2) by Cauchy-Schwarz inequality.
For simplicity, we adopt the notation ., which means the left hand side is bounded by a constant
times the right, where the constant does not affect related analysis. Below are three lemmas needed
for proving Theorems 4.2 to A.6. Their proofs are in Supplement F. Suppose that M and E are
matrices and k·k is a matrix norm and that {An } is a sequence of random d × d matrices and A a
deterministic d × d matrix, and denote Σ̂n = (1/n)Sn , the sample covariance matrix.
Lemma A.1 (Stewart (1969)). If kIk = 1 and kM −1 kkEk < 1, then

Lemma A.2. If kA−1 kF,d

k(M + E)−1 − M −1 k
kM −1 kkEk
≤
.
kM −1 k
1 − kM −1 kkEk
√
P
P
−1
is bounded, An −→ A, and rd ≥ 1/ d, then A−1
n −→ A , where the

convergence in probability is wrt rd k·kF .
P

Lemma A.3. If Assumption (E2) holds and rd2 d4 /n → 0, then Σ̂n −→ ΣX wrt rd k·kF .
Proof of Theorem 4.2. By the proof of Lemma 4.1, wpg1, the solution β̂ n to ϕn (β) = 0 on
BC (β ? ) is explicitly given by β̂ n = β ? + T0−1 (T1 + T2 + T3 − T4 ), where T0 = (1/n)Ss1 +1,n ,
?
? and T4 = (λ/n)SS ? . Then, rd kβ̂ − β k2 ≤
T1 = (1/n)SµS ? , T2 = (1/n)Ss1 +1,n , T3 = (λ/n)SS21
n
31
12
√
P
kT0−1 kF,d 4i=1 rd dkTi k2 . We will show that kT0−1 kF,d is bounded by a positive constant wpg1
√
P
and rd dkTi k2 −→ 0 for i = 1, 2, 3, 4. Then, rd kβ̂ n − β ? k2 = oP (1). Consider T0 . By Lemma
P

A.3, kT0 − ΣX kF,d −→ 0 under Assumption (E2) and the condition d3 /n → 0. Then, by Lemma
P

−1
A.2, together with Assumption (D), kT0−1 − Σ−1
X kF,d −→ 0. This implies that, wpg1, kT0 kF,d is
√
√
bounded by a positive constant. Consider T1 . Wpg1, rd dkT1 k2 ≤ rd ds2 κn γn /n = o(1) for s2 =
√
P
o(n/(rd dκn γn )). Consider T2 . For any δ > 0, P (kT2 k2 > δ) ≤ (1/δ 2 )P k(1/n) ni=s1 +1 X i i k22 ≤
√
Pd
2 /(nδ 2 ), where σ̄ 2 = (1/d)
2
2 2 2 2
2
dσ 2 σ̄X
j=1 σj . Thus, P (rd dkT2 k2 > δ) ≤ rd d σ σ̄X /(nδ ) → 0 by AsX
√
√
sumption (E1) and (rd d)2 /n → 0. Consider T3 and T4 . Wpg1, rd dkT3 k2 ≤ rd dλs1 κn /n = o(1)
√
√
for s1 = o(n/(rd dλκn )). Similarly, rd dkT4 k2 = oP (1).

The next lemma is needed for proving Theorem 4.3 and its proof is in Supplement F. Suppose
2
{ξ i } are i.i.d. copies of ξ 0 , a d-dimensional random vector with mean zero. Denote σξ,max
=
2
max1≤j≤d Var[ξ0j ], σξ,min
= min1≤j≤d Var[ξ0j ] and γξ,max = max1≤j≤d E|ξ0j |3 .

Lemma A.4. Suppose σξ,max and γξ,max are bounded from above and σξ,max is bounded from zero.
√
√
√ P
If d = o( n), then (1/ n) ni=1 ξ i = OP ( d log d) wrt k·k2 .
27

Proof of Theorem 4.3. We reuse the notations Ti ’s in the proof of Theorems 4.2, from which,
√
√
nAn (β̂ n − β ? ) = V1 + V2 + V3 − V4 , where Vi = B n Ti for i = 1, 2, 3, 4 and B n = nAn T0−1 .
d

It is sufficient to show that V2 −→ N (0, σ 2 GX ) and other Vi ’s are oP (1). Consider V1 . We
√
have kV1 k2 ≤ ndkAn kF kT0−1 kF,d kT1 k2 . By Assumption (F), kAn kF is bounded. By Lemmas
A.2 and A.3 and Assumption (D), for d = o(n1/3 ), wpg1, kT0−1 kF,d is bounded. We have, wpg1,
p
√ √
kT1 k2 ≤ s2 κn γn /n. Then, kV1 k2 . d/ns2 κn γn , Thus, kV1 k2 = oP (1) for s2 = o( n/( dκn γn )).
√
√
nAn (T0−1 − Σ−1
Consider V2 . We have V2 = V21 + V22 , where V21 = nAn Σ−1
X T2 and V22 =
X )T2 .
p
Pn
√
First, note that V21 = (n − s1 )/n i=s1 +1 Z n,i , where Z n,i = (1/ n − s1 )An Σ−1
X X i i . On one
Pn
hand, for every δ > 0, i=s1 +1 EkZ n,i k22 {kZ n,i k2 > δ} ≤ (n − s1 )EkZ n,0 k42 /δ 2 , and EkZ n,0 k42 =
2

−1
T
1
2
E40 E(X T0 Σ−1
X An An ΣX X 0 ) ,
(n−s1 )2

−2
d
4
2
which is ≤ (n−s
2 E0 λmax (Gn )λmin (ΣX )κX . Then, by As1)
P
√
sumptions (D’), (E) and (F) and for d = o( n), ni=s1 +1 EkZ n,i k22 {kZ n,i k2 > δ} → 0. On the other
P
T
2
hand, ni=s1 +1 Cov(Z n,i ) = σ 2 An Σ−1
X An → σ GX by Assumption (F). Thus, by central limit thed

orem (see Proposition 2.27 in van der Vaart (1998)), V21 −→ N (0, σ 2 GX ). Next, consider V22 . Note
−1/2 k√nT k . By Assumption (F), kA k
that kV22 k2 ≤ kAn kF (d log(d))1/2 kT0−1 −Σ−1
2 2
n F
X kF (d log(d))
5
is O(1); by Lemmas A.2 and A.3, (d log(d))1/2 kT0−1 −Σ−1
X kF is oP (1) for d log(d) = o(n); by Lemma
√
√
P
A.4, (d log(d))−1/2 k nT2 k2 = (d log(d))−1/2 k √1n Ss1 +1,n k2 is OP (1) for d = o( n). Then, V22 −→ 0.
d

Thus, by slutsky’s lemma, V2 −→ N (0, σ 2 GX ). Consider V3 and V4 . First consider V3 . By noting
√
√
√
√
√
that s1 = o( n/(λ dκn )), wpg1, kV3 k2 ≤ ndkAn kF kT0−1 kF,d kT3 k2 . dλs1 κn / n → 0. Thus,
kV3 k2 = oP (1). In the same way, kV4 k2 = oP (1).
T1
Proof of Theorem 4.4. By the definition of E, we have P (E) = T1 T2 T3 , where T1 = P ( si=1
{|µ?i +
T
T
X Ti (β ? −β̂)+i | > λ}), T2 = P ( si=s1 +1 {|µ?i +X Ti (β ? −β̂)+i | ≤ λ}) and T3 = P ( ni=s+1 {|X Ti (β ? −
β̂) + i | ≤ λ}). We will show that each Ti converges to one. Then, P (E) → 1. Denote C =
{rd kβ̂ − β ? k2 ≤ 1}. Then P (C) → 1 since β̂ is a consistent estimator of β ? wrt rd k·k2 . Consider
S
T1 . We have 1 − T1 ≤ T11 + T12 , where T11 = P ( i∈S ? {|µ?i + X Ti (β ? − β̂) + i | ≤ λ}) and
21
S
T12 = P ( i∈S ? {|µ?i +X Ti (β ? − β̂)+i | ≤ λ}). It is sufficient to show that both T11 and T12 converge
√31
S
to zero. By dκn  λ  µ? , T11 ≤ P ( i∈S ? {i ≤ λ − µ? + kX i k2 · kβ̂ − β ? k2 }, C) + P (C c ), which
21
√
S
?
c
is ≤ P ( i∈S ? {i ≤ λ − µ + dκn }) + P (C ) ≤ s1 P {0 ≤ −γn } + P (C c ) −→ 0. Similarly, T12 → 0.
21
√
T
Thus T1 → 1. Consider T2 and T3 . By αγn ≤ λ and dκn  λ, T2 ≥ P ( si=s1 {−λ−µ?i +(1/rd )κn ≤
√
√
T
i ≤ λ − µ?i − (1/rd )κn }, C), which is ≥ P ( si=s1 {−λ − µ?i + dκn ≤ i ≤ λ − µ?i − dκn }, C) ≥
T
P ( si=s1 {−γn ≤ i ≤ γn }, C) → 1. Then T2 → 1. Similarly, T3 → 1.
Proof of Theorem 4.5. Note that

√

1/2

nAn ΣX (β̃ − β ? ) = R̃1 + R̃2 + V1 + V2 , where R̃1 =
28

√

nAn R1 ,

R̃2 =

√

nAn R2 , R1 = (X TIˆ X Iˆ0 )−1 X TIˆ Y Iˆ0 {Iˆ0 6= I0 }, R2 = −(X TI0 X I0 )−1 X TI0 Y I0 {Iˆ0 6= I0 },
0

0

and Vi ’s are defined in the proof of Theorem 4.3. Since P (kR̃1 k2 = 0) ≥ P {Iˆ0 = I0 } → 1,
we have R̃1 = oP (1). Similarly, R̃2 = oP (1). By the proof of Theorem 4.3, V1 = oP (1) and
d

V2 −→ N (0, σ 2 GX ). Therefore, the desired result follows by Slutsky’s lemma.
Lemma A.5 (Consistency on σ̂). Suppose the assumptions and conditions of Theorem 4.2 hold
√
P
with rd ≥ d. If s2 = o(n/γn2 ), then σ̂ −→ σ.
√
Proof of Lemma A.5. Since the assumptions and conditions of Theorem 4.2 hold with rd ≥ d,
√
the penalized estimators β̂ and β̃ are consistent estimators of β ? wrt dk·k2 by Theorems 4.2 and
F.4 in Supplement F. Let A = {Iˆ0 = I0 }. Then A occurs wpg1 by Theorem 4.4.
Note that σ̂ 2 = T A + σ̂ 2 Ac , where T = (n − s1 )−1 kY I0 − X TI0 β̃k22 . It suffices to show that
Pn
P6
P
?
T
2
−1
T −→ σ 2 . Note that T =
i=s1 +1 [X i (β − β̃)] , T2 = (n −
i=1 Ti , where T1 = (n − s1 )
P
P
P
s1 )−1 ni=s1 +1 2i , T3 = 2(n − s1 )−1 ni=s1 +1 X Ti (β ? − β̃)i , T4 = (n − s1 )−1 si=s1 +1 µ?2
i , T5 = 2(n −
P
P
P
s1 )−1 si=s1 +1 µi X Ti (β ? − β̃) and T6 = 2(n−s1 )−1 si=s1 +1 µ?i i . It is clear that T2 −→ σ 2 . Thus, it
√
is sufficient to show other Ti ’s are oP (1). For every η > 0, wpg1, dkβ ? − β̃k2 ≤ η. By Assumption
√
T 2
1 Pn
2 . η 2 . For every
dkβ ? − β̃k2 )2 ≤ 2η 2 d1 EkX T0 k22 = 2η 2 σ̄X
(E1), wpg1, |T1 | ≤ d1 n−s
kX
k
(
i
2
i=s
+1
1
1
√
T
1 Pn
kX

k
η > 0, wpg1, |T3 | ≤ 2 √1d n−s
dkβ ? − β̃k2 ≤ 4η √1d EkX T0 0 k2 = 4σησ̄X . η. For
i
2
i
i=s
+1
1
1
√
√
1
s
γ
κ
dkβ ? −
s2 = o(n/γn2 ), |T4 | ≤ (n − s1 )−1 s2 γn2 → 0. For s2 = o( dn/(γn κn )), |T5 | ≤ 2 √1d n−s
2
n
n
1
P

1
1
β̃k2 ≤ 2η √1d n−s
s2 γn κn −→ 0. For s2 = o(n/γn ), wpg1, |T6 | ≤ 4 n−s
γn s2 E|0 | → 0.
1
1

Theorem A.6 (Asymptotic Distributions on β̂ and β̃ with ĜX,n ). Under the assumptions and
√ −1/2
d
conditions of Theorem 4.3, if d8 (log(d))2 = o(n), then nĜX,n An (β̂ − β ? ) −→ N (0, σ 2 I q ).
Similarly, under the assumptions and conditions of Theorem 4.5, If d8 (log(d))2 = o(n), then
√ −1/2
d
nĜX,n An (β̃ − β ? ) −→ N (0, σ 2 I q ).
−1/2

Note that a stronger requirement on d is required to handle ĜX,n in Theorem A.6.. Below is
a lemma needed for proving Theorem A.6.
Lemma A.7 (Wihler (2009)). Suppose A and B are m × m symmetric positive-semidefinite
matrices. Then, for p > 1, kA1/p − B 1/p kpF ≤ m(p−1)/2 kA − BkF . Specifically, for p = 2,
kA1/2 − B 1/2 kF ≤ (m1/2 kA − BkF )1/2 .
Proof of Theorem A.6. We only show the result on β̂. since the result on β̃ can be obtained in a sim√ −1/2
ilar way. We reuse the notations Ti ’s in the proof of Theorems 4.2, from which, nĜX,n An (β̂ n −
√ −1/2
√
−1/2
−1/2
β ? ) = M + R, where M = nGX,n An (β̂ n − β ? ) and R = n(ĜX,n − GX,n )An (β̂ n − β ? ). By
29

d

P

Theorem 4.3, M −→ N (0, σ 2 GX ). Then, it is sufficient to show that R −→ 0 wrt k·k2 . We have
√
−1/2
−1/2
R = R1 + R2 + R3 − R4 , where Ri = B n Ti for i = 1, 2, 3, 4 and B n = n(ĜX,n − GX,n )An T0−1 .
We will show each Ri converges to zero in probability, which finishes the proof. Before that,
−1/2

−1/2

−1/2

−1/2

we first establish an inequality for kĜX,n − GX,n kF . By Lemma A.7, kĜX,n − GX,n kF ≤
−1
√
P
1/2 . Note that, by Lemma A.3, kΣ̂ − Σ k −→
( qkĜX,n − G−1
0 for d4 = o(n). Then, by
n
X F
X,n kF )
−1

P

2
Lemma A.2, kĜX,n − GX,n kF ≤ kAn k2F kΣ̂n − Σ−1
X kF . kAn kF kΣ̂n − ΣX kF −→ 0. Thus, by
−1

2
Lemma A.2, kĜX,n − G−1
X,n kF . kĜX,n − GX,n kF . kAn kF kΣ̂n − ΣX kF . Since q is a fixed in−1/2
√
−1/2
teger, it follows kĜX,n − GX,n kF . kAn kF ( qkΣ̂n − ΣX kF )1/2 . kAn kF (kΣ̂n − ΣX kF )1/2 .
√ √
−1/2
−1/2
Consider R1 . Note that kR1 k2 ≤ n dkĜX,n − GX,n kF kAn kF kT0−1 kF,d kT1 k2 , which is .
√
n(dkΣ̂n − ΣX kF )1/2 kAn k2F kT0−1 kF,d kT1 k2 . By Lemmas A.2 and A.3, dkΣ̂n − ΣX kF = oP (1)

for d6 = o(n). By Assumption (F), kAn kF is bounded. By Lemmas A.2 and A.3 and Assumption (D), for d = o(n1/3 ), wpg1, kT0−1 kF,d is bounded. Also note that, wpg1, kT1 k2 ≤
√
√
s2 κn γn /n. Then, kR1 k2 . s2 κn γn / n. Thus, kR1 k2 = oP (1) for s2 = o( n/(κn γn )). Consider
√
−1/2
−1/2
R2 . Note that kR2 k2 ≤ kĜX,n − GX,n kF kAn kF kT0−1 kF k nT2 k2 , which is . (d2 log(d)kΣ̂n −
√
ΣX kF )1/2 kAn k2F kT0−1 kF,d (d log(d))−1/2 k nT2 k2 . By Lemmas A.2 and A.3, d2 log(d)kΣ̂n −ΣX kF is
oP (1) for d8 (log(d))2 = o(n). By Assumption (F), kAn kF is O(1). By Lemmas A.2 and A.3 and As√
sumption (D), for d = o(n1/3 ), wpg1, kT0−1 kF,d is bounded. By Lemma A.4, (d log(d))−1/2 k nT2 k2 =
√
P
(d log(d))−1/2 k √1n Ss1 +1,n k2 is OP (1) for d = o( n). Thus, R2 −→ 0. Consider R3 and R4 .
√
√
−1/2
−1/2
By s1 = o( n/(λκn )), wpg1, kR3 k2 ≤ nkĜX,n − GX,n kF kAn kF kT0−1 kF kT3 k2 , which is .
√
√
n(dkΣ̂n − ΣX kF )1/2 kAn k2F kT0−1 kF,d kT3 k2 . λs1 κn / n → 0. Thus, kR3 k2 = oP (1). In the same
way, kR4 k2 = oP (1).

B

Supplementary Materials

Additional materials for Sections 1 to 4 can be found in the file of supplementary materials.

References
Debabrata Basu. On the elimination of nuisance parameters. Journal of the American Statistical
Association, 72:355–366, 1977.
Derek Bean, Peter Bickel, Noureddine El Karoui, Chinghway Lim, and Bin Yu. Penalized robust
regression in high-dimension. 2012.

30

Louis H.Y. Chen, Larry Goldstein, and Qi-Man Shao. Normal Approximation by Stein’s Method.
Springer-Verlag, 2010.
Jianqing Fan and Runze Li. Variable selection via nonconcave penalized likelihood and its oracle
properties. Journal of the American Statistical Association, 96(456):1348–1360, Dec 2001. ISSN
0162-1459.
Jianqing Fan and Jinchi Lv. A selective overview of variable selection in high dimensional feature
space. Statistica Sinica, 20:101–148, 2010.
Jianqing Fan and Heng Peng. On non-concave penalized likelihood with diverging number of
parameters. The Annals of Statistics, 32:928–961, 2004.
Jianqing Fan, Yuan Liao, and Martina Mincheva. High dimensional covariance matrix estimation
in approximate factor models. 2011.
Jianqing Fan, Yingying Fan, and Emre Barut. Adaptive robust variable selection. 2012a.
Jianqing Fan, Yang Feng, and Xin Tong. A road to classification in high dimensional space: the
regularized optimal affine discriminant. Journal of the Royal Statistical Society Series B, 2012b.
Peter J. Huber. Robust estimation of a location parameter. The Annals of Mathematical Statistics,
35:73–101, 1964.
Peter J. Huber. Robust regression: Asymptotics, conjectures and monte carlo. The Annals of
Statistics, 1:799–821, 1973.
Johannes Jahn. Introduction to the Theory of Nonlinear Optimization. Springer Berlin Heidelberg,
2007.
J. Kiefer and J. Wolfowitz. Consistency of the maximum likelihood estimator in the presence of
infinitely many incidental parameters. The Annals of Mathematical Statistics, 27:887–906, 1956.
Michael R. Kosorok. Bootstrapping the grenander estimator. In Beyond Parametrics in Interdisciplinary Research: Festschrift in Honor of Professor Pranab K. Sen, pages 282–292. Institute
of Mathematical Statistics: Hayward, CA., 2008.
S. Lambert-Lacroix and L. Zwald. Robust regression through the huber’s criterion and adaptive
lasso penalty. Electronic Journal of Statistics, 5:1015–1053, 2011. ISSN 1935-7524.
31

Tony Lancaster. The incidental parameter problem since 1948. Journal of Econometrics, 95:
391–413, 2000.
Marcelo J. Moreira. A maximum likelihood method for the incidental parameter problem. The
Annals of Statistics, 37(6A):3660–3696, 2009. ISSN 0090-5364.
Jerzy Neyman and Elizabeth L. Scott. Consistent estimates based on partially consistent observations. Econometrica, 16:1–32, 1948.
Stephen Portnoy and Xuming He. A robust journey in the new millennium. Journal of the American
Statistical Association, 95:1331–1335, 2000.
Albert N. Shiryaev. Probability. Springer-Verlag, second edition, 1995.
G. W. Stewart. On the continuity of the generalized inverse. SIAM Journal on Applied Mathematics,
17:33–45, 1969.
Runlong Tang, Moulinath Banerjee, and Michael R. Kosorok. Likelihood based inference for current
status data on a grid: A boundary phenomenon and an adaptive inference procedure. The Annals
of Statistics, 40(1):45–72, 2012.
Robert Tibshirani. Regression shrinkage and selection via the lasso. J. Royal. Statist. Soc B, 58:
267–288, 1996.
Aad W. van der Vaart. Asymptotic Statistics. Cambridge University Press, 1998.
Aad W. van der Vaart and Jon A. Wellner. Weak Convergence and Empirical Processes. Springer,
1996.
Thomas P. Wihler. On the holder continuity of matrix functions for normal matrices. Journal of
inequalities in pure and applied mathematics, 10, 2009.
Peng Zhao and Bin Yu. On model selection consistency of lasso. The Journal of Machine Learning
Research, 7(Nov):2541–2563, 2006.

32

Supplementary Materials for Paper:
“Partial Consistency with Sparse Incidental Parameters”
by Jianqing Fan, Runlong Tang and Xiaofeng Shi
C

Supplement for Section 1

In this supplement, we first show the method proposed by Neyman and Scott (1948) does not work
for model (1.1) and then explain which assumptions or conditions for the consistent results of the
penalized methods in Zhao and Yu (2006), Fan and Peng (2004) and Fan and Lv (2011) are not
satisfied for model (1.1).
Although the modified equations of maximum likelihood method proposed by Neyman and Scott
(1948) could handle “a number of important cases” with incidental parameters, unfortunately, it
does not work for model (1.1). More specifically, consider the simplest case of model (1.1) with
d = 1:
Yi = µ?i + Xi β ? + i , for i = 1, 2, · · · , n,
where {i } are i.i.d. copies of N (0, σ 2 ). Using the notations of Neyman and Scott (1948), the
√
likelihood function for (Xi , Yi ) is pi = pi (β, σ, µi |Xi , Yi ) = ( 2πσ)−1 exp{−(2σ 2 )−1 (Yi −µ?i −Xi β)2 },
√
and the log-likelihood function is log pi = − log( 2πσ) − (2σ 2 )−1 (Yi − µ?i − Xi β)2 . Then, the score
functions are
1
∂ log pi
= 2 (Yi − µ?i − Xi β)Xi ,
∂β
σ
∂ log pi
1
1
φi2 =
= + 3 (Yi − µ?i − Xi β)2 ,
∂σ
σ σ
∂ log pi
1
ωi =
= 2 (Yi − µ?i − Xi β).
∂µi
σ
φi1 =

From the equation ωi = 0, we have µ̂i = Yi − Xi β. Plugging this µ̂i into φi1 and φi2 (replacing
µi with µ̂i ), we obtain φi1 = 0 and φi2 = 1/σ. Then, Ei1 = Eφi1 = 0 and Ei2 = Eφi1 = 1/σ.
Thus, Ei1 and Ei2 do only depend on the structural parameters (β ? and σ). However, we then
have Φi1 = φi1 − Ei1 = 0 and Φi2 = φi2 − Ei2 = 0. This means Fn1 = Fn2 = 0, independent of
structural parameters! Consequently, the estimation equations degenerate to two 0 = 0 equations,
which means the modified equation of maximum likelihood method does not work for model (1.1).

i

Next, we explicitly explain which assumptions or conditions for the consistent results of the
penalized methods in Zhao and Yu (2006), Fan and Peng (2004) and Fan and Lv (2011) are not
valid for model (1.1).
Zhao and Yu (2006) derive strong sign consistency for lasso estimator. However, their consistency results Theorems 3 and 4 do not apply to model (1.1), since the above specific design matrix
X does not satisfy their regularity condition (6) on page 2546. More specifically, with model (1.1),




0
0
I
X
1
s
1,s
a.s. 
n
 −→
,
C11
= 
n X T Pn X X T
0 Σ
1,s

i

i=1

X

i

where ΣX is the covariance matrix of the covariates. This means that some of the eigenvalues of
n goes to 0 as n → ∞. Then the regularity condition (6), which is
C11
n
αT C11
α ≥ a positive constant , for all α ∈ Rs+d such that kαk22 = 1,

does not hold any more. Thus the consistency results Theorems 3 and 4 in Zhao and Yu (2006) is
not applicable for model (1.1).
Fan and Peng (2004) show the consistency with Euclidean metric of a penalized likelihood
estimator when the dimension of the sparse parameter increases with the sample size in Theorem
1 on Page 935. Under their framework, the log-likelihood function of the data point Vi = (X i , Yi )
for each i from model (1.1) with random errors being i.i.d. copies of N (0, σ 2 ) is given by
log fn (Vi , µi , β) ∝ −

1
(Yi − µi − X Ti β)2 ,
2σ 2

where ∝ means “proportional to”. As we can see that log-likelihood functions with different i’s
might different since µi ’s might be different for different i’s. This violates a condition that all the
data points are i.i.d. from a structural density in Assumption (G) on Page 934.
This violation might not be essential, however, since we could consider the log-likelihood function for all the data directly. That is, we consider
Ln (µ, β) =

n
X

log fn (Vi , µi , β) ∝ −

n
1 X
(Yi − µi − X Ti β)2 .
2σ 2
i=1

i=1

Then, the Fisher information matrix for (µ, β) is given by


σ −2 I n
0
,
In+d (µ, β) = 
0
nσ −2 Σ2X

ii

where In is the n × n identity matrix. Then, the Fisher information for one data point is


−1 σ −2 I
n
0
1
n
.
In+d (µ, β) = 
n
0
σ −2 Σ2
X

It is clear that the minimal eigenvalue λmin (In+d (µ, β)/n) = n−1 σ 2 → 0 as n → ∞. This violates
the condition that the minimal eigenvalue should be lower bounded from 0 in Assumption (F) on
Page 934. Thus, the consistency result Theorem 1 in Fan and Peng (2004) can not be applied to
model (1.1).
Fan and Lv (2011) “consider the variable selection problem of nonpolynomial dimensionality
in the context of generalized linear models” by taking the penalized likelihood approach with
folded-concave penalties. Theorem 3 on page 5472 of Fan and Lv (2011) shows that there exists a
consistent estimator of the unknown parameters with the Euclidean metric under certain conditions.
In Condition 4 on page 5472, there is a condition on a minimal eigenvalue
min λmin [X TI Σ(XI δ)X I ] ≥ cn,

δ∈N0

where X I consists of the first s + d columns of the design matrix X. With model (1.1), this
condition becomes
λmin [X TI X I ] ≥ cn,
which is
n
λmin [(1/n)X TI X I ] = λmin [C11
] ≥ c,
n is the matrix defined in Zhao and Yu (2006) and c is a positive constant. Since the minwhere C11
n ] converges to 0, the above condition does not hold. Thus, the consistency
imal eigenvalue λmin [C11

result Theorem 3 of Fan and Lv (2011) is not applicable for model (1.1).

D

Supplement for Section 2

In this supplement, we provide Lemmas D.1 and D.3, Proposition D.2 and their proofs. Before
that, there are two graphs, Figures 5 and 6, illustrating the incidental parameters and the step of
updating the responses in the iteration algorithm with d = 1.

iii

n−s

γn

λ

µ∗

{µ i∗ : s1 + 1 ≤ i ≤ s}
s2

0

∞

{µi∗ : 1 ≤ i ≤ s1}
s1

−∞

− λ − µ∗

−γn

Figure 5: An illustration of three types of µ?i ’s, that is, large µ?1 , bounded µ?2 and zero µ?3 . The
negative half of the real line is folded at 0 under the positive half for convenience. For the penalized least square method with a soft penalty function and under the assumption of fixed d, the
√
specification of the regularization parameter λ is that κn  λ, αγn ≤ λ, and λ  min{µ? , n}.

Y

0
X

f ( x ) = xβ ( k −1)

Figure 6: An illustration for the updating of responses with d = 1. The solid black line is a fitted
regression line. The dashed black lines are the corresponding shifted regression lines. The circle
and diamond points are the original data points. The circle and triangle points are the updated
data points. That is, the diamond points are drawn onto the shifted regression lines.

iv

Lemma D.1. A necessary and sufficient condition for (µ̂, β̂) to be a minimizer of L(µ, β) is that
β̂ = (X T X)−1 X T (Y − µ̂),
for i ∈ Iˆ0c ,

Yi − µ̂i − X Ti β̂ = λSign(µ̂i ),
|Yi − X Ti β̂| ≤ λ,

for i ∈ Iˆ0 ,

where Sign(·) is a sign function and Iˆ0 = {1 ≤ i ≤ n : µ̂i = 0}.
Proof of Lemma D.1. By subdifferential calculus (see, for example, Theorem 3.27 in Jahn (2007)),
a necessary and sufficient condition for (µ̂, β̂) to be a minimizer of L(µ, β) is that zero is in the
subdifferential of L at (µ̂, β̂), which means that, for each i,
β̂ = (X T X)−1 X T (Y − µ̂),
Yi − µ̂i − X Ti β̂ = λSign(µ̂i ),
|Yi − X Ti β̂| ≤ λ,

if µ̂i 6= 0,

if µ̂i = 0.

Thus, the conclusion of Lemma D.1 follows.
Proposition D.2. Suppose Assumptions (A) and (B) hold and there exist positive constants C1
and C2 such that kβ ? k2 < C1 and kβ (0) k2 < C2 wpg1. If s1 λ/n = O(1) and s2 γn /n = o(1), then,
for every K ≥ 1 and k ≤ K, wpg1 as n → ∞,
kβ (K+1) − β (K) k2 ≤ O((s1 /n)K ),

and

√
kβ (k) k2 ≤ 2 dC1 + C2 .

Remark 2. For any prespecified critical value in the stopping rule, Proposition D.2 implies that
the algorithm stops at the second iteration wpg1. In practice, the sample size n might not be
large enough for the two-iteration estimator to have a decent performance so that more iterations
are usually needed to activate the stopping rule. By Proposition D.2, K iterations will make the
distance kβ (K+1) −β (K) k2 of the small order (s1 /n)K . When s1 /n is small, the algorithm converges
quickly, which has been verified by our simulations.
√
Proof of Proposition D.2. First, we show that, wpg1, kβ (1) k2 is bounded by 2 dC1 + C2 . For each
k ≥ 1, we have
Sβ (k) = SµS11 + SµS12 + SS1 β ? + SS1 + SS2 ∪S3 β (k−1) + λ(SS2 − SS3 ),

v

where Si = ∪3j=1 Sij (β (k−1) ) for i = 1, 2, 3 and Sij ’s are defined at the end of Section 2. Denote
Ak−1 as the event
?
?
?
?
?
{S11 (β (k−1) ) = ∅, S12 (β (k−1) ) = S12
, S1 (β (k−1) ) = S10
∪ S12
; S2 (β (k−1) ) = S21
; S3 (β (k−1) ) = S31
},
? ’s are defined at the beginning of Section 3.
where Sij

By Lemma 3.1, P (A0 ) → 1. Thus, wpg1,
β (1) = T0−1 T1 + T0−1 T2 + T0−1 T3 + T0−1 T4 (β (0) ) + T0−1 T5 ,
where T0 = S/n, T1 = SµS ? /n, T2 = Ss1 +1,n β ? /n, T3 = Ss1 +1,n /n, T4 (β (0) ) = S1,s1 β (0) /n and T5 =
12
√
−1
−1
−1
? − SS ? )λ/n. We will show that, wpg1, kT
(SS21
0 T1 k2 ≤ C2 /4, kT0 T2 k2 ≤ 2 dC1 , kT0 T3 k2 ≤
31
C2 /4, kT0−1 T4 (β (0) )k2 ≤ C2 /4 and kT0−1 T5 k2 ≤ C2 /4. Then, wpg1,
kβ (1) k2 ≤

5
X

√
kT0−1 Ti k2 ≤ 2 dC1 + C2 .

i=1

On T0−1 T1 . For s2 γn /n = o(1), wpg1,
1
s2
1
kT0−1 T1 k2 ≤ k( S)−1 kF k SµS ? k2 ≤ 4kΣ−1
γn → 0.
X kF EkX 0 k2
12
n
n
n
Thus, wpg1, kT0−1 T1 k2 ≤ C2 /4.
On T0−1 T2 . Wpg1,
√
1
1
kT0−1 T2 k2 ≤ k( S)−1 Ss1 +1,n kF kβ ? k2 ≤ 2kI d kF C1 = 2 dC1 .
n
n
On T0−1 T3 . Wpg1,
1 
P
kT0−1 T3 k2 ≤ 2kΣ−1
X kF k Ss1 +1,n k2 −→ 0.
n
Thus, wpg1, kT0−1 T3 k2 ≤ C2 /4.
On T0−1 T4 (β (0) ). For s1 /n = o(1),
kT0−1 T4 (β (0) )k2 ≤

s1 √
s1 1 −1 1
P
k( S)
S1,s1 kF kβ (0) k2 ≤ 2 dC2 −→ 0.
n n
s1
n

Thus, wpg1, kT0−1 T4 (β (0) )k2 ≤ C2 /4.
On T0−1 T5 . For s1 λ/n = O(1), wpg1,
kT0−1 T5 k2 ≤ 2kΣ−1
X kF

s1 λ 1
1
P
? k2 + k
? k2 ) −→ 0.
(k SS21
SS31
n
s1
s1

Thus, wpg1, kT0−1 T5 k2 ≤ C2 /4.
vi

Next, consider kβ 2 − β 1 k2 . Since β (1) is bounded wpg1, by Lemma 3.1, A1 occurs wpg1. Then,
β (2) = T0−1 T1 + T0−1 T2 + T0−1 T3 + T0−1 T4 (β (1) ) + T0−1 T5 ,
where T4 (β (1) ) = (1/n)S1,s1 β (1) . Thus, wpg1,
β (2) − β (1) = S−1 S1,s1 (β (1) − β (0) ).
It follows that, for s1 = o(n), wpg1,
√
√
kβ (2) − β (1) k2 ≤ kS−1 S1,s1 kF kβ (1) − β (0) k2 ≤ (2 ds1 /n)(4 dC1 + 2C2 ) → 0.
Then, wpg1, β (2) = β (1) , which means that, wpg1, the iteration algorithm stops at the second
iteration.
Finally, for any K ≥ 1, repeat the above arguments. Then, with at least probability pn,K =
TK
P ( k=0 Ak ), which increases to one by Lemma 3.1, we have
√
√
kβ (K+1) − β (K) k2 ≤ (2 ds1 /n)K (4 dC1 + 2C2 ) = O((s1 /n)K ) → 0,
√
and kβ (k) k2 ≤ 2 dC1 + C2 for all k ≤ K.
Lemma D.3. A necessary and sufficient condition for (µ̂, β̂) to be a minimizer of L(µ, β) is that
it is a solution to equations (2.5) and (2.6).
Proof of Lemma D.3. First, we show a solution of (2.5) and (2.6) satisfies the necessary and sufficient condition in Lemma D.1. Denote a solution of (2.5) and (2.6) as (µ̂, β̂). Then β̂ =
(X T X)−1 X T (Y − µ̂), which is exactly the first condition in Lemma D.1, and, for each i =
1, 2, · · · , n, (µ̂, β̂) satisfies one of three cases: |Yi − X Ti β̂| ≤ λ and µ̂i = 0; Yi − X Ti β̂ > λ and
µ̂i = Yi − X Ti β̂ − λ; Yi − X Ti β̂ < −λ and µ̂i = Yi − X Ti β̂ + λ. If (µ̂, β̂) satisfies the first case,
it satisfies the third condition in Lemma D.1. If (µ̂, β̂) satisfies the second case, then µ̂i > 0 and
Yi − µ̂i − X Ti β̂ = λ = λSign(µ̂i ), which means that the second case satisfies the second condition
in Lemma D.1. Similarly, the third case also satisfies the second condition in Lemma D.1. Thus
(µ̂, β̂) satisfies the necessary and sufficient condition in Lemma D.1.
In the other direction, suppose (µ̂, β̂) satisfies the necessary and sufficient condition in Lemma
D.1. Then, the first condition in Lemma D.1 exactly (2.5). For each i, (µ̂, β̂) satisfies one of three
cases: µ̂i = 0 and |Yi − X Ti β̂| ≤ λ; µ̂i > 0 and Yi − µ̂i − X Ti β̂ = λ; µ̂i < 0 and Yi − µ̂i − X Ti β̂ = −λ.
If (µ̂, β̂) satisfies the first case, it satisfies the first case in (2.6). If (µ̂, β̂) satisfies the second case,
vii

then µ̂i = Yi − X Ti β̂ − λ and Yi − X Ti β̂ > λ, which means that (µ̂, β̂) satisfies the second case
of (2.6). Similarly, If (µ̂, β̂) satisfies the third case, then it satisfies the third case of (2.6). Thus,
(µ̂, β̂) satisfies (2.5) and (2.6).

E

Supplement for Section 3

In this supplement, we provide the proofs of the theoretical results in Section 3. Before that, we
point out that those two different sufficient conditions in Theorem 3.2 come from the different
analysis on the term SµS ? . Each of the two different sufficient conditions does not imply the other.
12

Specifically, on one hand, suppose the absolute values of µ?i ’s are all equal for i = s1 +1, s2 +2, · · · , s.
P
(2+δ)/2 ? 2+δ
Then, kµ?2 k2+δ
= s2
|µs |
and si=s1 +1 |µ?i |2+δ = s2 |µ?s |2+δ . Thus Assumption (C) holds
2
automatically since s2 → ∞. This means that Assumption (C) holds at least when the absolute
magnitudes of µ?i ’s are similar to each other. For this case, there still exists a consistent estimator
even if n/(κn γn )  s2  n. On the other hand, suppose µ?s = γn and the other µ?i ’s are all equal to
P
= [γn2 +(s2 −1)c2 ](2+δ)/2 and si=s1 +1 |µ?i |2+δ = γn2+δ +(s2 −1)c2+δ .
a constant c > 0. Then, kµ?2 k2+δ
2
If s2  γn2  n/(κn γn ), the previous two terms are both asymptotically equivalent to γn2+δ . Thus
Assumption (C) fails but the other sufficient condition holds.
Proof of Lemma 3.1. The proof is the similar to that of Lemma 4.1 and omitted.
Proof of Theorems 3.2. By Lemma 3.1, wpg1, the solution β̂ n to ϕn (β) = 0 on BC (β ? ) is explicitly
given by
β̂ n = β ? + T0−1 (T1 + T2 + T3 − T4 ),
? and T4 = (λ/n)SS ? .
where T0 = (1/n)Ss1 +1,n , T1 = (1/n)SµS ? , T2 = (1/n)Ss1 +1,n , T3 = (λ/n)SS21
31
12

P

P

We will show that T0 −→ Σ−1
X > 0 with the Frobenius norm and Ti −→ 0 with the Euclidean norm
for i = 1, 2, 3, 4. Thus, by Slutsky’s lemma (see, for example, Lemma 2.8 on page 11 of van der
Vaart (1998)), β̂ n is a consistent estimator of β ? .
P

On T0−1 . By law of large number, T0 −→ ΣX > 0. Then, by continuous mapping theorem,
P

T0−1 −→ Σ−1
X > 0.
On T1 : Approach One. Suppose s2 = o(n/(κn γn )). Then,
kT1 k2 ≤

s
s
1 X
1 X
kX i µ?i k2 =
kX i k2 · |µ?i | ≤ s2 κn γn /n = o(1).
n
n
i=s1 +1

i=s1 +1

viii

On T1 : Approach Two. Under Assumption (C), it follows (ΣX

?2 −1/2 Sµ
?
i=s1 +1 µi )
S12

Ps

d

−→

N (0, Id ). In fact, Assumption (C) implies the Lyapunov condition for sequence of random vectors
(see, e.g. Proposition 2.27 on page 332 of (van der Vaart, 1998)). More specifically, recall the
Lyapunov condition is that there exists some constant δ > 0 such that
s
X

s
X

Ek(ΣX

i=s1 +1

−1/2
µ?2
X i µ?i k2+δ
→ 0.
j )
2

j=s1 +1

Then, by Assumption (C),
s
X

Ek(ΣX

i=s1 +1

s
X

1

−2
µ?2
X i µ?i k2+δ
≤(
j )
2

j=s1 +1

s
X

−
µ?2
j )

2+δ
2

j=s1 +1

s
X

− 2+δ

−→ 0,
|µ?i |2+δ λmin2 EkX0 k2+δ
2

i=s1 +1

where λmin > 0 is the minimum eigenvalue of ΣX . Then,
s
s
X
X
1 µ
1
?2 1/2
−1/2 µ
kT1 k2 = k SS ? k2 ≤ k(ΣX
µi ) kF k(ΣX
µ?2
SS ? k2
i )
12
n 12
n
i=s1 +1

=

1
(
n

s
X

i=s1 +1

1/2

1/2
µ?2
kΣX kF OP (1) ≤
i )

i=s1 +1

1
1
(s2 γn2 )1/2 OP (1) ≤ √ γn OP (1) = oP (1).
n
n

On T2 . By law of large number, T2 = oP (1).
√
On T3 and T4 . By noting λ  n,
√
s1 1
1
λ
?
? k2 ≤ √ OP (1) = oP (1).
kT3 k2 = kλ SS21 k2 = λ
k √ SS21
n
n
s1
n
Thus T3 = oP (1). In the same way, we can show that T4 = oP (1) holds.
In Theorem, one condition is Dn /n = o(1). It turns out we can consider other conditions on
Dn and derive more possible asymptotic distributions for β̂ n .
Theorem E.1 (Asymptotic Distributions on β̂ n : more cases). Under Assumptions (A), (B) and
(C), for all constants b, c ∈ R+ ,
√
d
(1) when s1  n/λ2 and Dn2 /n = o(1), n(β̂ n − β ? ) −→ N (0, σ 2 Σ−1
X ); [main case]
√
d
(2) when s1  n/λ2 and Dn2 /n ∼ c, n(β̂ n − β ? ) −→ N (0, (c + σ 2 )Σ−1
X );
d

√

d

√

(3) when s1  n/λ2 and Dn2 /n → ∞, rn (β̂ n − β ? ) −→ N (0, Σ−1
X ), where rn ∼ n/Dn 
√
d
(4) when s1 ∼ bn/λ2 and Dn2 /n = o(1), n(β̂ n − β ? ) −→ N (0, (b + σ 2 )Σ−1
X );
√
d
(5) when s1 ∼ bn/λ2 and Dn2 /n ∼ c, n(β̂ n − β ? ) −→ N (0, (b + c + σ 2 )Σ−1
X );
(6) when s1 ∼ bn/λ2 and Dn2 /n → ∞, rn (β̂ n − β ? ) −→ N (0, Σ−1
X ), where rn ∼ n/Dn 

ix

n;

n;

d

(7) when s1  n/λ2 and Dn2 /n = o(1) or Dn2 /n ∼ c, rn (β̂ n − β ? ) −→ N (0, Σ−1
X ), where rn ∼
√
√
n/(λ s1 )  n;
√
√
√
(8) when s1  n/λ2 and Dn2 /n → ∞, letting rn ∼ min{ bn/(λ s1 ), n/Dn }  n,
√

√
d
bn/(λ s1 )  n/Dn , then rn (β̂ n − β ? ) −→ N (0, Σ−1
X );
√
√
d
(8b) if bn/(λ s1 ) ∼ n/Dn , then rn (β̂ n − β ? ) −→ N (0, (1 + b)Σ−1
X );
√
√
d
(8c) if bn/(λ s1 )  n/Dn , then rn (β̂ n − β ? ) −→ N (0, bΣ−1
X ).

(8a) if

Theorem 3.3 groups the results according to the asymptotic magnitude of s1 given an upper
bound of the diverging speed of s2 . Alternatively, Theorem E.1 groups the results according to
the asymptotic magnitudes of s1 and Dn2 . Since both s1 and Dn2 have three cases, Theorem E.1
basically contains nine cases. For the last case, there are further three cases on the relationship
√
√
between bn/(λ s1 ) and n/Dn . As in Theorem 3.3, the first case of Theorem E.1 is denoted as
the main case since for this case the incidental parameters are sparse in the sense that the size
and magnitude of the nonzero incidental parameters µ?1 and µ?2 are well controlled. Note that
√
s2 = o( n/(κn γn )) implies Dn2 /n = o(1). which means that, under Assumption (C), the cases (1),
(4) and (7) of Theorem E.1 actually imply the three results of Theorem 3.3. As in Theorem 3.3,
√
the convergence rate of β̂ n becomes less than n when s1  n/λ2 or Dn2 /n → ∞, that is, when
the size and magnitude of the nonzero incidental parameters are large; the boundary phenomenon
also appears.
Proof of Theorems 3.3 and E.1. It is sufficient to provide the proof for the case where the sizes of
? = {1 ≤ i ≤ s : µ? > 0} and S ? = {1 ≤ i ≤ s : µ? < 0} are both asymptotically
index sets S21
1
1
31
i
i

s1 /2 and b = 2.
From the proof of Theorems 3.2, wpg1,
µ

?
?
β̂ = β ? + S−1
s1 +1,n [SS ? + Ss1 +1,n + λ(SS21 − SS31 )].
12

Let rn be a sequence going to infinity. Then, rn (β̂ n −β ? ) = T0−1 (V1 +V2 +V3 −V4 ), where V1 = rn T1 ,
V2 = rn T2 , V3 = rn T3 , V4 = rn T4 and Ti ’s are defined in the proof of Theorem 3.2. Next we derive
the asymptotic properties of T0 and Vi ’s, from which the desired results follow by Slutsky’s lemma.
P

On T0 . By the proof of Theorem 3.2, T0−1 −→ Σ−1
X
√
√
On V1 : Approach One. If rn = n and s2 = o( n/(κn γn )), then
s
1 µ
1 X
1
1
kT1 k2 = krn SS ? k2 ≤ rn
kX i k2 · |µ?i | ≤ rn s2 κn γn = √ s2 κn γn = o(1).
12
n
n
n
n
i=s1 +1

x

√
n and s2 = o( n/(κn γn )), then T1 = oP (1).
√
On V1 : Approach Two. If rn = n, then

Thus, if rn =

√

n or rn 

√

1
Dn 1 µ
Dn 1 µ
T1 = rn SµS ? = rn
SS ? = √
S ?,
n 12
n Dn 12
n Dn S12
P
1/2 . There are three cases on D /√n or D 2 /n. If D 2 /n → 0,
where Dn = kµ?2 k2 = ( si=s1 +1 µ?2
n
n
n
i )
√
d
P
2
2
then T1 −→ 0. If Dn /n → 1, then T1 −→ N (0, ΣX ). If Dn /n → ∞, it means that rn = n is too
√ p
√
d
fast. Let rn ∼ n/Dn = n n/Dn2  n. Then T1 −→ N (0, ΣX );
√
√
√
d
P
On V2 . If rn = n, then T2 −→ N (0, σ 2 ΣX ). Thus, if rn  n, T2 −→ 0; if rn  n;
P

T2 −→ ∞.
On V3 and V4 . First consider T3 . Denote #(·) as the size function. If rn =
r
s1 /2
1
1
? = λ
p
S ?.
T3 = λrn SS21
? ) S21
n
n
#(S21

√
n, then

p
p
P
? ) = s /2. There are three cases on λ s /(2n). If λ s /(2n) → 0, then T −→
Note that #(S21
0.
1
1
1
3
p
p
d
Note that λ s1 /(2n) → 0 is equivalent to s1 = o(2n/λ2 ). If λ s1 /(2n) → 1, then T3 −→
p
p
N (0, ΣX ). Note that λ s1 /(2n) → 1 is equivalent to s1 ∼ 2n/λ2 . If λ s1 /(2n) → ∞, it means
p
√
√ √
√
√
rn = n is too large. Let rn ∼ n/(λ (s1 /2)) = n 2n/(λ s1 )  n. With this rate rn ,
p
d
T3 −→ N (0, ΣX ). Note that λ s1 /2n → ∞ is equivalent to s1  O(2n/λ2 ). In the same way, T4
can be analyzed and parallel results can be obtained.
Proof of Theorem 3.4. The proof is similar to that of Theorem 4.4 and omitted.

E.1

Supplement for Subsection 3.1

The following Theorem implies Theorem 3.5 since it contains more details.
Theorem E.2 (Consistency and Asymptotic Normality on β̃). Suppose Assumptions (A) and (B)
√
P
hold. If either s2 = o(n/(κn γn )) or Assumption (C) holds, then β̃ −→ β ? . If s2 = o( n/(κn γn )),
√
d
then n(β̃ − β ? ) −→ N (0, σ 2 Σ−1
X ). On the other hand, under Assumption (C),
(1) if Dn2 /n = o(1), then
(2) if Dn2 /n ∼ c, then

√

d

n(β̃ − β ? ) −→ N (0, σ 2 Σ−1
X ); [main case]

√
d
+
n(β̃ − β ? ) −→ N (0, (c + σ 2 )Σ−1
X ), for every constant c ∈ R ;
d

(3) if Dn2 /n → ∞, then rn (β̃ − β ? ) −→ N (0, Σ−1
X ) where rn ∼ n/Dn 

xi

√

n.

Proof of Theorem E.2. Denote I0 = {s1 + 1, s1 + 2, · · · , s = s1 + s2 , s + 1, · · · , n}. Note that
√
s2 = o( n/(κn γn )) ensures that β̂ is consistent by Theorem 3.2. By Theorem 3.4, P {Iˆ0 = I0 } goes
to 1. Then,
β̃ = R1 + R2 + T0−1 (T1 + T2 ),
where R1 = (X TIˆ X Iˆ0 )−1 X TIˆ Y Iˆ0 {Iˆ0 6= I0 } and R2 = −(X TI0 X I0 )−1 X TI0 Y I0 {Iˆ0 6= I0 } and Ti ’s are
0

0

defined in the proof of Theorem 3.2. The proof for the consistency is similar to that of Theorem
3.2 and is omitted. Next we show the asymptotic normality. We have,
rn (β̃ − β ? ) = rn R1 + rn R2 + T0−1 (V1 + V2 ),
√
where Vi ’s are defined in the proof of Theorem E.1. Since P ( nR1 = 0) ≥ P {Iˆ0 = I0 } → 1, we
√
√
have nR1 = oP (1). Similarly, nR2 = oP (1). From the analysis on Vi ’s in the proof of Theorem
E.1, the asymptotic distributions follows by Slutsky’s lemma.
Lemma E.3 (Consistency on σ̂). Suppose Assumptions (A) and (B) hold and either s2 = o(n/(κn γn ))
P

or Assumption (C) holds. If s2 = o(n/γn2 ), then σ̂ −→ σ.
Proof of Lemma E.3. When Assumption (C) or s2 = o(n/(κn γn )) holds, the penalized estimators β̂
and β̃ are consistent estimators of β ? by Theorems 3.2 and 3.5. Denote C = {Iˆ0 = I0 }. By Theorem
3.4, C occurs wpg1. Then, σ̂ 2 = T C + σ̂ 2 C c , where T = an kY I0 − X TI0 β̃k22 and an = 1/(n − s1 ). It
P
P
P
is sufficient to show T −→ σ 2 . We have T = 6i=1 Ti , where T1 = an ni=s1 +1 [X Ti (β ? − β̃)]2 , T2 =
Ps
P
P
P
?
T
an ni=s1 +1 2i , T3 = 2an ni=s1 +1 X Ti (β ? − β̃)i , T4 = an si=s1 +1 µ?2
i , T5 = 2an
i=s1 +1 µi X i (β −
P
P
P
β̃) and T6 = 2an si=s1 +1 µ?i i . It is straightforward to show that T2 −→ σ 2 and each other Ti −→ 0
P

under the condition s2 = o(n/γn2 ) and by noting that β̃ −→ β ? . Then σ̂ is a consistent estimator
of σ.

E.2

Supplement for Subsection 3.2

In this supplement, we consider a special case with exponentially tailed covariates and errors. For
convenience, we first introduce the definition of Orlicz norm and related inequalities. For a strictly
increasing and convex function ψ with ψ(0) = 0, the Orlicz norm of a random variable Z with
respect to ψ is defined as
kZkψ = inf{C > 0 : Eψ(|Z|/C) ≤ 1}.
Then, for each x > 0,
P (|Z| > x) ≤ 1/ψ(x/kZkψ ).
xii

(E.1)

(See Page 96 of van der Vaart and Wellner (1996)). Next, we introduce a lemma on Orlicz norm with
ψ1 . Suppose {Zi }ni=1 is a sequence of random variables and {Z i }ni=1 is a sequence of d-dimensional
random vectors with Z i = (Zi1 , Zi2 , · · · , Zid )T . From Lemma 8.3 on Page 131 of Kosorok (2008),
we have the following extension.
Lemma E.4. If for each 1 ≤ i ≤ n and 1 ≤ j ≤ d,
x2
1
x2
1
} and P (|Zij | > x) ≤ c exp{− ·
},
P (|Zi | > x) ≤ c exp{− ·
2 ax + b
2 ax + b
with a, b ≥ 0 and c > 0, then
k max |Zi k|ψ1 ≤ K{a(1 + c) log(1 + n) +
1≤i≤n

p
p
b(1 + c) log(1 + n)},

p
p
√
k max kZ i k2 kψ1 ≤ K{a d(1 + cd) log(1 + n) + bd(1 + cd) log(1 + n)}.
1≤i≤n

where K is a universal constant which is independent of a, b, c, {Zi } and {Z i }.
Proof of Lemma E.4. The proof for random variables {Zi } is the same to the proof of Lemma 8.3
on Page 131 of Kosorok (2008). For random vectors {Z i },
d
X
√
√
1 x2
P (kZ i k2 ≥ x) ≤ P ( max |Zij | > x/ d) ≤
P (|Zij | > x/ d) ≤ c0 exp{− 0
},
1≤j≤d
2 a x + b0
j=1

√

where a0 = a d, b0 = bd and c0 = cd. Then, by the result on random variables, the desired result
on random vectors follows.
Now, suppose, for every x > 0,
1
x2
1
x2
P (|i | > x) ≤ c1 exp{− ·
} and P (|Xij | > x) ≤ c2 exp{− ·
},
2 a1 x + b1
2 a2 x + b2

(E.2)

with ai , bi ≥ 0 and ci > 0 for i = 1, 2. By Lemma E.4, it follows
k max |i k|ψ1 ≤ K{a1 (1 + c2 ) log(1 + n) +
1≤i≤n

p

b1 (1 + c1 )

p
log(1 + n)},

p
p
√
k max kX i k2 kψ1 ≤ K{a2 d(1 + c2 d) log(1 + n) + b2 d(1 + c2 d) log(1 + n)}.
1≤i≤n

p
Thus, from the inequality (E.1), if a1 > 0, let γn  log(n); otherwise, let γn  log(n).
p
Similarly, if a2 > 0, let κn  log(n); otherwise, let κn  log(n). Then, such γn and κn satisfy
the condition (2.2). Suppose both a1 and a2 are positive, which means both i and Xij ’s have
exponential tails. As before, set κn = γn = log(n)τn . For this case, the regularization parameter
√
specification (2.4) becomes log(n)τn  λ  min{µ? , n}.
xiii

At the end of this supplement, we simply list explicit expressions of κn under different assumptions on the covariates for the case with a diverging number of covariates, which are the extension
of the results in Section 3.2. The magnitude of κn becomes larger than that for the case with d
√
fixed while γn keeps the same. Specifically, if X 0 is bounded with CX > 0, then κn = dCX .
q
2 [(3/2) log(d) + log(n)]. If the
If X 0 follows a Gaussian distribution N (0, ΣX ), then κn = 2dσX
P
Orlicz norm kX0j kψ exists for 1 ≤ j ≤ d and their average (1/d) dj=1 kX0j kψ is bounded, then
κn  dψ −1 (n); for instance, if ψ = ψp with p ≥ 1, then κn  d(log(n))1/p . Finally, if the data {X i }
satisfies the right inequality of (E.2) with a2 > 0, that is, each component of X i is sub-exponentially
tailed, then κn  d3/2 log(n). It is worthwhile to note that these expressions of κn depend on a
factor involving the diverging number of covariates d, which will influence the specification of the
regularization parameter and the sufficient conditions of all the theoretical results in Section 4.

F

Supplement for Section 4

In this supplement, we provide Proposition D.2 and its proof, the proofs of the lemmas in the
appendix and some additional results.
We first extend Proposition D.2 to the case with d → ∞ and d  n. Before that, we list two
simple lemmas for a diverging d. Suppose {ξ i } is a sequence of i.i.d. copies of ξ 0 , a d-dimensional
P
random vector with mean zero. Denote σ̄ξ2 = (1/d) dj=1 Var[ξ0j ].
Lemma F.1. Suppose σ̄ξ2 is bounded. If d/n = o(1), then
n

k

1X
P
ξ i k2 −→ 0.
n
i=1

Lemma F.2. Suppose σ̄ξ2 is bounded. If d/n = o(1), then
n

1X
P
kξ i k2 − P kξ 0 k2 −→ 0.
n
i=1

Suppose the specification of the regularization parameter is given by
dκn  λ, αγn ≤ λ, and λ  µ? ,

(F.1)

where α is a constant greater than 2.
Proposition F.3. Suppose assumptions (D) and (G) hold and the regularization parameter satisfies
√
√
(F.1). Suppose there exist constants C1 and C2 such that kβ ? k2 < C1 d and kβ (0) k2 < C2 d
xiv

√
√
wpg1. If the regularization parameter satisfies (2.4), s1 λκn /(n d) = o(1) and s2 κn γn /(n d) =
o(1), then, for every K ≥ 1, with at least probability pn,K which increases to one as n → ∞,
√
kβ (K+1) − β (K) k2 ≤ O(( ds1 κ2n /n)K d) and kβ (k) k2 ≤ (2C1 + C2 )d for all k ≤ K. Specifically,
wpg1, the iterative algorithm stops at the second iteration.
Proof of Proposition F.3. Reuse the notations in the proof of Lemma D.2. First, we show that,
wpg1, kβ (1) k2 ≤ (2C1 + C2 )d. For each k ≥ 1,
Sβ (k) = SµS11 + SµS12 + SS1 β ? + SS1 + SS2 ∪S3 β (k−1) + λ(SS2 − SS3 ),
Since the regularization parameter satisfies (F.1), it is easy to check that the conclusion of Lemma
4.1 continues to hold, which implies P (A0 ) → 1.
Thus, wpg1,
β (1) = T0−1 T1 + T0−1 T2 + T0−1 T3 + T0−1 T4 (β (0) ) + T0−1 T5 .
We will show that, wpg1,
kT0−1 T1 k2 ≤ (C2 /4)d,
kT0−1 T2 k2 ≤ 2C1 d,
kT0−1 T3 k2 ≤ (C2 /4)d,
kT0−1 T4 (β (0) )k2 ≤ (C2 /4)d,
kT0−1 T5 k2 ≤ (C2 /4)d.
Thus, wpg1,
kβ (1) k2 ≤

5
X

kT0−1 Ti k2 ≤ (2C1 + C2 )d.

i=1

√
On T0−1 T1 . Under Assumption (D), for s2 κn γn /(n d) = o(1), wpg1,
1
1
s2
kT0−1 T1 k2 ≤ k( S)−1 kF k SµS ? k2 ≤ 2kΣ−1
X kF,d √ κn γn d → 0.
n
n 12
n d
Thus, wpg1, kT0−1 T1 k2 ≤ C2 d/4.
On T0−1 T2 . Wpg1,
1
1
kT0−1 T2 k2 ≤ k( S)−1 Ss1 +1,n kF kβ ? k2
n
n
√
√
1
1
≤ kI d kF C1 d + k( S)−1 S1,s1 kF C1 d
n
n
√
1
1 −1
≤ C1 d + k( S) kF k S1,s1 kF C1 d,
n
n
xv

and
s

1
1X
s1
1
kX i k22 ≤ κ2n .
k S1,s1 kF =
n
n
n

i=1

Thus, Under Assumption (D), for s1 κ2n /n = o(1), wpg1,
√
ds1 2 √
−1
−1
kT0 T2 k2 ≤ C1 d + 2kΣX kF,d
κn C1 d ≤ 2C1 d.
n
On T0−1 T3 . Under assumptions (D) and (G), for log(d)/n = o(1), wpg1,
√

1 p
1
1
d√
d log(d)k( S)−1 kF,d (d log(d))−1/2 k √ Ss1 +1,n k2
n
n
n
p
d log(d)
P
√
≤
2kΣ−1
X kF,d OP (1) −→ 0.
n

kT0−1 T3 k2 =

Thus, wpg1, kT0−1 T3 k2 ≤ C2 d/4.
On T0−1 T4 (β (0) ). Under Assumption (D), for s1 κ2n /n, wpg1,
√

1
1
dk( S)−1 kF,d k S1,s1 kF kβ (0) k2
n
n
√
s1 2 √ P
−1
≤ d2kΣX kF,d κn C2 d −→ 0.
n

kT0−1 T4 (β (0) )k2 ≤

Thus, wpg1, kT0−1 T4 (β (0) )k2 ≤ C2 d/4.

√
On T0−1 T5 . Under Assumption (D), for s1 κn λ/(n d) = o(1), wpg1,
√

λ
1
? k2 + kSS ? k2 )
dk( S)−1 kF,d (kSS21
31
n
n
√
λ
≤ d2kΣ−1
X kF,d s1 κn ≤ C2 d/4.
n

kT0−1 T5 k2 ≤

Next, consider kβ 2 − β 1 k2 . Since β (1) ≤ (2C1 + C2 )d wpg1, the conclusion of Lemma 4.1 holds,
which implies A1 occurs wpg1.
Then,
β (2) = T0−1 T1 + T0−1 T2 + T0−1 T3 + T0−1 T4 (β (1) ) + T0−1 T5 ,
where
T4 (β (1) ) =

1
S1,s1 β (1) .
n

Thus, wpg1,
β (2) − β (1) = S−1 S1,s1 (β (1) − β (0) ).

xvi

Thus, for d3/2 s1 κ2n /n = o(1), wpg1,
√

1
1
dk S−1 kF,d k S1,s1 kF kβ (1) − β (0) k2
n
n
√ s1 2
−1
≤ 2kΣX kF,d d κn (2C1 + C2 )d . d3/2 s1 κ2n /n → 0.
n

kβ (2) − β (1) k2 ≤

Thus, wpg1, β (2) = β (1) , which means that, wpg1, the iteration algorithm stops at the second
iteration.
T
For any K ≥ 1, repeating the above arguments, with at least probability pn,K = P ( K
k=0 Ak ),
which increases to one, we have β (k) ≤ (2C1 + C2 )d for k ≤ K and
√ s1 2 K
√
kβ (K+1) − β (K) k2 ≤ (2kΣ−1
κn ) (2C1 + C2 )d . ( ds1 κ2n /n)K d → 0.
X kF,d d
n
This completes the proof.
Next, we provide the proofs of Lemmas A.2, A.3 and A.4 in the appendix.
√
P
Proof of Lemma A.2. Let E = An − A. Note that rd ≥ 1/ d. Then, rd kEkF −→ 0 implies
P

kEkF,d −→ 0. Thus, wpg1, kEkF,d is bounded by a constant C > 0. By Lemma A.1,
−1
−1
kA−1
n − A kF,d ≤ kA kF,d

kEkF,d
kA−1 kF,d kEkF,d
≤ C2
.
−1
1 − CkEkF,d
1 − kA kF,d kEkF,d

Therefore,
rd kEkF
P
−→ 0.
1 − CkEkF,d

−1
2
rd kA−1
n − A kF ≤ C

This completes the proof.
Proof of Lemma A.3. For any δ > 0, we have
P (kΣ̂n − ΣX kF > δ) ≤

d X
d
X
d2
k=1 l=1

δ

n

P(
2

1X
d4 1 2
Xik Xil − σkl )2 ≤
σ̄ .
n
n δ 2 XX
i=1

2 r 2 d4 /(nδ 2 ) = o(1) by Assumption (E2) and for r 2 d4 /n → 0.
Thus, P (rd kΣ̂n − ΣX kF > δ) ≤ σ̄XX
d
d

Thus, Σ̂n is a consistent estimator of ΣX wrt rd k·kF .
Proof of Lemma A.4. Let αd =
n

√

d log d and C1 ≥

√

d

2σξ,max . Then
n

X
1 X
1 X ξij
αd C1
P (k √
ξ i k2 > αd C1 ) ≤
P (| √
| > √ ),
σ
n
n
σj d
i=1
j=1
i=1 j

xvii

where σj is the standard deviation of ξ0j . By Berry and Esseen Theorem (see, for example, P375
√ P
in Shiryaev (1995)), there exists a constant C2 > 0 such that P (k(1/ n) ni=1 ξ i k2 > αd C1 ) ≤
T1 + 2T2 , where
T1 =

d
X

P (|N (0, 1)| >

j=1

αd C1
√ ),
σj d

T2 =

d
X
C2 E|ξ0j |3
√ .
σj3 n
j=1

By noting d2 = o(n),
√
σξ,max d
αd C1
αd C1
√ ) < 2d
√ ) → 0,
φ(
T1 ≤
P (|N (0, 1)| >
α d C1
σξ,max d
σξ,max d
j=1
d
X

T2 ≤

d
X
C2 γξ,max
C2 γξ,max
√ = d 3 √ → 0.
3
σ
n
σmin n
j=1 ξ,min

√ P
Therefore, k(1/ n) ni=1 ξ i k2 = OP (αd ).
Next result is on the consistency of the penalized two-step estimator β̃.
Theorem F.4 (Consistency on β̃). Suppose the assumptions and conditions of Theorem 4.2 hold.
√
P
If rd ≥ 1/ d, then β̃ −→ β ? wrt rd k·k2 .
P
Proof of Theorems F.4. By Theorem 4.2, β̂ −→ β ? wrt rd k·k2 . By Theorem 4.4, P {Iˆ0 = I0 } → 1
√
for rd ≥ 1/ d, where I0 = {s1 + 1, s1 + 2, · · · , s = s1 + s2 , s + 1, · · · , n}. Then, wpg1,

β̃ − β ? = R1 + R2 + T0−1 T1 + T0−1 T2 ,
where R1 = (X TIˆ X Iˆ0 )−1 X TIˆ Y Iˆ0 {Iˆ0 6= I0 }, R2 = −(X TI0 X I0 )−1 X TI0 Y I0 {Iˆ0 6= I0 } and Ti ’s are
0

0

defined in the proof of Theorem 4.2. Then,
√
√
rd kβ̃ − β ? k2 ≤ rd kR1 k2 + rd kR2 k2 + kT0−1 kF,d rd dkT1 k2 + kT0−1 kF,d rd dkT2 k2 .
Since P (kR1 k2,d = 0) ≥ P {Iˆ0 = I0 } → 1, we have R1 = oP (1). Similarly, R2 = oP (1). By the proof
√
P
P
of Theorem 4.2, kT0−1 kF,d is bounded and rd dkTi k2 −→ 0 for i = 1, 2. Thus, β̃ −→ β ? wrt rd k·k2
√
and rd ≥ 1/ d.
Finally, we provide some additional results on the asymptotic distributions of β̂ and β̃ with a
√
different scaling. Specifically, the scaling in Section 4 is nAn . Next, we consider another natural
√
1/2
scaling nAn ΣX .

xviii

Theorem F.5 (Asymptotic Distribution on β̂). Suppose assumptions (D’), (D”), (E), (F) and
√
√
(G) hold. If d6 log d = o(n), s1 = o( n/(λdκn )) and s2 = o( n/(dκn γn )), then
√

d

1/2

nAn ΣX (β̂ n − β ? ) −→ N (0, σ 2 G).

Theorem F.6 (Asymptotic Distribution on β̃). Suppose the assumptions and conditions of Theo√
rem F.5 hold except the condition s1 = o( n/(λdκn )). Then
√

d

1/2

nAn ΣX (β̃ − β ? ) −→ N (0, σ 2 G).

By Theorems F.5 and F.6, Wald-type confidence regions can be constructed. In order to validate
these confidence regions with estimated σ and ΣX , we need Lemma A.5 and the following result.
Theorem F.7 (Asymptotic Distributions on β̂ and β̃ with Σ̂n ). Suppose the assumptions and
conditions of Theorem F.5 hold. If d9 (log(d))2 = o(n), then
√

1/2

d

nAn Σ̂n (β̂ − β ? ) −→ N (0, σ 2 G).

Similarly, suppose the assumptions and conditions of Theorem F.6 hold. If d9 (log(d))2 = o(n), then
√

1/2

d

nAn Σ̂n (β̃ − β ? ) −→ N (0, σ 2 G).

Remark 3. A comparison of the assumptions and conditions of Theorem F.7 with those of Theorems F.5 and F.6 reveals that a much stronger requirement on d is needed to ensure Σ̂n is a good estimator of ΣX . Precisely, the former require that d9 (log(d))2 = o(n) and the latter d6 log(d) = o(n).
This stronger requirement on d is a price paid for estimating ΣX .
Remark 4. The condition on the dimension d in Theorems 4.3 and 4.5 is d5 log(d) = o(n), slightly
weaker than the condition d6 log(d) = o(n) in Theorems F.5 and F.6. Accordingly, The condition
on the dimension d in Theorem A.6 is d8 (log(d))2 = o(n), slightly weaker than the condition
√
d9 (log(d))2 = o(n) in Theorem F.7. This means that the scaling nAn is slightly better than the
√
1/2
scaling nAn ΣX in terms of the condition on d. Further, the former scaling is more suitable for
constructing confidence regions for some entries of β ? .
At the end of this supplement, we provide the proofs of the above theorems.
Proof of Theorems F.5. Reuse the notations Ti ’s in the proof of Theorems 4.2, from which,
√

1/2

nAn ΣX (β̂ n − β ? ) = V1 + V2 + V3 − V4 ,
xix

where Vi = B n Ti for i = 1, 2, 3, 4 and B n =

√

d

1/2

nAn ΣX T0−1 . We will show V2 −→ N (0, σ 2 G) and

other Vi ’s are oP (1), from which the desired result follows by applying Slutsky’s lemma.
√
1/2
On V1 . We have kV1 k2 ≤ ndkAn kF kΣX kF,d kT0−1 kF,d kT1 k2 . By Assumption (F), kAn kF is
1/2

bounded. By Assumption (D”), kΣX kF,d is bounded. By Lemmas A.2 and A.3 and Assumption
(D), for d = o(n1/3 ), wpg1, kT0−1 kF,d is bounded. Further, wpg1, kT1 k2 ≤
kV1 k2 .

1
n s2 κn γn .

Then,

√1 s2 dκn γn ,
n

where . means that the left side is bounded by a constant times the right
√
side, as noted at the beginning of the appendix. Thus, kV1 k2 = oP (1) for s2 = o( n/(dκn γn )).
On V2 . We have V2 = V21 + V22 , where
V21 =

√

−1/2

nAn ΣX

First, consider V21 . We have V21 =

T2 ,

V22 =

1/2

nAn ΣX (T0−1 − Σ−1
X )T2 .

p
P
(n − s1 )/n ni=s1 +1 Z n,i , where

Z n,i = √
On one hand, for every δ > 0,

√

1
−1/2
An ΣX X i i .
n − s1

Pn

2
i=s1 +1 EkZ n,i k2 {kZ n,i k2

> δ} ≤ (n − s1 )EkZ n,0 k42 /δ 2 and

1
−1/2
−1/2
E4 E(X T0 ΣX ATn An ΣX X 0 )2
(n − s1 )2 0
1
≤
E4 λmax (Gn )λmin (ΣX )−1 E(X T0 X 0 )2
(n − s1 )2 0

EkZ n,0 k42 =

d

≤

1X
d2
4 1/2 2
E40 λmax (Gn )λmin (ΣX )−1 (
(EX0j
) ) .
2
(n − s1 )
d
j=1

P
√
Thus, by assumptions (D’), (E) and (F), ni=s1 +1 EkZ n,i k22 {kZ n,i k2 > δ} → 0 for d = o( n). On
P
the other hand, ni=s1 +1 Cov(Z n,i ) = σ 2 An ATn → σ 2 G. Thus, by central limit theorem (see, for
d

example, Proposition 2.27 in van der Vaart (1998)), V21 −→ N (0, σ 2 G). Next, consider V22 . We
have
1/2
−1/2 √
kV22 k2 ≤ kAn kF kΣX kF,d d(log(d))1/2 kT0−1 − Σ−1
k nT2 k2 .
X kF (d log(d))
1/2

By Assumption (F), kAn kF is O(1); By Assumption (D”), kΣX kF,d is O(1); by Lemmas A.2
6
and A.3, d(log(d))1/2 kT0−1 − Σ−1
X kF is oP (1) for d log(d) = o(n); By Lemma A.4, together with
√
√
Assumption (G), (d log(d))−1/2 k nT2 k2 = (d log(d))−1/2 k √1n Ss1 +1,n k2 is OP (1) for d = o( n).
P

d

Thus, V22 −→ 0. By slutsky’s lemma, V2 −→ N (0, σ 2 G).
√
On V3 and V4 . First consider V3 . By noting that s1 = o( n/(λdκn )), wpg1, kV3 k2 ≤
√
√
1/2
d nkAn kF kΣX kF,d kT0−1 kF,d kT3 k2 . dλs1 κn / n → 0. Thus, kV3 k2 = oP (1). In the same way,
kV4 k2 = oP (1). This completes the proof.
xx

√
1/2
Proof of Theorem F.6. From the proof of Theorem F.4, we have nAn ΣX (β̃ − β ? ) = R̃1 + R̃2 +
√
√
1/2
1/2
V1 + V2 , where R̃1 = nAn ΣX R1 , R̃2 = nAn ΣX R2 , and Ri ’s and Vi ’s are defined in the
proofs of Theorems F.4 and F.5. Since P (kR̃1 k2 = 0) ≥ P {Iˆ0 = I0 } → 1, we have R̃1 = oP (1).
d

Similarly, R̃2 = oP (1). By the proof of Theorem F.5, V1 = oP (1) and V2 −→ N (0, σ 2 G). Thus, the
asymptotic distribution of β̃ is Gaussian by Slutsky’s lemma.
Proof of Theorem F.7. We only show the result on β̂. since the result on β̃ can be obtained in a
similar way. We reuse the definitions of Ti ’s in the proof of Theorems 4.2, from which,
√

where M =

√

1/2

nAn Σ̂n (β̂ n − β ? ) = M + R,

1/2

nAn ΣX (β̂ n − β ? ) and R =

√

1/2

nAn (Σ̂n

d

1/2

− ΣX )(β̂ n − β ? ). By Theorem F.5,

P

M −→ N (0, σ 2 G). Then, it is sufficient to show that R −→ 0 wrt k·k2 . We have
R = R1 + R2 + R3 − R4 ,
where Ri = B n Ti for i = 1, 2, 3, 4 and B n =

√

1/2

nAn (Σ̂n

1/2

− ΣX )T0−1 . We will show each Ri

converges to zero in probability, which finishes the proof.
1/2

1/2

On R1 . By Lemma A.7, kΣ̂n − ΣX kF ≤ (d1/2 kΣ̂n − ΣX kF )1/2 . Then,
kR1 k2 ≤
≤

√
√

1/2

1/2

nkAn kF kΣ̂n − ΣX kF kT0−1 kF kT1 k2
ndkAn kF (kΣ̂n − ΣX kF,d )1/2 kT0−1 kF,d kT1 k2 .

By Assumption (F), kAn kF is bounded. By Lemma A.3, kΣ̂n − ΣX kF,d = oP (1) for d = o(n1/3 ).
By Lemmas A.2 and A.3 and Assumption (D), for d = o(n1/3 ), wpg1, kT0−1 kF,d is bounded.
We have, wpg1, kT1 k2 ≤
√
o( n/(dκn γn )).

1
n s2 κn γn .

Then, kR1 k2 .

√1 s2 dκn γn .
n

Thus, kR1 k2 = oP (1) for s2 =

On R2 . We have
√
1/2
1/2
kR2 k2 ≤ kAn kF d(log(d))1/2 kΣ̂n − ΣX kF kT0−1 kF,d (d log(d))−1/2 k nT2 k2 ,
and
1/2

1/2

d(log(d))1/2 kΣ̂n − ΣX kF ≤ (d5/2 log(d)kΣ̂n − ΣX kF )1/2 .
By Assumption (F), kAn kF is O(1); by Lemma A.3, d5/2 log(d)kΣ̂n −ΣX kF = oP (1) for d9 (log(d))2 =
6
o(n); by Lemmas A.2 and A.3, d(log(d))1/2 kT0−1 − Σ−1
X kF is oP (1) for d log(d) = o(n); by Lemma
√
√
P
A.4, (d log(d))−1/2 k nT2 k2 = (d log(d))−1/2 k √1n Ss1 +1,n k2 is OP (1) for d = o( n). Thus, R2 −→ 0.

xxi

√
On R3 and R4 . First consider R3 . By noting that s1 = o( n/(λdκn )), wpg1,
√
√
1/2
1/2
kR3 k2 ≤ d nkAn kF (kΣ̂n − ΣX kF,d )1/2 kT0−1 kF,d kT3 k2 . dλs1 κn / n → 0.
Thus, kR3 k2 = oP (1). In the same way, kR4 k2 = oP (1).

References
Jianqing Fan and Jinchi Lv. Non-concave penalized likelihood with np-dimensionality. IEEE
Transactions On Information Theory, 57:5467–5484, 2011.
Jianqing Fan and Heng Peng. On non-concave penalized likelihood with diverging number of
parameters. The Annals of Statistics, 32:928–961, 2004.
Michael R. Kosorok. Introduction to Empirical Processes and Semiparametric Inference. Springer
New York, 2008.
Jerzy Neyman and Elizabeth L. Scott. Consistent estimates based on partially consistent observations. Econometrica, 16:1–32, 1948.
Albert N. Shiryaev. Probability. Springer-Verlag, second edition, 1995.
Aad W. van der Vaart. Asymptotic Statistics. Cambridge University Press, 1998.
Peng Zhao and Bin Yu. On model selection consistency of lasso. The Journal of Machine Learning
Research, 7(Nov):2541–2563, 2006.

xxii

