Goodness of fit tests for high-dimensional linear models
Rajen D. Shah∗
University of Cambridge

Peter Bühlmann
ETH Zürich

arXiv:1511.03334v4 [stat.ME] 8 Apr 2017

April 11, 2017

Abstract
In this work we propose a framework for constructing goodness of fit tests in both low
and high-dimensional linear models. We advocate applying regression methods to the scaled
residuals following either an ordinary least squares or Lasso fit to the data, and using some
proxy for prediction error as the final test statistic. We call this family Residual Prediction
(RP) tests. We show that simulation can be used to obtain the critical values for such tests
in the low-dimensional setting, and demonstrate using both theoretical results and extensive
numerical studies that some form of the parametric bootstrap can do the same when the highdimensional linear model is under consideration. We show that RP tests can be used to test for
significance of groups or individual variables as special cases, and here they compare favourably
with state of the art methods, but we also argue that they can be designed to test for as diverse
model misspecifications as heteroscedasticity and nonlinearity.

1

Introduction

High-dimensional data, where the number of variables may greatly exceed the number of observations, has become increasingly more prevalent across a variety of disciplines. While such data
pose many challenges to statisticians, we now have a variety of methods for fitting models to highdimensional data, many of which are based on the Lasso [Tibshirani, 1996]; see Bühlmann and van
de Geer [2011] for a review of some of the developments.
More recently, huge strides have been made in quantifying uncertainty about parameter estimates. For the important special case of the high-dimensional linear model, frequentist p-values for
individual parameters or groups of parameters can now be obtained through an array of different
techniques [Wasserman and Roeder, 2009, Meinshausen et al., 2009, Bühlmann, 2013, Zhang and
Zhang, 2014, Lockhart et al., 2014, van de Geer et al., 2014, Javanmard and Montanari, 2014,
Meinshausen, 2015, Ning and Liu, 2014, Voorman et al., 2014, Zhou, 2015]—see Dezeure et al.
[2015] for an overview of some of these methods. Subsampling techniques such as Stability Selection [Meinshausen and Bühlmann, 2010] and its variant Complementary Pairs Stability Selection
(CPSS) [Shah and Samworth, 2013] can also be used to select important variables whilst preserving
error control in a wider variety of settings.
Despite these advances, something still lacking from the practitioner’s toolbox is a corresponding
set of diagnostic checks to help assess the validity of, for example, the high-dimensional linear
∗

Supported in part by the Forschungsintitut für Mathematik (FIM) at ETH Zürich.

1

model. For instance, there are no well-established methods for detecting heteroscedasticity in
high-dimensional linear models, or whether a nonlinear model may be more appropriate.
In this paper, we introduce an approach for creating diagnostic measures or goodness of fit
tests that are sensitive to different sorts of departures from the ‘standard’ high-dimensional linear
model. As the measures are derived from examining the residuals following e.g. a Lasso fit to the
data, we use the name Residual Prediction (RP) tests. To the best of our knowledge, it is the
first methodology for deriving confirmatory statistical conclusions, in terms of p-values, to test for
a broad range of deviations from a high-dimensional linear model. In Section 1.2 we give a brief
overview of the idea, but first we discuss what we mean by goodness of fit in a high-dimensional
setting.

1.1

Model misspecification in high-dimensional linear models

Consider the Gaussian linear model
y = Xβ + σε,

(1)

where y ∈ Rn is a response vector, X ∈ Rn×p is the fixed design matrix, β ∈ Rp is the unknown
vector of coefficients, ε ∼ Nn (0, I) is a vector of uncorrelated Gaussian errors, and σ 2 > 0 is
the variance of the noise. In the low-dimensional situation where p < n, we may speak of (1)
being misspecified such that E(y) 6= Xβ. When X has full row rank however, any vector in Rn
can be expressed as Xβ for some β ∈ Rp , leaving in general no room for nonlinear alternatives.
When restricting to sparse linear models specified by (1), the situation is different though and
misspecification can happen [Bühlmann and van de Geer, 2015]; we will take a sparse Gaussian
linear model as our null hypothesis (see also Theorems 3 and 4). We discuss an approach to handle
a relaxation of the Gaussian error assumption in Section B of the supplementary material.
When there is no good sparse approximation to Xβ, a high-dimensional linear model may
not be an appropriate model for the data-generating process: a sparse nonlinear model might
be more interpretable and may generalise better, for example. Moreover, the Lasso and other
sparse estimation procedures may have poor performance, undermining the various different highdimensional inference methods mentioned above that make use of them. Our proposed RP tests
investigate whether the Lasso is a good estimator of the signal.

1.2

Overview of Residual Prediction (RP) tests and main contributions

Let R̂ be the residuals following a Lasso fit to X. If Xβ is such that it can be well-estimated by
the Lasso, then the residuals should contain very little signal and instead should behave roughly
like the noise term σε. On the other hand, if the signal is such that the Lasso performs poorly
and instead a nonlinear model were more appropriate, for example, some of the (nonlinear) signal
should be present in the residuals, as the Lasso would be incapable of fitting to it.
Now if we use a regression procedure that is well-suited to predicting the nonlinear signal (an
example may be Random Forest [Breiman, 2001]), applying this to the residuals and computing
the resulting mean residual sum of squares (RSS) or any other proxy for prediction error will give
us a test statistic that under the null hypothesis of a sparse linear model we expect to be relatively
large, and under the alternative we expect to be relatively small. Different regression procedures
applied to the residuals can be used to test for different sorts of departures from the Gaussian linear
model. Thus RP tests consist of three components.

2

1. An initial procedure that regresses y on X to give a set of residuals; this is typically the Lasso
if p > n or could be ordinary least squares if X is low-dimensional.
2. A residual prediction method (RP method) that is suited to predicting the particular signal
expected in the residuals under the alternative(s) under consideration.
3. Some measure of the predictive capability of the RP method. Typically this would be the
residual sum of squares (RSS), but in certain situations a cross-validated estimate of prediction
error may be more appropriate, for example.
We will refer to the composition of a residual prediction method and an estimator of prediction
error as a residual prediction function. This must be a (measurable) function f of the residuals and
all available predictors, pall of them in total, to the reals f : Rn × Rn×pall → R. For example, if
rather than testing for nonlinearity, we wanted to ascertain whether any additional variables were
significant after accounting for those in X, we could consider the mean RSS after regressing the
residuals on a matrix of predictors containing both X and the additional variables, using the Lasso.
If the residuals can be predicted better than one would expect under the null hypothesis with a
model as in (1), this provides evidence against the null.
Clearly in order to use RP tests to perform formal hypothesis tests, one needs knowledge of
the distribution of the test statistic under the null, in order to calculate p-values. Closed form
expressions are difficult if not impossible to come by, particularly when the residual prediction
method is something as intractable as Random Forest.
In this work, we show that under certain conditions, the parametric bootstrap [Efron and Tibshirani, 1994] can be used, with some modifications, to calibrate any RP test. Thus the RP method
can be as exotic as needed in order to detect the particular departure from the null hypothesis that
is of interest, and there are no restrictions requiring it to be a smooth function of the data, for example. In order to obtain such a general result, the conditions are necessarily strong; nevertheless,
we demonstrate empirically that for a variety of interesting RP tests, bootstrap calibration tends
to be rather accurate even when the conditions cannot be expected to be met. As well as providing
a way of calibrating RP tests, we also introduce a framework for combining several RP tests in
order to have power against a diverse set of alternatives.
Although formally the null hypothesis tested by our approach is that of the sparse Gaussian
linear model (1), an RP test geared towards nonlinearity is unlikely to reject purely due to nonGaussianity of the errors, and so the effective null hypothesis typically allows for more general error
distributions. By using the nonparametric bootstrap rather than the parametric bootstrap, we can
allow for non-Gaussian error distributions more explicitly. We discuss this approach in Section B of
the supplementary material, where we see that type I error is very well controlled even in settings
with t3 and exponential errors.
Some work related to ours here is that of Chatterjee and Lahiri [2010], Chatterjee and Lahiri
[2011], Camponovo [2014] and Zhou [2014] who study the use of the bootstrap with the (adaptive)
Lasso for constructing confidence sets for the regression coefficients. Work that is more closely
aligned to our aim of creating diagnostic measures for high-dimensional models is that of Nan and
Yang [2014], though their approach is specifically geared towards variable selection and they do not
provide theoretical guarantees within a hypothesis testing framework as we do.

3

1.3

Organisation of the paper

Simulating the residuals under the null is particularly simple when rather than using the Lasso
residuals, ordinary least squares residuals are used. We study this simple situation in Section 2 not
only to help motivate our approach in the high-dimensional setting, but also to present what we
believe is a useful method in its own right. In Section 3 we explain how several RP tests can be
aggregated into a single test that combines the powers of each of the tests. In Section 4 we describe
the use of RP tests in the high-dimensional setting, and prove the validity of a calibration procedure
based on the parametric bootstrap. We give several applications of RP tests in Section 5 along
with the results of extensive numerical experiments, and conclude with a discussion in Section 6.
The supplementary material contains further discussion of the power of RP tests; a proposal for
how to test null hypotheses of the form (1) allowing for more general error distributions; additional
numerical results; a short comment concerning the interpretation of p-values; and all of the proofs.
The R [R Development Core Team, 2005] package RPtests provides an implementation of the
methodology.

2

Ordinary least squares RP tests

A simple but nevertheless important version of RP tests uses residuals from ordinary least squares
(OLS) in the first stage. For this, we require p < n in the set-up of (1). Let P denote the
orthogonal projection on to the column space of X. Then under the null hypothesis that the model
(1) is correct, the scaled residuals R̂ are
R̂ :=

(I − P)y
(I − P)ε
=
,
k(I − P)yk2
k(I − P)εk2

and so their distribution does not depend on any unknown parameters: they form an ancillary
statistic. Note that the scaling of the residuals eliminates the dependence on σ 2 . It is thus simple
to simulate from the distribution of any function of the scaled residuals, and this allows critical
values to be calculated for tests using any RP method.
We note that by using OLS applied to a larger set of variables as the RP method, and the RSS
from the resulting fit as the estimate of prediction error, the overall test is equivalent to a partial
F -test for the significance of the additional group of variables. To see this let us write Z ∈ Rn×q for
an additional group of variables. Let Pall be the orthogonal projection on to all available predictors,
that is projection on to Xall = (X, Z) ∈ Rn×pall , where pall = p + q. When the RP method is OLS
regression of the scaled residuals on to Xall , the resulting RSS is
k(I − Pall )R̂k22 =

k(I − Pall )(I − P)yk22
k(I − Pall )yk22
=
,
k(I − P)yk22
k(I − P)yk22

since (I − Pall )P = 0. We reject for small values of the quantity above, or equivalently large values
of
k(Pall − P)yk22 n − pall
×
,
pall − p
k(I − Pall )yk22
which is precisely the F -statistic for testing the hypothesis in question.
An alternative way to arrive at the F -test is to first residualise Z with respect to X and define
new variables Z̃ = (I − P)Z. Let us write P̃ for the orthogonal projection on to Z̃. Now if our RP
4

method is OLS regression of the scaled residuals on to Z̃, we may write our RSS as
k(I − P̃)R̂k22 =

k(I − P̃)(I − P)yk22
k{I − (P + P̃)}yk22
=
,
k(I − P)yk22
k(I − P)yk22

the final equality following from the fact that the column spaces of X and Z̃ and hence P and P̃
are orthogonal. It is easy to see that P + P̃ = Pall , and so we arrive at the F -test once more.
We can use each of the two versions of the F -test above as starting points for generalisation,
where rather than using OLS as a prediction method, we use other RP methods more tailored to
specific alternatives of interest. The distribution of the output of an RP method under the null
hypothesis of a linear model can be computed via simulation as follows. For a given B > 1 we
generate independent n-vectors with i.i.d. standard normal components ζ (1) , . . . , ζ (B) . From these
we form scaled residuals
(I − P)ζ (b)
R̂(b) =
.
(2)
k(I − P)ζ (b) k2
Let Xall be the full matrix of predictors. Writing the original scaled residuals as R̂ we apply our
chosen RP function f to all of the scaled residuals to obtain a p-value


B
X
1
1+
1{f (R̂(b) ,Xall )≤f (R̂,Xall )} .
B+1

(3)

b=1

See also Section 4 for the extension to the case using Lasso residuals.
Even in situations where the usual F -test may seem the natural choice, an RP test with a
carefully chosen RP method can often be more powerful against alternatives of interest. This is
particularly true when we aggregate the results of various different RP methods to gain power over
a diverse set of alternatives, as we describe in the next section.

3

Aggregating RP tests

