Discussion Paper No.92
IPAD: Stable Interpretable Forecasting with
Knockoffs Inference
Yingying Fan, Jinchi Lv, Mahrad Sharifvaghefi
and Yoshimasa Uematsu

January 29, 2019

Data Science and Service Research
Discussion Paper

Center for Data Science and Service Research
Graduate School of Economic and Management
Tohoku University27-1 Kawauchi, Aobaku
Sendai 980-8576, JAPANIPAD: Stable Interpretable Forecasting with Knockoffs
Inference ∗Yingying Fan1 , Jinchi Lv1 , Mahrad Sharifvaghefi1 and Yoshimasa Uematsu2
University of Southern California1 and Tohoku University2
January 29, 2019

Interpretability and stability are two important features that are desired in many
contemporary big data applications arising in economics and finance. While the former
is enjoyed to some extent by many existing forecasting approaches, the latter in the senseof controlling the fraction of wrongly discovered features which can enhance greatly theinterpretability is still largely underdeveloped in the econometric settings. To this end, in
this paper we exploit the general framework of model-X knockoffs introduced recently in
Candès, Fan, Janson and Lv (2018), which is nonconventional for reproducible large-scaleinference in that the framework is completely free of the use of p-values for significancetesting, and suggest a new method of intertwined probabilistic factors decoupling (IPAD)
for stable interpretable forecasting with knockoffs inference in high-dimensional models.
The recipe of the method is constructing the knockoff variables by assuming a latent factor model that is exploited widely in economics and finance for the association structureof covariates. Our method and work are distinct from the existing literature in that we
estimate the covariate distribution from data instead of assuming that it is known whenconstructing the knockoff variables, our procedure does not require any sample splitting,we provide theoretical justifications on the asymptotic false discovery rate control, andthe theory for the power analysis is also established. Several simulation examples and the∗

Yingying Fan is Dean’s Associate Professor in Business Administration, Data Sciences and Operations
Department, Marshall School of Business, University of Southern California, Los Angeles, CA 90089, USA(E-mail: fanyingy@marshall.usc.edu). Jinchi Lv is Professor and McAlister Associate Professor in Business
Administration, Data Sciences and Operations Department, Marshall School of Business, University of Southern California, Los Angeles, CA 90089, USA (E-mail: jinchilv@marshall.usc.edu). Mahrad Sharifvaghefi is
Ph.D. candidate, Department of Economics, University of Southern California, Los Angeles, CA 90089, USA(E-mail: sharifva@usc.edu). Yoshimasa Uematsu is Assistant Professor, Department of Economics and Management, Tohoku University, Sendai 980-8576, Japan (E-mail: yoshimasa.uematsu.e7@tohoku.ac.jp). Mostof this work was completed while Uematsu visited USC Marshall as JSPS Overseas Research Fellow and
Postdoctoral Scholar. This work was supported by NIH Grant 1R01GM131407-01, NSF CAREER Award
DMS-1150318, a grant from the Simons Foundation, Adobe Data Science Research Award, and a Grantin-Aid for JSPS Overseas Research Fellowship 29-60. We are grateful for comments from participants ofseminars at USC, the 2018 Institute of Mathematical Statistics Asia Pacific Rim Meeting in Singapore, the2018 International Symposium on Financial Engineering and Risk Management in Shanghai, and the 2018ICSA China Statistics Conference in Qingdao.

1real data analysis further demonstrate that the newly suggested method has appealing
finite-sample performance with desired interpretability and stability compared to some
popularly used forecasting methods.

Running title: IPAD
Key words: Reproducibility; Power; Big data; Interpretable forecasting; Stability; Latent factors; Model-X knockoffs; Large-scale inference and FDR; Scalability; Intertwined
probabilistic factors decoupling; Lasso and random forest

1

Introduction

Forecasting is a fundamental problem that arises in economics and finance. With the availability of big data, many machine learning algorithms such as the Lasso and random forest
can be resorted to for such a purpose by exploring a large pool of potential features. Many ofthese existing procedures provide a certain measure of feature importance which can then be
utilized to judge the relative importance of selected features for the goal of interpretability.
Yet the issue of stability in the sense of controlling the fraction of wrongly discovered features
is still largely underdeveloped in the econometric settings. As argued in [20], it is difficult
to obtain interpretability and stability simultaneously even in simple Lasso forecasting. A
natural question is how to ensure both interpretability and stability for flexible forecasting.
Naturally stability is related to statistical inference. The recent years have witnesseda growing body of work on high-dimensional inference in the econometrics and statistics
literature. For example, [42] proposed a simple procedure for inference of the average partial effects based on a debiased `1 -regularized method in approximately sparse panel probit
models. [38] used the de-sparsified estimator for constructing pointwise and group confidence sets. [43] conducted simultaneous inference for high-dimensional sparse linear models
based on a bootstrap and desparsifying Lasso estimator. They also applied their procedure
for the family-wise error rate control. [16] provided a double/debiased machine learning(DML) method for estimation and inference of treatment effects, which utilizes the Neyman orthogonal scores and cross-fitting. [18] then extended this idea to linear functionals.
[17] considered debiased simultaneous inference in a system of high-dimensional regression
equations with temporal and cross-sectional dependency based on a uniform robust postselection procedure. [36] proposed Lasso residual-based tests for checking goodness-of-fit in(low- and) high-dimensional linear models. [29] presented a method for estimating the effectof the treatment on the outcome by using instrumental variables where the instruments are
not necessarily valid.
Most existing work on high-dimensional inference for interpretable models has focused
primarily on the aspects of post-selection inference known as selective inference and debiasing for regularization and machine learning methods. In real applications, one is often
interested in conducting global inference relative to the full model as opposed to local inference conditional on the selected model. Moreover, many statistical inferences are basedon p-values form significance testing. However, oftentimes obtaining valid p-values even for2the Lasso in relatively complicated high-dimensional nonlinear models also remains largely
unresolved, not to mention for the case of more complicated model fitting procedures suchas random forest. Indeed high-dimensional inference is intrinsically challenging even in theparametric settings [27].
The desired property of stability for interpretable forecasting in this paper concentrates
on global inference by controlling precisely the fraction of wrongly discovered features inhigh-dimensional models, which is also known as reproducible large-scale inference. Such aproblem involves testing the joint significance of a large number of features simultaneously,
which is known widely as the problem of multiple testing in statistical inference. For thisproblem, the null hypothesis for each feature states that the feature is unimportant in the
joint model which can be understood as the property that this individual feature and the
response are independent conditional on all the remaining features, while the corresponding
alternative hypothesis states the opposite. Conventionally p-values from the hypothesis
testing are used to decide whether or not to reject each null hypothesis with a significancelevel to control the probability of false discovery in a single hypothesis test, meaning rejecting
the null hypothesis when it is true. When performing multiple hyothesis tests, the probabilityof making at least one false discovery which is known as the family-wise error rate can beinflated compared to that for the case of a single hypothesis test. The work on controllingsuch an error rate for multiple testing dates back to [13], where a simple, useful idea islowering the significance level for each individual test as the target level divided by the
total number of tests to be performed. The Bonferroni correction procedure is, however,well known to be conservative with relatively low power. Later on, [30] proposed a stepdown procedure which is less conservative than the Bonferroni procedure. More recently,[35] suggested a procedure in which the critical values of individual tests are constructedsequentially.
A more powerful and extremely popular approach to multiple testing is the Benjamini–
Hochberg (BH) procedure for controlling the false discovery rate (FDR) which was originatedin [9], where the FDR is defined as the expectation of the fraction of falsely rejected null
hypotheses known as the false discovery proportion. Given the p-values from the multiple
hypothesis tests, this procedure sorts the p-values from low to high and chooses a simple,
intuitive cutoff point, which can be viewed as an adaptive extension of the Bonferroni correction for multiple comparisons, of the p-values for rejecting the null hypotheses. The BH
procedure was shown to be capable of controlling the FDR at the desired level for independent test statistics in [9] and for positive regression dependency among the test statistics in[10], where it was shown that a simple modification of the procedure can control the FDR
under other forms of dependency but such a modification is generally conservative. Thereis a huge literature on the theory, applications, and various extensions of the original BH
procedure for FDR control. See, for instance, [8] for a review of related developments, [24]for a factor model approach to FDR control under arbitrary covariance dependence, and [7]for a review of key results on estimation and inference including multiple testing with FDR
control in high-dimensional models.
It is worth mentioning that [19] recently introduced a one covariate a time, multiple3

testing procedure for high-dimensional variable selection in linear regression models. In particular, their method was shown to have asymptotic FDR equal to the ratio of the number ofpseudo signals and the total number of pseudo signals and true signals, where the true signals
have nonzero regression coefficients and the pseudo signals have zero regression coefficients
but nonzero marginal correlations with the response. Unlike [19], the main interest of ourpaper is the FDR control with respect to only the set of true signals.
The aforementioned econometric and statistical inference methods including the BH-type
procedures for FDR control are all rooted on the availability and validity of computable pvalues for evaluating variable importance. As mentioned before, such a prerequisite canbecome a luxury that is largely unclear how to obtain in high dimensions even for the caseof Lasso in general nonlinear models and random forest. In contrast, [4] proposed a novel
procedure named the knockoff filter for FDR control that bypasses the use of p-values inGaussian linear model with deterministic design matrix, where the dimensionality is no larger
than the sample size, and [5] generalized the method to high-dimensional linear models asa two-step procedure based on sample splitting, where a feature screening approach is used
to reduce the dimensionality to below sample size (see, e.g., [23] and [25]) and then theknockoff filter is applied to the set of selected features after the screening step for selective
inference. The key ingredient of the knockoff filter is constructing the so-called knockoff
variables in a geometrical way that mimic perfectly the correlation structure among theoriginal covariates and can be used as control variables to evaluate the importance of originalvariables. Recently, [15] extended the work of [4] by introducing the framework of model-X
knockoffs for FDR control in general high-dimensional nonlinear models. A crucial distinctionis that the knockoff variables are constructed in a probabilistic fashion such that the joint
dependency structure of the original variables and their knockoff copies is invariant to theswapping of any set of original variables and their knockoff counterparts, which enables us
to go beyond linear models and handle high dimensionality. As a result, model-X knockoffs
enjoys exact finite-sample FDR control at the target level. However, a major assumption in[15] is that the joint distribution of all the covariates needs to be known for the valid FDR
control.
Motivated by applications in economics and finance, in this paper we model the association structure of the covariates using the latent factor model, which reduces effectively
the dimensionality and enables reliable estimation of the unknown joint distribution of all
the covariates. By taking into account the latent factor model structure, we first estimate
the association structure of covariates and then construct empirical knockoffs matrix using
the estimated dependency structure. Our empirical knockoffs matrix can be regarded as an
approximation to the oracle knockoffs matrix in [15] that requires the knowledge of the truecovariate distribution. Exploiting the general framework of model-X knockoffs in [15], wesuggest the new method of intertwined probabilistic factors decoupling (IPAD) for stableinterpretable forecasting with knockoffs inference in high-dimensional models. The innovations of our method and work are fourfold. First, we estimate the covariate distribution fromdata instead of assuming that it is known when constructing the knockoff variables. Second,
our procedure does not require any sample splitting and is thus more practical when the4sample size is limited. Third, we provide theoretical justifications on the asymptotic false
discovery rate control when the estimated dependency structure is employed. Fourth, thetheory for power analysis is also established which reveals that there can be asymptoticallyno power loss in applying the knockoffs procedure compared to the underlying variable selection method. Therefore, FDR control by knockoffs can be a pure gain. Compared to earlier
work, an additional challenge of our study is that knowing the true underlying distribution
does not lead to the most efficient construction of the oracle knockoffs matrix due to thepresence of latent factors. The appealing interpretability and stability of our new method
compared to some popularly used forecasting methods are confirmed with several simulation
and real data examples.
The rest of the paper is organized as follows. Section 2 introduces the model settingwith a review of the model-X knockoffs inference framework and presents the new IPAD
procedure. We establish the asymptotic properties of IPAD in Section 3. Sections 4 and 5
present several simulation and real data examples to showcase the finite-sample performanceand the advantages of our newly suggested procedure compared to some popularly used ones.
We discuss some implications and extensions of our work in Section 6. The proofs of themain results and additional technical details are relegated to the Appendix.

2

Intertwined probabilistic factors decoupling

To facilitate the technical presentation, we will introduce the model setting for the highdimensional FDR control problem in Section 2.1 with a review of the model-X knockoffs
inference framework in Section 2.2, and present the new IPAD procedure in Section 2.3.

2.1

Model setting

Consider the high-dimensional linear regression modely = Xβ + ε,(1)where y ∈ Rn is the response vector, X ∈ Rn×p is the random matrix of a large number ofpotential regressors, β = (β1 , · · · , βp )0 ∈ Rp is the regression coefficient vector, ε ∈ Rn is thevector of model errors, and n and p denote the sample size and dimensionality, respectively.
Here without loss of generality, we assume that both the response and the covariates are
centered with mean zero and thus there is no intercept. Motivated by many applications ineconomics and finance, we further assume that the design matrix X follows the exact factor
model
0

X = F0 Λ0 + E = C0 + E,

(2)

where F0 = (f10 , . . . , fn0 )0 ∈ Rn×r is a random matrix of latent factors, Λ0 = (λ01 , . . . , λ0p )0 ∈Rp×r is a matrix of deterministic factor loadings, and E ∈ Rn×p captures the remaining
variation that cannot be explained by these latent factors. We assume that the number offactors r is fixed but unknown and the components of E are independent and identically5distributed (i.i.d.) from some unknown parametric distribution with cumulative distributionfunction G(·; η 0 ), where η 0 ∈ Rm is a finite-dimensional parameter vector. For simplicity,models (1) and (2) are assumed to have no endogeneity.
In this paper, we focus on the high-dimensional scenario when the dimensionality p canbe much larger than sample size n. Therefore, to ensure model identifiability we impose thesparsity assumption that the true regression coefficient vector β has only a small portionof nonzeros; specifically, β takes nonzero values only on some (unknown) index set S 0 ⊂{1, . . . , p} and βj = 0 for all j ∈ S 1 := {1, . . . , p}\S 0 . Denote by s = |S 0 | the size of S 0 . We
assume that s = o(n) throughout the paper.
We are interested in identifying the index set S 0 with a theoretically guaranteed error
rate. To be more precise, we try to select variables in S 0 while keeping the false discovery
rate (FDR) under some prespecified desired level q ∈ (0, 1), where the FDR is defined as
FDR := E [FDP]with

FDP :=|Sb ∩ S 1 |
.
b ∨1
|S|

(3)

Here the FDP stands for the false discovery proportion and Sb represents the set of variables
selected by some procedure using observed data (X, y). A slightly modified version of FDR
is defined as"#|Sb ∩ S 1 |mFDR := E.
(4)b + q −1|S|
Clearly, FDR is more conservative than mFDR in that the latter is always under control if
the former is.
It is easy to see that FDR is a measurement of type I error for variable selection. The
other important aspect of variable selection is power, which is defined as"#"#|Sb ∩ S 0 ||Sb ∩ S 0 |Power := E
=E.
(5)
|S 0 |s
It is well known that FDR and power are two sides of the same coin. We aim at developinga variable selection procedure with theoretically guaranteed FDR control and meanwhileachieving high power.

2.2

Review of model-X knockoffs framework

The key idea of the model-X knockoffs framework is to construct the so-called model-X
knockoff variables, which were introduced originally in [15] and whose definition is statedformally as follows for completeness.
Definition 1 (Model-X knockoff variables [15]) For a set of random variables x =e1 , · · · , Xep ) is called a set of model-Xe = (X(X1 , . . . , Xp ), a new set of random variables x
knockoff variables if it satisfies the following properties:e]swap(S) = [x, xe] in distribution, where the1) For any subset S ⊂ {1, . . . , p}, we have [x, xee]swap(S) is obtained by swapping Xj and Xj for each j ∈ S.
vector [x, x6e is independent of response Y .
2) Conditional on x, the knockoffs vector x
An important consequence is that the null regressors {Xj : j ∈ S 1 } can be swapped with theirknockoffs without changing the joint distribution of the original variables x, their knockoffse, and response Y . That is, we can obtain for any S ⊂ S 1 ,xde]swap(S) , Y ) = ([x, xe], Y ),([x, x(6)dwhere = denotes equal in distribution. Such a property is known as the exchangeability
property using the terminology in [15]. For more details, see Lemma 3.2 therein. Followinge ∈ Rn×p given observed design matrix X.
[15], one can obtain a knockoffs matrix X
e and response vector y constructed by stacking
Using the augmented design matrix [X, X]e y), j ∈the n observations, [15] suggested constructing knockoff statistics Wj = wj ([X, X],{1, . . . , p}, for measuring the importance of the jth variable, where wj is some functionthat satisfies the property that swapping xj ∈ Rn with its corresponding knockoff variableej ∈ Rn changes the sign of Wj ; that is,xw ([X, X],
e y),
j∈
/ S,
j
e swap(S) , y) =wj ([X, X](7)
−w ([X, X],e y), j ∈ S.
je y) satisfy the so-called sign-flip
The knockoff statistics constructed above Wj = wj ([X, X],
property; that is, conditional on |Wj |’s the signs of the null Wj ’s with j 6∈ S 0 are i.i.d. coinflips (with equal chance 1/2). For the examples on valid constructions of knockoff statistics,
see [15].
Let t > 0 be a fixed threshold and define Sb = {j : Wj ≥ t} as the set of discovered
variables. Then intuitively, the sign-flip property entailsd
Sb ∩ S 1 = {j : Wj ≤ −t} ∩ S 1 ≤ |{j : Wj ≤ −t}| .

Therefore, the FDP function can be estimated (conservatively) as
FDP =|{j : Wj ≤ −t}|
|Sb ∩ S 1 |[≤=: FDPbb ∨1|S| ∨ 1|S|for each t. In light of this observation, [15] proposed to choose the threshold by resorting to[ Their results are summarized formally as follows.
the above FDP.
Result 1 ([15]) Let q ∈ (0, 1) denote the target FDR level. Assume that we choose athreshold T1 > 0 such that|{j : Wj ≤ −t}|≤q
T1 = min t > 0 :|{j : Wj ≥ t}| ∨ 1
or T1 = +∞ if the set is empty. Then the procedure selecting the variables Sb = {j : Wj ≥ T1 }controls the mFDR in (4) to no larger than q. Moreover, assume that we choose a slightly
more conservative threshold T2 > 0 such that

1 + |{j : Wj ≤ −t}|
T2 = min t > 0 :≤q|{j : Wj ≥ t}| ∨ 1
7

or T2 = +∞ if the set is empty. Then the procedure selecting the variables Sb = {j : Wj ≥ T2 }controls the FDR in (3) to no larger than q.
It is worth noting that Result 1 was derived under the assumption that the joint distribution of the p covariates is known. In our model setting (1) and (2), however there exist
unknown parameters that need to be estimated from data. In such case, it is natural toconstruct the knockoff variables and knockoff statistics with estimated distribution of the pcovariates. Such a plug-in principle usually leads to breakdown of the exchangeability property in Definition 1, preventing us from using directly Result 1. To address this challenging
issue, we will introduce our new method in the next section and provide detailed theoreticalanalysis for it.
It is also worth mentioning that recently, [6] provided an elegant new line of theory which
ensures FDR control of model-X knockoffs procedure under the approximate exchangeability
assumption, which is weaker than the exact exchangeability condition required in Definition
1. However, the conditions they need on estimation error of the joint distribution of x
is difficult to be satisfied in high dimensions. [26] investigated the robustness of model-X
knockoffs procedure with respect to unknown covariate distribution when covariates x follow
a joint Gaussian distribution. Their procedure needs data splitting and their proofs rely
heavily on the Gaussian distribution assumption, and thus their development may not besuitable for economic data with limited sample size and heavy-tailed distribution. For thesereasons, our results complement substantially those in [15], [26], and [6].

2.3

IPAD

It has been seen from the previous section that the key for the model-X knockoffs frameworkis the construction of valid knockoff variables. We begin with the ideal situation where thethe factor model structure (2) is fully available to us; that is, we know the realization C0 andthe distribution G(·; η 0 ) for the error matrix E. In such case, the oracle knockoffs matrixe 0 ) can be constructed asX(θe 0 ) = C0 + Eη0 ,
X(θ

(8)

where Eη0 is an i.i.d. copy of E and θ 0 = (C0 , η 0 ) is the augmented parameter vector. Notethat Eη0 itself is not a function of η0 , but we slightly abuse the notation to emphasize thee 0)dependence of the distribution function on parameter η0 . It is easy to check that X(θconstructed above is a valid knockoffs matrix and satisfies the properties in Definition 1.
e 0 ) is generally unavailable to us, it plays an important role in our theoretical
Although X(θdevelopments.
We remark that in the construction above, we slightly misuse the concept and call C0 a
parameter. This is because although C0 is a random matrix, for the construction of valid
knockoff variables it is the particular realization C0 leading to the observed data matrix
X that matters. In other words, a valid construction of knockoff variables requires theknowledge of the specific realization C0 instead of the distribution of C0 . To understandthis, consider the scenario where the underlying parameter η 0 and the exact distribution of8

C0 are fully available to us. If we independently generate random variables from this known
distribution and form a new data matrix X1 , because of the independence between X1 andX, the exchangeability assumption in Definition 1 will be violated and thus X1 cannot bea valid knockoffs matrix. On the other hand, as long as we know the realization C0 ande 0 ) can be constructed using (8) regardless ofparameter η 0 , a valid knockoffs matrix X(θwhether the exact distribution of C0 is available to us or not.
e 0 ) is inaccessible. To
In practice, however θ 0 is unavailable to us and consequently, X(θovercome this difficulty, we next introduce our new method IPAD. We start with introducing
the knockoff generating function – for each given parameter vector θ = (C, η), defineeX(θ)= C + Eη ,(9)where Eη is a matrix composed of i.i.d. random samples from the distribution G(·; η). Lete θ̂) as the empiricalting θ̂ denote an estimator (obtained using data X) of θ 0 , we name X(e 0 ) as the oracle (ideal) knockoffs matrix.
knockoffs matrix while X(θWith the aid of empirical knockoffs matrix, we suggest the following IPAD procedure for
FDR control with knockoffs inference.
Procedure 1 (IPAD)
1) (Estimation of parameters) Estimate the unknown parame0b η̂) the resulting estimatedters in θ using the design matrix X. Denote by θ̂ = (C,parameter vector.
2) (Construction of empirical knockoffs matrix) Construct the empirical knockoffs matrix
by applying the knockoff generating function in (9) to the estimated parameter θ̂; thatis,e θ̂) = Cb + Eη̂ ,X((10)

where Eη̂ ∈ Rn×p is a matrix composed of i.i.d. random variables from G(·; η̂), and isindependent of (X, y) conditional on η̂.
3) (Application of knockoffs inference) Calculate knockoff statistics Wj (θ̂) using datae θ̂)], y) and then construct Sb by applying knockoffs inference to Wj (θ̂).
([X, X(
Intuitively, the accuracy of the estimator θ̂ in Step 1 will affect the performance of our
IPAD procedure. In fact, as shown later in our Theorem 1, the consistency rate of θ̂ isindeed reflected in the asymptotic FDR control of the IPAD procedure. For the specificcase when the error distribution is N (0, σ 2 ), the parameter σ 2 can be estimated naturallyP
b = X − C.
b In Step 3, various methods
as (np)−1 i,j ê2ij , where êij ’s are the entries of E
can be used to construct knockoff statistics. For the illustration purpose, we use the Lasso
coefficient difference (LCD) statistic as in [15]. Specifically, with y the response vector ande θ̂)]) the augmented design matrix we consider the variable selection procedure Lasso([X, X([39] which solves the following optimization problemnoe θ̂)]bk2 + λkbk1 ,β̂ aug (θ̂; λ) = arg min ky − [X, X(
(11)
2
b∈R2p

9

where λ ≥ 0 is the regularization parameter and k · km with m ≥ 1 denotes the vector`m -norm. Then for each variable xj , the knockoff statistic can be constructed asaug
Wj (θ̂; λ) = |β̂jaug (θ̂; λ)| − |β̂p+j(θ̂; λ)|,(12)where β̂`aug (θ̂; λ) is the `th component of the Lasso regression coefficient vector β̂ aug (θ̂; λ).
It is seen that intuitively the LCD knockoff statistics evaluate the relative importance of thejth original variable by comparing its Lasso coefficient β̂jaug (θ̂; λ) with that of its knockoffauge 0 ) is used instead of
copy β̂j+p(θ̂; λ). In the ideal case when the oracle knockoffs matrix X(θe θ̂) in (11), it is easy to verify that the LCD is a valid construction of knockoff statistics
X(and satisfies the sign-flip property in (7). Consequently, the general theory in [15] can beapplied to show that the FDR is controlled in finite sample. We next show that even with theempirical knockoffs matrix employed in (11), the FDR can still be asymptotically controlledwith delicate technical analyses.

3

Asymptotic properties of IPAD

We now provide theoretical justifications for our IPAD procedure suggested in Section 2e θ̂)], y; λ) defined in (12). We willwith the LCD knockoff statistics Wj (θ̂; λ) = wj ([X, X(
first present some technical conditions in Section 3.1, then prove in Section 3.2 that the
FDR is asymptotically under control at desired target level q, and finally in Section 3.3 showthat asymptotically IPAD has no power loss compared to the Lasso under some regularity
conditions.

3.1

Technical conditions

We first introduce some notation and definitions which will be used later on. We use X ∼subG(Cx2 ) to denote that X is a sub-Gaussian random variable with variance proxy Cx2 > 0
if E[X] = 0 and its tail probability satisfies P(|X| > u) ≤ 2 exp(u2 /Cx2 ) for each u ≥ 0. In
all technical assumptions below, we use M > 1 to denote a large enough generic constant.
Throughout the paper, for any vector v = (vi ) let us denote by kvk1 , kvk2 , and kvkmax the `1 P
Pnorm, `2 -norm, and max-norm defined as kvk1 = i |vi |, kvk2 = ( i vi2 )1/2 , and kvkmax =maxi |vi |, respectively. For any matrix M = (mij ), we denote by kMkF , kMk1 , kMk2 ,and kMkmax the Frobenius norm, entrywise `1 -norm, spectral norm, and entrywise `∞ -norm
defined as kMkF = k vec(M)k2 , kMk1 = k vec(M)k1 , kMk2 = supv6=0 kMvk2 /kvk2 , andkMkmax = k vec(M)kmax , respectively, where vec(M) represents the vectorization of matrix
M. For a symmetric matrix M, vech(M) stands for the vectorization of the lower triangularpart of M.
Condition 1 (Regression errors) The model error vector ε has i.i.d. components from
subG(Cε2 ).
Condition 2 (Latent factors) The rows of F0 consist of mean zero i.i.d. random vectors
fi0 ∈ Rr such that kF0 kmax ≤ M almost surely (a.s.) and kΣf k2 + kΣ−1
f k2 ≤ M , where000Σf := E[fi fi ].
10

Condition 3 (Factor loadings) The rows of Λ0 consist of deterministic vectors λ0j ∈ Rr00such that kΛ0 kmax ≤ M and kp−1 Λ0 Λ0 k2 + k(p−1 Λ0 Λ0 )−1 k2 ≤ M .
Condition 4 (Factor errors) The entries of matrix Eη0 are i.i.d. copies of eη0 ∼ subG(Ce2 )with continuous distribution function G(·; η 0 ). For each 1 ≤ ` ≤ m, the `th element of η 0m
is specified as η`0 = h` (E[eη0 ], . . . , E[em
η 0 ]) with h` : R → R some local Lipschitz continuousfunction in the sense thath` (t1 , . . . , tm ) − h` (E[eη0 ], . . . , E[em
η 0 ]) ≤ M

maxk∈{1,...,m}tk − E[ekη0 ]for each tk ∈ {t : |t − E[ekη0 ]| ≤ M cnp } and 1 ≤ k ≤ m, where cnp := (p−1 log n)1/2 +(n−1 log p)1/2 . Moreover, there exists some stochastic process (eη )η such thati) for each η ∈ {η ∈ Rm : kη − η 0 kmax ≤ M cnp }, the entries of Eη in (9) have identical
distribution to eη ,ii) for some sub-Gaussian random variable Z ∼ subG c2e with some positive constant ce ,supη: kη−η 0 kmax ≤M cnp

1/2
|Z|.
|eη − eη0 | ≤ M 1/2 cnp

(13)

0

Condition 5 (Eigenseparation) The r eigenvalues of p−1 Λ0 Λ0 Σf are distinct for all p.
The number of factors r is assumed to be known for developing the theory with simplification, but in practice it can be estimated consistently using methods such as information
criteria [3] and test statistics [1]. The sub-Gaussian assumptions in Conditions 1 and 4 canbe replaced with some other tail conditions as long as similar concentration inequalities hold.
Condition 3 is standard in the analysis of factor models. Stochastic loadings can be assumed
in Condition 3 with some appropriate distributional assumption, such as sub-Gaussianity, atthe cost of much more tedious technical arguments. The boundedness of the eigenvalues ofΣf in Condition 2 is standard while the i.i.d. assumption and boundedness of fi0 are stronger
compared to the existing literature (e.g., [3] and [2]). However, these conditions are imposed
mostly for technical simplicity. In fact, the boundedness condition on fi0 can be replacedwith (unbounded) sub-Gaussian or other heavier-tail assumption whenever concentration
inequalities are available at the cost of slower convergence rates and stronger sample size
requirement. Our theory on FDR control is based on that in [15], which applies only to thecase of i.i.d. rows of design matrix X. This is the main reason for imposing the i.i.d. assumption on εi and fi in Conditions 1 and 2. However, we conjecture that similar results can also
hold in the presence of some sufficiently weak serial dependence in εi and fi . Condition 4
introduces a sub-Gaussian process eη with respect to η. The norm in (13) can be replacedwith any other norm since η is finite dimensional. In the specific case when the componentsof E have Gaussian distribution such that η is a scalar parameter representing variance, by
the the reflection principle for the Wiener process ([12], p.511), eη can be constructed asa Wiener process and the inequality (13) can be satisfied. For more information on subGaussian processes, see, e.g., [41]. Condition 5 guarantees that F̂0 F0 /n is asymptoticallynonsingular, which has been proved in [2] and is used in the proof of Lemma 6 in Appendix.
11

Recall that in the IPAD procedure, we first obtain the augmented Lasso estimator
Denote by Aaug (θ; λ) = supp(β̂ aug (θ; λ)) ⊂{1, . . . , 2p} the active set of the augmented Lasso regression coefficient vector. Throughout
this section, we content ourselves with sparse estimates satisfyingeβ̂ aug (θ; λ) ∈ R2p by regressing y on [X, X(θ)].

|Aaug (θ; λ)| ≤ k/2(14)

for some positive integer k which may diverge with n at an order to be specified later; see,e.g., [28] and [32] for a similar constraint and justifications therein. This can always beachieved since users have the freedom to choose the size of the Lasso model.

3.2

FDR control

b for theTo develop the theory for IPAD, we consider the principle component estimator C
realization C0 . More specifically, we first conduct the singular value decomposition (SVD)
X = UDV0 with U and V the left and right singular matrices and D a diagonal matrixof singular values, and then threshold the diagonal matrix D by setting the smallest n − rsingular values to zero. Let us denote the thresholded matrix as Dr . Then matrix C0 can beb = UDr V0 . Denote by E
b = (êij ) = X − C.
b The estimator η̂ = (η̂1 , · · · , η̂m )0
estimated as C
is constructed as η̂` = h` (Enp ê, . . . , Enp êm ) with h` , 1 ≤ ` ≤ m, introduced in ConditionP
4 and Enp êk = (np)−1 i,j êkij the empirical moments of êij . Throughout our theoreticalanalysis, we consider the regularization parameter fixed at λ = C0 n−1/2 log p with C0 somelarge enough constant for all the Lasso procedures. Therefore, we will drop the dependence ofvarious quantities on λ whenever there is no confusion. For example, we will write Aaug (θ; λ)and β̂ aug (θ; λ) as Aaug (θ) and β̂ aug (θ), respectively.
0 [X, X(θ)]
0 y and defineeeeDenote by U(θ) := n−1 [X, X(θ)]
and v(θ) := n−1 [X, X(θ)]T(θ) := vec(vech U(θ), v(θ)) ∈ RP with P := p(2p + 3). The following lemma statesthat the statistic T(θ) plays a crucial role in our procedure.
Lemma 1 The set of variables Sb selected by Procedure 1 depends only on T(θ).
aug
For any given θ, define the active set A∗ (θ) := Aaug
1 (θ) ∪ A2 (θ) ⊂ {1, . . . , p}, whereaug (θ)} and Aaug (θ) := {j−p : j ∈ {p+1, . . . , 2p}∩Aaug (θ)}.
Aaug
1 (θ) := {j : j ∈ {1, . . . , p}∩A
2
That is, A∗ (θ) is equal to the support of knockoff statistics (W1 (θ), · · · , Wp (θ))0 if there areno ties on the magnitudes of the augmented Lasso coefficient vector β̂ aug (θ).
We next focus on the low-dimensional structure of T(θ) inherited from the augmented
Lasso because it will be made clear that this is the key to controlling the FDR without sample
splitting. For any subset A ⊂ {1, . . . , p}, define a lower-dimensional expression of the vector
as TA (θ) := vec(vech UA (θ), vA (θ)) with UA (θ) the principle submatrix of U(θ) formed bycolumns and rows in A and vA (θ) the subvector of v(θ) formed by components in A. Then ite A (θ)]0 [XA , X
e A (θ)] and vA (θ) = n−1 [XA , X
e A (θ)]0 y.
is easy to see that UA (θ) = n−1 [XA , X
Motivated by Lemma 1, we define a family of mappings indexed by A that describes thee A (θ)], y) that forms TA (θ).
selection algorithm of Procedure 1 with given data set ([XA , X
Formally, define a mapping SA : R|A|(2|A|+3) → 2A as tA 7→ SA (tA ) for given TA (θ) = tA ,12where 2A refers to the power set of A. That is, SA (tA ) represents the outcome of first
restricting ourselves to the smaller set of variables A and then applying IPAD to TA (θ) = tA
to further select variables from set A.
Lemma 2 Under Conditions 1–4, for any subset A ⊃ A∗ (θ) we have S{1,...,p} (T(θ)) =SA (TA (θ)).
When restricting on set A, we can apply Procedure 1 to a lower-dimensional data set
e A (θ)], y) that forms TA (θ) to further select variables from A. The previous two([XA , Xlemmas ensure that this gives us a subset of A that is identical to S{1,...,p} (T(θ)). Note thatthe lower-dimensional problem based on TA (θ) can be easier compared to the original one.
We also would like to emphasize that the dimensionality reduction to a smaller model A is
only for assisting the theoretical analysis and our Procedure 1 does not need any knowledgeof such set A.
It is convenient to define t0 = E T(θ 0 ) ∈ RP . Denote byno
I := t ∈ RP : kt − t0 kmax ≤ anp := C1 (k 1/2 + s3/2 )c̃np ,(15)
where C1 is some positive constant and c̃np = p−1/2 log n + n−1/2 log p. For any subset
A ⊂ {1, · · · , p}, let IA be the subspace of I when taking out the coordinates corresponding
to E TA (θ 0 ). Thus IA ⊂ R|A|(2|A|+3) . In addition to Conditions 1–5, we need an assumption
on the algorithmic stability of Procedure 1.
Condition 6 (Algorithmic stability) For any subset A ⊂ {1, . . . , p} that satisfies |A| ≤k ≤ n ∧ p, there exists a positive sequence ρnp → 0 as n ∧ p → ∞ such that
SA (t2 )4SA (t1 )= O(ρnp ),|S
A (t1 )| ∧ |SA (t2 )||A|≤k t1 ,t2 ∈IA
sup

sup

where 4 stands for the symmetric difference between two sets.
Intuitively the above condition assumes that the knockoffs procedure is stable with respect
to a small perturbation to the input t in any lower-dimensional subspace IA . Under theseregularity conditions, the asymptotic FDR control of our IPAD procedure can be established.
Theorem 1 (Robust FDR control) Assume that Conditions 1–6 hold. Fix an arbitrary
positive constant ν. If (s, k, n, p) satisfies s ∨ k ≤ n ∧ p, cnp ≤ c/[r2 M 2 C(ν + 2)]1/2 , and(k 1/2 + s3/2 )c̃np → 0 as n ∧ p → ∞ with c and C some positive constants defined in Lemma
7 in Appendix, then the set of variables Sb obtained by Procedure 1 (IPAD) with the LCD
knockoff statistics controls the FDR (3) to be no larger than q + O (ρnp + n−ν + p−ν ).
Recall that by definition, the FDR is a function of T(θ̂) and can be written as E FDP(T(θ̂))
while the FDR computed with the oracle knockoffs, E FDP(T(θ 0 )), is perfectly controlled to
be no larger than q. This observation motivates us to first establish asymptotic equivalenceof T(θ̂) and T(θ 0 ) with large probability. Then a natural idea is to show that E FDP(T(θ̂))
13

converges to E FDP(T(θ 0 )) in probability, which turns out to be highly nontrivial because of
the discontinuity of FDP(·) (the convergence would be straightforward via the Portmanteau
lemma if FDP(·) were continuous). Condition 6 above provides a remedy to this issue by
imposing the algorithmic stability assumption.

3.3

Power analysis

We have established the asymptotic FDR control for our IPAD procedure in Section 3.2. We
now look at the other side of the coin – the power (5). Recall that in IPAD, we apply theknockoffs inference procedure to the knockoff statistics LCD, which are constructed using
the augmented Lasso in (11). Therefore the final set of variables selected by IPAD is a subset
of variables picked by the augmented Lasso. For this reason, the power of IPAD is always
upper bounded by that of Lasso. We will show in this section that there is in fact no power
loss relative to the augmented Lasso in the asymptotic sense.
Condition 7 (Signal strength I) For any subset A ⊂ S 0 that satisfies |A|/s > 1 − γ forsome γ ∈ (0, 1], it holds that kβA k1 > bnp sn−1/2 log p for some positive sequence bnp → ∞.
Condition 8 (Signal strength II) There exists some constant C2 ∈ (2(qs)−1 , 1) such that
|S2 | ≥ C2 s with S2 = {j : |βj |  (s/n)1/2 log p}.
Condition 7 requires that the overall signal is not too weak, but is weaker than theconventional beta-min condition minj∈S 0 |βj |  n−1/2 log p. Under Condition 8, we can showb ≥ C2 s with probability at least 1 − O(p−ν + n−ν ) using similar techniques to thosethat |S|of Lemma 6 in [26]. The intuition is that given s → ∞, for a variable selection procedure tohave high power it should select at least a reasonably large number of variables. The resultb ≥ C2 s will be used to derive the asymptotic order of threshold T , which is in turn crucial|S|
to establish the theorem below on power.
Theorem 2 (Power guarantee) Assume that Conditions 1–5 and 7–8 hold. Fix an arbitrary positive constant ν. If (s, k, n, p) satisfies 2s ≤ k ≤ n ∧ p, cnp ≤ c/(r2 M 2 C(ν + 2))1/2 ,and sk 1/2 c̃np → 0 as n ∧ p → ∞ with c and C some positive constants defined in Lemma 7,
then both the Lasso procedure based on (X, y) and our IPAD procedure (Procedure 1) havepower bounded from below by γ − o(1) as n ∧ p → ∞. In particular, if γ = 1 IPAD has no
power loss compared to Lasso asymptotically.

4

Simulation studies

We have shown in Section 3 that IPAD can asymptotically control the FDR in high-dimensionalsetting and there can be no power loss in applying the procedure. We next move on to numerically investigate the finite-sample performance of IPAD using synthetic data sets. We
will compare IPAD with the knockoff filter in [4] (BCKnockoff) and the high-dimensional
knockoff filter in [5] (HD-BCKnockoff). In what follows, we will first explain in detail themodel setups and simulation settings, then discuss the implementation of the aforementioned
methods, and finally summarize the comparison results.
14

4.1

Simulation designs and settings

In all simulations, the design matrix X ∈ Rn×p is generated from the factor model
X = F0 (Λ0 )0 +√rθE = C0 +

√
rθE,

(16)

where F0 = (f10 , · · · , fn0 )0 ∈ Rn×r is the matrix of latent factors, Λ0 = (λ01 , . . . , λ0p )0 ∈ Rp×r
is the matrix of factor loadings, E ∈ Rn×p is the matrix of model errors, and θ is a constant√controlling the signal-to-noise ratio. The term r is used to single out the effect of the numberof factors in calculating the signal-to-noise ratio in factor model (16). We then rescale eachcolumn of X to have `2 -norm one and simulate the response vector y = (y1 , · · · , yn )0 fromthe following model√(17)
yi = f (xi ) + cεi , i = 1, · · · , n,where f : Rp → R is the link function which can be linear or nonlinear, c > 0 is a constantcontrolling the signal-to-noise ratio, and ε = (ε1 , · · · , εn )0 is the vector of model error. We
next explain the four different designs of our simulation studies.
4.1.1

Design 1: linear model with normal factor design matrix

The elements of F0 , Λ0 , E, and ε are drawn independently from N (0, 1). The link function
takes a linear form, that is,√
y = Xβ + cε,
where the coefficient vector β = (β1 , · · · , βp )0 ∈ Rp is generated by first choosing s random
locations for the true signals and then setting βj at each location to be either A or −A
randomly with A some positive value. The remaining p − s components of β are set to zero.
4.1.2

Design 2: linear model with fat-tail factor matrix and serial dependence

The elements of E are generated as
eij =

ν−2
χ2ν,j

!
uij ,

(18)

where uij ∼ i.i.d. N (0, 1) for all i = 1, · · · , n and j = 1, · · · , p, and χ2ν,j , j = 1, · · · , p are
i.i.d. random variables from chi-square distribution with ν = 8 degrees of freedom. The
rest of the design is the same as in Design 1. It is worth mentioning that in this case, theentries of matrix E have fat-tail distribution with serial dependence in each column becauseof the common factor χ2ν,j . This design is used to check the robustness of IPAD method withrespect to the serial dependence and the fat-tail distribution of E.
4.1.3

Design 3: linear model with misspecified design matrix

To evaluate the robustness of IPAD procedure to the misspecification of the factor model
structure (16), we set Λ = 0, rθ = 1 and simulate the rows of matrix E independently fromN (0, Σ) with Σ = (σij ), σij = ρ|i−j| for ≤ i, j ≤ p. The remaining design is the same as in15

Design 1. It is seen that our assumption on the independence of the entries of E is violated.
This design is used to test the robustness of IPAD to misspecification of the factor model
structure of X.
4.1.4

Design 4: nonlinear model with normal factor design matrix

Our last design is used to evaluate the performance of IPAD method when the link function
f is nonlinear. To be more specific, we assume the following nonlinear model between the
response and covariates√y = sin(Xβ) exp(Xβ) + cε,
where the coefficient vector β, design matrix X, and model error ε are generated similarly
as in Design 1.
4.1.5

Simulation settings

The target FDR level is set to be q = 0.2 in all simulations. For Design 1 and Design 2, we
set n = 2000, p = 2000, A = 4, s = 50, c = 0.2, r = 3, and θ = 1. In order to evaluate thesensitivity of our method to the dimensionality p and the model sparsity s, we also explore
the settings of p = 1000, 3000 and s = 100, 150. In Design 3, we set r = 0 and ρ = 0, 0.5.
In Design 4, since the model is nonlinear, we use nonparametric method to fit the modeland consider lower-dimensional settings of p = 50, 250, 500. We also decrease the number ofobservations to n = 1000 and number of true variables to s = 10. Moreover, we set θ = 1, 2
and c = 0.1, 0.2, 0.3 to test the effects of signal-to-noise ratio on the performance of IPAD
procedure in Design 4.

4.2

Estimation procedure

In implementing the IPAD algorithm suggested in Section 2, we use the P Cp1 criterion
proposed in [3] to estimate the number of factors r. With an estimated number of factorsbr̂, we use the principle component method discussed in Section 3.2 to obtain an estimate C0
b = (êij ) = X − C.
b Recall that in the construction of knockoffof matrix C . Denote by E
variables, the distribution of E needs to be estimated. Throughout our simulation studies,
we misspecify the model and treat the entries of E as i.i.d. Gaussian random variables.
Under this working model assumption, the only unknown parameter is the variance which
can be estimated by the following maximum likelihood estimator2−1σ̂ = (np)pn XXê2ij .

i=1 j=1b is constructed using (10) with the entries of Eη̂ drawn indepenThen the knockoffs matrix X
dently from N (0, σ̂ 2 ). For the two comparison methods BCKnockoff and HD-BCKnockoff,
we follow the implementation in [4] and [5], respectively. Thus it is seen that neither BCKnockoff nor HD-BCKnockoff uses the factor structure in X when constructing the knockoff
variables.
16

b we apply the Lasso
In Designs 1–3, with the constructed empirical knockoffs matrix X
b the augmented design matrix.
method to fit the model with y the response vector and [X, X]
Then the LCD discussed in Section 2.3 is used in the construction of knockoff statistics. In
Design 4, we assume the nonlinear relationship between the response and the covariates. In
this case, random forest is used for estimation of the model. To construct the knockoff statistics, we use the variable importance measure of mean decrease accuracy (MDA) introducedin [14]. This measure is based on the idea that if a variable is unimportant, then rearranging
its values should not degrade the prediction accuracy. The MDA for the jth variable, de\ j , measures the amount of increase in prediction error when the values of thenoted as MDA\j
jth variable in the out-of-sample prediction are permuted randomly. Then intuitively, MDA
will be small and around zero if the jth variable is unimportant in predicting the response.
\ j | − |MDA
\ j+p |, j = 1, · · · , p.
For each original variable xj , we compute Wj statistic as |MDA

4.3

Simulation results

For each method, we use 100 simulated data sets to calculate its empirical FDR and power,
which are the average FDP and TDP (true discovery proportion as in (5)) over 100 repetitions, respectively. Two different thresholds, knockoff and knockoff+ (T1 and T2 in Result
1, respectively), are used in the knockoffs inference implementation. It is worth mentioningthat as shown in [15] and summarized in Result 1, knockoff+ controls FDR (3) exactly whileknockoff controls only the modified FDR (4).
Tables 1 and 2 summarize the results from Designs 1 and 2, respectively. As shownin Table 1, all approaches can control empirical FDR at the target level (q = 0.2) andknockoff+, which is more conservative, reduces power negligibly. It is worth mentioning thateven for Design 2, in which the design matrix X is drawn from fat-tail distribution with
serial dependence, we still have FDR under control with decent level of power. This suggeststhat the no serial correlation assumption in our theoretical analysis could just be technical.
Compared to the results by BCKnockoff and HD-BCKnockoff, we see that using the extra
information from the factor structure in constructing knockoff variables can help with both
FDR and power. Table 2 also shows the effects of model sparsity on the performance ofvarious approaches. It can be seen that when the number of true signals is increased from50 to 150, the FDR is still under control and the empirical power of IPAD remains steady.
Table 3 is devoted to the case of Design 3, where the rows of matrix X are generated
independently from multivariate normal distribution with AR(1) correlation structure. Thisis a setting where the factor model structure in X is misspecified. Since BCknockoff and
HD-BCknockoff make no use of the factor structure in generating knockoff variables, in bothlow- and high-dimensional examples both methods control FDR exactly at the target level.
IPAD based methods have empirical FDR slightly over the target level, which may be causedby the misspecification of the factor structure. On the other hand, IPAD based approaches
have much higher empirical power than comparison methods.
Table 4 corresponds to Design 4 in which response y is related to X nonlinearly. Since
BCKnockoff and HD-BCKnockoff are designed for linear models, only the results from IPAD

17

Table 1: Simulation results for Designs 1 and 2 of Section 4.1 with different values of dimensionality p
Design 1
FDR

Power

FDR+

Design 2
Power+

R2

FDR

Power

FDR+

Power+

R2

0.961
0.887

0.180
0.152

0.960
0.885

0.652
0.653

0.935
0.607

0.183
0.113

0.933
0.581

0.656
0.644

0.913
0.559

0.171
0.099

0.911
0.524

0.658
0.651

p = 1000
IPAD
BCKnockoff

0.195
0.207

0.991
0.942

0.180
0.192

0.990
0.938

0.659
0.659

0.199
0.172

p = 2000IPAD
HD-BCKnockoff

0.194
0.142

0.979
0.706

0.179
0.127

0.979
0.691

0.649
0.649

0.199
0.136

p = 3000IPAD
HD-BCKnockoff

0.191
0.172

0.964
0.668

0.176
0.149

0.963
0.658

0.652
0.652

0.188
0.125

Note that FDR+ and Power+ are the values of FDR and Power corresponding to the knockoff+ threshold
T2 .

Table 2: Simulation results for Designs 1 and 2 of Section 4.1 with different sparsity level s
Design 1
FDR

Power

FDR+

Design 2
Power+

R2

FDR

Power

FDR+

Power+

R2

0.935
0.607

0.183
0.113

0.933
0.581

0.656
0.644

0.937
0.583

0.174
0.097

0.936
0.573

0.789
0.778

0.935
0.578

0.182
0.106

0.935
0.570

0.848
0.843

s = 50IPAD
HD-BCKnockoff

0.194
0.142

0.979
0.706

0.179
0.127

0.979
0.691

0.649
0.649

0.199
0.136

s = 100
IPAD
HD-BCKnockoff

0.191
0.152

0.978
0.703

0.183
0.140

0.977
0.698

0.783
0.787

0.181
0.106

s = 150IPAD
HD-BCKnockoff

0.183
0.139

0.973
0.660

0.178
0.130

0.972
0.654

0.842
0.858

0.188
0.115

Table 3: Simulation results for Design 3 of Section 4.1

ρ=0
FDR

Power

FDR+

ρ = 0.5
Power+

R2

FDR

Power

FDR+

Power+

R2

0.984
0.827

0.216
0.117

0.984
0.821

0.446
0.445

0.982
0.522

0.202
0.090

0.980
0.467

0.445
0.442

0.979
0.439

0.206
0.064

0.978
0.381

0.443
0.447

p = 1000
IPAD
BCKnockoff

0.204
0.188

0.995
0.919

0.189
0.172

0.995
0.917

0.444
0.444

0.226
0.137

p = 2000IPAD
HD-BCKnockoff

0.203
0.151

0.993
0.630

0.189
0.126

0.993
0.603

0.447
0.449

0.220
0.115

p = 3000IPAD
HD-BCKnockoff

0.225
0.150

0.988
0.589

0.205
0.126

0.987
0.560

0.445
0.446

18

0.219
0.092

method are reported. It can be seen form Table 4 that IPAD approach can control FDR
with reasonably high power even in the nonlinear setting. We also observe that in nonlinear
setting, the power of IPAD deteriorates faster as dimensionality p increases compared to thelinear setting.
Table 4: Simulation results for Design 4 of Section 4.1
θ=1

c = 0.1
c = 0.2
c = 0.3

θ=2

FDR

Power

FDR+

Power+

R2

0.109
0.137
0.137

0.839
0.847
0.765

0.081
0.068
0.091

0.720
0.726
0.582

0.707
0.547
0.451

FDR
p = 50
0.110
0.097
0.123

Power

FDR+

Power+

R2

0.943
0.920
0.907

0.061
0.061
0.076

0.858
0.837
0.774

0.707
0.547
0.451

0.876
0.831
0.766

0.139
0.118
0.141

0.788
0.660
0.599

0.702
0.552
0.451

0.831
0.749
0.691

0.173
0.126
0.156

0.740
0.594
0.550

0.702
0.543
0.451

p = 250c = 0.1
c = 0.2
c = 0.3

0.189
0.218
0.200

0.740
0.666
0.569

0.104
0.131
0.101

0.504
0.522
0.361

0.702
0.552
0.451

0.174
0.209
0.224

p = 500
c = 0.1
c = 0.2
c = 0.3

5

0.243
0.204
0.247

0.661
0.507
0.478

0.169
0.111
0.128

0.497
0.266
0.299

0.702
0.543
0.451

0.223
0.216
0.241

Empirical analysis

Our simulation results in Section 4 suggest that IPAD is a powerful approach with asymptotic
FDR control. We further examine the application of IPAD to the quarterly data on 109
macroeconomic variables from the third quarter of year 1960 (1960Q3) to the fourth quarterof year 2008 (2008Q4) in the United States discussed in [37]. These variables are transformedby taking logarithms and/or differencing following [37]. Our real data analysis consists of twoparts. In the first part, we focus on the performance of IPAD method in terms of empirical
FDR and power. In the second part, the forecasting performance of IPAD method will be
evaluated.

5.1

Simulation study

To evaluate the performance of IPAD approach in terms of empirical FDR and power withreal economic data, we set up one additional Monte Carlo simulation study. In this design,
we use the transformed macroeconomic variables described above as the design matrix X,
but simulate response y from the model in Design 1 in Section 4.1. We set the numberof true signals, the amplitude of signals, and the target FDR level to s = 10, A = 4, andq = 0.2, respectively.
Table 5 shows the results for IPAD and HD-BCKnockoff approaches. As expected, HDBCKnockoff can control FDR but suffers from lack of power. On the other hand, IPAD has
empirical FDR slightly higher than the target level (q = 0.2) while its power is reasonably
high. These results are consistent with our theory in Section 3 because IPAD only controls
19

FDR asymptotically. Additional reason for having slightly higher FDR than the target
level can be deviation of the design matrix from our factor model assumption. Overall
this simulation study indicates that IPAD can control FDR at around the target level withreasonably high power when we use the macroeconomic data set. In the next section, using
the same data set we will compare the forecasting performance of IPAD with that of somecommonly used forecasting methods in the literature.
Table 5: Real data simulation results with (n, p) = (195, 109)
FDR

Power

FDR+

Power+

R2

0.796
0.002

0.747
0.758

0.723
0.036

0.665
0.678

0.571
0.001

0.560
0.554

c = 0.2
IPAD
HD-BCKnockoff

0.278
0.096

0.812
0.009

0.223
0.010
c = 0.3

IPAD
HD-BCKnockoff

0.280
0.149

0.757
0.121

0.221
0.027
c = 0.5

IPAD
HD-BCKnockoff

5.2

0.286
0.119

0.661
0.009

0.215
0.008

Forecasting results

In this section, we apply the IPAD approach to the real economic data set for forecasting.
One-step ahead prediction is conducted using rolling window of size 120. More specifically,
one of the 109 variables is chosen as the response and the remaining 108 variables are treatedas predictors. For each quarter between 1990Q3 and 2008Q4, we use the previous 120 periods
for model fitting and then one-step ahead prediction is conducted based on the fitted model.
We compare the following different methods, where each method is implemented in a same
way as IPAD for one-step ahead prediction.
1) Autoregression of order one (AR(1)). Assume thatyt = α0 + ρyt−1 + εt ,where yt is regressed on yt−1 , and α0 and ρ are the AR(1) coefficients that need tobe estimated. With the ordinary least squares estimates α̂0 and ρ̂, the one-step ahead
prediction based on this model is ŷT +1 = α̂0 + ρ̂yT .
2) Factor augmented AR(1) (FAR). We first extract m factors f1 , · · · , fm form the 109transformed macroeconomic variables by principal component analysis (PCA). Denoteby f̃t ∈ Rm the factor vector at time t extracted from the rows of matrix [f1 , · · · , fm ] ∈Rn×m . Then we regress yt on yt−1 and f̃t−1 and fit the following modelyt = α0 + ρyt−1 + γ 0 f̃t−1 + εt20

with γ ∈ Rm . The number of factors m is determined using the P Cp1 criterion in [3].
Similar to AR(1) model, one-step ahead forecast of yt at time T isŷT +1 = α̂0 + ρ̂yT + γ̂ 0 f̃T .
3) Lasso method. The yt is regressed on yt−1 , f̃t−1 , and the 108 transformed macroeconomic variables zt−1 ∈ R108 at time t − 1
yt = α0 + ρyt−1 + γ 0 f̃t−1 + δ 0 zt−1 + εt ,where f̃t is the same as in the FAR(1) model, and α0 , ρ, and δ ∈ R108 are regression
coefficients that need to be estimated. The coefficients are estimated by Lasso methodwith regularization parameter chosen by the cross-validation. With the estimated Lasso
coefficient vector β̂Lasso , one-step ahead forecast of yt at time T is0ŷT +1 = β̂LassoxT ,where xT is the augmented predictor vector at time T .
4) IPAD method. We regress yt on the augmented vector (yt−1 , z0t−1 )0 . The lagged variableyt−1 is assumed to be always in the model. To account for this, we implement IPAD
in three steps. First, we regress yt on yt−1 and obtain the residuals ey,t . Second, we
regress each of the 108 variables in zt−1 on yt−1 and obtain the residual vector ez,t−1 .
At last, we fit model (1)–(2) using the IPAD approach by treating ey,t as the response
and ez,t−1 as predictors, which returns us a set of selected variables (a subset of the
108 macroeconomic variables). With the set of variables Sb selected by IPAD, we fit
the following model by the least-squares regression
yt = α0 + ρyt−1 + δ 0 zt−1,Sb + εt ,(19)