In many situations, we would like to try a variety of different RP methods, in order to have
power against various different alternatives. A key example is when an RP method involves a
tuning parameter such as the Lasso. Each different value of the tuning parameter effectively gives
a different RP method. One could also aim to create a generic omnibus test to test for, say,
nonlinearity, heteroscedasticity and correlation between the errors, simultaneously.
To motivate our approach for combining the results of multiple RP tests, we consider the famous
diabetes dataset of Efron et al. [2004]. This has p = 10 predictors measured for n = 442 diabetes
patients and includes a response that is a quantitative measure of disease progression one year after
baseline. Given the null hypothesis of a Gaussian linear model, we wish to test for the presence of
interactions and quadratic effects. In order to have power against alternatives composed of sparse
coefficients for these effects, we consider as RP methods the Lasso applied to quadratic effects
residualised with respect to the linear terms via OLS. We regress the OLS scaled residuals onto the
transformed quadratic effects using the Lasso with tuning parameters on a grid of λ values, giving
a family of RP tests.
We plot the residual sums of squares from the Lasso fits to the scaled residuals in Figure 1, as
a function of λ. Also shown are the residual sums of squares from Lasso fits to scaled residuals
5

30
10

0.94

0

0

RSS

0.96

0.98

1.00

RSS

RSS

0.80

0.85

0.90

0.95

1.00

0.80

Density

50

50
30

Density

10

0.90
0.85

0.00

0.05

0.10

0.15

λ

Figure 1: Bottom plot: the residual sums of squares from Lasso fits to the original scaled residuals
(blue) and simulated residuals (red), as well as the mean of the latter (black) and the mean displaced
by one and two standard deviations (black, dotted). Top plots: kernel density plots for the simulated
residual sums of squares at λ = 0 (left) and λ = 0.8 (right) with the original residual sums of squares
in blue.
simulated under the null hypothesis of a Gaussian linear model, as simulation under the null
hypothesis is the general principle which we use for deriving p-values.
At the point λ = 0, the observed RSS is not drastically smaller than those of the simulated
residuals, as the top left density plot shows. Indeed, were we to calculate a p-value just based
on the λ = 0 results corresponding to the F -test, we would obtain roughly 10%. The output at
λ = 0.8, however, does provide compelling evidence against the null hypothesis, as the top right
density plot shows. Here the observed RSS is far to the left of the support of the simulated residual
sums of squares. In order to create a p-value for the presence of interactions based on all of the
output, we need a measure of how ‘extreme’ the entire blue curve is, with respect to the red curves,
in terms of carrying evidence against the null. Forming a p-value based on such a test statistic is
straightforward, as we now explain.
Suppose we have residual prediction functions fl , l = 1, . . . , L (in our example these would be
the RSS when using the Lasso with tuning parameter λl ) and their evaluations on the true scaled
(b)
residuals R̂(0) := R̂ and simulated scaled residuals {R̂(b) }B
= fl (R̂(b) , Xall ), let
b=1 . Writing fl
(b)
f (b) = {fl }L
l=1 be the curve or vector of RP function evaluations at the bth scaled residuals, and
0
(−b)
denote by f
= {f (b ) }b0 6=b the entire collection of curves, potentially including the curve for the

6

true scaled residuals R(0) , but excluding the bth curve. Let
Q̃ : RL × RL×B → R
(f (b) , f (−b) ) 7→ Q̃(f (b) , f (−b) )
be any measure of how extreme the curve f (b) is compared to the rest of the curves f (−b) (larger
values indicating more extreme). Here Q̃ can be any function such that Q̃b := Q̃(f (b) , f (−b) ) does
not depend on the particular ordering of the curves in f (−b) ; we will give a concrete example below.
We can use the {Q̃b }b6=0 to calibrate our test statistic Q̃0 as detailed in the following proposition.
Proposition 1. Suppose the simulated scaled residuals are constructed as in (2). Setting


B
X
1
1+
Q=
1{Q̃b ≥Q̃0 } ,
B+1
b=1

we have that under the null hypothesis (1), P(Q ≤ x) ≤ x for all x ∈ [0, 1], so Q constitutes a valid
p-value.
The result above is a straightforward consequence of the fact that under the null {Q̃b }B
b=0 form
an exchangeable sequence, and standard results on Monte Carlo testing (see Davison and Hinkley
[1997] Ch. 4 for example). Under an alternative, we expect Q̃0 to be smaller and Q̃b for b ≥ 1 to
be larger on average, than under the null. Thus this approach will have more power than directly
comparing Q̃0 to a sample from its null distribution.
(−b)
(−b)
respectively be the empirical
We recommend constructing Q̃ as follows. Let µ̂l
and σ̂l
(b0 )

mean and standard deviation of {fl

}b0 6=b . We then set
(−b)

Q̃b = max {(µ̂l
l

(b)

(−b)

− fl )/σ̂l

},

(4)

the number of standard deviations by which the bth curve lies below the rest of the curves, max(1)
imised along the curve. The intuition is that were fl to have a Gaussian distribution under the
(0)
(−0)
null for each l, Φ{(fl − µ̂l )/σ̂ (−0) } would be an approximate p-value based on the lth RP
function, whence Φ(Q̃0 ) would be the minimum of these p-values. Though it would be impossible
to match the power of the most powerful test for the alternative in question (perhaps that corresponding to λ = 0.8 in our diabetes example) among the L tests considered, one would hope
to come close. We stress however that this choice of Q̃ (4) yields valid p-values regardless of the
(1)
distribution of fl under the null.
Using this approach with a grid of L = 100 λ values, we obtain a p-value of under 1% for the
diabetes example. As discussed in Section 1.2, this low p-value is unlikely to be due to a deviation
from Gaussian errors, and indeed when we take our simulated errors ζ (b) to be resamples from the
vector of residuals (see Section B of the supplementary material), we also obtain a p-value under
1%; clear evidence that a model including only main effects is inappropriate for the data. Further
simulations demonstrating the power of this approach are presented in Section 5.

4

Lasso RP tests

When the null hypothesis is itself high-dimensional, we can use Lasso residuals in the first stage
of the RP testing procedure. Although unlike scaled OLS residuals, scaled Lasso residuals are not
7

ancillary, we will see that under certain conditions, the distribution of scaled Lasso residuals are
not wholly sensitive to the parameters β and σ in (1).
Let us write R̂λ (β, σε) for the scaled Lasso residuals when the tuning parameter is λ (in squareroot parametrisation, see below):
√
β̂ λ (β, σε) ∈ arg min {kX(β − b) + σεk2 / n + λkbk1 }
(5)
b∈Rp

R̂λ (β, σε) =

X{β − β̂(β, σε)} + σε
kX{β − β̂(β, σε)} + σεk2

.

√
Note that under (1) β̂ λ (β, σε) ∈ arg min {ky − Xbk2 / n + λkbk1 } and R̂λ (β, σε) = {y −
b∈Rp

Xβ̂(β, σε)}/ky − Xβ̂(β, σε)k2 . Sometimes we will omit the first argument of β̂ for convenience
in which case it will always be the true parameter value under the null, β. Here we are using
the Lasso in the square-root parametrisation [Belloni et al., 2011, Sun and Zhang, 2012] rather
than the conventional version where the term in the objective assessing the model fit would be
kX(β − b) + σεk22 . We note that the two versions of the Lasso have identical solution paths but
these will simply be parametrised differently. For this reason, we will simply refer to (5) as the
Lasso solution. Note that while the Lasso solution may potentially be non-unique, the residuals
are always uniquely defined as the fitted values from a Lasso fit are unique (see Tibshirani [2013],
for example). Throughout we will assume that the columns of X have been scaled to have `2 -norm
√
n.
We set out our proposal for calibrating RP tests based on Lasso residuals using the parametric
bootstrap in Algorithm 1 below. In the following section we aim to justify the use of the parametric
Algorithm 1 Lasso RP tests
1. Let β̌ be an estimate of β, typically a Lasso estimate selected by cross-validation.
√
2. Set σ̌ = ky − Xβ̌k2 / n.
(b)
are i.i.d. draws from
3. Form B scaled simulated residuals {R̂λ (β̌, σ̌ζ (b) )}B
b=1 where the ζ
Nn (0, I), and λ chosen according to the proposal of Sun and Zhang [2013].

4. Based on the scaled simulated residuals {R̂λ (β̌, σ̌ζ (b) )}B
b=1 , compute a p-value (3) or use these
to form an aggregated p-value as described in Section 3.

bootstrap from a theoretical perspective and also discuss the particular choices β̌, σ̌ and λ used
above.

4.1

Justification of Lasso RP tests

Given b ∈ Rp and a set A ⊆ {1, . . . , p}, let bA be the subvector of b with components consisting
of those indexed by A. Also for a matrix M, let MA be the submatrix of M containing those
columns indexed by A, and let Mk = M{k} , the kth column. The following result shows that
if sgn(β̌) = sgn(β), with the sign function understood as being applied componentwise, we have

8

partial ancillarity of the scaled residuals. In the following we let S = {j : βj 6= 0} be the support
set of β.
Theorem 2. Suppose β̌ is such that sgn(β̌) = sgn(β). For t ∈ [0, 1) and λ > 0, consider the
deterministic set

Λλ,t = {ζ ∈ Rn : sgn β̂ λ,S (β, σζ) = sgn(β S ) and min β̂j (β, σζ)/βj > t}.
j∈S

Then we have that for all ζ ∈ Λλ,t , R̂λ (β, σζ) = R̂λ (β̌, σ̌ζ) provided 0 < σ̌/σ < minj∈S β̌j /{(1 −
t)βj }.
In words, provided the error ζ is in the set Λλ,t and conditions for β̌ and σ̌ are met, the scaled
residuals from a Lasso fit to y = Xβ + σζ are precisely equal to the scaled residuals from a Lasso
fit to Xβ̌ + σ̌ζ. Note that all of the quantities in the result are deterministic. Under reasonable
conditions and for a sensible choice of λ (see Theorem 3), when ζ ∼ Nn (0, I), we can expect the
event ζ ∈ Λλ,t to have large probability. Thus Theorem 2 shows that the scaled residuals are not
very sensitive to the parameter σ or to the magnitudes of the components of β, but instead depend
largely on the signs of the latter. It is the square-root parametrisation that allows the result to
hold for a large range of values of σ̌, and in particular for all σ̌ sufficiently small.
Theorem 2 does not directly justify a way to simulate from the distribution of the scaled Lasso
residuals as in Algorithm 1 since the sign pattern of β̌ must equal that of β. Accurate estimation
of the sign pattern of β using the Lasso requires a strong irrepresentable or neighbourhood stability
condition [Meinshausen and Bühlmann, 2006, Zhao and Yu, 2006]. Nevertheless, we now show that
we can modify Algorithm 1 to yield provable error control under more reasonable conditions. In
Section 4.1.2 we argue heuristically that the same error control should hold for Algorithm 1 in a
wide range of settings.
4.1.1

Modified Lasso RP tests

Under a so-called beta-min condition (see Theorem 3 below), with high probability we can arrive
at an initial estimate of β, β 0 via the Lasso for which sgn(β 0S ) = sgn(β S ), and where minj∈S |βj0 | >
maxj∈S c |βj0 |. With such a β 0 , we can aim to seek a threshold τ for which the Lasso applied only
on the subset of variables XSτ where Sτ := {j : |βj0 | > τ } yields an estimate that has the necessary
sign agreement with β. This then motivates Algorithm 2 based on maximising over the candidate
p-values obtained through different β estimates derived from applying the Lasso to different subsets
of the initial active set (see also Chatterjee and Lahiri [2011] which introduces a related scheme).
Note we do not recommend the use of Algorithm 2 in practice; we only introduce it to facilitate
theoretical analysis which sheds light on our proposed procedure Algorithm 1. Let s = |S| and
s0 = |{j : βj0 6= 0}|. The theorem below gives conditions under which with high probability, s0 ≥ s
(s)

and residuals from responses generated around β̌ will equal the true residuals. This then shows
that the maximum p-value Q will in general be a conservative p-value as it will always be at least
as large as Qs , on an event with high probability.
As well as a beta-min condition, the result requires some relatively mild assumptions on the
design matrix. Let C (ξ, T ) = {u : kuT c k1 ≤ ξkuT k1 , u 6= 0}. The restricted eigenvalue [Bickel
et al., 2009, Koltchinskii, 2009] is defined by


√
kXuk2 / n
φ(ξ) = inf
: u ∈ C (ξ, S) .
(6)
kuk2
9

Algorithm 2 Modified Lasso RP tests (only used for Theorem 3)
1. Let β 0 = β̂ λ (σε) be the Lasso estimate of β.
2. Let s0 = |{j : βj0 6= 0}| and suppose 0 < |βj0 s0 | ≤ · · · ≤ |βj0 1 | are the non-zero components of β 0
arranged in order of non-decreasing magnitude. Define Ŝ (k) = {j1 , j2 , . . . , jk }.
3. For k = 1, . . . , s0 let β̌
β̌

(0)

(k)

be the Lasso estimate from regressing y on XŜ (k) . Further set

= 0.

√
(k)
(k)
4. Using each of the β̌
in turn and σ̌ (k) = ky − Xβ̌ k2 / n, generate sets of residuals
(k)
(b)
{R̂λ (β̌ , σ̌ (k) ζ (b) )}B
are i.i.d. draws from Nn (0, I). Use these to create
b=1 where the ζ
corresponding p-values Qk for RP tests based on (3) or the method introduced in Section 3.
5. Output Q = maxk=0,...,s0 Qk as the final approximate p-value.

For a matrix M ∈ Rn×p , T ⊂ {1, . . . , p} and ξ > 1, the compatibility factor κ(ξ, T, M) [van de
Geer and Bühlmann, 2009] is given by


√
kMuk2 / n
κ(ξ, T, M) = inf
: u ∈ C (ξ, T ) .
(7)
kuT k1 /|T |
When either of the final two arguments are omitted, we shall take them to be S and X respectively;
the more general form is required in Section 4.2. The sizes of κ(ξ) and φ(ξ) quantify the ill-posedness
of the the design matrix X; we will require κ( ξ), φ(ξ) > 0 for some ξ > 1. Note that in the random
design setting where the rows of X are i.i.d. multivariate normal with the minimum eigenvalue of
the covariance matrix bounded away from zero, the factors (6) and (7) can be thought of as positive
constants in asymptotic regimes where s log(p)/n → 0. We refer the reader to van de Geer and
Bühlmann [2009] and Zhang and Zhang [2012] for further details.
p
Theorem√3. Suppose the data follows the Gaussian linear model (1). Let λ = A 2 log(p/η)/n
with A > 2 and pe−s−2 > η > 0. Suppose for ξ > 1 that
√


s log(p/η)
1
2(ξ + 1) 1
≤ 2
min 1 −
,
.
(8)
nκ2 (ξ)
A (ξ + 1)
A(ξ − 1) 5
Assume a beta-min condition
√
σ
min |βj | > 10 2Aξ
j∈S

p
s log(p/η)
√
.
φ2 (ξ) n

(9)

Then for all x ∈ [0, 1],
2(1 + rn−s )η
P(Q ≤ x) ≤ x + p
+ e−n/8
π log(p/η)
where rm → 0 as m → ∞.

10

(10)

p Although the beta-min condition, which is of the form is of the form minj∈S |βj | ≥ const. ×
s log(p)/n, may be regarded as somewhat strong, the conclusion is correspondingly strong: any
RP method or collection of RP methods with arbitrary Q̃ for combining tests can be applied to
the residuals and the result remains valid. It is also worth noting however that the conditions are
only required under the null. For example, if the alternative of interest was that an additional
variable z was related to the response after accounting for those in the original design matrix X,
no conditions on the relationship between X and z are required for the test to be valid.
More importantly though, the conditions are certainly not necessary for the conclusion to hold.
The scaled residuals are a function of the fitted values and the response, and do not involve Lasso
parameter estimates directly. Thus whilst duplicated columns in X could be problematic for inferential procedures relying directly on Lasso estimates such as the debiased Lasso [Zhang and Zhang,
2014], they pose no problem for RP tests. In addition, given a particular RP method, exact equality
of the residuals would not be needed to guarantee a result of the form (10).
4.1.2

Relevance of Theorem 3 to Algorithm 1

In the special case of testing for the significance of a single predictor described above, we have a
much stronger result than Theorem 3 (see Theorems 4 and 5) which shows that neither the betamin condition nor the maximisation over candidate p-values of Algorithm 2 is necessary for error
control to hold. More generally, in our experiments we have found Qs0 is usually equal to or close to
the maximum Q for large B across a variety of settings. Thus selecting Qs0 rather than performing
the maximisation (which amounts to Algorithm 1) is able to deliver conservative error control as
evidenced by the simulations in Section 5.
A heuristic explanation for why the error is controlled is that typically the amount of signal
(k)
(k)
remaining in R̂λ (β̌ , σ̌ζ) increases with k, simply because typically kXβ̌ k2 also increases with
k. This can result in the prediction error of a procedure applied to the various residuals decreasing
with k because the signal-to-noise ratios tend to be increasing; thus the p-values tend to increase
with k.
In addition, when the Lasso performs well, we would expect residuals to contain very little
signal, and any differences in the signals contained in R̂λ (β, σζ) and R̂λ (β̌, σ̌ε) to be smaller still,
particularly when Xβ and Xβ̌ are close. Typically the RP function will be insensitive to such small
differences since they are unlikely to be too close to directions against which power is desired. We
now discuss the choices of β̌, λ and σ̌ in Algorithm 1.
4.1.3

Practical considerations

Choice of β̌. In view of the preceding discussion, it suffices for β̌ to satisfy a screening-type
property: we would like the support of β̌ to contain that of β. Though Theorem 3 suggests a
fixed λ, since β̌ only needs to be computed once, we can use cross-validation. This is the perhaps
the most standard way of producing an estimate that performs well for screening (see for example
Section 2.5.1 of Bühlmann and van de Geer [2011]).
If the folds for cross-validation are chosen at random, the estimate will have undesirable randomness beyond that of the data. We thus suggest taking many random partitions into folds and
using an estimate based on a λ that minimises the cross-validation error curve based on all of the
folds used. In our simulations in Section 5 we partition the observations into 10 random folds a
total of 8 times.
11

Choice of σ̌. The normalised RSS is perhaps the most natural choice for σ̌ 2 (see also Reid et al.
[2016]), though as Theorem 2 suggests, the results are essentially unchanged when this is doubled
or halved, for example.
Choice of λ for the Lasso residuals. The choice of λ should be such that with high probability,
the resulting estimate
contains the support
of β (see Theorem 2). Though Theorem 3 suggests
p
√
taking λ = A 2 log(p)/n for A > 2, the restriction on A is an artefact of basing our result
on oracle inequalities from Sun and Zhang [2012], which place relatively simple conditions on the
design. Sun and Zhang [2013] has a more involved theory which suggests a slightly smaller λ. We
therefore use their method, the default in the R package Sun [2013], as a convenient fixed choice of
λ.

4.2

Testing the significance of individual predictors

Here we consider the collection of null hypotheses Hk : βk = 0 and their corresponding alternatives
that βk 6= 0. Note that for this setting there are many other approaches that can perform the
required tests. Our aim here is to show that RP tests can be valid under weaker assumptions than
those laid out in Theorem 3, and moreover that the simpler approach of Algorithm 1 can control
type I error.
We begin with some notation. For Ak := {1, . . . , p} \ {k} and b ∈ Rp let b−k = bAk and
X−k = XAk . For each variable k, our RP method will be a least squares regression onto a version
of Xk that has been residualised with respect to X−k . Since in the high-dimensional setting X−k
will typically have full row rank, an OLS regression of Xk on X−k will return the 0-vector as
residuals. Hence we will residualise Xk using the square-root Lasso:
√
Ψk = arg min {kXk − X−k bk2 / n + γkbk1 }.
b∈Rp−1

This RP method is closely related to the pioneering idea by Zhang and Zhang [2014] and similar
to that of Ning and Liu [2014], who consider using the regular Lasso (without the square-root
parametrisation) at each stage. If Xk were not residualised with respect to X−k , and the regular
Lasso were used, the resulting RP method would be similar to that of Voorman et al. [2014]. The
work of Ren et al. [2015] studies an analogous procedure in the context of the Gaussian graphical
model.
Let Wk be the residual Xk − X−k Ψk . Note for each k we may write
y = X−k Θk + βk Wk + σε
where Θk = β −k + βk Ψk ∈ Rp−1 . Let Θ̂k be the square-root Lasso regression of y on to X−k with
tuning parameter λ. Our RP function will be the RSS from OLS regression of the scaled Lasso
residuals (y − X−k Θ̂k )/ky − X−k Θ̂k k2 on to Wk . Note this is an RP function even though it
involves the residualised version of Xk , Wk ; the latter is simply a function of X. Equivalently, we
can consider the test statistic Tk2 with Tk defined by
Tk =

WkT (y − X−k Θ̂k )

√ .
kWk k2 ky − X−k Θ̂k k2 / n

12

Note that Tk is simply a regularised partial correlation between y and Xk given X−k . The bootstrap
version is
∗

Tk∗ =

WkT (yk∗ − X−k Θ̂k )

√ ,
∗
kWk k2 kyk∗ − X−k Θ̂k k2 / n
∗

where yk∗ = X−k Θ̂k + σ̌ε∗ , ε∗ ∼ Nn (0, I) and Θ̂k is the Lasso regression of y∗ on X−k . Here we
√
will consider taking σ̌ = ky − Xβ̂k2 / n where β̂ is the square-root Lasso regression of y on the
full design matrix X.
As before, let S be the support of β, which without loss of generality we will take to be {1, . . . , s},
and also let N = {1, . . . , p} \ S be the set of true nulls. Assume ε ∼ Nn (0, I). The following result
shows that only a relatively mild compatibility condition is needed in order to ensure that the type
I error is controlled. We consider an asymptotic regime with n → ∞ where β, X and p are all
allowed to vary with n though we suppress this in the notation. In the following we denote the
cumulative distribution function of the standard normal by Φ.
p
p
Theorem 4. Let λ = A1 2 log(p)/n for some constant
A1 > 1 and suppose that s log(p)2 /n/κ2 (ξ, S) →
p
0 for some ξ > (A1 + 1)/(A1 − 1). Let γ = A2 2 log(p)/n for some constant A2 > 0. Define
B = {b ∈ Rp : bN = 0}. Then
|P(Tk ≤ x) − Φ(x)| → 0,

sup
k∈N, β∈B, x∈R

p

sup
k∈N, β∈B, x∈R

|P(Tk∗ ≤ x|ε) − Φ(x)| → 0.

We see that a bootstrap approach can control the type I error uniformly across the noise
variables and β ∈ B. We note that this result is for a fixed design X and does not require any
sparsity assumptions on the inverse covariance matrix of a distribution that could have generated
the rows of X [Ning and Liu, 2014], for example.
Theorem 5. Let λ and γ be as in Theorem 4. Assume that for some ξ p
with ξ > (A1 + 1)/(A1 − 1)
there ispa sequence of sets {1, . . . , s−1} ⊆ T ⊂ {1, . . . , p−1} such that |T | log(p)2 /n/κ2 (ξ, T, X−k ) →
√
0 and log(p)kΘk,T c k1 → 0 where T c = {1, . . . , p−1}\T . Further assume that βk kWk k2 / n → 0.
Define Bk = {b ∈ B : bk = βk }. Then
sup
β∈Bk ,x∈R


q
P(Tk ≤ x) − Φ x − βk kWk k2
σ 2 + βk2 kWk k22 /n → 0,


sup
β∈Bk ,x∈R

p

|P(Tk∗ ≤ x|ε) − Φ(x)| → 0.

If Ψk and hence Θk were sparse, we could take T as the set of nonzeroes and the second
condition involving kΘk,T c k1 would be vacuous. This would be the case with high probability in
the random design setting where X has i.i.d. Gaussian rows with sparse inverse covariance matrix
[van de Geer et al., 2014]. However, Θk,T c can also have many small coefficients provided they
have small `1 -norm. The result above shows that the power of our method is comparable to the
proposals of Zhang and Zhang [2014] and van de Geer et al. [2014] based on the debaised Lasso.
√
If kWk k2 = O( n) as would typically be the case in the random design setting discussed above,

13

√
we would have power tending to 1 if βk → 0 but n|βk | → ∞. Further results on power to detect
nonlinearities are given in Section A of the supplementary material.
The theoretical results do not suggest any real benefit from using the bootstrap as to test
hypotheses we can simply compare Tk to a standard normal distribution. However our experience
has been that this can be slightly anti-conservative in certain settings. Instead, we propose to use
the bootstrap to estimate the mean and standard deviation of the null distribution of the Tk by
computing the empirical mean m̂k and standard deviation v̂k of B samples of Tk∗ . Then we take as
our p-values 2[1 − Φ{|(Tk − m̂k )/v̂k |}].
This construction of p-values appears to yield tests that very rarely have size exceeding their
nominal level. Indeed in all our numerical experiments we found no evidence of this violation
occurring. An additional advantage is that only a modest number of bootstrap samples is needed
to yield the sort of low p-values that could fall below the threshold of a typical multiple testing
procedure. We recommend choosing B between 50 and 100.
4.2.1

Computational considerations

Using our bootstrap approach for calibration presents a significant computational burden when it is
applied to test for the significance of each of a large number of variables in turn. Some modifications
to Algorithm 1 can help to overcome this issue and allow this form of RP tests to be applied to
typical high-dimensional data with large p.
Firstly rather than using cross-validation to choose λ for computation of Θ̂k , we recommend
using the fixed λ of Sun and Zhang [2013] (see also Section 4.1.3). The tuning parameter γ required
to compute Wk can be chosen in the same way, and we also note that these nodewise regressions
only need to be done once rather than for each bootstrap sample. Great computational savings can
be realised by first regressing y on X to yield coefficients β̂. Writing Ŝ = {k : β̂k 6= 0}, we know
that for each k ∈
/ Ŝ, Θ̂k = β̂ −k , so we only need to compute Θ̂k for those k in Ŝ. The same logic
∗
can be applied to computation of Θ̂k for the bootstrap replicates.
We also remark that approaches for directly simulating Lasso estimates [Zhou, 2014] may be used
to produce simulated residuals. These have the potential to substantially reduce the computational
burden; not just in the case of testing significance of individual predictors but for RP tests in
general.

5
5.1

Applications
Low-dimensional nulls