where zt,Sb stands for the subvector of zt corresponding to the set of variables Sb selectedby IPAD at time t. Since Sb from IPAD is random due to the randomness in generatingknockoff variables, we apply the IPAD procedure 100 times and compute the averageof these 100 one-step ahead predictions based on (19) and use the mean value as the
final predicted value of yT +1 .
Table 6 shows the root mean-squared prediction error (RMSE) of these methods. As canbe seen, the RMSE of IPAD is very close to those of comparison methods. To statisticallycompare the relative prediction accuracy of IPAD versus other approaches, we have used theDiebold–Mariano test [21], where the square of one-step ahead prediction error is used as
the loss function. Table 7 reports the test results. The results indicate that one-step ahead
prediction accuracy of IPAD is comparable to other approaches.
It is worth mentioning that one main advantage of IPAD is its interpretability and stability. Using IPAD for forecasting, we not only enjoy the same level of accuracy as othermethods but also obtain the information on variable importance with stability. Recall that21

Table 6: Root mean-squared error of one-period ahead forecast of various macroeconomic variables

RGDP
CPI-ALL
Imports
IP: cons dble
Emp: TTU
U: mean duration
HStarts: SouthNAPM new ordrs
PCED-NDUR-ENERGY
Emp. Hours
FedFunds
Cons credit
EX rate: Canada
DJIA
Consumer expect