Here we return to the problem of testing for quadratic effects in the diabetes dataset used in the
example of Figure 1. In order to further investigate the power of the aggregate RP test constructed
through Lasso regressions on a grid of 100 λ values as described in Section 3, we created artificial
signals from which we simulated responses. The signals (mean responses) were constructed by
selecting at random s of the quadratic terms and giving these coefficients generated using i.i.d.
Unif[−1, 1] random variables. The remaining coefficients for the variables were set to 0, so s
determined the sparsity level of the signal. Responses were generated by adding i.i.d. Gaussian
noise to the signals, with variance chosen such that the F -test for the presence of quadratic effects
has power 0.5 when the size is fixed at 0.05. We created 25 artificial signals at each sparsity level
s ∈ {1, 4, 10, 20, 35, 54}. Note that the total number of possible quadratic effects was 54 (as one
14

0.6
0.4

Power

0.8

1.0

of the variables was binary), so the final sparsity level represents fully dense alternatives where we
might expect the F -test to have good power. We note however that the average power of the F -test
in the dense case rests critically on the form of the covariance between the generated quadratic
coefficients, with optimality guarantees only in special circumstances (see Section 8 of Goeman
et al. [2006]). For the RP tests, we set the number of bootstrap samples B to be 249.
We also compare the power of RP tests to the global test procedure of Goeman et al. [2006].
The results, shown in Figure 2, suggest that RP tests can outperform the F -test in a variety of
settings, most notably when the alternative is sparse, but also in dense settings. When there are
small effects spread out across many variables (s ∈ {35, 54}), the global test tends to do best;
indeed in such settings it is optimal. In the sparser settings, RP tests perform better.

1

4

10

20

35

54

Sparsity

Figure 2: Boxplots of the power of RP tests (grey) and the global test (white) across the 25 signals
estimated through 100 repetitions, for each of the sparsity levels s; the power of the F -test is fixed
at 0.5 and shown as a dotted line.

5.2

High-dimensional nulls

In this section we report the results of using RP tests (Algorithm 1) tailored to detect particular
alternatives on a variety of simulated examples where the null hypothesis is high-dimensional. We
investigate both control of the type I error and the powers of the procedures.
Our examples are inspired by Dezeure et al. [2015]. We use n × p simulated design matrices
with p = 500 and n = 100 except for the setting where we test for heteroscedasticity in which we
increase n to 300 in order to have reasonable power against these alternatives. The rows of the
matrices are distributed as Np (0, Σ) with Σ given by the three types described in Table 1.
Table 1: Generation of Σ.
Toeplitz:
Exponential decay:
Equal correlation:

Σjk = 0.9|j−k|
(Σ−1 )jk = 0.4|j−k|/5
Σjk = 0.8 if j 6= k and 1 otherwise.

In addition to the randomly generated design matrices, we also used a publicly available real
design matrix from gene expression data of Bacillus Subtilis with n = 71 observations and p =
15

4088 predictors [Bühlmann et al., 2014]. Similarly to Dezeure et al. [2015], in order to keep the
computational burden of the simulations manageable, we reduced the number of variables to p = 500
by selecting only those with the highest empirical variance. For each of the four design settings,
we generated 25 design matrices (those from the real data were all the same). The columns of the
√
design matrices were mean-centred and scaled to have `2 -norm n.
In order to create responses under the null hypothesis, for each of these 100 design matrices, we
randomly generated a vector of coefficients β as follows. We selected a set S of 12 variables from
{1, . . . , p}. We then assigned β S c = 0 and each βk with k ∈ S was generated according to Unif[−2, 2]
independently of other coefficients. This form of signal is similar to the most challenging signal
settings considered in Dezeure et al. [2015] and also resembles the estimated signal from regression
of the true response associated with the gene expression data on to the predictors using the Lasso
or MCP [Zhang, 2010]. Other constructions for generating the non-zero regression coefficients are
considered in Section C in the supplementary material. Given X and β, we generated r = 100
responses according to the linear model (1) with σ = 1. Thus in total, here we evaluate the type
I error control of our procedures on over 100 data-generating processes. The number of bootstrap
samples B used was 100 when testing for significance of individual predictors and fixed at 249 in
all other settings.
We now explain interpretation of the plots in Figures 3–5; a description of Figure 6 is given
in Section 5.2.4. The top and bottom rows of each of Figures 3–5 concern settings under null and
alternative hypotheses respectively. Thin red curves trace the empirical cumulative distribution
functions (CDFs) of the p-values obtained using RP tests, whilst thin blue curves, where shown,
represent the same for debiased Lasso-based approaches. In all plots, thickened coloured curves are
averages of their respective thin coloured curves; note these are averages over different simulation
settings.
The black dashed line is the CDF of the uniform distribution; thus we would hope for the
empirical CDFs to be close to this in the null settings (top rows), and rise above it in the bottom
rows indicating good power. Of course, even if all of the p-value distributions were stochastically
larger than uniform so the type I error was always controlled, we would not expect their estimated
distributions i.e. the empirical CDFs to always lie below the dashed line. The black dotted curve
allows us to assess type I error control across the simulation settings more easily. It is constructed
such that in each of the plots, were the type I error to be controlled exactly, we would expect on
average 1 out of the 25 empirical CDFs for RP tests to escape above the region the line encloses.
Thus several curves not completely enclosed under the dotted line in a given plot would indicate
poor control of type I error. More precisely, the line is computed as follows. Let qα (x) be the upper α
quantile of a Bin(B +1, x)/(B +1) distribution. Note this is the marginal distribution of Û (x) where
Û is the empirical CDF of B samples from the uniform distribution on {1/(B +1), 2/(B +1), . . . , 1}.
The curve then traces qα (x) with α chosen such that

P max (Û (x) − qα (x)) > 0 = 1/25.
x∈[0,0.1]

We see that across all of the data-generating processes and for each of the three RP testing methods,
it appears the size never exceeds the nominal level by a significant amount. Moreover the same
holds for the additional 100 data-generating processes whose results presented in the supplementary
material: the type I error is controlled well uniformly across all settings considered.
We now describe the particular RP tests used in Figures 3–5, and the alternatives investigated,
as well as the results shown in Figure 6 concerning testing for the significance of individual predictors
16

False positive rate
True positive rate

Figure 3: Testing significance of groups: the empirical distribution functions of the p-values from
RP tests (red) and the debiased Lasso (blue) under the null (top row) and alternative (bottom row)
respectively. The dashed line equals the 45 degree line corresponding to the Unif[0, 1] distribution
function, and the dotted curve is explained in the main text.
as detailed in Section 4.2.
5.2.1

Groups

We consider the problem of testing the null hypothesis β G = 0 within linear model (1). One
approach is to regress each column of XG on to XGc in turn using the square-root Lasso (c.f.
Section 4.2), and consider a matrix of residuals X̃ ∈ Rn×|G| . We may then use Lasso regression on
to X̃ as our family of RP methods and combine the resulting test statistics as in Section 3.
We use this approach on our simulated data and the results are displayed in red in Figure 3. For
the null settings (top row) we took Gc to be a randomly selected set of size p/2 containing S. Thus
under the null, β Gc had 12 non-zero components whilst β G = 0. The alternatives, corresponding to
the bottom row, also modify the signal such that β A is non-zero (in addition to β S being non-zero
as was the case under the null) with coefficients generated in exactly the same way as for β S and
A being a randomly selected set of 12 variables chosen from G.
The blue lines trace the empirical CDFs of p-values constructed using the debiased Lasso proposal of van de Geer et al. [2014] and implemented in the hdi package Dezeure et al. [2015] for
R. More specifically, we use the minimum of the p-values associated with each of the coefficients
in G (see Section 2.3 of van de Geer et al. [2014]) as our test statistic, and calibrate this using
the Westfall–Young procedure [Westfall and Young, 1993] as explained in Bühlmann [2013]. This
ensures that no power is lost due to correlations among the individual p-values, as would be the
case with Bonferroni correction, for example. Remaining parameters were set to the defaults in the
hdi package.
17

False positive rate
True positive rate
True positive rate

False positive rate

Figure 4: Testing for nonlinearity; the interpretation is similar to that of Figure 3.

Figure 5: Testing for heteroscedasticity; the interpretation is similar to that of Figure 3.

18

False positive rate
True positive rate

Figure 6: Testing individual variables: the plots give the proportion of |S c | = 488 null (top row)
and |S| = 12 true variables (bottom row) selected at various threshold levels with RP tests (red)
and the debiased Lasso (blue).
Although the sizes of the debiased Lasso-based tests averaged over the equal correlation design
examples are very close to the nominal level, this is due to the several settings where the size exceeds
the desired level being compensated for by other examples where the tests are more conservative.
On the other hand, RP tests have slightly conservative type I error control across all the examples,
and greater power among the Toeplitz and Exponential decay settings.
5.2.2

Nonlinearity

In order to test for nonlinearity, we consider an RP method based on Random Forest [Breiman,
2001]. We used the default settings for Random Forest as implemented by Liaw and Wiener [2002],
but rather than using a direct application to the residuals we apply it to the equicorrelation set:
the set of variables with maximum absolute correlation with the residuals. This is invariably the
set of variables selected in the initial Lasso fit, though in situations where the Lasso solution is not
unique this will in general be a superset of the support of any Lasso solution. Using this smaller
set of variables reduces the computational burden of a Random Forest fit, and also gives the test
greater power in situations where the variables contributing to the nonlinear signal also feature in
sparse linear approximations to the truth. Applying a Random Forest to the entire set of variables
may have slightly greater power when this is not the case, but would have greatly diminished power
in the more natural situations where this this holds. Rather than using the RSS from the Random
Forest fits as our proxy for prediction error, we use the out of bag error. This has the advantage of
being more insensitive to the size of the equicorrelation set and tends to result in greater power.
To create the nonlinear signal for the alternative settings, we randomly divided S into four
groups of three. Each variable x was transformed via a sigmoid composed with a random affine
19

mapping as below:
x 7→ [1 + exp{−5(a + bx)}]−1 .
Here a, b ∈ N (0, 1) independently. The transformed variables in each group were multiplied together, and a linear combination of these resulting products with Unif[−1, 1] generated coefficients
formed the nonlinear component of the signal. This nonlinear signal was then scaled such that the
residuals from an OLS fit to the variables in S had an empirical variance of 2, and finally added to
the linear signal.
The results displayed in Figure 4 show that RP tests are able to deliver reasonable power in many
of the settings considered, though the real design examples appear to be particularly challenging.
5.2.3

Heteroscedasticity

As testing for heteroscedasticity in a high-dimensional setting is rather challenging, here we increase
the number of observations for the simulated design settings to n = 300 in order to have reasonable
power against the alternative. The data-generation procedure under the null was left unchanged. In
order to generate vectors of variances for the alternative settings, we randomly selected 3 variables
from S and formed a linear combination of these variables with Unif[−2, 2] coefficients. A constant
was then added, so the minimum component was 0.01, and finally the vector was scaled so the
average of its components was 1. This vector then determined the variance of normal errors added
to the signal.
To detect this heteroscedasticity, we used a family of RP methods given by Lasso regression of
the absolute values of the residuals onto the equicorrelation set. The results are shown in Figure 5.
RP tests are able to deliver reasonable power in the the simulated design settings, but do struggle
to detect the heteroscedasticity with the real design which has a lower number of observations
(n = 71).
5.2.4

Testing significance of individual predictors

Figure 6 shows the results of using RP tests as described in Section 4.2 to test hypotheses Hk :
βk = 0. The red curves give the average proportions of false (top row) and true positives (bottom
row) that would be selected given p-value thresholds varying along the x-axis. Thus for example
in order to obtain the expected number of false positives selected at a given threshold, the y values
should be multiplied by p − |S| = 488. The blue curves display the same results for the debaised
Lasso as implemented in the hdi package. The dashed 45 degree line gives the expected proportion
of false positives that would be incurred by an exact test.
We see that even at the low p-value thresholds particularly relevant for multiple testing correction, RP tests give consistent error control whilst also delivering superior or equal power. Such
error control effectively requires accurate knowledge of the extreme tails of the null distribution of
the test statistics. We see here that the debiased Lasso approach is not always able to achieve this
in the Toeplitz and Exponential decay settings, and indeed error control for multiple testing is rare
among the currently available methods [Dezeure et al., 2015].

6

Discussion

The RP testing methodology introduced in this work allows us to treat model checking as a prediction problem: that of fitting any (prediction) function to the scaled residuals from OLS or Lasso.
20

This makes the problem of testing goodness of fit amenable to the entire range of prediction methods that have been developed across statistics and machine learning. We have investigated here RP
tests for detecting significant single or groups of variables, heteroscedasticity, or deviations from
linearity, and we expect that effective RP methods can also be found for testing for correlated errors, heterogeneity and other sorts of departures from the standard Gaussian linear model. Related
ideas should be applicable to test for model misspecification in high-dimensional generalised linear
models, for example.

References
A. Belloni, V. Chernozhukov, and L. Wang. Square-root lasso: pivotal recovery of sparse signals
via conic programming. Biometrika, 98(4):791–806, 2011.
P. Bickel, Y. Ritov, and A. Tsybakov. Simultaneous Analysis of Lasso and Dantzig selector. Ann.
Statist., 37:1705–1732, 2009.
L. Breiman. Random Forests. Machine Learning, 45:5–32, 2001.
P. Bühlmann. Statistical significance in high-dimensional linear models. Bernoulli, 19:1212–1242,
2013.
P. Bühlmann and S. van de Geer. Statistics for High-Dimensional data: Methods, Theory and
Applications. Springer, 2011.
P. Bühlmann and S. van de Geer. High-dimensional inference in misspecified linear models. Electron.
J. Statist., 9:1449–1473, 2015.
P. Bühlmann, M. Kalisch, and L. Meier. High-dimensional statistics with a view toward applications
in biology. Annual Review of Statistics and Its Application, 1:255–278, 2014.
L. Camponovo. On the validity of the pairs bootstrap for lasso estimators. Biometrika, to appear,
2014.
A. Chatterjee and S. Lahiri. Asymptotic properties of the residual bootstrap for lasso estimators.
Proc. Am. Math. Soc., 138(12):4497–4509, 2010.
A. Chatterjee and S. N. Lahiri. Bootstrapping lasso estimators. J. Am. Statist. Ass., 106(494):
608–625, 2011.
A. C. Davison and D. V. Hinkley. Bootstrap methods and their application, volume 1. Cambridge
university press, 1997.
R. Dezeure, P. Bühlmann, L. Meier, and N. Meinshausen. High-dimensional Inference: Confidence
intervals, p-values and R-Software hdi. Statistical Science, 30:533–558, 2015.
B. Efron and R. J. Tibshirani. An introduction to the bootstrap. CRC press, 1994.
B. Efron, T. Hastie, I. Johnstone, and R. Tibshirani. Least Angle Regression. Ann. Statist., 32:
407–451, 2004.

21

J. J. Goeman, S. A. van de Geer, and H. C. van Houwelingen. Testing against a high-dimensional
alternative. J. R. Statist. Soc. B, 68:477–493, 2006.
A. Javanmard and A. Montanari. Confidence intervals and hypothesis testing for high-dimensional
regression. J. Mach. Learn. Res., 15:2869–2909, 2014.
V. Koltchinskii. The dantzig selector and sparsity oracle inequalities. Bernoulli, 15(3):799–828,
2009.
A. Liaw and M. Wiener. Classification and regression by randomforest. R News, 2(3):18–22, 2002.
URL http://CRAN.R-project.org/doc/Rnews/.
R. Lockhart, J. Taylor, R. J. Tibshirani, and R. Tibshirani. A significance test for the lasso. Ann.
Statist., 42:413–468, 2014.
N. Meinshausen. Group bound: confidence intervals for groups of variables in sparse high dimensional regression without assumptions on the design. J. R. Statist. Soc. B, 77(5):923–945,
2015.
N. Meinshausen and P. Bühlmann. High dimensional graphs and variable selection with the Lasso.
Ann. Statist., 34:1436–1462, 2006.
N. Meinshausen and P. Bühlmann. Stability selection (with discussion). J. R. Statist. Soc. B, 72:
417–473, 2010.
N. Meinshausen, L. Meier, and P. Bühlmann. P-values for high-dimensional regression. J. Am.
Statist. Ass., 104:1671–1681, 2009.
Y. Nan and Y. Yang. Variable selection diagnostics measures for high-dimensional regression. J.
Computnl Graph. Statist., 23(3):636–656, 2014.
Y. Ning and H. Liu. A general theory of hypothesis tests and confidence regions for sparse high
dimensional models. arXiv preprint arXiv:1412.8765, 2014.
R Development Core Team. R: A language and environment for statistical computing. R Foundation
for Statistical Computing, Vienna, Austria, 2005. URL http://www.r-project.org.
S. Reid, R. Tibshirani, and J. Friedman. A study of error variance estimation in lasso regression.
Statistica Sinica, to appear, 2016.
Z. Ren, T. Sun, C.-H. Zhang, and H. H. Zhou. Asymptotic normality and optimalities in estimation
of large Gaussian graphical models. Ann. Statist., 43:991–1026, 2015.
R. D. Shah and R. J. Samworth. Variable selection with error control: another look at stability
selection. J. R. Statist. Soc. B, 75(1):55–80, 2013.
T. Sun. scalreg: Scaled sparse linear regression, 2013. URL https://CRAN.R-project.org/
package=scalreg. R package version 1.0.
T. Sun and C.-H. Zhang. Scaled sparse linear regression. Biometrika, 99(4):879–898, 2012.

22

T. Sun and C.-H. Zhang. Sparse matrix inversion with scaled lasso. The J. Mach. Learn. Res., 14
(1):3385–3418, 2013.
R. Tibshirani. Regression shrinkage and selection via the lasso. J. R. Statist. Soc. B, 58:267–288,
1996.
R. J. Tibshirani. The lasso problem and uniqueness. Electron. J. Statist., 7:1456–1490, 2013.
S. van de Geer and P. Bühlmann. On the conditions used to prove oracle results for the lasso.
Electron. J. Statist., 3:1360–1392, 2009.
S. van de Geer, P. Bühlmann, Y. Ritov, and R. Dezeure. On asymptotically optimal confidence
regions and tests for high-dimensional models. Ann. Statist., 42:1166–1202, 2014.
A. Voorman, A. Shojaie, and D. Witten. Inference in high dimensions with the penalized score
test. arXiv preprint arXiv:1401.2678, 2014.
L. Wasserman and K. Roeder. High dimensional variable selection. Ann. Statist., 37:2178–2201,
2009.
P. Westfall and S. Young. Resampling-based multiple testing: Examples and methods for p-value
adjustment. John Wiley & Sons, 1993.
C.-H. Zhang. Nearly unbiased variable selection under minimax concave penalty. Ann. Statist., 38:
894–942, 2009.
C.-H. Zhang and S. S. Zhang. Confidence intervals for low dimensional parameters in high dimensional linear models. J. R. Statist. Soc. B, 76:217–242, 2014.
C.-H. Zhang and T. Zhang. A general theory of concave regularization for high-dimensional sparse
estimation problems. Statistical Science, 27(4):576–593, 2012.
P. Zhao and B. Yu. On Model Selection Consistency of Lasso. J. Mach. Learn. Res., 7:2541–2563,
2006.
Q. Zhou. Monte carlo simulation for lasso-type problems by estimator augmentation. J. Am.
Statist. Ass., 109(508):1495–1516, 2014.
Q. Zhou. Uncertainty quantification under group sparsity. arXiv preprint arXiv:1507.01296, 2015.

23

Supplementary material
This supplementary material is organised as follows. Section A contains results on the power of
the RP tests approach for detecting nonlinearity.
In Section B we discuss how the RP tests methodology can be extended to test for null hypotheses of linear models with non-Gaussian errors, and present numerical results in support of our
proposed scheme. Additional numerical results to complement those of Section 5 in the main paper
are presented in Section C. In Section D we provide some brief comments on the interpretation of
p-values derived from RP tests. Finally the proofs of all of the results in the main paper, we well
as those stated in Section A, are collected in Section E. Note that all equations numbered 1–10 are
in the main paper.

A

Power of Lasso RP tests

In this section we briefly discuss the power of RP tests for detecting nonlinearity. Suppose the
response is generated according to
y = Xβ + f + σε,
where β ∈ Rp is a sparse vector with S = {j : βj 6= 0}, s = |S| and as before ε ∼ Nn (0, I). The
nonlinear term f is to be thought of as a vector of function evaluations of some nonlinear function:
fi = f (xi,S ) where f : R|S| → R, though this is not assumed in the sequel.
As in Section 4 of the main paper, here we require that the columns of X have been scaled
√
to have `2 -norm
p n. To facilitate theoretical analysis, we will assume β̌ is a Lasso estimate with
fixed λ = A1 2 log(p)/n and A1 > 1, rather than with the tuning parameter selected by crossvalidation as in Algorithm 1. Furthermore, we will also take this to be the tuning parameter used
in the construction of the Lasso scaled residuals R̂ = R̂λ (β, f + σε). Let the bootstrap scaled
residuals be R̂∗ = R̂λ (β̌, σ̌ζ). We will also take this λ to be the tuning parameter used in the
construction of the Lasso scaled residuals R̂ = R̂λ (β, f + σε). Let the bootstrap scaled residuals
√
derived from β̌ and σ̌ := ky − Xβ̌k2 / n be R̂∗ := R̂λ (β̌, σ̌ζ) where ζ ∈ Nn (0, I).
To quantify the potential power of RP tests, we define
√
ψ γ = arg min {kf − Xbk2 / n + γkbk1 }
b∈Rp