AR

FAR

Lasso

IPAD

2.245
1.526
7.549
9.683
1.112
0.573
0.074
4.800
31.927
2.102
0.421
2.573
10.13223.117
6.496

1.929
1.552
5.871
8.353
0.989
0.487
0.071
4.378
32.121
1.899
0.396
2.537
10.13923.997
6.888

2.070
1.579
6.595
8.424
1.167
0.502
0.076
4.659
33.546
2.080
0.406
2.648
10.12224.585
6.681

2.106
1.571
6.993
9.175
1.100
0.494
0.074
4.673
32.164
1.944
0.392
2.580
10.11323.398
6.661

Table 7: Diebold–Mariano test for comparing prediction accuracy of IPAD against other procedures

RGDP
CPI-ALL
Imports
IP: cons dble
Emp: TTU
U: mean duration
HStarts: SouthNAPM new ordrs
PCED-NDUR-ENERGY
Emp. Hours
FedFunds
Cons credit
EX rate: Canada
DJIA
Consumer expect

IPAD vs. AR

IPAD vs. FAR

IPAD vs. Lasso-0.7800.521-0.976-1.026-0.140-3.383∗∗∗0.096-0.5170.753-1.200-0.9710.207-0.466
0.585
1.212

1.160
0.394
2.631∗∗
1.567
1.692
0.672
0.821
1.814
0.049
0.297
-0.1340.359-0.138-0.959-1.0380.462-0.2181.464
2.487∗-1.845-0.505-0.7660.076-1.759-2.063∗-0.625-0.661-0.037-1.428-0.27722

for each one-step ahead prediction, we apply IPAD 100 times and obtain 100 sets of selectedvariables. Thus we can calculate the selection frequency of each variable in each one-stepahead prediction. Figure 1 depicts the frequencies of top five selected variables in predicting
real GDP growth before and after year 2000, where the variable importance is ranked according to the aggregated frequencies over the entire time period before or after 2000. We haveexperimented with different cutoff years around year 2000, and the top five ranked variables
stay the same so only the results corresponding to cutoff year 2000 are reported. Changes
in index of help wanted advertising in newspapers, percentages of changes in real personal
consumption of services, and percentage of changes in real gross private domestic investment
in residential sector were the top three important variables in predicting real GDP growthduring the whole period. It is interesting to see that percentage of changes in residential
price index was among top five important variables in predicting GDP growth during the90s, and then starting from year 2000 it was replaced by changes in index of consumer expectations about stability of economy. Moreover, it is also seen that the percentage of changes
in industrial production of fuels was of great importance for predicting real GDP growthduring some periods but not the others.
As a comparison, it is very difficult to interpret the results of FAR. As for Lasso based
method, there is no theoretical guarantee on FDR control and in addition, Lasso usually
gives us models with much larger size. For instance, in predicting real GDP growth, IPAD on
average selects 5.42 macroeconomic variables while Lasso on average selects 13.32 variables.
To summarize, our real data analysis indicates that IPAD is an applicable approach forcontrolling FDR with competitive prediction power and high interpretability and stability.

6

Discussions

We have suggested in this paper a new procedure IPAD for feature selection in highdimensional linear models that achieves asymptotic FDR control while retaining high power.
Our model setting involves a latent factor model that is motivated by applications in economics and finance. Our method falls in the general model-X knockoffs framework in [15],but allows the unknown covariate distribution for the knockoff variable construction. With
the LCD knockoff statistics, we have shown that the FDR of IPAD can be asymptotically under control while the power can be asymptotically the same as that of Lasso. Our simulation
study and empirical analysis also suggest that IPAD has highly competitively performancecompared to many widely used forecasting methods such as Lasso and FAR, but with muchhigher interpretability and stability.
Our work has focused on the scenario of static models. It would be interesting to extend
the IPAD procedure to high-dimensional dynamic models with time series data. It is alsointeresting to consider nonlinear models and more flexible machine learning methods forforecasting as well as more refined factor model structures on the covariates for the knockoffs
inference with IPAD, and develop theoretical guarantees for the IPAD framework in thesemore general model settings. These extensions are beyond the scope of the current paper
and are interesting topics for future research.
23

(a) 1990-1999(b) 2000-2008

Figure 1: Frequencies of top selected variables in predicting real GDP growth. The set of
selected variables are index of help-wanted advertising in newspapers (Help wanted indx),
real personal consumption expenditures - services (Cons-Serv), real gross private domestic
investment - residential (Res.Inv), residential price index (PFI-RES), industrial production
index - fuels (IP:fuels), and University of Michigan index of consumer expectations (Consumer expect).

References[1] Ahn, S. C. and A. R. Horenstein (2013). Eigenvalue ratio test for the number of factors.
Econometrica 81, 1203–1227.
[2] Bai, J. (2003). Inferential theory for factor models of large dimensions. Econometrica 71,24

135–171.
[3] Bai, J. and S. Ng (2002). Determining the number of factors in approximate factor
models. Econometrica 70, 191–221.
[4] Barber, R. F. and E. J. Candès (2015). Controlling the false discovery rate via knockoffs.
The Annals of Statistics 43, 2055–2085.
[5] Barber, R. F. and E. J. Candès (2016). A knockoff filter for high-dimensional selective
inference. arXiv preprint arXiv:1602.03574 .
[6] Barber, R. F., E. J. Candès, and R. J. Samworth (2018). Robust inference with knockoffs.
arXiv preprint arXiv:1801.03896 .
[7] Belloni, A., V. Chernozhukov, D. Chetverikov, C. Hansen, and K. Kato (2018). Highdimensional econometrics and regularized GMM. arXiv preprint arXiv:1806.01888 .
[8] Benjamini, Y. (2010). Discovering the false discovery rate. Journal of the Royal StatisticalSociety Series B 72, 405–416.
[9] Benjamini, Y. and Y. Hochberg (1995). Controlling the false discovery rate: a practicaland powerful approach to multiple testing. Journal of the Royal Statistical Society Series
B 57, 289–300.
[10] Benjamini, Y. and D. Yekutieli (2001). The control of the false discovery rate in multiple
testing under dependency. The Annals of Statistics 29, 1165–1188.
[11] Bercu, B., B. Delyon, and E. Rio (2015). Concentration Inequalities for Sums andMartingales (1st ed.). Springer.
[12] Billingsley, P. (1995). Probability and Measure (3rd ed.). Wiley-Interscience.
[13] Bonferroni, C. E. (1935). Il calcolo delle assicurazioni su gruppi di teste. Studi in Onoredel Professore Salvatore Ortu Carboni , 13–60.
[14] Breiman, L. (2001). Random forests. Machine Learning 45, 5–32.
[15] Candès, E. J., Y. Fan, L. Janson, and J. Lv (2018). Panning for gold: ‘modelX’ knockoffs
for high dimensional controlled variable selection. Journal of the Royal Statistical Society
Series B 80, 551–577.
[16] Chernozhukov, V., D. Chetverikov, M. Demirer, E. Duflo, C. Hansen, W. Newey, andJ. Robins (2018). Double/debiased machine learning for treatment and structural parameters. Econometrics Journal 21, C1–C68.
[17] Chernozhukov, V., W. K. Härdle, C. Huang, and W. Wang (2018). Lasso-driven inference in time and space. arXiv preprint arXiv:1806.05081 .
[18] Chernozhukov, V., W. Newey, and J. Robins (2018). Double/de-biased machine learning
using regularized Riesz representers. arXiv preprint arXiv:1802.08667 .
25