and let wγ = f − Xψ γ be the nonlinear signal f residualised with respect to X. As in Section 4.2
we consider an asymptotic regime where p, X, β, S and f can all change as n → ∞, though we
suppress this in the notation. Also, as in Theorem 4, let B = {b ∈ Rp : bS c = 0}.
The result, which follows from Theorem 5 and its proof, shows that whilst the true residuals
are positively correlated with the residualised signal wγ , the bootstrap residuals are not.
√
√
Corollary 6. Suppose kf k2 / n → 0 and for some γ we have nγkψ γ,S c k1 → 0 and kf k2 γ =
p
p
o( log(p)). Assume there exists ξ > (A1 + 1)/(A1 − 1) with sγ log(p)/κ2 (ξ, S) → 0. We have
sup
β∈B,x∈R

P(wγT R̂/kwγ k2



q
≤ x) − Φ x − kwγ k2
σ 2 + kwγ k22 /n → 0,

sup
β∈B,x∈R

p

|P(wγT R̂∗ /kwγ k2 ≤ x|ε) − Φ(x)| → 0.

1

An interesting application of the result above is quantification of the power to detect interactions, as we now discuss. Consider a random design setting where X is a scaled version of a
matrix Z whose rows zi are independent with zi ∼ Np (0, Σ) and Σjj = 1 for all j. That is we have
√
Xk = nZk /kZk k2 . Let fi = zTi,S Θzi,S where without loss of generality, Θ ∈ Rs×s is a symmetric
matrix. Thus the nonlinear component of the signal is a quadratic function of the variables in
S. As before, we will consider asymptotics where n, p → ∞ and now also Σ ∈ Rp×p will change
as p → ∞, though we suppress this in the notation. We will assume a restricted eigenvalue-type
condition on the sequence of covariance matrices Σ: let
 1/2

kΣ uk2
φ0 (ξ) = inf
: u ∈ C (ξ, S)
kuk2
and assume that for ξ > (A1 + 1)/(A1 − 1), we have φ0 (ξ) > φ > 0 as n → ∞. Note that this is
weaker than assuming the minimum eigenvalue of Σ is bounded away from zero, for example.
Theorem 7. Suppose n1/3 E(f12 ) → 0 and s log(p)/n1/3 → 0. We have
q
 p
sup P(f T R̂/kf k2 ≤ x|Z) − Φ x − kf k2 / σ 2 + kf k22 /n → 0,
β∈B,x∈R

sup

p

|P(f T R̂∗ /kf k2 ≤ x|ε, Z) − Φ(x)| → 0.

β∈B,x∈R

Note that the theorem allows for E(kf k22 ) = nE(f12 ) → ∞, though we do need E(f12 ) → 0. We
see that the nonlinear signal is positively correlated with the true Lasso residuals, but not with
the bootstrap residuals. Thus the nonlinear signal is present in the true Lasso residuals, and in
principle can be detected by a suitable RP method.

B

Non-Gaussian errors

Although the null hypothesis that the Gaussian linear model (1) is correct is often of interest,
one may wish to consider a larger null hypothesis that allows for non-Gaussian errors. Theorem 3
cannot easily be extended to this setting as it allows for arbitrary (collections of) RP functions to
be used, including those that might directly test for normality. We do not pursue this further here
but note that knowing errors are non-Gaussian can be helpful for designing a different objective
function to use with `1 penalisation that may be more efficient for estimation. We also note that
one could in principle extend the results of Theorems 4 and 5 to allow for non-Gaussian error
distributions under the null. The result for a single variable follows via the central limit theorem,
but the uniformity of variables in S c requires the deep results of Chernozhukov et al. [2014].
Nevertheless, it is desirable that a test for e.g. nonlinearity should not reject more often when a
sparse linear model with non-normal errors holds. When non-Gaussian errors must be included in
the null hypothesis, we recommend taking the simulated errors ζ (b) to be a sample with replacement
from the original scaled residuals R̂.
Figures 7 and 8 are identical to Figures 3 and 4 but with exponential errors rather than Gaussian
errors used in all simulations. Similarly Figures 9 and 10 use t-distributed errors with 3 degrees of
freedom scaled to have variance 1. We use the nonparametric bootstrap approach described above.
We see that type I error is very well controlled for RP tests across all the settings with the power
also competitive.
2

False positive rate
True positive rate
True positive rate

False positive rate

Figure 7: Testing significance of groups with exponential errors; the interpretation is similar to
that of Figure 3.

Figure 8: Testing for nonlinearity with exponential errors; the interpretation is similar to that of
Figure 4.

3

False positive rate
False positive rate
False positive rate

False positive rate

Figure 9: Testing significance of groups with t-distributed errors with 3 degrees of freedom; the
interpretation is similar to that of Figure 3.

Figure 10: Testing for nonlinearity with t-distributed errors with 3 degrees of freedom; the interpretation is similar to that of Figure 4.

4

False positive rate
True positive rate

Figure 11: Testing significance of groups when coefficients are chosen according to (11); the interpretation is similar to that of Figure 3.

C

Additional numerical results

In this section we present additional numerical results of the same format as those in Section 5
in the main paper, but where the nonzero coefficients given active variables S = {j1 , . . . , j12 } are
chosen as follows:
1
βjk ∝ √ .
(11)
k
The coefficients are then scaled to have an `1 -norm of 12. Modulo these modifications, Figures 11
and 12 are exactly analogous to Figures 3 and 4 respectively.
We see the results are very much in line with those of Section 5. Note that the reduced power
of RP tests compared to the debiased Lasso in the equal correlation and real design settings of
Figure 11 is due to the poor calibration of the debiased approach, which here tends to greatly
exceeds its nominal level.