[19] Chudik, A., G. Kapetanios, and H. Pesaran (2018). A one covariate at a time, multiple testing approach to variable selection in high-dimensional linear regression models.
Econometrica, to appear .
[20] De Mol, C., D. Giannone, and L. Reichlin (2008). Forecasting using a large number ofpredictors: Is Bayesian shrinkage a valid alternative to principal components? Journal of
Econometrics 146, 318–328.
[21] Diebold, F. X. and R. S. Mariano (1995). Comparing predictive accuracy. Journal ofBusiness & Economic Statistics 20, 134–144.
[22] Durrett, R. (2010). Probability: Theory and Examples (4th ed.). Cambridge University
Press.
[23] Fan, J. and Y. Fan (2008). High-dimensional classification using features annealed
independence rules. The Annals of Statistics 36, 2605–2637.
[24] Fan, J., X. Han, and W. Gu (2012). Estimating false discovery proportion under arbitrary covariance dependence (with discussion). Journal of American Statistical Association 107, 1019–1045.
[25] Fan, J. and J. Lv (2008). Sure independence screening for ultrahigh dimensional feature
space (with discussion). Journal of the Royal Statistical Society Series B 70, 849–911.
[26] Fan, Y., E. Demirkaya, G. Li, and J. Lv (2018). RANK: large-scale inference withgraphical nonlinear knockoffs. Journal of the American Statistical Association, to appear .
[27] Fan, Y., E. Demirkaya, and J. Lv (2017). Nonuniformity of p-values can occur early indiverging dimensions. arXiv preprint arXiv:1705.03604 .
[28] Fan, Y. and J. Lv (2013). Asymptotic equivalence of regularization methods in thresholded parameter space. Journal of the American Statistical Association 108, 1044–1061.
[29] Guo, Z., H. Kang, T. T. Cai, and D. S. Small (2018). Confidence intervals for causal
effects with invalid instruments by using two-stage hard thresholding with voting. Journalof the Royal Statistical Society Series B 80, 793–815.
[30] Holm, S. (1979). A simple sequentially rejective multiple test procedure. Scandinavian
Journal of Statistics 6, 65–70.
[31] Horn, R. A. and C. R. Johnson (2012). Matrix Analysis (2nd ed.). Cambridge University
Press.
[32] Lv, J. (2013). Impacts of high dimensionality in finite samples. The Annals of Statistics 41, 2236–2262.
[33] Negahban, S. N., P. Ravikumar, M. J. Wainwright, and B. Yu (2012). A unified framework for high-dimensional analysis of M -estimators with decomposable regularizers. Statistical Science 27, 538–557.
26

[34] Rigollet, P. and J.-C. Hütter (2017). High Dimensional Statistics. Massachusetts Institute of Technology, MIT Open CourseWare.
[35] Romano, J. P. and M. Wolf (2005). Exact and approximate stepdown methods for
multiple hypothesis testing. Journal of the American Statistical Association 100, 94–108.
[36] Shah, R. D. and P. Bühlmann (2018). Goodness-of-fit tests for high dimensional linear
models. Journal of the Royal Statistical Society Series B 80, 113–135.
[37] Stock, J. H. and M. W. Watson (2012). Generalized shrinkage methods for forecasting
using many predictors. Journal of Business & Economic Statistics 30, 481–493.
[38] Stucky, B. and S. van de Geer (2018). Asymptotic confidence regions for highdimensional structured sparsity. IEEE Transactions on Signal Processing 66, 2178–2190.
[39] Tibshirani, R. (1996). Regression shrinkage and selection via the lasso. J. Roy. Statist.
Soc. Ser. B 58, 267–288.
[40] Vershynin, R. (2012). Introduction to the non-asymptotic analysis of random matrices.
In Y. C. Eldar and G. Kutyniok (Eds.), Compressed Sensing: Theory and Practice, pp.
210–268. Cambridge University Press.
[41] Vizcarra, A. B. and F. G. Viens (2007). Some applications of the Malliavin calculus
to sub-Gaussian and non-sub-Gaussian random fields. In R. C. Dalang, M. Dozzi, andF. Russo (Eds.), Seminar on Stochastic Analysis, Random Fields and Applications V, pp.
363–395. Springer Science & Business Media.
[42] Wooldridge, J. M. and Y. Zhu (2018). Inference in approximately sparse correlated
random effects probit models. Journal of Bussiness & Economic Statistics, to appear .
[43] Zhang, X. and G. Cheng (2017). Simultaneous inference for high-dimensional linear
models. Journal of the American Statistical Association 112, 757–768.

27

Appendix
This appendix contains all the proofs and technical details for the theoretical results ofthe paper. In particular, Section A details the proofs of Lemmas 1–2 and Theorems 1–2,
Section B presents some key lemmas and their proofs, and Section C provides some additional
technical lemmas and their proofs.
To ease the technical presentation, let us introduce some notation. We denote by .
e θ̂)the inequality up to some positive constant factor. Restricting the columns of X and X(
to the variables in index set A such that |A| ≤ k, we obtain the n × k submatrices XA
e A (θ̂), respectively. Moreover, we define TA (θ̂) := vec(vech UA (θ̂), vA (θ̂)) ∈ Rk(2k+3)
and X
with UA (θ̂) the principle submatrix of U(θ̂) formed by columns and rows in set A, andvA (θ̂) the subvector of v(θ̂) formed by components in set A. Then it is easy to see thate A (θ̂)]0 [XA , X
e A (θ̂)] and vA (θ̂) = n−1 [XA , X
e A (θ̂)]0 y. For the oracle
UA (θ̂) = n−1 [XA , X
factor loading matrix Λ0 , with a slight abuse of notation we use Λ0A to denote the row
restricted to the variables in A for notational convenience. Recall that ν > 0 is a fixed
positive number, cnp = (p−1 log n)1/2 + (n−1 log p)1/2 , and c̃np = p−1/2 log n + n−1/2 log p. We
define πnp = n−ν + p−ν . Since λ is fixed at C0 n−1/2 log p, in all the proofs we will drop the
dependence of various quantities on λ whenever there is no confusion.

A

Proofs of main results

A.1

Proof of Lemma 1e
For λ fixed at C0 n−1/2 log p and each given θ, Wj (θ) = wj ([X, X(θ)],y) depends only onaugaugβ̂ (θ) by the LCD construction. Moreover, the Lasso solution β̂ (θ) satisfies the Karush–
Kuhn–Tucker (KKT) conditions:v(θ) − U(θ)β̂ aug (θ) = n−1 λz,

(A.1)


{sgn(β̂ )} if β̂ 6= 0,jjwhere z = (z1 , · · · , z2p )T with zj ∈[−1, 1]if β̂j = 0,

for j = 1, . . . , 2p.

(A.2)

e
This means that β̂ aug (θ) depends on the data ([X, X(θ)],y)) only through U(θ) and v(θ).
Thus using notation T(θ) = vec(vech U(θ), v(θ)) with the fact that U(θ) is symmetric,e
we can reparametrize wj ([X, X(θ)],y) as wj (T(θ)) with a slight abuse of notation. Furthermore, note that the thresholds T1 and T2 are both completely determined by wj (T(θ)).
Consequently, by the construction of Sb we can see that Sb depends only on T(θ), which
completes the proof of Lemma 1.

A.2

Proof of Lemma 2

We continue to use the same λ and θ as in Lemma 1 and its proof. Recall that SA (tA )represents the outcome of first restricting ourselves to the smaller set of variables A and28then applying IPAD to TA (θ) = tA to further select variables from A. Also recall thatA∗ (θ) is the support of knockoff statistics Wj (θ). Thus the knockoff threshold T1 or T2
depends only on Wj (θ) with j ∈ A∗ (θ).
On the other hand, when we restrict ourselves to A ⊃ A∗ (θ) we solve the following KKT
conditions with respect to β̃ := (β̃1 , · · · , β̃2|A| )T ∈ R2|A| to get the Lasso solution:β̃ = (UA (θ)0 UA (θ))−1 (vA (θ) − n−1 λz̃),{sgn(β̃ )} if β̃ 6= 0,jjwhere z̃ = (z̃1 , · · · , z̃2|A| )T with z̃j ∈[−1, 1]if β̃j = 0,

(A.3)
for j = 1, . . . , 2|A|. (A.4)

Since λ is always fixed at the same value C0 n−1/2 log p, it is seen that the solution to theaugabove KKT conditions is identical to β̂AA(θ), where the latter denotes the subvector of
aug
aug
aug
(θ), j2 ∈ A all together. Therefore,β̂ (θ) formed by stacking β̂j1 (θ), j1 ∈ A and β̂p+j2
the Lasso solution to (A.3)–(A.4) and the Lasso solution to (A.1)–(A.2) have the identicalsupport (when viewed in the original 2p-dimensional space) and in addition, identical values
on the support. This guarantees that S{1,...,p} (T(θ)) and SA (TA (θ)) are identical and thus
concludes the proof of Lemma 2.

A.3

Proof of Theorem 1

Recall that for a given θ, A∗ (θ) is the support of knockoff statistics (W1 (θ), · · · , Wp (θ))0 .
Define setb θ̂) := A∗ (θ̂) ∪ A∗ (θ 0 ).
A(b θ̂) is bounded by k. Hereafter we write A(b θ̂)
It follows from (14) that the cardinality of A(
as Ab for notational simplicity.
By Lemmas 1–2 and the definition of the FDP, we know that S{1,...,p} (T(θ̂)) = SAb(TAb(θ̂))
and thus the resulting FDR’s are the same. Therefore, we can restrict ourselves to the smaller
model Ab when studying the FDR of IPAD. The same arguments as above also hold for theoracle knockoffs; that is, the FDR of IPAD applied to T(θ 0 ) is the same as that applied
to TAb(θ 0 ). Note that all the FDR’s we discuss here are with respect to the full model{1, · · · , p}. For this reason, in what follows we will abuse the notation and use FDRAb(TAb(θ̂))
and FDRAb(TAb(θ 0 )) to denote the FDR of IPAD based on TAb(θ) and TAb(θ 0 ), respectively.
We want to emphasize that although we put a subscript Ab in FDR’s, their values are still
deterministic as argued above. Summarizing the facts, we obtain
FDRAb(TAb(θ̂)) = FDR{1,··· ,p} (T(θ̂)),
FDRAb(TAb(θ 0 )) = FDR{1,··· ,p} (T(θ 0 )).
e 0 ) satisfies the two properties in Definition 1 and is a
Meanwhile, by construction X(θvalid model-X knockoffs matrix. Therefore, for any value of the regularization parameter,e 0 )], y) together with Result 1 ensure the exact
the LCD statistics Wj (θ 0 ) based on ([X, X(θ29FDR control at some target level q ∈ (0, 1). Summarizing this, we obtain that the FDR of
IPAD applied to T(θ 0 ) is controlled at target level q.
Combining the arguments in the previous two paragraphs, we deduce
FDRAb(TAb(θ 0 )) = FDR{1,··· ,p} (T(θ 0 )) ≤ q.
Thus the desired results follow automatically if we can prove that FDRAb(TAb(θ̂)) is asymptotically close to FDRAb(TAb(θ 0 )). We next proceed to prove it.
Recall the definitions of I and IA as in (15). Define the eventno 
Enp = TAb(θ̂) ∈ IAb ∩ TAb(θ 0 ) ∈ IAb .
Lemma 3 in Section B.1 establishes θ̂ ∈ Θnp with probability at least 1 − O(πnp ) andθ 0 ∈ Θnp . Hence, Lemma 4 in Section B.2 guarantees that
!

c
(A.5)
P Enp
≤ 2P
sup
TA (θ) − E[TA (θ 0 )] max > anp = O(πnp ),|A|≤k, θ∈Θnp

where anp = C1 (k 1/2 + s3/2 )c̃np for some constant C1 > 0.
For a given deterministic set A ⊂ {1, · · · , p}, let FDPA (·) be the FDP function corresponding to FDRA (·). By the definition of FDP function, we have for any t1 , t2 ∈ R|A|(2|A|+3) ,|S 1 ∩ SA (t2 )| |S 1 ∩ SA (t1 )|−|SA (t2 )||SA (t1 )|1|S ∩ SA (t2 )| · (|SA (t1 )| − |SA (t2 )|) |S 1 ∩ SA (t2 )| − |S 1 ∩ SA (t1 )|=+.
|SA (t1 )| · |SA (t2 )||SA (t1 )|FDPA (t2 ) − FDPA (t1 ) =Further, note that|S 1 ∩ SA (t2 )|/|SA (t2 )| ≤ 1,||SA (t2 )| − |SA (t1 )|| ≤ |SA (t2 )4SA (t1 )| ,|S ∩ SA (t2 )| − |S ∩ SA (t1 )| ≤ {SA (t2 )4SA (t1 )} ∩ S 1 .
1

1Combining the results above yields| FDPA (t1 ) − FDPA (t2 )|≤|SA (t1 )| − |SA (t2 )|{SA (t2 )4SA (t1 )} ∩ S 1
SA (t2 )4SA (t1 )+≤2.
|SA (t1 )||SA (t1 )||SA (t1 )|Similarly we have
| FDPA (t1 ) − FDPA (t2 )| ≤ 2

SA (t2 )4SA (t1 ).
|SA (t2 )|Thus it holds thatsupSA (t2 )4SA (t1 )|A|≤k t1 ,t2 ∈IA |SA (t1 )| ∧ |SA (t2 )|sup | FDPA (t1 ) − FDPA (t2 )| ≤ sup|A|≤k t1 ,t2 ∈IA

sup= O(ρnp ),
30

(A.6)

where the last two steps are due to Condition 6. Therefore, (A.5) and (A.6) together withthe fact that FDP(·) ∈ [0, 1] entail that
FDRAb(TAb(θ̂)) − FDRAb(TAb(θ 0 )) = E FDPAb(TAb(θ̂)) − E FDPAb(TAb(θ 0 ))≤ E FDPAb(TAb(θ̂)) − FDPAb(TAb(θ 0 ))ihc
≤ E FDPAb(TAb(θ̂)) − FDPAb(TAb(θ 0 )) | Enp P (Enp ) + 2 P Enp≤ sup

sup |FDPA (t1 ) − FDPA (t2 )| + O(πnp )|A|≤k t1 ,t2 ∈IA= O(ρnp ) + O(πnp ).
This completes the proof of Theorem 1.

A.4

Proof of Theorem 2

By the definition of the LCD statistics, we construct the augmented Lasso estimator for eachθ ∈ Θnp , which is defined asβ̂ aug (θ) = arg min

b∈R2p

e
y − [X, X(θ)]b2
2+ λkbk1 .

(A.7)

The Lasso estimator of regressing y on only X is also given byβ̂ = arg minp ky − Xbk22 + λkbk1 ,
b∈R

(A.8)

where λ = O(n−1/2 log p). According to the true model S 0 , the underlying true parameter vector corresponding to β̂ aug (θ) should be given by β aug := (β 0 , 00 )0 ∈ R2p with
β = (βS0 0 , 00 )0 ∈ Rp and |S 0 | = s for any θ ∈ Θnp . By Lemma 5 in Section B.3, withprobability at least 1 − O(πnp ) the Lasso estimators satisfyβ̂ aug (θ) − β aug

sup
θ∈Θnp

β̂ − β11= O(sλ),= O(sλ),where λ = O(n−1/2 log p).
We now prove that under Condition 7, the power of the augmented Lasso (A.7) is boundedfrom below by γ ∈ [0, 1]; that is,
E Sbauglasso ∩ S 0 /s ≥ γ,

(A.9)

where Sbauglasso = {j : β̂jaug (θ) 6= 0}. To this end, we first show that with asymptotic
probability one,c
|Sbauglasso∩ S 0 |/s ≤ 1 − γ.

31

(A.10)

c
The key is to use proof by contradiction. Suppose |Sbauglasso∩ S 0 |/s > 1 − γ. Then we cansee that

sup
θ∈Θnp

β̂ aug (θ) − β aug

1

≥ sup
θ∈Θnp

= β augbc

β̂ augbc

Sauglasso

Sauglasso1(θ) − β aug
bcSauglasso≥ β augbc

Sauglasso ∩S 0

1

> bnp sn−1/2 log p,1

where the last step is by Condition 7. However, by Lemma 5 with probability at least
1−O(πnp ), the left hand side above is bounded from above by O(sλ) with λ = O(n−1/2 log p).
These two results contradict with each other since bnp → ∞. Hence (A.10) is proved.
cTherefore, the result in (A.9) follows immediately since |Sbauglasso ∩ S 0 | = s − |Sbauglasso∩ S 0|andE Sbauglasso ∩ S 0 /s ≥ γ P |Sbauglasso ∩ S 0 |/s > γc= γ P |Sbauglasso∩ S 0 |/s ≤ 1 − γ = γ(1 − O(πnp )).
Using the same argument, we can show that the power of the Lasso (A.8) is also boundedfrom below by γ(1 − O(πnp )) under Condition 7. That is, we have
E Sblasso ∩ S 0 /s ≥ γ(1 − O(πnp )),where Sblasso = {j : β̂j 6= 0}.
Next we show that our knockoffs procedure has at least the same power as the augmented
Lasso and hence the Lasso itself. Namely, we prove
E Sb ∩ S 0 /s ≥ γ(A.11)

with threshold T2 . Note that the same argument is still valid for T1 . Let |W(1) | ≥ · · · ≥ |W(p) |and define j ∗ as |W(j ∗ ) | = T2 . Then by the definition of T2 , it holds that −T2 < Wj ∗ +1 ≤ 0.
Here we have assumed that there are no ties on the magnitudes of Wj ’s which should be areasonable assumption considering the continuity of the Lasso solution. As in the proof of
Theorem 3 in [26], it is sufficient to consider the following two cases.
Case 1. Consider the case of −T2 < W(j ∗ +1) < 0. In this case, from the definition ofthreshold T2 we have
2 + |{j : W(j) ≤ −T2 }|> q.
|{j : W(j) ≥ T2 }|Using the same argument as in Lemma 6 of [26] together with Lemma 5, we can prove fromb ≥ C2 s with probability at least 1 − O(πnp ). This leads to |{j : W(j) ≤Condition 8 that |S|−T2 }| > C2 qs − 2 with the same probability. Now from the same argument as in A.5 of [26],32we can obtain T2 = O(λ). On the other hand, Lemma 5 and some algebra establish thatO(sλ) = kβ̂aug(θ̂) − β

aug

k1 =

p
X

|β̂jaug (θ̂) − βj | +j=1|β̂jaug (θ̂) − βj | +

X

=

p
X

aug
|β̂j+p
(θ̂)|

j=1

X

|β̂jaug (θ̂)|

j∈S 1

b 0
j∈S∩S

X

+

|β̂jaug (θ̂) − βj | +

p
X

aug
|β̂j+p(θ̂)|.

(A.12)

j=1

j∈Sbc ∩S 0

We then consider the lower bound of the last term in (A.12). For any j ∈ Sbc , it holds thataug|β̂j+p(θ̂)| > |β̂jaug (θ̂)| − T2 . Hence we obtainpXaug|β̂j+p(θ̂)| ≥

j=1

X

aug
|β̂j+p(θ̂)| +≥aug|β̂j+p(θ̂)|

j∈Sbc ∩S 0

b 0
j∈S∩S

X

X

aug
|β̂j+p(θ̂)| +

X|β̂jaug (θ̂)| − T2 |Sbc ∩ S 0 |.

(A.13)

j∈Sbc ∩S 0

b 0
j∈S∩S

Plugging (A.13) into (A.12) and applying the triangle inequality yield
X
X aug
X
O(sλ) ≥
|β̂jaug (θ̂) − βj | +|β̂j (θ̂)| +|β̂jaug (θ̂) − βj |
j∈S 1

b 0
j∈S∩S

+

X

aug
|β̂j+p(θ̂)| +

b 0
j∈S∩S

≥

X

|β̂jaug (θ̂)| − T2 |Sbc ∩ S 0 |j∈Sbc ∩S 0|β̂jaug (θ̂) − βj | +X|β̂jaug (θ̂)|

j∈S 1

b 0
j∈S∩S

+

X
j∈Sbc ∩S 0≥

j∈Sbc ∩S 0X

X

|βj | +Xaug
|β̂j+p(θ̂)| − T2 |Sbc ∩ S 0 |

b 0
j∈S∩S

|βj | − T2 |Sbc ∩ S 0 | = kβSbc ∩S 0 k1 − T2 |Sbc ∩ S 0 |.

j∈Sbc ∩S 0

Since T2 |Sbc ∩S 0 | = O(sλ) for λ = O(n−1/2 log p) due to the discussion above, we consequently
obtainkβSbc ∩S 0 k1 = O(sn−1/2 log p).

(A.14)

Suppose |Sbc ∩ S 0 |/s > 1 − γ. Then Condition 7 gives kβSbc ∩S 0 k1 > bnp sn−1/2 log p for somepositive diverging sequence bnp ; this contradicts with (A.14). Thus we obtain |Sbc ∩ S 0 |/s ≤1 − γ with asymptotic probability one, which leads to (A.11) by taking expectation.
Case 2. Consider the case of W(j ∗ +1) = 0. In this case, by the definition of threshold T2
1 + |{j : W(j) < 0}|≤ q.
|{j : W(j) > 0}|(A.15)

If |{j : W(j) < 0}| > C3 s for some constant C3 > 0, then from the same argument as in A.5
of [26], we can obtain T2 = O(λ), and the rest of the proof is the same as in Case 1. On the33other hand, if |{j : W(j) < 0}| ≤ o(s) we have|{j : W(j) 6= 0} ∩ S 0 | = |{j : W(j) > 0} ∩ S 0 | + |{j : W(j) < 0} ∩ S 0 |≤ |Sb ∩ S 0 | + o(s).
Now note that |{j : W(j) 6= 0}| ≥ |{j : |β̂jaug | =6 0, j = 1, · · · , p}|. Then we can see that withasymptotic probability one,|{j : W(j) 6= 0} ∩ S 0 | ≥ |{j : β̂jaug 6= 0, j = 1, · · · , p} ∩ S 0 |= |Sbauglasso ∩ S 0 |.
≥ γs(1 − o(1)).
Consequently, we obtain |Sb∩S 0 |/s ≥ γ(1−o(1)), which leads to (A.11) by taking expectation.
Combining these two cases concludes the proof of Theorem 2.

B

Some key lemmas and their proofs

B.1

Lemma 3 and its proof

Lemma 3 Assume that Conditions 2–5 hold. Then with probability at least 1 − O(πnp ), theb 0 , η̂ 0 )0 lies in the shrinking set given by
estimator θ̂ = (vec(C)Θnp = θ = (vec(C)0 , η 0 )0 : C − C0 max + η − η 0 max ≤ O(cnp ) ,where cnp = (n−1 log p)1/2 + (p−1 log n)1/2 and πnp = p−ν + n−ν .
b − C0 kmax in Part 1
Proof. We divide the proof into two parts. We prove the bound for kC
and then for kη̂ − η 0 kmax in Part 2.
b − C0 kmax = maxi,j |ĉij − c0 |, where the maximum is taken overPart 1. Note that kC
iji ∈ {1, . . . , n} and j ∈ {1, . . . , p}. We write fi∗ = H0 fi0 and λ∗j = H−1 λ0j with rotation matrix
H defined in Lemma 6 in Section C.1. From the definition of cij , it holds thatĉij − c0ij = (f̂i − fi∗ )0 λ∗j + f̂i0 (λ̂j − λ∗j ).
From Lemma 6, we can assume kHk2 + kH−1 k2 + kVk2 + kV−1 k2 . 1, which occurs withprobability at least 1 − O(p−ν ). We also have maxi∈{1,...,n} kf̂i k22 . 1 a.s. by the assumed
restriction F̂0 F̂/n = Ir as mentioned on p.213 of [3]. Hence, the triangle and Cauchy–Schwarz
inequalities with Conditions 2 and 3 givemax |ĉij − cij | ≤ max kf̂i − fi∗ k2 max kλ∗j k2 + max kf̂i k2 max kλ̂j − λ∗j k2
i,j

i

j

i

. max kf̂i − fi∗ k2 + max kλ̂j − λ∗j k2 .
i
j

j

(B.1)

Then it is sufficient to derive upper bounds for maxi kf̂i −fi∗ k2 and maxj kλ̂j −λ∗j k2 that holdwith high probability. Using the decomposition of A.1 in [2] along with taking maximum34

over i, ` ∈ {1, . . . , n}, we can deduce
max kf̂i − fi∗ k2i≤ kV−1 k2 max (σe2 /n)kf̂i k2 + n−1
i

+ n−1

n
X

0

kf̂` f`0 k2 p−1

p
X

i,`

p
X

(e`j eij − E[e`j eij ])

j=1

+ n−1

λ0j eijj=1. O(n−1 ) + max p−1
. O(n

p
X

kf̂` k2 p−1

`=1

`=1

−1

n
X

n
X

0

kf̂` fi0 k2 p−1(e`j eij − E[e`j eij ]) + max p−1
i

j=1


λ0j e`j

j=1

`=1

2

p
X

p
X


2

λ0j eij

j=1

2) + R1 + R2 ,(B.2)

where we have used the boundedness of kf̂` k2 discussed above and kf`0 k2 ≤ r1/2 kf`0 kmax . 1 inPCondition 2 for the second inequality, and defined R1 = maxi,` p−1 pj=1 (e`j eij − E[e`j eij ])Pand R2 = maxi,k p−1 pj=1 λ0jk eij . Similarly, the expression on p.165 of [2] with taking
maximum over i ∈ {1, . . . , n} and j ∈ {1, . . . , p} leads to
max kλ̂j − λ∗j k2
j≤ kHk2 max n−1
j+ max n−1j. max n
j

−1

n
X

n
X

+ n−1

fi0 eij

i=1
n
X

n
X
i=1

2

H−1 2 max λ0j 2

f̂i (f̂i − fi∗ )0j

2(f̂i − fi∗ )eij

i=1

2

fi0 eijf̂i − fi∗

i=1

+ max
i

2

+ max
i

f̂i − fi∗

2

max n
j

2

= R3 + max kf̂i − fi∗ k2 (1 + R4 ),
i

−1

n
X

!1/2
e2ij

i=1

(B.3)


P
Pn 2 1/2
0e−1
where R3 = maxj,k n−1 ni=1 fik
and
R
=
max
n
, and the Cauchy–ij 24ji=1 eij
Schwarz inequality has been used to obtain the second inequality. To evaluate R4 , we notethat
R42 ≤ max E e2ij + max n−1
j
j

n
X

e2ij − E e2ij



.

i=1

The first term is bounded by 2Ce2 . For the second term, Lemma 7(a) in Section C.2 with p
replaced by n and the union bound give
!
!
n
n
X
X


−1
2
2
−1
2
2
eij − E eij > u
eij − E eij > u ≤ p max P n
P max n
j

j

i=1

i=1

2

≤ 2p exp(−nu /C)for all 0 ≤ u ≤ c. Thus putting u = (C(ν + 1)n−1 log p)1/2 and using condition cnp ≤c/(r2 M 2 C(ν + 2))1/2 , we obtain R42 = O(1) + O((n−1 log p)1/2 ) = O(1) with probability at
35

least 1 − O(p−ν ). This together with the observation from (B.1)–(B.3) yields
max |ĉij − c0ij | . R3 + R1 + R2 + O(n−1 ) (1 + R4 )i,j

. R1 + R2 + R3 + O(n−1 ).
Hence the convergence rate of maxi,j |ĉij − c0ij | is determined by the slowest term out of R1 ,R2 , R3 , and O(n−1 ). We evaluate these terms by Lemma 7 in Section C.2 and the union
bound with condition cnp ≤ c/(r2 M 2 C(ν + 2))1/2 as above. First for R1 , Lemma 7(a) byletting u1 = (C(ν + 2)p−1 log n)1/2 results in
P (R1 > u1 ) ≤ 2n2 exp −p(ν + 2)p−1 log n = O(n−ν ).
Next for R2 , Lemma 7(c) with u2 = (2(ν + 1)p−1 log n)1/2 givesP (R2 > u2 ) ≤ 2rn exp −p(ν + 1)p−1 log n = O(n−ν ).
Finally for R3 , Lemma 7(b) with putting u3 = (C(ν + 1)n−1 log p)1/2 leads toP (R3 > u3 ) ≤ 2rp exp −n(ν + 1)n−1 log p = O(p−ν ).
b − C0 kmax = O(cnp ), which holds with probability
Consequently, we obtain the first result kC
at least 1 − O(πnp ).
Part 2. Next we derive the convergence rate of η̂. It is sufficient to prove only thecase when η 0 is a scalar (so that we write η 0 = η10 ) since dimensionality m is fixed and ηk0 ’s
Pshare the identical property thanks to Condition 4. Recall notation Enp ek = (np)−1 i,j ekij .
Letting δij = c0ij −ĉij , we have êij = xij −ĉij = eij +δij . For an arbitrary fixed k ∈ {1, . . . , m},
the binomial expansion entails
Enp êk − E ek = Enp (e + δ)k − E ek
k−1  Xk ` k−`= Enp (e − E e ) + Enpeδ``=0
k−1  Xk
≤ Enp (ek − E ek ) +max |δij |k−` Enp |e|`
` i,j
`=0

X
k−1
k
k
. Enp (e − E e ) + O max |δij |Enp |e|` .
k

k

i,j

(B.4)

`=0

For all k ∈ {1, . . . , m}, the strong law of large numbers with Theorem 2.5.7 in [22] entails|Enp ek − E ek | = o((np)−1/2 log(np)) a.s. under Condition 4. Furthermore, the second termof (B.4) is O(cnp ) with probability at least 1 − O(πnp ) from Part 1 and the same law of large
numbers. Consequently, we obtain
Enp êk − E ek . cnp .

36

Therefore by the construction of η̂1 and local Lipschitz continuity of h1 in Condition 4, we
see thatη̂1 − η10 = |h1 (Enp ê, . . . , Enp êm ) − h1 (E e, . . . , E em )|.

max
k∈{1,...,m}Enp êk − E ek

with probability at least 1 − O(πnp ). This completes the proof of Lemma 3.

B.2

Lemma 4 and its proof

Lemma 4 Assume that Conditions 1–4 hold. Then with probability at least 1 − O(πnp ), thefollowing statements holdUA (θ) − E[UA (θ 0 )] max = O k 1/2 c̃np ,(a)sup|A|≤k, θ∈Θnp

(b)

sup
|A|≤k, θ∈ΘnpvA (θ) − E[vA (θ 0 )] max = O s3/2 c̃np ,where Θnp was defined in Lemma 3 and c̃np = n−1/2 log p + p−1/2 log n. Consequently, we
have



sup
TA (θ) − E[TA (θ 0 )] max = O k 1/2 + s3/2 c̃np .
|A|≤k, θ∈Θnp

Proof. To complete the proof of (a), we verify the following
(a–i)

sup
|A|≤k, θ∈Θnp

UA (θ) − UA (θ 0 ) max . k 1/2 c̃np ,U(θ 0 ) − E[U(θ 0 )] max . (n−1 log p)1/2 .

(a–ii)

From (a–i) and (a–ii), we can conclude thatsup|A|≤k, θ∈Θnp≤

UA (θ) − E[UA (θ 0 )] max
sup|A|≤k, θ∈Θnp≤

sup|A|≤k, θ∈Θnp

UA (θ) − UA (θ 0 ) max + sup|A|≤k

UA (θ 0 ) − E[UA (θ 0 )] max

UA (θ) − UA (θ 0 ) max + U(θ 0 ) − E[U(θ 0 )] max. k 1/2 c̃np ,which yields result (a).
We begin with showing (a–i); this is the uniform extension of Lemma 8(a) in Section
C.3 over |A| ≤ k. In fact, the proof is almost the same, with the only difference that bound(B.15) should be replaced with the bound derived in Lemma 9(c); that is,
max n−1/2 EA

|A|≤k

2

. 1 ∨ kn−1 log p

1/2

,

(B.5)