D

Uses of RP tests and interpretation of p-values

One use of RP tests is as a form of reassurance that Lasso-based inference is safe to use on the
data at hand: a large p-value indicates a lack of evidence for the Lasso having poor performance.
A small value on the other hand suggests a sparse linear model is inappropriate.
As we explain in this work however, the RP testing framework is general enough to include
tests for significance of groups or individual predictors. It is important to note though that as with
all p-values, low values in these cases can only indicate inadequacy of the null models considered;
extra assumptions are required to draw conclusions, such as that a variable is significant, on the
basis of having observed a low p-value.
5

False positive rate
True positive rate

Figure 12: Testing for nonlinearity; the interpretation is similar to that of Figure 4.

E

Proofs

In the proofs that follow, we will let c1 , . . . , c4 ≥ 0 denote constants, which may change from line
to line.

E.1

Proof of Theorem 2

In the following we suppress the dependence of β̂ λ (β, σζ) on λ for notational simplicity. We know
that every Lasso solution β̂(β, σζ) is characterised by the KKT conditions
1 XT [X{β − β̂(β, σζ)} + σζ]
√
= λν̂,
n kX{β − β̂(β, σζ)} + σζk2
where kν̂k∞ ≤ 1 and ν̂ Ŝ = sgn(β̂ Ŝ (β, σζ)) with Ŝ = {j : β̂j (β, σζ) 6= 0}.
Now picking a particular Lasso solution β̂(β, σζ) in the case where it is not unique, let
β̃(ζ) = β̌ +

σ̌
{β̂(β, σζ) − β}.
σ

Note that when ζ ∈ Λλ,t , the upper bound on σ̌ ensures that we have sgn(β̃(ζ)) = sgn(β̂(β, σζ)).
Next observe that
σ̌
σ̌
X{β̌ − β̃(ζ)} + σ̌ζ = − X{β̂(β, σζ) − β} + σ̌ζ = [X{β − β̂(β, σζ)} + σζ],
σ
σ
so

1 XT [X{β̌ − β̃(ζ)} + σ̌ζ]
1 XT [X{β − β̂(β, σζ)} + σζ]
√
=√
= λν̂(β, ζ).
n kX{β̌ − β̃(ζ)} + σ̌ζk2
n kX{β − β̂(β, σζ)} + σζk2
6

(12)

But this shows that β̃(ζ) satisfies the KKT conditions for β̂(β̌, σ̌ζ). Since Lasso fitted values
are unique, we must have Xβ̃(ζ) = Xβ̂(β̌, σ̌ζ). Now substituting into (12) finally shows that
R̂λ (β̌, σ̌ζ) = R̂λ (β, σζ) as required.

E.2

Results from Sun and Zhang [2012]

The proofs of Theorems 3–7 make use of Theorem 2 and Corollary 1 in Sun and Zhang [2012].
We re-state a subset of these results here for convenience. We have modified the notation in Sun
and Zhang [2012] in order to avoid clashes with our own notation. Furthermore, we have replaced
the sign-restricted cone invertibility factor F2 (ξ, S) (equation 21 in Sun and Zhang [2012]) with its
lower bound F2 (ξ, S) ≥ φ2 (ξ)/(1 + ξ) [Zhang and Zhang, 2012].
Consider the linear model setup of (1) though without any assumptions on the distribution of
ε initially. For ξ > 1 and λ > 0, define


kβ T c k1
λ|T |/{2(1 − ν)}
µ(λ, ξ, β, X) = (ξ + 1) min inf max
,
(13)
, 2
T 0<ν<1
ν
κ {(ξ + ν)/(1 − ν), T, X}
√
where the minimum is over T ⊂ {1, . . . , p} and T c = {1, . . . , p} \ T . Further let σ̃ = ky − Xβk2 / n
√
and τ 2 = τ 2 (σ̃, λ, ξ, β, X) = λµ(σ̃λ, ξ, β, X)/σ̃. Writing β̂ = β̂ λ (β, σε), define σ̂ = ky − Xβ̂k2 / n.
Theorem 8 (Theorem 2 and Corollary 1 in Sun and Zhang [2012]). Let


T
ξ−1
n kX ζk∞
2
Λ1 = ζ ∈ R : √
≤ (1 − τ )λ
.
ξ+1
nkζk2

(14)

When ε ∈ Λ1 ,
max(1 − σ̂/σ̃, 1 − σ̃/σ̂) ≤ τ 2 ,
kβ̂ − βk1 ≤ µ(σ̃λ, ξ, β, X)/(1 − τ 2 )
√
2σ̃ξ sλ
kβ̂ − βk2 ≤
.
(1 − τ 2 )φ2 (ξ)
Lemma 9. Let n, p, m ∈ N with n ≥ m ≥ 3. Let ε ∼ Nm (0, I). Let a ∈ Rm and suppose
0 < η < pe−(n−m)−2 . Then
P(aT ε/kεk2 >

p

2 log(p/η)/n kak2 ) ≤

1 (1 + rm )η
p
p π log(p/η)

where rm → 0 as m → ∞.
T
Proof. We
p follow the proof of part (ii) of Theorem 2 of Sun and Zhang [2012]. Let z = a ε/(kak2 kεk2 ).
Then z/ (1 − z 2 )/(m − 1) follows a t-distribution with m − 1 degrees of freedom. The only change
we need to make in the aforementioned proof is to note that

η < pe−(n−m)−2 ⇒ m − 2 − log(p/η) ≥ n − 2 log(p/η).
and modify (A8) in Sun and Zhang [2012] appropriately.

7

E.3

Proof of Theorem 3

p
Without loss of generality assume S = {1, . . . , s}. Let λ0 = 2 log(p/η)/n and let P denote the
orthogonal projection on to XS . Define the following subsets of Rn :
√
|XTj (I − P)ζ|
≤ λ0 k(I − P)Xj k2 / 2
Λ1 = ζ : max
j
k(I − P)ζk2


√
|XTj Pζ|
≤ λ0 kPXj k2 / 2
Λ2 = ζ : max
j
kζk2
√
√
Λ3 = {ζ : kζk2 / n ≤ 2}.




Let Λ = Λ1 ∩ Λ2 ∩ Λ3 and let Ω = {ε ∈ Λ}. Lemma 9 shows that on the event Ω we have the
following properties.
√
(i) sgn(β 0S ) = sgn(β S ) and minj∈S βj0 /βj > 1 − 1/(2 2).
(s)

(ii) Ŝ (s) = S and sgn(β̌ ) = sgn(β S ).
√
(s)
(iii) 2 2σ minj∈S β̌j /βj > σ̌ (s) .
(s)

In the following, we suppress dependence on λ. Let β̌ ∈ Rp be β̌ with p − s zeroes added:
(s)
β̌ = (β̌ , 0, . . . , 0) and let σ̌ = σ̌ (s) . Note that by Theorem 2, on Ω we know that R̂(β̌, σ̌ζ) =
R̂(β, σζ) for all ζ ∈ Λ. Moreover, Lemma 11 shows that conditional on Ω, R̂(β, σε) and Xβ̌/σ̌ are
independent. Write R̂(b) = R̂λ (β̌, σ̌ζ (b) ) for b = 1, . . . , B and let R̂(0) = R̂λ (β, σε). We see that
(0)
(b)
B
conditional on Ω, {R̂(b) }B
b=0 are independent. Also R̂ 1{ε∈Λ} , {R̂ 1{ζ (b) ∈Λ} }b=1 are independent
and identically distributed.
Let Q̃b , b = 0, . . . , B be derived as in Section 3 by applying
the function Q̃ to appropriate
PB
1
.
Recall
that
Q
=
functions of scaled residuals {R̂(b) }B
b=0 {Q̃b ≥Q̃0 } /(B + 1). Let R = {b :
b=0
ζ (b) ∈ Λ} and note that letting δ = P(Ωc ) we have |R| ∼ Bin(1 − δ, B). We have
P(Q ≤ x) ≤ P(Q ≤ x, Ω) + δ
= (1 − δ)E{P(Q ≤ x|R, Ω)} + δ.
Lemma 10 gives the required bound on δ; it only remains to show the first term on the RHS is at
most x. From the above, conditional on R and the event Ω, {Q̃P
0 , {Q̃b }b∈R } are exchangeable. Now
there can be at most bx(B + 1)c values b ∈ {0, . . . , B} with b0 6=b 1{Q̃ 0 ≥Q̃b } /(B + 1) ≤ x. This
b
entails that
(|R| + 1)P(Q ≤ x|R, Ω) ≤ bx(B + 1)c
x(B + 1)
P(Q ≤ x|R, Ω) ≤
.
|R| + 1

8

Therefore


B+1
(1 − δ)E{P(Q ≤ x|R, Ω)} ≤ x(1 − δ)E
|R| + 1
 
B
X
r+1 B−r B B + 1
=x
(1 − δ) δ
r r+1
r=0


B+1
X
r B+1−r B + 1
=x
(1 − δ) δ
r
r=1

= x(1 − δ B+1 ) ≤ x.

E.4

Proofs of Theorems 4 and 5

The proofs of Theorems 4 and 5 rest on the following decomposition of Tk and an analogous one
for Tk∗ :
Tk =
where

σ
1
1
kWk k2 βk + Zk + δk
σ̂k
σ̂k
σ̂k

(15)

√
σ̂k = ky − X−k Θ̂k k2 / n.

The first term in (15) is zero for k ∈ N and
WkT
X−k (Θk − Θ̂k ),
kWk k2
WkT ε
Zk =
∼ N (0, 1).
kWk k2
δk =

(16)
(17)

The main term we have to control is δk . The result of Sun and Zhang [2012] shows that kΘk − Θ̂k k1
is small with high probability. Next appealing to the KKT conditions for the square-root Lasso,
we obtain
√
kXT−k Wk k∞ /kWk k2 ≤ nγ.
(18)
Hölder’s inequality then gives
|δk | ≤

√

nγkΘk − Θ̂k k1 ,

(19)

which forms the basis of the proofs.
E.4.1

Proof of Theorem 4

In view of Lemma 14 applied with Fn simply a constant, for the first part we need only show that
p
p
supβ∈B,k∈N δk → 0 and supβ∈B,k∈N |σ̂k − σ| → 0. First note that
sup
β∈B,k∈N

µ(λ, ξ, β −k , X−k ) ≤ sup µ(λ, ξ, β, X) → 0
β∈B

in view of Lemma 12 and consequently (as clearly λ → 0)
sup
β∈B,k∈N

τ 2 (σ̃, λ, ξ, β −k , X−k ) ≤ sup τ 2 (σ̃, λ, ξ, β, X) → 0.
β∈B

9

Now writing τ 2 = supβ∈B τ 2 (σ̃, λ, ξ, β, X) for convenience, we see that for n sufficiently large it
p
must be the case that (1 − τ 2 )λ(ξ − 1)/(ξ + 1) > 2 log(p)/n. By Theorem 8 there is a sequence
of events with probability tending to 1 on which we have


σ̂k
σ̃
sup max 1 − , 1 −
≤ τ 2,
(20)
σ̃
σ̂k
β∈B,k∈N
(ξ + 1)σ̃λs
,
(1 − τ 2 )κ2 (ξ, S)
β∈B,k∈N


σ̌
σ̃
≤ τ 2,
sup max 1 − , 1 −
σ̃
σ̌
β∈B
sup

kΘ̂k − Θk k1 ≤

(21)
(22)

√
p
where σ̃ = σkεk2 / n. From Lemma 13 we know that σ̃ → σ, so in particular we have supβ∈B,,k∈N |σ̂k −
p

σ| → 0. Thus also, on a sequence of events with probability tending to 1, applying (19) to (21) we
have
p
log(p)λs
→ 0,
sup |δk | ≤ c1 σ 2
κ (ξ, S)
β∈B,k∈N
which completes the proof of the first part.
Turning to the bootstrap results, we know that that on a sequence of events of the form Ωn =
{ε ∈ ∆n } with probability tending to 1, we have
√
√
sup
nγkΘ̂k − Θk k1 ≥ sup
nγkΘ̂k,Nk k1 → 0 and sup |σ̌ − σ| → 0.
β∈B,k∈N

β∈B,k∈N

β∈B

Here Nk corresponds to the set of noise components of Θk . Thus on Ωn , by Lemma 12, we have
√
sup
nγµ(c1 λ, ξ ∗ , Θ̂k , X−k ) → 0
(23)
β∈B,k∈N

for any fixed c1 > 0 and some ξ ∗ > (A1 + 1)/(A1 − 1). Let δk∗ and σ̂k∗ be the bootstrap equivalents
of δk and σ̂k respectively. By applying Lemma 14 now with Fn = ε, we see that it is enough to
show that on Ωn , for all η > 0 we have
P(

sup
β∈B,k∈N

|δk∗ | > η|ε) → 0 and P(

sup
β∈B,k∈N

|σ̂k∗ − σ| > η|ε) → 0.

For this it is sufficient to exhibit a sequence Ω∗n = {ε∗ ∈ ∆∗n } whose probability tends to 1 such
that on Ωn ∩ Ω∗n we have
sup

|δk∗ | → 0

(24)

σ| → 0.

(25)

β∈B,k∈N
sup |σ̂k∗ −
β∈B,k∈N

√
Let σ̃ ∗ = σ̌kε∗ kn / n. Define
τ∗2 =

sup

τ 2 (σ̃ ∗ , λ, ξ ∗ , Θ̂k , X−k ).

β∈B,k∈N

10

p
Provided τ∗2 → 0 we have (1 − τ∗2 )λ(ξ ∗ − 1)/(ξ ∗ + 1) > 2 log(p)/n for n sufficiently large. This
gives us the equivalent of (20) and (21) for n sufficiently large:


σ̂k∗
σ̃ ∗
(26)
sup max 1 − ∗ , 1 − ∗ ≤ τ∗2 ,
σ̃
σ̂k
β∈B,k∈N
∗

sup
β∈B,k∈N

kΘ̂k − Θ∗k k1 ≤ c2

sup

µ(σ̃ ∗ λ, ξ ∗ , Θ̂k , X−k ).

(27)

β∈B,k∈N
p

By Lemma 13, conditional on ε, supβ∈B |σ̃ ∗ − σ̌| → 0. Thus τ∗2 tending to 0 and therefore also (26)
and (27) occur on a sequence of events with probability tending to 1, Ω∗n ∩ Ωn , where Ω∗n is of the
form {ε∗ ∈ ∆∗n }. As on Ω∗n ∩ Ωn , supβ∈B |σ̌ − σ| → 0, (22) gives (25). Applying (19) to (27) and
(23) then gives (24), which completes the proof.
E.4.2

Proof of Theorem 5

The proof proceeds similarly to that of Theorem 4. For the first result, we use Lemma 14 with
p
p
Fn simply constant. Thus it suffices to show supβ∈Bk |σ̂k − σ| → 0, supβ∈Bk |δk | → 0 and
q
p
supβ∈Bk |βk |kWk k2 |σk0 −1 − σ̂k−1 | → 0 where σk0 =
σ 2 + kWk k22 βk2 /n. To this end, first note
that by Lemma 12, for some ξ 0 > (A1 + 1)/(A1 − 1) we have log(p) supβ∈Bk µ(λ, ξ 0 , Θk , X−k ) :=
log(p)µk → 0. Let
1
1
σ̃k = √ ky − X−k Θk k2 = √ kσε + βk Wk k2 .
n
n
Then
σ̃k2



1
2
2
2
2
=
σ kεk2 + βk kWk k2 + 2σβk kWk k2 Zk
n

√
p
where Zk is defined as in (17). Since βk kWk k2 / n → 0, we have that σ̃k → σ by Lemma 13. For
later use we also note that
√
√
√
2
n(σ̃k2 − σk0 ) = n(σ 2 kεk22 /n − 1) + 2σβk kWk k2 Zk / n = OP (1).
(28)
√
p
by the central limit theorem and as βk kWk k2 / n → 0. Thus we have τk2 := supβ∈Bk τ 2 (σ̃k , λ, ξ 0 , Θk , X−k ) →
0. Note that by (18),
r
1
1
log(p)
T
βk kX−k Wk k∞ ≤ √ βk γkWk k2 = o
.
n
n
n
Therefore
 T

kX−k (σε + βk Wk )k∞ /n
ξ−1
2
lim P
≥ (1 − τk )λ
n→∞
σ̃k
ξ+1
 T

kX−k εk∞ /n p
√
≤ 2 log(p)/n = 1.
≥ lim P
n→∞
kεk2 / n
Thus by Theorem 8 we have that on a sequence of events Ωn = {ε ∈ ∆n } with probability tending
√
to 1, |σ̃k − σ| → 0, n supβ∈Bk |σ̂k − σ̃k | → 0, and log(p) supβ∈Bk kΘ̂k − Θk k1 → 0. The latter in
11

p

p

conjunction with (19) shows supβ∈Bk |δk | → 0. That supβ∈Bk |βk |kWk k2 |σk0 −1 − σ̂k−1 | → 0 follows
from (28).
Now we derive the result concerning the bootstrap test statistic. Note that on Ωn ,
√
√
nγ sup kΘ̂k,T c k1 ≤ nγ{kΘk,T c k1 + sup kΘ̂k − Θk k1 } → 0,
(29)
β∈Bk

β∈Bk

√
and so also nγ supβ∈Bk µ(c1 λ, ξ ∗ , Θ̂k , X−k ) → 0 for any fixed c1 > 0 and some ξ ∗ > (A1 +
1)/(A1 − 1). The rest of the proof then proceeds exactly as the proof for the bootstrap statistic in
Theorem 4.

E.5

Proof of Corollary 6

The proof is essentially identical to that of Theorem 5, but with βk Wk replaced by wγ and Θk
replaced by β + ψ γ .

E.6

Proof of Theorem 7