1/2
which holds with probability at least 1 − O(p−ν ). Notice that kn−1 log p≤ log1/2 p.
Therefore, even if we use (B.5) instead of (B.15) in the proof of Lemma 8(a) we can still37derive the same convergence rate k 1/2 c̃np as in Lemma 8(a), and hence (a–i) holds withprobability at least 1 − O(πnp ).
For (a–ii), we see thatU(θ 0 ) − E[U(θ 0 )] max ≤ n−1 X0 X − E[n−1 X0 X] maxe 0 )0 X(θe 0 ) − E[n−1 X(θe 0 )0 X(θe 0 )]+ n−1 X(θ
+2 n

−1

0e

0

X X(θ ) − E[n

−1

0e

0

X X(θ )]max

max=: W1 + W2 + 2W3 .

(B.6)

We derive the bounds for each of these terms. First, W1 is bounded as
00

W1 ≤ n−1 C0 C0 − E[n−1 C0 C0 ]max

+ n−1 E0 E − E n−1 E0 E max + 2 n−1 E0 C0 max=: W1,1 + W1,2 + W1,3 .
Under Condition 3, we deduce
W1,1 =rX

maxj,`∈{1,...,p}≤ rM 2

λ0jk λ0`m n−1

0 0
0 0
fik
fim − E fik
fim



i=1

k,m=1

n−1

max

n
X

j,`∈{1,...,p}

n
X

0 0
0 0
fim
fim − E fik
fik



.

i=1

From Lemma 7(d) with Condition 2 and the union bound, we have
!
n
X

−1
0 0
0 0
P
max
n
fik fim − E fik fim > uj,`∈{1,...,p}

≤ p2

i=1

max
j,`∈{1,...,p}

n−1

P

n
X

!

0

0 0
0
fik
fim − E fik
fim

>u


≤ 2p2 exp −nu2 /C .

i=1

Hence, letting u = (C(ν + 2)n−1 log p)1/2 above yields the bound W1,1 . (n−1 log p)1/2 withprobability at least 1 − O(p−ν ). Next for W1,2 , we can find from Lemma 7(a) with p replacedby n and the union bound that!
nX
P n−1 E0 E − E n−1 E0 E max > u ≤ p2 max P n−1(eij ei` − E eij ei` ) > u
j,`

2

2

i=1



≤ 2p exp −nu /C .
Letting u = (C(ν + 2)n−1 log p)1/2 and using n−1 log p ≤ c2 /(C(ν + 2)), we obtain W1,2 .
(n−1 log p)1/2 with probability at least 1 − O(p−ν ). Next for W1,3 , the union bound gives
!
r X
n


X
0
0 0
P n−1 E0 F0 Λ0
>u =P
max
n−1
eij fikλ`k > u
max

≤P r
≤ rp

max

j,`∈{1,...,p}

max

n

−1

n
X

−1

n
X

j,`∈{1,...,p} k∈{1,...,r}

max

max P

k∈{1,...,r} j∈{1,...,p}k=1 i=1!
0
eij fik

λ0`k

>u

i=1

n

i=1

38

!
0
eij fik> u/(rM ) .

Lemma 7(b) states that for all 0 ≤ u/(rM ) ≤ c/(rM ) it holds that
!
n
X

−1
0
P n
eij fik > u/(rM ) ≤ 2 exp −nu2 /(Cr2 M 2 ) .
i=1

Therefore, if we put u = rM (C(ν + 1)n−1 log p)1/2 using n−1 log p ≤ c2 /(r2 M 2 C(ν + 1)),the upper bound of the probability is further bounded by 2rp−ν . Thus we obtain W13 .
(n−1 log p)1/2 with probability at least 1 − O(p−ν ). Consequently, the bound of W1 is
W1 ≤ W1,1 + W1,2 + W1,3 . (n−1 log p)1/2
with probability at least 1 − O(p−ν ). Note that we have the same result for W2 since it has
the same distribution as W1 . Finally, W3 is bounded as
00W3 ≤ n−1 C0 C0 − E[n−1 C0 C0 ]max

+ n−1 E0 C0 max + n−1 E0η0 C0+ n−1 E0 Eη0 max
max=: W1,1 + W3,1 + W1,3 + W3,2 .

The upper bound of W3,1 turns out to be O((n−1 log p)1/2 ) that holds with probability atleast 1 − O(p−ν ). We check this claim. Using the union bound and the inequality of Lemma
7(a) with p replaced by n and putting u = (C(ν + 2)n−1 log p)1/2 yield
!
n
X


P n−1 E0 Eη0 max > u ≤ p2 max P n−1
eij eη0 ,i` > u ≤ 2p−ν .
j,`

i=1

Finally, W3,2 is found to have the same bound as W1,3 because Eη0 is an independent copyof E. Consequently, with probability at least 1 − O(p−ν ), we obtainU(θ 0 ) − E[U(θ 0 )] max . (n−1 log p)1/2 .
This completes the proof of (a) since p−ν /πnp = O(1).
Next we show (b) by verifying the following
(b–i)

sup
|A|≤k, θ∈Θnp(b–ii)vA (θ) − vA (θ 0 ) max . s3/2 c̃np ,v(θ 0 ) − E[v(θ 0 )] max . s(n−1 log p)1/2 .

Similar to the proof of (a), we need to modify the proof of Lemma 8(b) in Section C.3 forobtaining the uniform bound with respect to A, but the obtained result is already uniform
over the choice of A. Thus the same upper bound holds and (b–i) follows. Next we show(b–ii). It holds thatv(θ 0 ) − E v(θ 0 ) max
e 0 )0 y − E n−1 X(θe 0 )0 y≤ n−1 X0 y − E n−1 X0 y max + n−1 X(θ
max

−1 0
−1 0
−1 0
−1 0
≤ n X X − E n X X β max + n X ε − E n X ε maxe 0 )0 X − E[n−1 X(θe 0 )0 X] βe 0 )0 ε − E[n−1 X(θe 0 )0 ε]+ n−1 X(θ+ n−1 X(θmax=: Z1 + Z2 + Z3 + Z4 .
39

max

These terms can be bounded by the results obtained in the proof of (a–ii). We see that
Z1 ≤ s1/2 n−1 X0 XS 0 − E n−1 X0 XS 0 max kβS 0 k2 . sW1 . s(n−1 log p)1/2
with probability at least 1 − O(p−ν ). Next we deduce0Z2 ≤ n−1 Λ0 F0 εmax

+ n−1 E0 ε max .

The first and second terms can be bounded by the same ways as W1,3 and W3,1 in the proofof (a) above with E and Eη0 replaced by ε, respectively. Then the first term dominates thesecond and hence Z2 . (n−1 log p)1/2 with probability at least 1 − O(p−ν ). Similarly, we canobtaine 0 )0 XS 0 − E n−1 X(θe 0 )0 XS 0
Z3 ≤ s1/2 n−1 X(θmax

kβS 0 k2 . sW3 . s(n−1 log p)1/2

with probability at least 1 − O(p−ν ). Note that Z4 has the same bound as Z2 . Consequently,
collecting terms leads to the result, Z1 + · · · + Z4 . s(n−1 log p)1/2 with probability at least
1 − O(p−ν ). This proves (b–ii) and concludes the proof of Lemma 4.

B.3

Lemma 5 and its proof

Lemma 5 Assume that all the conditions of Theorem 2 hold. Then with probability at least
1 − O(πnp ), the Lasso solution in (A.7) satisfiesβ̂ aug (θ) − β aug

sup
θ∈Θnp

β̂ aug (θ) − β aug

sup
θ∈Θnp

2

1

= O(s1/2 λ),= O(sλ),where λ = c1 n1/2 log p with c1 some positive constant.
Proof. Let δ(:= δ(θ)) := β̂ aug (θ) − β aug . We start with introducing two inequalities
sup0en−1 [X, X(θ)]
ε

θ∈Θnp

inf

θ∈Θnp , δ∈Vmax≤ 2−1 λ,(B.7)

δ 0 U(θ)δ/kδk22 ≥ σe2 (1 + o(1)),(B.8)

where λ = c1 n−1/2 log p for some positive constant c1 and
V = δ ∈ R2p : kδS 1 k1 ≤ 3kδS 0 k1 , kδk0 ≤ k .

(B.9)

It is well known that the rate of convergence of the Lasso estimator can be obtained provided
that (B.7) and (B.8) hold. Thus we show that these two inequalities actually hold with highprobability in Step 1, and then derive the convergence rate using (B.7) and (B.8) in Step 2.
Step 1. We check whether (B.7) and (B.8) actually hold with high probability. We first
verify (B.7). By the proofs of Lemmas 8 and 4, we have0en−1 [X, X(θ)]
ε

sup
θ∈Θnp

max

≤ n−1 X0 ε max + sup

0ee 0 )0 εn−1 X(θ)ε − n−1 X(θ

θ∈Θnp

40

max

e 0 )0 ε+ n−1 X(θmax.

The first and third terms can both be upper bounded by O(n−1/2 log p) with probability atleast 1 − O(p−ν ), following the same lines for deriving bound for Z2 in the proof of Lemma 4.
To evaluate the second term, we can use the argument about V2 and its upper bound (B.16)
in the proof of Lemma 8. That bound still holds with the same rate O(n−1/2 log p) even if
we take A = {1, . . . , p}. Thus we conclude that (B.7) is true for the given λ by choosing anappropriate positive large constant c1 , with probability at least 1 − O(πnp ).
Next to verify (B.8), we derive the population lower bound first and then show that thedifference is negligible. From the construction, we have00e 0 ) X(θe 0 )] = E[n−1 X0 X] = Λ0 Σf Λ0 + σ 2 Ip ,E[n−1 X(θee 0 )0 X] = E[n−1 X0 X(θe 0 )] = Λ0 Σf Λ0 0 .
E[n−1 X(θUsing these equations, we obtain the lower bound0!
00
Λ0 Σf Λ0 + σe2 Ip
Λ0 Σf Λ0δ/kδk2200δ∈VΛ0 Σf Λ0
Λ0 Σf Λ0 + σe2 Ip
(
!
)
1 1
0
00
2
0
⊗ Λ Σf Λ + σe I2p δ/kδk22= inf δδ∈V1 1

inf δ E U(θ ) δ/kδk22 = inf δ 0


0

δ∈V

≥ σe2 .

(B.10)

Because δ ∈ V is sparse and satisfies |B| ≤ k for B := supp(δ), it holds that δ 0 U(θ 0 )δ =
δB0 UB (θ 0 )δB and δ 0 E U(θ 0 ) δ = δB0 E UB (θ 0 ) δB . Hence from Lemma 4 together withthe condition on dimensionality, we obtainsup
UB (θ) − E[UB (θ 0 )] max = O k 1/2 c̃np
|B|≤k, θ∈Θnp= o(s−1 )(B.11)

with probability at least 1 − O(πnp ). Thus using (B.11), we have for any δ ∈ V,
δ 0 E[U(θ 0 )]δ − δ 0 U(θ)δ = δB0 E[UB (θ 0 )] − UB (θ) δB≤ kδk21sup|B|≤k, θ∈Θnp

UB (θ) − E[UB (θ 0 )] max = (kδS 0 k1 + kδS 1 k1 )2 o(s−1 ). kδS 0 k21 o(s−1 ) ≤ kδS 0 k22 o(1) ≤ kδk22 o(1).
Rearranging the terms with (B.10) yields
infθ∈Θnp , δ∈V

δ 0 U(θ)δ/kδk22 ≥ inf δ 0 E[U(θ 0 )]δ/kδk22 − |o(1)| ≥ σe2 − |o(1)|,δ∈V

resulting in (B.8). In consequence, two inequalities (B.7) and (B.8) hold with probability atleast 1 − O(πnp ).
Step 2. This part is well known in the literature (e.g., [33]) so we briefly give the proof
omitting the details. Because the objective function is given byeβ̂ aug (θ) = arg min n−1 y − [X, X(θ)]b
b∈R2p

41

2
2

+ λkbk1 ,the global optimality of the Lasso estimator implieseβ̂ aug (θ)(2n)−1 y − [X, X(θ)]

22+ λ β̂ aug (θ)

auge≤ (2n)−1 y − [X, X(θ)]β

2
21+ λ kβ aug k1 ,where the true parameter vector β aug was defined in the proof of Theorem 2. Note thatsupθ∈Θnp kδ(θ)k0 ≤ k by the assumption. Expanding the inequality and collecting termswith (B.7) yielde2−1 δ 0 U(θ)δ ≤ n−1 ε0 [X, X(θ)]

maxkδk1 + λkδk1 ≤ (3/2)λkδk1 .

(B.12)On the other hand, applying Lemma 1 of [33] to our model reveals that δ ∈ V. Thus we canuse (B.8), (B.12), and (B.9) to getkδk22 (σe2 + o(1)) ≤ 3λkδk1 = 3λ (kδS 1 k1 + kδS 0 k1 ) ≤ 12λkδS 0 k1 .
Since |S 0 | = s and kδS 0 k1 ≤ s1/2 kδS 0 k2 , it holds that kδk2 ≤ 12s1/2 λ/(σe2 + o(1)). SincekδS 0 k2 ≤ kδk2 , we obtain the desired bound kδk1 ≤ 48sλ/(σe2 + o(1)). This bound holds
uniformly over θ ∈ Θnp , which completes the proof of Lemma 5.

C

Additional technical lemmas and their proofs

C.1

Lemma 6 and its proof

Lemma 6 Denote by V ∈ Rr×r a diagonal matrix with its entries the r largest eigenval00
ues of (np)−1 XX0 and define H = (Λ0 Λ0 /p)(F0 F̂/n)V−1 . Assume that Conditions 2–5
hold. Then kHk2 + kH−1 k2 + kVk2 + kV−1 k2 is bounded from above by some constant withprobability at least 1 − O(p−ν ).
Proof. Let λk [A] denote the kth largest eigenvalue of square matrix A throughout the proof.
0
Because kΛ0 Λ0 /pk2 ≤ M and0

kF0 F̂/nk2 ≤ kn−1/2 F0 k2 kn−1/2 F̂k2

1/2
≤ (rn)1/2 kn−1/2 F0 kmax λ1 [n−1 F̂0 F̂]≤ r1/2 M
by Conditions 2–3, and F̂0 F̂/n = Ir , we have0

kHk2 ≤ Λ0 Λ0 /p

02

F0 F̂/n2

V−1 2 . V−1 2 ,where kV−1 k2 is equal to the reciprocal of the rth largest eigenvalue of (np)−1 XX0 . Similarly,
under Conditions 2–3 we also have0

H−1 2 ≤ kVk2 (F0 F̂/n)−102(Λ0 Λ0 /p)−102. kVk2 (F0 F̂/n)−12,where kVk2 is equal to the largest eigenvalue of (np)−1 XX0 and the inverse matrix in the0upper bound is well defined by [2]. To see if k(F0 F̂/n)−1 k2 is bounded from above, it suffices42

0

to bound the minimum eigenvalue of F0 F̂F̂0 F0 /n2 away from zero uniformly in n. Regardingr eigenvalues of the matrix, Sylvester’s law of inertia (e.g., [31], Theorem 4.5.8) entails that
all the r eigenvalues are positive for all n. Moreover, by Proposition 1 of [2] we know thatthe limiting matrix of F̂0 F0 /n is nonsingular under Conditions 2 and 5. Therefore, we can0conclude that lim inf n→∞ λr [F0 F̂F̂0 F0 /n2 ] > 0 a.s., and hence kH−1 k2 . kVk2 follows.
To complete the proof, it is sufficient to show that the maximum and rth largest eigenvalues of (np)−1 XX0 are bounded from above and away from zero, respectively, for all largen and p. By the definition of the spectral norm and triangle inequality, we have
 1
 1/2
λ (np)−1 XX0 )= (np)−1/2 X
≤ n

−1/2

F

0
2

2

p

≤ (np)−1/2 F0 Λ0
−1/2

0

Λ

2

+ (np)02+ (np)−1/2 E

−1/2

E

2

2

.

By Conditions 2 and 3, the first term is a.s. bounded by a constant as discussed above. The
second term is O((n ∧ p)−1/2 ) = o(1) with probability at least 1 − 2 exp(−|O(n ∨ p)|) byLemma 9(a) under Condition 4. Therefore, the largest eigenvalue of (np)−1 XX0 is boundedfrom above by some constant with probability at least 1 − 2 exp(−|O(n ∨ p)|).
Next we bound the rth largest eigenvalue of (np)−1 XX0 away from zero. Since the matrix
is symmetric, Weyl’s inequality (e.g., [31], Theorem 4.3.1) yields
oi

n

h


0
0
0
0
λr (np)−1 XX0 = λr (np)−1 F0 Λ0 Λ0 F0 + EΛ0 F0 + F0 Λ0 E0 + EE0
i

h
i
h


0
0
0
0
≥ λr (np)−1 F0 Λ0 Λ0 F0 + λn (np)−1 EΛ0 F0 + F0 Λ0 E0 + λn (np)−1 EE0 .
(B.13)The third term of lower bound (B.13) is obviously nonnegative. For the first term of00lower bound (B.13), let V denote a subspace of Rn . Because F0 Λ0 Λ0 F0 is symmetric,
the Courant–Fischer min-max Theorem (e.g., [31], Theorem 4.2.6) yields()
ih0 F0 Λ0 0 Λ0 F0 0 v
v
00(np)−1min
λr (np)−1 F0 Λ0 Λ0 F0 = maxv0 v
V:dim(V)=r v∈V\{0}!
!
0 F0 F0 0 v
0 F0 Λ0 0 Λ0 F0 0 v
v
v
≥ max
min
n−1
min
p−1
v0 v
V:dim(V)=r v∈V\{0}v0 F0 F0 0 v
F0 0 v∈Rr \{0}
i
i h
i
h
h
i h
0
0
0
0
= λr n−1 F0 F0 λr p−1 Λ0 Λ0 = λr n−1 F0 F0 λr p−1 Λ0 Λ0hi
00≥ λr [Σf ] λr p−1 Λ0 Λ0 − n−1 F0 F0 − Σf
2
h
i
r
r
−1 0 0 0−1 0 0 0≥ λ [Σf ] λ p Λ Λ − r n F F − Σf.
max

In this lower bound, the first term is bounded away from zero by Conditions 2–3. Meanwhile,
to evaluate the second term we use Lemma 7(d) in Section C.2, which together with the union
bound establishes
!
n


X

−1 0 0 0
2
−1
0 0
0 0
P n F F − Σf>u ≤rmax P nfik fi` − E fik fi` > umax

2k,`∈{1,...,r}2≤ 2r exp(−nu /C)
43

i=1

for any 0 ≤ u ≤ c. Thus the second one turns out to be O((n−1 log p)1/2 ) = o(1) withprobability at least 1 − O(p−ν ) once we set u = (Cνn−1 log p)1/2 and assume n−1 log p ≤c2 /(Cν) without loss of generality. Therefore, the first term of lower bound (B.13) is bounded
away from zero eventually. For the second term of (B.13), since the spectral norm gives theupper bound of the spectral radius we have


h

i
0
0
0
0
≤ (np)−1 EΛ0 F0 + F0 Λ0 E0
λn (np)−1 EΛ0 F0 + F0 Λ0 E0
2

−1/2

−1/2

−1/2

0

≤ 2 (np)
E
p
Λ
n
F
2
2


= O (n ∧ p)−1/2 O(1)O(1) = o(1),02which holds with probability at least 1 − 2 exp(−|O(n ∨ p)|) by Lemma 9(a) in Section C.4.
As a consequence, the desired result holds with probability at least 1 − O(p−ν ) and thisconcludes the proof of Lemma 6.

C.2

Lemma 7 and its proof

Lemma 7 Assume that Conditions 2–4 hold. Then there exist some positive constants c
and C such that the following inequalities hold(a) For all `, i ∈ {1, . . . , n} and 0 ≤ u ≤ c, we have


p
X

P  p−1(e`j eij − E[e`j eij ]) > u ≤ 2 exp −pu2 /C .
j=1(b) For all k ∈ {1, . . . , r}, j ∈ {1, . . . , p}, and 0 ≤ u ≤ c, we have
!
n
X

0
eij > u ≤ 2 exp −nu2 /C .
P n−1
fik
i=1

(c) For all k ∈ {1, . . . , r}, i ∈ {1, . . . , n}, and u ≥ 0, we have


p
X

P  p−1λ0jk eij > u ≤ 2 exp −pu2 /C .
j=1(d) For all k, ` ∈ {1, . . . , r} and 0 ≤ u ≤ c, we have
P

n−1

n
X

!


0 0
0 0
fik
fi` − E[fikfi` ] > u ≤ 2 exp −nu2 /C .

i=1

Proof. (a) To obtain the first result, we rely on the Hanson–Wright inequality. Let ξ =(ξ1 , . . . , ξm )0 ∈ Rm denote a random vector whose components are independent copies ofe ∼ subG(Ce2 ). Then the inequality states that for any (nonrandom) matrix A ∈ Rm×m ,



2

u
u
0
0
eH minP ξ Aξ − E ξ Aξ > u ≤ 2 exp −C,,(B.14)
K 4 kAk2F K 2 kAk2
44

eH is a positive
where K is a positive constant such that supk≥1 k −1/2 (E |e|k )1/k ≤ K and C
2
constant. In our setting, we can take K = 3Ce (e.g., Lemma 1.4 of [34]). Using thisinequality, we first prove the case when ` = i. If we set m = p and A = diag(p−1 , . . . , p−1 ),then we have
ξ 0 Aξ − E ξ 0 Aξ = p−1

pXd(ξj2 − E ξj2 ) = p−1

p
X

j=1

e2ij − E[e2ij ]j=1

for all i. Moreover, we obtain kAk2F = p−1 and kAk2 = p−1 in this case. The assumed
condition 0 < u ≤ 9Ce2 = K 2 entails that u2 /K 4 ≤ u/K 2 so the result follows from (B.14)
eH replaced by CH = 81C 4 /CeH .
with Ce
Similarly, we prove the case when ` 6= i. We set m = p + 1 and A = (a1 , . . . , ap+1 ), where
a1 = (0, p−1 , . . . , p−1 )0 and aj = 0 for j = 2, . . . , p + 1. That is, the entries of A are all zero
except that the second to (p + 1)th components in the first column vector are p−1 . Under
this setting, we observe that
ξ 0 Aξ − E ξ 0 Aξ = p−1

p+1
X
j=2

d

ξ1 ξj = p−1

pXe`j eij

j=1

for all ` 6= i. Moreover, we obtain kAk2F = kAk2 = p−1 in this case. Therefore, the same
bound holds as in the case of ` = i from (B.14) again. Consequently, for any 0 ≤ u ≤ 9Ce2
we havepXP  p−1(e`j eij − E[e`j eij ]) > u ≤ 2 exp −pu2 /CH .
j=1(b) We prove the second assertion by Bernstein’s inequality for the sum of a martingale
difference sequence (e.g., Theorem 3.14 in [11]). Fix k = 1 and j = 1. Define Fi−1 as
0 : ` = i, i − 1, . . . }. Then (f 0 e , F ) forms a martingale
the σ-field generated from {f`1
i
i1 i1
0
0
difference sequence because E |fi1 ei1 | < ∞ and E[fi1 ei1 |Fi−1 ] = 0 under Conditions 2 and
4. Sincehthe sub-Gaussianityof ei1 implies E e2i1 ≤ 4Ce2 (e.g., Lemma 1.4 of [34]), we have
i
Pn
2
2
2
2
0 2 e2 | F
Vi := E fik
i−1 ≤ 4Ce M , and henceij
i=1 Vi ≤ 4nCe M a.s. due to boundedness
0 | ≤ M a.s. On the other hand, by the sub-Gaussianity of e and boundedness of |f 0 |
|fi1
ij
i1
again we observe that for all p ≥ 3 and i ∈ {1, . . . , n},0
E (0 ∨ fi1ei1 )p | Fi−1 ≤ M p (2Ce2 )p/2 pΓ(p/2) ≤ p!(2Ce M )p−2 Vi /2,
where Γ denotes the Gamma function and we have used the estimates pΓ(p/2) ≤ p! and2p/2−2 ≤ 2p−2 /2 for p ≥ 3 in the last inequality. Then an application of Theorem 3.14 in [11]by putting x = u, y = 4M 2 Ce2 , and c = 2M Ce in their notation gives the one-sided result.
Making twice the bound yields
!


n
X
nu2
−1
0
.
P n
fik eij > u ≤ 2 exp −8M 2 Ce2 + 4M Ce u
i=145For all 0 ≤ u ≤ M Ce2 , the upper bound is further bounded by 2 exp(nu2 /(12M 2 Ce2 )). We
set CI = 12M 2 Ce2 . Consequently, for any 0 ≤ u ≤ M Ce2 we have
!
n
X

0
P n−1
fik
eij > u ≤ 2 exp −nu2 /CI .
i=1(c) We prove the third inequality. Note that
(
)


2

u
u2
0
P λjk eij > u ≤ 2 exp − 02 2 ≤ 2 exp −.
2M 2 Ce2
2λjk Ce
This implies that λ0jk eij is a sequence of i.i.d. subG(M 2 Ce2 ). Thus the result is obtained
directly by Bernstein’s inequality for the sum of independent sub-Gaussian random variables.
Consequently, for any u ≥ 0 putting CJ = M 2 Ce2 leads to


p
X

P  p−1λ0jk eij > u ≤ 2 exp −pu2 /CJ .
j=1(d) We show the last inequality. Note that for each k, (fik )i ∼ i.i.d. subG(M 2 ) since
0 | ≤ M a.s. by Lemma 1.8 of [34] under Condition 2. Thus the remaining is the same as|fikeH here. Then for any 0 ≤ u ≤ 9M 2 , we have(a). Set CK = 81M 4 /C
!
n
X


−1
0 0
0 0
P n
fik fi` − E[fik fi` ] > u ≤ 2 exp −nu2 /CK .
i=1

Finally the obtained inequalities hold even if the constant in the upper bound is replacedwith arbitrary fixed constant C such that C ≥ max{CH , CI , CJ , CK }. Similarly, we can also
restrict the range of u for each inequality to be 0 ≤ u ≤ c for arbitrary fixed constant c thatsatisfies 0 < c ≤ min(9Ce2 , M Ce2 , 9M 2 ). This completes the proof of Lemma 7.

C.3

Lemma 8 and its proof

Lemma 8 Assume that Conditions 1–4 hold. Then for any set A satisfying |A| ≤ k, thefollowing statements hold with probability at least 1 − O(πnp )(a) sup UA (θ) − UA (θ 0 ) max = O k 1/2 c̃np ,
θ∈Θnp

(b)



vA (θ) − vA (θ 0 ) max = O s3/2 c̃np ,

sup
θ∈Θnp

where Θnp was defined in Lemma 3 and c̃np = n−1/2 log p + p−1/2 log n. Consequently, we
have



sup TA (θ) − TA (θ 0 ) max = O k 1/2 + s3/2 c̃np .
θ∈Θnp

Proof. We first state some results that are useful in the proof. Since kn−1/2 F0 k2 = O(1) a.s.
by Condition 2 and kk −1/2 Λ0A k2 = O(1) for any A such that |A| ≤ k under Condition 3, we
first haven−1/2 C0A2

≤ n−1/2 F0

2

46

k 1/2 k −1/2 Λ0A

2. k 1/2 .

Next Lemma 9(b) in Section C.4 gives directlyn−1/2 Eη0 A.1

2(B.15)

with probability at least 1 − O(p−ν ). By Condition 4, we also deduce!
PEη − Eη0 max > u

sup
η∈Θnp

≤ np max Peηij − eη0 ij > u

sup

i,j

!

η∈Θnp



≤ np max P |Z| > u/(M 1/2 c1/2
)
np
i,j

≤ 2np exp −u2 / c2e M cnp
1/2

for any u ≥ 0. Thus setting u = 2ce M 1/2 cnp log1/2 (np) with some large enough positive
constant M , we obtain that with probability at least 1 − O((np)−ν ),
sup
η∈Θnp

Eη − Eη0 max . cnp log1/2 (np) = O(c̃np ).

We will use these results and Lemma 10 in Section C.5 in the proofs below.
To prove (a), we havee A (θ)0 X
e A (θ) − n−1 X
e A (θ 0 )0 X
e A (θ 0 )UA (θ) − UA (θ 0 ) max ≤ n−1 X
+2 n

−1

0
e A (θ) − n−1 X0 Xe
X0A X
A A (θ )max

max=: U1 + U2 .

Observe that U1 is further bounded as
0

U1 ≤ n−1 C0A CA − n−1 C0A C0A
+2 n

−1

+ n−1 E0ηA EηA − n−1 E0η0 A Eη0 A

max

E0ηA CA − n−1 E0η0 A C0A

max

max

=: U11 + U12 + U13 .

By Lemma 10, it is easy to see thatU11 ≤ n−1 CA − C0A
≤n

−1/2

0

CA − C0A

CA − C0A max

0


max

+ 2 n−1 C0A CA − C0A

CA − C0A 2 + 2

n

−1/2

C0A

2


max

CA − C0A max

2. k 1/2 C − C0 max + k 1/2 C − C0 max
= O k 1/2 c2np + k 1/2 cnp = O k 1/2 cnp ,where the last estimate follows from Lemma 3. Similarly, we deduceU12 ≤ n−1 EηA − Eη0 A
≤n

−1/2

0

EηA − Eη0 A

max+ 2 n−1 E0η0 A EηA − Eη0 A

EηA − Eη0 A max EηA − Eη0 A 2 + 2 n
2. k 1/2 Eη − Eη0 max + Eη − Eη0 max
= O k 1/2 c̃2np + c̃np

47

−1/2

Eη 0 A2max

EηA − Eη0 A max

andU13 ≤ n−1 EηA − Eη0 A

0
CA − C0A
max



0
−1 0
0
+ n Eη0 A CA − CA
+ n−1 C0A EηA − E0η0 A

≤k

1/2

Eη − Eη0 max

max
C − C0 max

max

Eη − E0η0+ n−1/2 Eη0 A
C − C0 max + n−1/2 C0A
2
max
2




1/2
1/2
1/2
= O k c̃np cnp + cnp + k c̃np = O k c̃np .
Combining these bounds of U11 –U13 , we haveU1 ≤ U11 + U12 + U13 . k 1/2 c̃np .
This holds uniformly in θ ∈ Θnp with probability at least 1 − O(πnp ) by Lemma 3 and thediscussion above. Next we obtain0

U2 ≤ n−1 C0A (CA − C0A )+0+ n−1 C0A (EηA − Eη0 A )

max
−1 0
0
n EA (CA − CA ) max +≤ n−1/2 C0A

2

n

−1

max
0
EA (EηA − Eη0 A ) max

CA − C0A max + n−1/2 C0A

2EηA − Eη0 A max

Eη − Eη0 max
CA − C0A max + n−1/2 Eη0 A+ n−1/2 Eη0 A
22
= O k 1/2 cnp + k 1/2 c̃np + cnp + c̃np
= O k 1/2 c̃np .
This also holds uniformly in θ ∈ Θnp with probability at least 1 − O(πnp ) by Lemma 3 andthe discussion above. Consequently, it holds that
sup
θ∈Θnp

UA (θ) − UA (θ 0 ) max . k 1/2 c̃npwith probability at least 1 − O(πnp ).
To prove (b), we havee A (θ)0 y − n−1 X
e A (θ 0 )0 yvA (θ) − vA (θ 0 ) max ≤ n−1 X
≤ n

−1 e

0

XA (θ) Xβ − n

−1 e

0 0

XA (θ ) Xβ

max

+ n

max

−1 e

e A (θ 0 )0 ε
XA (θ)0 ε − n−1 X=: V1 + V2 .
First, because Xβ = XS 0 βS 0 we see thate A (θ)0 XS 0 − n−1 X
e A (θ 0 )0 XS 0
V1 ≤ s1/2 n−1 X
e A (θ)0 XS 0 − n−1 X
e A (θ 0 )0 XS 0
. s n−1 X

48

max

max

.

kβS 0 k2max

Recall that |S 0 | = s and s ≤ n ∧ p. By a similar bound of U2 , the norm just above can bebounded further as
n−1/2 C0S 0

2CA − C0A max + n−1/2 C0S 0+ n−1/2 Eη0 S 02

2EηA − Eη0 A max

CA − C0A max + n−1/2 Eη0 S 02

EηA − Eη0 A max. s1/2 C − C0 max + s1/2 Eη − Eη0 max + C − C0 max + Eη − Eη0 max
= O s1/2 cnp + s1/2 c̃np + cnp + c̃np = O s1/2 c̃np .
Thus we have
V1 . ss1/2 c̃np = s3/2 c̃npwith probability at least 1 − O(πnp ). Next the same procedure yieldse A (θ) − X
e A (θ 0 )V2 ≤ X
0

. C−C

max

max

n1/2 ε2+ Eη − Eη0 max . c̃np ,(B.16)

where kn1/2 εk2 = (E ε2 )1/2 + o(1) a.s. by the law of large numbers for independent randomvariables. Since the results hold uniformly in θ ∈ Θnp , combining them leads to
sup
θ∈Θnp

vA (θ) − vA (θ 0 ) max . s3/2 c̃npwith probability at least 1 − O(πnp ). This concludes the proof of Lemma 8.

C.4

Lemma 9 and its proof

Lemma 9 Assume that Condition 4 holds. Then the following statements hold(a) We haveP(n ∨ p)−1/2 E

2. 1 ≥ 1 − 2 exp(−|O(n ∨ p)|);(b) For any fixed set A with |A| ≤ k ≤ n, we have
P n−1/2 EA . 1 ≥ 1 − 2p−ν ;2(c) For all k ≤ n, we have
P max n−1/2 EA
|A|≤k

2

.1∨ n

−1

k log p1/2



≥ 1 − 2p−ν ,where ν > 0 is a predetermined constant.
Proof. Result (a) is obtained by Theorem 5.39 of [40]. Moreover, by the same theorem there
exist some positive constants c and C such that for any A with |A| ≤ k ≤ n and every t ≥ 0,
P σe−1 kn−1/2 EA k2 > 1 + C + n−1/2 t ≤ 2 exp −ct2 ,49

where σe2 = E e2 . Therefore, result (b) is immediately obtained by putting t2 = c−1 ν log psince n−1/2 t = o(1) and exp −ct2 = p−ν in this case.
For (c), taking the union bound leads to


−1
−1/2
−1/2
P σe max kn
E A k2 > 1 + C + n
t
|A|≤k
 



p
≤
max P σe−1 kn−1/2 EA k2 > 1 + C + n−1/2 t ≤ 2pk exp −ct2 .
k |A|≤k
Set t2 = c−1 (ν + k) log p in this inequality. Then we have n−1/2 t = O (n−1 k log p)1/2 and
2pk exp −ct2 ≤ 2pk exp (−(ν + k) log p) = 2p−ν ,which gives result (c) and completes the proof of Lemma 9.

C.5

Lemma 10 and its proof

Lemma 10 For matrices A ∈ Rk1 ×n and B ∈ Rn×k2 , we have kABkmax ≤ n1/2 kAk2 kBkmax
and kABkmax ≤ n1/2 kAkmax kBk2 .
Proof. For any matrix M = (mij ) ∈ Rk×n , let kMk∞,∞ denote the induced `∞ -norm. First,
we have
kMk∞,∞ :=supv∈Rn \{0}kMvk2 kvk2
kMvkmax≤ sup≤ n1/2 kMk2 .
kvkmax
kvk
kvk
n
2
max
v∈R \{0}Therefore, by a simple calculation we see thatkABkmax = k vec(AB)kmax = k(Ik2 ⊗ A) vec(B)kmax=

k(Ik2 ⊗ A) vec(B)kmax
k vec(B)kmax
k vec(B)kmax≤ kIk2 ⊗ Ak∞,∞ k vec(B)kmax = kAk∞,∞ kBkmax ≤ n1/2 kAk2 kBkmax .
The second assertion follows from applying this inequality to B0 A0 . This concludes the proof
of Lemma 10.

50