By Corollary 1 of Raskutti et al. [2010],
p with probability tending to 1 we have κ(ξ) ≥ φ(ξ) > φ/8 >
E(f12 ) = o(n−1/6 ) we see that kf k2 = oP (n1/3 ). Indeed, we
0. By Lemma 15 and the
fact
that
p
have P(kf k2 /n1/3 ≤ n1/6 2E(f12 )) → 1 as n → ∞. Also, with probability tending to 1,
p
p
p
log(p) 1
log(p)
1
T
−1/6
√
kX f k∞ ≤ c1
log(p)/n),
·
kf
k
≤
c
·
o(n
)
=
o(
2
1
n
n
n1/3
n1/3
p
i.e. kXT f k∞ /n = oP ( log(p)/n). We now temporarily make the dependence of Z and p on n
explicit by writing Z(n) and pn respectively, in order to explain the structure of the argument to
(n) ∈ Λ ) → 1,
follow. From the above, we have a sequence of sets Λn ⊆ Rn×pn for
n
p which P(Z
1
1/3
T
and on the respective events kf k2 = o(n ) and n kX f k∞ = o( log(p)/n) (uniformly). Let
p
√
σ̂ = ky − Xβ 0 k2 / n. In relation to Corollary 6, we will take γ = c1 log(p)/n1/3 sufficiently large
such that the Lasso regression of f on X produces the zero vector. By Lemma 14, it suffices to
show that for each η > 0,

sup P sup |σ̂ − σ| > η | Z(n) → 0
Z(n) ∈Λn

β∈B


sup P sup kX f k∞ kβ − β 0 k1 /kf k2 > η | Z(n) → 0
T

Z(n) ∈Λn

β∈B


sup P sup kf k2 |σ 0−1 − σ̂ −1 | > η | Z(n) → 0.
Z(n) ∈Λn

(30)

β∈B

p
Here σ 0 := σ 2 + kf k22 /n. The remainder of the argument to arrive at the first result is essentially
identical to that in Theorem 5, with f playing the role of βk Wk , and with the probabilities being
conditional on Z(n) . The only difference is that in place of (28) (which leads to the equivalent of
(30)), we have
n1/3 |σ̃ 2 − σ 02 | = |n1/3 σ 2 (kεk22 /n − 1) + 2σf T ε/n2/3 |,

12

√
where σ̃ = kσε + f k2 / n. Since for all Z(n) ∈ Λn , kf k2 = o(n1/3 ) uniformly, it is straightforward
to show that

sup P sup kf ||2 |σ 0 − σ̃| > η | Z(n) → 0.
Z(n) ∈Λn

β∈B

for any η > 0, which then leads to (30).
The bootstrap result is simpler. We know there is a sequence of events depending only on (Z, ε)
on which
p
sup kβ − β̌k1 ≤ c2 σs log(p)/n/φ2 (ξ).
β∈B

Thus on the same sequence of events
p
√
sup n · c1 log(p)/n1/3 · µ(λ, ξ, β̌, X) → 0,
β∈B

from which the result follows by arguing along the lines of the second part of the proof of Theorem 4.

E.7

Technical lemmas

Lemma 10. Consider the setup of Theorem 3 and its proof. Recall that


√
|XTj (I − P)ζ|
≤ λ0 k(I − P)Xj k2 / 2 ,
Λ1 = ζ : max
j
k(I − P)ζk2


√
|XTj Pζ|
Λ2 = ζ : max
≤ λ0 kPXj k2 / 2 ,
j
kζk2
√
√
Λ3 = {ζ : kζk2 / n ≤ 2},
Λ = Λ1 ∩ Λ2 ∩ Λ3 ,
p
with λ0 = 2 log(p/η)/n. On the event Ω = {ε ∈ Λ} we have the following properties.
√
(i) sgn(β 0S ) = sgn(β S ) and minj∈S βj0 /βj > 1 − 1/(2 2).
(s)

(ii) Ŝ (s) = S and sgn(β̌ ) = sgn(β S ).
√
(s)
(iii) 2 2σ minj∈S β̌j /βj > σ̌ (s) .
Furthermore

2(1 + rn−s )η
P(Ω) ≥ 1 − p
− e−n/8 .
π log(p/η)

(31)

Proof.pFirst we bound P(Ω). From Lemma 9 and the union bound we have P(ε ∈ Λ2 ) ≤ (1 +
rn )η/ π log(p/η). Next note that
XTj (I − P)ε
aT ζ
d
=
k(I − P)εk2 k(I − P)Xj k2
kζk2 kak2
d

n−s and = denotes equality in distribution. Thus Lemma 9 gives P(ε ∈ Λ ) ≤
where a ∈ R
1
p
(1 + rn−s )η/ π log(p/η). Finally, Lemma 13 gives P(ε ∈ Λ3 ) ≥ 1 − e−n/8 .

13

Now turning to (i), observe that as k(I − P)ζk2 ≤ kζk2 and
√
√
(kPXj k2 + k(I − P)Xj k2 )/ 2 ≤ kXj k2 = n,
we have


Ω⊆

Since


kXT εk∞
√
≤ λ0 .
nkεk2

1 (ξ + 1)λs
≥ µ(λ, ξ, β, X),
2 κ2 (ξ, S)

(8) ensures that (1 − τ 2 )λ(ξ − 1)/(ξ + 1) ≥ λ0 . Note also that (8) gives τ 2 ≤ 1/5. Thus by
Theorem 8, the fact that ε ∈ Λ3 and (9) we have
√
2σ̃ξ sλ
σ̃
1
0
kβ − βk2 ≤
≤
min |βj | ≤ √ min |βj |,
(32)
2
2
(1 − τ )φ (ξ)
4σ j∈S
2 2 j∈S
√
where σ̃ = kεk2 / n. This shows sgn(β 0S ) = sgn(β S ). Next
min
j∈S

βj0
kβ 0 − βk2
1
≥1−
≥1− √ ,
βj
minj∈S |βj |
2 2

which shows (i). Now
min |βj0 | − maxc |βj0 | ≥ min |βj | − max |βj − βj0 | − maxc |βj0 |
j∈S
j∈S
j∈S
j∈S
j∈S
√
> min |βj | − 2kβ − β 0 k2 > 0.
j∈S

Thus Ŝ (s) = S. By Theorem 8, the `2 bound (32) is also satisfied by β̌
Also note that σ̌ (s) ≤ σ̃/(1 − τ 2 ) by Theorem 8. From (32) we have

(s)

, which then shows (ii).



(s)
√ σ
β̌j
σ
σ̃
2 2 (s) min
≥ 2 2 (s) 1 −
4σ
σ̌ j∈S βj
σ̌
σ̃
σ̃
≥ 2 (s) − √
.
σ̌
2σ̌ (s)
√
√
2
We
√ see that√the RHS is at least 1 when τ ≤ ( 2 − 1)/(2 2 − 1), which then shows (iii) as
( 2 − 1)/(2 2 − 1) < 1/5.
√

Lemma 11. Consider the setup of Theorem 3 and its proof. Conditional on the event Ω, R̂(β, σε)
(s)
and (β̌ , σ̌ (s) ) are independent.
Proof. Write vN (ζ) = N k(I − P)ζk2 . First we claim that on Ω, R̂(β,
√ σε) depends only on (I −
P)ε/v1 (ε). To this end, we argue that for ζ ∈ Λ with vN (ζ) > 1/(2 2),
R̂(β, σζ) = R̂(β + σδ(ζ), σ(I − P)ζ)
= R̂(β, σ(I − P)ζ)

(33)

= R̂(β, σ(I − P)ζ/vN (ζ)),

(34)

14

where δ(ζ) = ((XTS XS )−1 XTS ζ, 0, . . . , 0) ∈ Rp . Note the validity
of the above inequalities would
√
prove the initial claim √
for all ζ ∈ Λ with vN (ζ) > 1/(2 2) for each fixed N . However, since
∪N >0 {ζ : vN (ζ) > 1/(2 2)} = Rn the initial claim would then have to be true for all ζ ∈ Λ. We
now set about proving these equalities. The first equality is clear by definition of R̂. Turning to
the second (33), let ζ ∈ Λ and denote by k · k the operator norm. From Lemma 10 we have
σkδ(ζ)k∞ = σk(XTS XS )−1 XTS ζk∞
√
√
√
≤ σ(kζk2 / n)k(XTS XS )−1 k skXTS Pζk∞ /(kζk2 / n)
p
2σ s log(p/η)/n
≤
≤ min |βj |/10,
j∈S
φ2 (ξ)
√
using the facts that ξ > 1 and A > 2 and (9) in the final √
line. Note that by Lemma 10,
Λ ⊆ Λλ,t with the latter defined in Theorem 2 and t = 1 − 1/(2 2). Now clearly if ζ ∈ Λ then
also (I − P)ζ ∈ Λ. An application of Theorem 2 then gives the desired equality (33). Equality
√
(34) then follows from a further application of Theorem 2 noting that vN (ζ) < 1/(1 − t) = 2 2 by
assumption.
(s)
Next we examine β̌ and σ̌ (s) . Note that on Ω, the former is simply the Lasso estimate from
(s)
regressing on XS and the latter is the resulting normalised root-RSS. Write β̌ = β̌ and σ̌ = σ̌ (s) .
The least squares part of the Lasso objective decomposes as
kXS β S + ε − XS bk2 = {kXS β S + Pε − XS bk22 + k(I − P)εk22 }1/2 .
Thus it is clear that the fitted values XS β̌ do not depend on (I − P)ε/v1 (ε). This then implies
that σ̌ does not depend on (I − P)ε/v1 (ε) since it is determined by kεk22 = kPεk22 + k(I − P)εk22
and Pε.
Now observe that

Pε ⊥
⊥ k(I − P)εk2 , (I − P)ε/k(I − P)εk2
k(I − P)εk2 ⊥
⊥ (I − P)ε/k(I − P)εk2 ,
so Pε, k(I − P)εk2 , (I − P)ε/k(I − P)εk2 are jointly independent. Let E1 = {R̂(β, σε) ∈ B1 },
E2 = {(XS β̌, σ̌) ∈ B2 }, where B1 ⊆ Rn and B2 ⊆ Rn+1 are arbitrary Borel sets. Let Ωk = {ε ∈
Λk }, k = 1, 2, 3. From the above
P(E1 , E2 |Ω1 , Ω2 , Ω3 ) =

P(E1 , Ω1 ) P(E2 , Ω2 , Ω3 )
= P(E1 |Ω1 )P(E2 |Ω2 , Ω3 ).
P(Ω1 )
P(Ω2 , Ω3 )

The conditional probabilities on the RHS remain unchanged if we modify the conditioning event
to be (Ω1 , Ω2 , Ω3 ) since E1 ⊥
⊥ (Ω2 , Ω3 ) and E2 ⊥⊥ Ω1 . This completes the proof.
Lemma 12. Given a sequence of collections of matrices Mk,n ∈ Rn×pn , k ∈ Nn , n = 1, 2, . . .,
suppose
there exists a sequence of collections of sets Sk,n ⊆ {1, . . . , pn }, tuning parameters λn =
p
A log(p)/n and ξ > c for constants A, c > 0 such that the follow holds:
p
|Sk,n | log(pn )2 /n
sup 2
→ 0.
k∈Nn κ (ξ, Sk,n , Mk,n )
15

Moreover, suppose the collection of sequences β k,n ∈ Rpn is such that
p
sup log(pn )kβ n,S c k1 → 0.
k,n

k∈Nn

Then there exists a ξ 0 > c such that
function µ is defined in (13)).

p
log(pn )µ(c1 λn , ξ 0 , β k,n , Mk,n ) → 0 for any c1 > 0 (where the

Proof. Let ξ 0 = (ξ + c)/2 > c and let ν 0 > 0 be given by (ξ 0 + ν 0 )/(1 − ν 0 ) = ξ. Then we have
p
sup log(pn )µ(c1 λn , ξ 0 , β k,n , Mk,n )
k∈Nn

 kβ

c k1
p
k,n,Sk,n
c1 λn |Sk,n |/{2(1 − ν 0 )}
≤(ξ + 1) log(pn ) max max
→ 0.
,
k∈Nn
ν0
κ2 (ξ, Sk,n , Mk,n )
0

Lemma 13. Let Zn ∼ χ2n . We have the following tail bounds [Boucheron et al., 2013, pg. 29]:
√
P(Zn > n + 2 nγ + 2γ) ≤ e−γ ,
whence taking γ = n/8,
p
√
√
P( Zn /n > 2) < P[Zn > n{1 + 2(1/ 8 + 1/8)}] ≤ e−n/8 .
Lemma 14. Let Tk,n ∈ R, k ∈ Nn , n = 1, 2, . . . be a collection of random variables which we may
decompose as
Tk,n = ak,n Zk,n + bk,n
where each Zk,n is identically distributed with continuous distribution function F . Suppose further
that each Zk,n is independent of the random elements Fn , and for all δ > 0, supk∈Nn P(|bk,n −dk,n | >
p
δ|Fn ) → 0 for some random variables dk,n ∈ R that are functions of Fn , and supk∈Nn P(|ak,n − 1| >
p
δ|Fn ) → 0. Then
p
sup sup |P(Tk,n ≤ x|Fn ) − F (x − dk,n )| → 0.
k∈Nn x∈R

Proof. First note that for all  > 0, there exists δ > 0 such that for all c1 , c2 ∈ [−δ, δ],


x + c1
sup F
− F (x) < .
1 + c2
x∈R

(35)

Indeed, the function G(x, c1 , c2 ) := F {(x + c1 )/(1 + c2 )} is uniformly continuous on R × [−δ 0 , δ 0 ]2
for 0 < δ 0 < 1 sufficiently small (since F is uniformly continuous and compositions of uniformly
continuous functions are continuous) and the LHS of (35) is supx∈R |G(x, c1 , c2 ) − G(x, 0, 0)|.
Hence given  > 0, let δ > 0 be such that the LHS of (35) is at most /2. Then
sup |P(Tk,n ≤ x|Fn ) − F (x − dk,n )| = sup |P(Zk,n ≤ (x + dk,n − bk,n )/ak,n |Fn ) − F (x)|
x∈R

x∈R

≤ sup

sup

|P(Zk,n ≤ (x + c1 )/(1 + c2 )|Fn ) − F (x)|

x∈R c1 ,c2 ∈[−δ,δ]

+ P(|bk,n − dk,n | > δ|Fn ) + P(|ak,n − 1| > δ|Fn )
≤ /2 + P(|bk,n − dk,n | > δ|Fn ) + P(|ak,n − 1| > δ|Fn ).
16

Thus
P{ sup sup |P(Tk,n ≤ x|Fn ) − F (x − dk,n )| > }
k∈Nn x∈R

≤ P{ sup P(|bk,n − dk,n | > δ|Fn ) > /4} + P{ sup P(|ak,n − 1| > δ|Fn ) > /4}
k∈Nn

→0

k∈Nn

as n → ∞.

Lemma 15. Consider the setup of Theorem 7. There exists constants c1 , c2 > 0 such that for all
t < n1/4
2

P(1 − tn−1/4 ≤ kf k22 /E(kf k22 ) ≤ 1 + tn−1/4 ) ≤ c1 e−c2 t ,
p


log(p)
1 kXT f k∞
≤ c1
→ 0.
P √
n kf k2
n1/3

(36)
(37)

Proof. Although in Theorem 7 certain assumptions are placed on E(f12 ), since the probabilities
above do not depend on E(f12 ) here we may assume E(f12 ) = 1. Fix j ∈ {1, . . . , p}. Using the
1/2
eigendecomposition PDPT = ΣS,S S, S 1/2 (where D is diagonal and P is orthogonal) and writing
d

= for equality in distribution, we have
−1/2

−1/2

d

(f1 , Z1j ) = (zT1,S Bz1,S , Z1j ) = (zT1,S ΣS,S PDPT ΣS,S z1,S , Z1j ) = (uT Du, v)
where (u, v) is multivariate Gaussian with u ∼ Ns (0, I) and Var(v) = 1. Let the diagonal entries
of D be θ.
First we show (36). Note that
n X
2 o X
X
E(f12 ) = E
u2j θ
=
θj2 E(u4j ) +
θj θk E(u2j )E(u2k )
j

j

= 2kθk22 +

X

j6=k
2
θj

= 1.

(38)

j

√
P
Thus in particular, kθk∞ ≤ 1/ 2, kθk22 ≤ 1/2 and
j θj ≤ 1.
Now with a view to applying Lemma 16 below, consider
E exp

n1  X
4

u2j θj

2

1/2 o

≤ E exp

X


 X

u2j θj /4 + E exp −
u2j θj /4

j

j

=

Y

(1 − θj /2)−1/2 +

j

j

Y

(1 + θj /2)−1/2

j

using the fact that u2j ∼ χ21 . Note that
(1 − θj /2)−1 = 1 +

θj2
θj
θj
+
≤1+
+ θj2 .
2
4(1 − θj /2)
2

17

Thus, by the AM–GM inequality we have
Y
Y
(1 − θj /2)−1 ≤
(1 + θj /2 + θj2 )
j

j

 X
s
1
2
≤
(1 + θj /2 + θj )
s
j

≤ (1 + 1/s)s < e,
using (38). Similarly,

Q

j (1

+ θj /2)−1 < e. Putting things together, we see that

√
E exp(|f12 − E(f12 )|1/2 /4) ≤ e1/4 E exp(|f1 |/4) < e1/4 2e.

(39)

Lemma 16 then immediately gives (36).
To show (37), we first obtain a tail bound for |f T Zj |/n. Observe that if (w1 , w2 , w3 ) is multid

variate Gaussian with zero-mean, then since (w1 , w2 , w3 ) = −(w1 , w2 , w3 ), we have
E(w1 w2 w3 ) = E{(−w1 )(−w2 )(−w3 )} = −E(w1 w2 w3 ) = 0.
Then, by Hölder’s inequality
E exp(|f1 Z1j |2/3 /4) = E exp

 X

vu2j θj

2/3 


4

j

 X
∞
X
E
≤
u2j θj /4
≤

≤

r=0
∞
X

r=0

2r/3

|v/2|

.
r!

j
2

r

1/3

[E{(v /4) }]

r=0
∞
X

2r/3

  X

r  2/3 .
2
E
uj θj /4
r!
j

∞
 X
X
r .
E
u2j θj /4
r! +
E{(v 2 /4)r }/r!
r=0

j
2

= E exp(|f1 |/4) + E(v /4) ≤

√

2e +

√
2,

√ p
√ p
using (39) in the final line. Note kf k2 / n → 1 by (36), and kZj k/ n → 1 by Lemma 13. Thus
√
by Lemma 16, P(|f T Xj |/( nkf k2 ) ≥ t) ≤ c1 exp(−c2 n2/3 t2 ) for t ∈ [0, 1] and some constants
c1 , c2 > 0. Thus for c3 with c2 c23 − 1 > 0,
p
√
2
P(kZT f k∞ /( nkf k2 ) ≥ c3 log(p)/n1/3 ) ≤ c1 p exp(−c2 c23 log(p)) = c1 p−(c2 c3 −1) → 0
as p → ∞ and log(p)/n2/3 → 0.
Lemma 16 (Lemma B.4 of Hao and Zhang [2014]). Let W1 , . . . , Wn be independent random variables with zero mean such that E(exp(c1 |Wi |α )) ≤ c2 for constants c1 , c2 > 0 and α ∈ (0, 1]. Then
there exist constants c3 , c4 > 0 such that for t ∈ [0, 1],

n
1 X
P
Wi ≥ t ≤ c3 exp(−c4 nα t2 ).
n


i=1

18

References
S. Boucheron, G. Lugosi, and P. Massart. Concentration inequalities: A nonasymptotic theory of
independence. OUP Oxford, 2013.
V. Chernozhukov, D. Chetverikov, and K. Kato. Central limit theorems and bootstrap in high
dimensions. arXiv preprint arXiv:1412.3661, 2014.
N. Hao and H. H. Zhang. Interaction screening for ultrahigh-dimensional data. J. Am. Statist.
Ass., 109(507):1285–1301, 2014.
G. Raskutti, M. J. Wainwright, and B. Yu. Restricted eigenvalue properties for correlated gaussian
designs. The J. Mach. Learn. Res., 11:2241–2259, 2010.
T. Sun and C.-H. Zhang. Scaled sparse linear regression. Biometrika, 99(4):879–898, 2012.
C.-H. Zhang and T. Zhang. A general theory of concave regularization for high-dimensional sparse
estimation problems. Statistical Science, 27(4):576–593, 2012.

19

