Electronic Journal of Statistics
Vol. 11 (2017) 2519–2546
ISSN: 1935-7524

Optimal Two-Step
Prediction in Regression

arXiv:1410.5014v3 [stat.ME] 6 Jun 2017

Didier Chételat
Department of Decision Sciences
HEC Montréal
3000, chemin de la Côte-Sainte-Catherine
Montréal, Canada
e-mail: didier.chetelat@hec.ca

Johannes Lederer
Departments of Statistics and Biostatistics
University of Washington
Box 354322
Seattle, WA 98195-4322
e-mail: ledererj@uw.edu

Joseph Salmon
LTCI, CNRS, Télécom ParisTech,
Université Paris-Saclay,
75013, Paris, France
e-mail: joseph.salmon@telecom-paristech.fr
Abstract: High-dimensional prediction typically comprises two steps: variable selection and subsequent least-squares refitting on the selected variables. However, the standard variable selection procedures, such as the
lasso, hinge on tuning parameters that need to be calibrated. Cross-validation,
the most popular calibration scheme, is computationally costly and lacks
finite sample guarantees. In this paper, we introduce an alternative scheme,
easy to implement and both computationally and theoretically efficient.
MSC 2010 subject classifications: Primary 62G08; secondary 62J07.
Keywords and phrases: High-Dimensional Prediction, Tuning Parameter Selection, Lasso.

Contents
1
2

Introduction . . . . . . . .
AVPr and Its Properties .
2.1 The AVPr Algorithm
2.2 Assumption 2.1 . . .
2.3 Oracle Inequality . .
3 Experiments . . . . . . . .
3.1 General Setup . . . .
3.2 Practical Choice of a
3.3 Choice of the Tuning

. . . . . . . . . .
. . . . . . . . . .
. . . . . . . . . .
. . . . . . . . . .
. . . . . . . . . .
. . . . . . . . . .
. . . . . . . . . .
. . . . . . . . . .
Parameter Grids
2519

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

2520
2523
2523
2523
2525
2527
2527
2528
2529

Chételat, Lederer, Salmon/Two-Step Prediction

3.4 Computational and Statistical Performance
4 Discussion . . . . . . . . . . . . . . . . . . . . . .
Acknowledgements . . . . . . . . . . . . . . . . . . .
References . . . . . . . . . . . . . . . . . . . . . . . .
A Proofs . . . . . . . . . . . . . . . . . . . . . . . .
A.1 Definitions . . . . . . . . . . . . . . . . . . .
A.2 Lemmas . . . . . . . . . . . . . . . . . . . .
A.3 Proofs of Results From Section 2 . . . . . .
A.4 Proof of Theorem 2.2 . . . . . . . . . . . . .
B Description of the lslassoBIC . . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

2520

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

2529
2529
2530
2531
2534
2534
2534
2539
2541
2545

1. Introduction
Variable selection has become a basic tool for estimating linear models on large
data sets. The most popular method for variable selection is the lasso [37], which
minimizes the sum of squares errors under an `1 -penalty. Although efficient at
selecting variables when properly tuned, the lasso has the disadvantage that all
coefficients are shrunk towards zero. To mitigate this bias, practitioners typically
rely on a two-step estimation of the coefficients by computing a least-squares
estimate on the variables selected by the lasso.
For illustration, consider the leukemia micro-array data set of [16], which
consists of n = 38 bone marrow samples analyzed with p = 7129 probes from
several thousand human genes. A particular interest is to predict the type of
leukemia (AML or ALL) present in a patient. The data set also contains an
independent test set of 34 observations that are used for assessment of the
predictive performance.
In this problem, there are many more variables (7129 features) than available observations (38 samples), and in such a context, a least-squares fitting
is not appropriate. A standard solution is to perform variable selection using
the lasso, with tuning parameter chosen by cross-validation on the prediction
loss. However, since the lasso is known to involve a bias, practitioners commonly
refit a least-squares estimate on the selected variables. If the lasso tuning parameter is chosen using 10-fold cross-validation, this approach, called lassoCV
in the following, yields a prediction risk of 0.36 on the test set, computed in 463
seconds.
Although common among practitioners, this approach is suboptimal, because
the cross-validation does not take into account the least-squares refitting. Another alternative is to tune the cross-validation for the entire two-step procedure.
On the test set, this approach with 10-fold cross-validation, called lslassoCV in
the following, yields a prediction risk of 0.45 computed in 499 seconds.
This adjusted approach is natural, yet suffers from two drawbacks. First,
every cross-validation fold must fit a least-squares on each subset selected on
the lasso path, which becomes computationally intensive once larger data sets
are considered. Second, the method does not come with theoretical guarantees,
an issue shared by most cross-validation procedures.

Chételat, Lederer, Salmon/Two-Step Prediction
Lasso on the Leukemia dataset, p = 7129, n = 38, r = 100

0.5

500
Timing in seconds

0.4
0.3
0.2
0.1

400
300
200
100

C
BI
so

=
10

las

V
so
as
lsl

las

so

(C

(C

V

as

so

=

BI

10

)

)

C

r
lsl

BI
so

las

(C
so
as
lsl

las

so

(C

lsl

V

V

=

=
10

10

)

)

C
BI
so
as

AV
p

C

0

r

0.0

Lasso on the Leukemia dataset, p = 7129, n = 38, r = 100

AV
p

Prediction Error

0.6

2521

Fig 1. Prediction error and computing times of the AVPr , lslassoBIC [2], lassoCV, lslassoCV, and lassoBIC procedures. The bars represent the prediction error on the 34 left out
observations. Note that the grid for the tuning parameter contains the same 50 values for all
methods.

To address these problems, we propose Adaptive Validation for Prediction,
(AVPr ), a novel variable selection scheme. A pseudo-code description of the
algorithm is given as Algorithm 1. Our proposal is closely related to the recently
introduced `∞ -Adaptive Validation (AV∞ ) scheme [12], which is based on tests
inspired by isotropic versions of Lepski’s method [11, 24, 25]. This approach has
been shown to provide fast and optimal calibration of the lasso for (one-step)
estimation with respect to `∞ -loss. For the two-step prediction considered in this
paper, however, a considerably different and more technical approach inspired
by non-isotropic tests is required.
As a practical example, Figure 1 compares AVPr and standard methods on
the Leukemia dataset discussed above. The methods under consideration to
select the lasso tuning parameter are AVPr , 10-fold cross-validation, Bayesian
Information Criterion (lassoBIC), and an estimator obtained by selecting with
BIC a least-square estimator over the supports geneterated by a the lasso path
(lslassoBIC) following [2] (see the Appendix for further information about the
implementation of the latter approach). As can be seen, AVPr is faster (39 seconds) to compute than cross-validation, and it is nearly as fast as the lassoBIC
(40 seconds) and lslassoBIC (56 seconds). At the same time, it rivals the predictive performance of all competing approaches. Note at this point that the
lassoBIC is a variable selection method rather than a predictive method; two
goals that can be considerably different from each other.
The organization of this article is as follows. In the next section, we introduce
the algorithm and prove that AVPr predictions satisfy an oracle inequality, that
is, are optimal up to a constant factor. In Section 3, we show that on simulations,
AVPr is substantially faster than cross-validation while being comparable in
accuracy.
All proofs are deferred to the Supplementary Material.

Chételat, Lederer, Salmon/Two-Step Prediction

2522

Framework and Notation
Let us describe the framework and the notation. We are interested in linear
regression models of the form
Y = Xβ + ε,

(1.1)

where Y ∈ Rn is the data, X ∈ Rn×p the design matrix, β ∈ Rp the regression
vector, and ε ∈ Rn the random noise. For ease of exposition, we assume that
the noise is Gaussian with unknown variance σ 2 , that is ε ∼ N (0, σ 2 ). We
n
assume that the columns of the design matrix X
√1 , . . . , Xp ∈ R have been
standardized to have Euclidean norm kXj k2 = n, but we otherwise allow
for arbitrary correlations between the columns and noise distributions. We are
mainly motivated by (but not limited to) high-dimensional settings with sparse
regression vectors, where the number of parameters p can rival or even exceed the
number of samples n. We finally denote by S := supp[β] := {j ∈ [p] : βj 6= 0}
the true support, whose cardinality is usually smaller than n and p (where
throughout the paper [d] stands for the set {1, . . . , d}).
A standard approach to find a vector β̂ with small prediction loss kX β̂ −
Xβk22 /n is performing a least-squares refitting to the lasso. After reducing the
initially large set of variables to a small number of relevant ones, the subsequent
refitting aims to lessen the bias associated with the lasso. For a fixed tuning
parameter λ, the lasso β̂ λ is defined via the minimization of objective function
n
o
2
β̂ λ ∈ arg min kY − Xθk2 + 2λkθk1 .
(1.2)
θ∈Rp

For simplicity, we will assume that the support of the minimizer equals the
equicorrelation set (see Supplementary Material for details). The subsequent
least-squares refitting is defined as a minimizer of
β̄ λ ∈

2

arg min
λ

kY − Xθk2 .

(1.3)

supp[θ]= supp[β̂ ]

We call this estimator the least-squares lasso (lslasso). This two-step procedure
is very popular as it has smaller bias than the lasso for a range of models [3, 21].
Our goal is to find optimal tuning parameters for the lslasso (1.3) in terms of
prediction. In practice, only finitely many estimators can be computed. Therefore, we consider finite sets of tuning parameters Λ = {λ1 , . . . , λr }, r ∈ N and
the associated supports (Ŝ 1 , . . . , Ŝ r ), Ŝ i := supp[β̂ λi ]. We denote the collection of supports by S := {Ŝ i : i ∈ [r]}. Finally, we introduce surrogate sets
Ŝ i,j := Ŝ i ∪ Ŝ j and corresponding estimators
β

i,j

∈ arg min kY − Xξk22 .

(1.4)

supp[ξ]⊂Ŝ i,j
i

In the special case i = j, it holds that Ŝ i,j = Ŝ i , and hence, β := β

i,i

λi

=β .

Chételat, Lederer, Salmon/Two-Step Prediction

2523

The lasso is only one out of many variable selection procedures. Our algorithms and derivations can be easily adapted to other procedures, such as
the square-root lasso [1, 4, 9], scaled-lasso variants [28, 34, 35] or thresholded
ridge regression [33], combined with subsequent least-squares refitting. Note for
instance that by one-to-one correspondence, our results also hold for the squareroot lasso. However, due to its popularity, we focus here only on the lasso.
Related Literature
Besides the references to the papers that are most closely connected with our
study, we provide some additional pointers to related literature. A discussion of
multi-stage methods for regression can be found in [41]. Approaches to tuning
parameter calibration in the single-stage setting include [8, 10, 15, 22, 27, 30, 32].
Related papers that appeared recently include [40], which contains an alternative
to least-squares refitting, and [2], which discusses BIC-type selection as well as
Q-aggregation approaches to model selection over the lasso path. The latter
contains sparse oracle inequalities as well as prediction and estimation bounds
under the standard restricted eigenvalue condition [5] - both for (a refitted)
BIC-type procedure and for a Q-aggregation procedure. These methods enjoy
similar theoretical guarantees as the ones we provide for AVPr , and they are
also subject to the same issue, namely, that a preliminary estimate of the noise
level is required.
2. AVPr and Its Properties
2.1. The AVPr Algorithm
The AVPr scheme is summarized in Algorithm 1. As inputs, it takes the data
(Y, X), a set of tuning parameters Λ, and a constant a > 0 specified in the
following section. It then conducts simple tests along the tuning parameter path
of the lasso until a stopping criterion is met. It returns the index of the current
i

tuning parameter i as well as the corresponding two-stage estimator β .
The algorithm requires the computation of a single lasso path and leastsquares estimators along this path. The computation of the paths can be conducted with readily available, easy-to-use, and highly efficient software such as
glmnet (in R) or scikit-learn(in Python) [14, 29]. For the computation of the
least-squares estimators, off-the-shelf solvers can be used since the number of
active variables of the second step is typically small.
In Section 2.3, AVPr is shown to satisfy an optimal finite sample prediction
bound, and the practical performance of AVPr is illustrated in Section 3.
2.2. Assumption 2.1
Let us first introduce and motivate an assumption that ensures a certain stability
of lasso solution. In general, if an estimator is unstable for data very close to

Chételat, Lederer, Salmon/Two-Step Prediction

2524

Data: Y, X, Λ = {λ1 , . . . , λr }, a
Result: ī ∈ [r], β ∈ Rp
Initialize index: i ← 1
1
r
Compute β , . . . , β
If needed, re-sort the estimators such that |Ŝ 1 | ≤ · · · ≤ |Ŝ r |
while i ≤ r − 1 do
Initialize stopping criterion: T estF ailure ← F alse
Initialize comparisons: j ← i + 1
while (j ≤ r) and (TestFailure==False) do
i,j
Compute Ŝ i,j and β
i
i,j 2
if kXβ − Xβ k2 ≤ a|Ŝ i | + a|Ŝ i,j | then
j ←j+1
else
T estF ailure ← T rue
if T estF ailure == T rue then
i←i+1
else
break
ī
Set output: ī ← i and β ← β

Algorithm 1: AVPr

the (noiseless) underlying truth, accurate estimation and inference hardly seem
realistic. For the goal of refitting, we thus introduce an assumption that ensures
the stability of supports. In the case of the lasso, this means that we restrict Xβ
from being too close to hyperplanes generated by the geometric arrangement
of the columns in X. Figure 2 contains a schematic picture of this: Xβ needs
to lie outside of small neighborhoods (depicted in orange) around the black
boundaries that represent the geometry of X. Most importantly, we stress the
assumption does not imply restrictions on the correlations of the design, and
does not require estimated supports to be accurate.
While the assumption concerns the model, it is most convenient to put the
precise formulation in terms of the lasso itself. For this, recall that for a fixed X,
the support of the lasso evaluated at a vector z ∈ Rn is determined by which
“zone” of Rn the vector z falls into [18, 38]. These zones exactly correspond to
the zones in Figure 2 that are separated by the black lines. Importantly, note
that we do not require additional variable selection guarantees for the lasso, but
merely that the selection is unambiguous. We now define
D : Rn → [0, ∞)

√
D(z) := inf{kz − z 0 k∞ / n : z 0 ∈ Rn s.t. for some λ ∈ Λ,
0

supp[β̂ λ (z)] 6= supp[β̂ λ (z 0 )] for all λ0 ∈ Λ} .
The function √
D quantifies how far away a signal can be from the zone boundaries.
The factor 1/ n in the definition reflects our normalization of the design matrix.
We also stress that the function involves lasso solutions only at fixed, nonrandom vectors z, z 0 ; in particular, D is independent of ε and Y .

Chételat, Lederer, Salmon/Two-Step Prediction

2525

Fig 2. An illustration of Assumption 2.1: Xβ needs to be separated from the boundaries of
the zones that determine the active set of the lasso.

Assumption 2.1. There is an integer N such that for all n ≥ N , it holds that
r
6σ 2 log n
D(Xβ) >
.
n
This assumption now ensures that Xβ is sufficiently far from the zone boundaries. Note that the assumption is very different from restricted eigenvalues [7]
or similar hypothesis in the theory for the lasso [13, 39]. While the latter assumptions need to be strict to ensure a good performance of the lasso, our assumption only requires that the estimates are unambiguous. In the specific case
where X = In×n , some insight can be obtained, since the quantity D(β) can be
computed. Indeed, with the convention that |β(s) | is the s-th largest amplitude
of the vector |β|, D(β) represents the smallest difference |β(s) | − |β(s+1) |, where
s is a support size of a Lasso solution applied on β (i.e., a soft-thresholded
version of β) for a threshold λ ∈ Λ. In this case, the assumption represents
bounding by below such differences, which makes it clear to be an assumption
on the underlying signal itself.
To motivate this assumption further, we finally show that a slightly weaker
version of Assumption 2.1 automatically holds for all Xβ up to a set of measure
zero.
Theorem 2.2. For all Xβ ∈ Rn up to a set of Lebesgue measure zero, the lasso
satisfies
D(Xβ) > 0 .
Theorem 2.2 does not completely exclude cases that violate Assumption 2.1.
However, together with the above discussion, the result indicates that these
cases are hardly generic and of limited relevance in applications.
2.3. Oracle Inequality
Oracle inequalities are bounds for the risk of an estimator. More precisely, they
compare the risk of an estimator with the risk of an oracle estimator, an estimator that has knowledge of the best model [6, 20].

Chételat, Lederer, Salmon/Two-Step Prediction

2526

In this section, we show that our estimator AVPr satisfies such an oracle
inequality. To this end, we first introduce the oracle set.
∗

Definition 2.1 (Oracle). The oracle set S ∗ ∈ S is the set S ∗ := Ŝ i with index

i∗ := min i ∈ [r] : Ŝ i ⊃ S .
i∗

and the associated oracle estimator is β ∗ := β .
In other words, the oracle set contains the true support S and has minimal
cardinality among all such sets. The oracle set can therefore be viewed as the
best possible approximation of S in S; in particular, S ∗ = S whenever S ∈ S.
We implicitly assume that the oracle set exists, that is, the true support set
is a subset of an estimated support along the path. However, one can easily
generalize the definition to avoid this assumption. Let S ∗ be an arbitrary set
and P ∗ the projection onto the space spanned by the columns with indexes
in S ∗ . Adding this projection in our proofs (cf. (A.2) for example) yields the
same results as below except for an additional term k(I −P ∗ )βk22 in the bounds.
However, as the above definition exists in generic cases (since the lasso supports
tend to be very exhaustive for small tuning parameters), and as it provides a
concise formulation of the results, we do not consider the extended version in
the following.
Now, given the oracle, we can state a bound for the two-step lasso procedure
with the optimal tuning parameter, that is, the tuning parameter that leads to
the oracle set. Throughout this section we invoke Assumption 2.1, which helps
us rule out ambiguous design settings.
Proposition 2.1. Under Assumption 2.1, for any α > 0, there exist constants
t, N, R > 0 such that for all n ≥ N and r ≥ R, the oracle estimator satisfies
with probability at least 1 − α the bound
σ 2 |S ∗ |
kXβ ∗ − Xβk22
≤ (1 + t log r)
.
n
n
This is a bound for the lasso with refitting - under the assumption that the oracle
set S ∗ is known and incorporated in the selection of the tuning parameter. The
constants t, N, and R are specified in the proof section.
In practice, we do not have access to the oracle set S ∗ . Therefore, we hope
to find a procedure that does not require its knowledge and still satisfies the
bound (up to constants) stated in Proposition 2.1. The following result shows
that AVPr provides this.
Theorem 2.3 (Oracle inequality for AVPr ). If Assumption 2.1 holds, for any
α > 0, there exist constants t, N, R > 0 such that for all n ≥ N and all r > R,
our estimator AVPr with a ≥ 2σ 2 (1 + t log r) satisfies with probability at least
1 − α the bounds

Chételat, Lederer, Salmon/Two-Step Prediction

|Ŝ| ≤ |S ∗ |
and

kXβ −
n

Xβk22

≤

2527

(i)
h

i |S ∗ |
6a + 4σ 2 (1 + t log r)
.
n

(ii)

This proves optimality of AVPr : indeed, if a & 2σ 2 (1 + t log r), AVPr satisfies
the same bound (up to constants) as the two-step approach that is based on the
knowledge of the oracle set S ∗ . Explicit constants can be found in the proofs
section, though we did not attempt to optimize them.
Theorem 2.3 holds for any sufficiently large a. The question is now how to
choose a in practice. Theorem 2.3 entails precise guidance for this choice. In view
of the bounds, one should select the smallest a that is still allowed, that is, a =
2σ 2 (1+t log r). However, since σ 2 is typically unknown in practice, we suggest to
replace it with a rough estimate σ̂. Moreover, we argue that the term 2(1+t log r)
is an artifact of our proof technique rather than a fundamental aspect of the
bound. We thus suggest the simple choice a = σ̂ 2 , see the empirical section
below. Consequently, the bounds above provide a solid theoretical foundation
for AVPr ; however, there is still a gap between theory and practice that deserves
to be studied further.
We note that our approach is very different from just replacing the unknown
noise variance in the existing theoretical tuning parameters. Standard oracle
inequalities for
pthe lasso hold true with probability t for tuning parameters of the
form constt σ (log p)/n, where constt is a factor involving the level t, see [6] and
references herein. Thus, one might be tempted to use these tuning parameters
with an estimate of σ. However, the above form is valid only for Gaussian noise,
while we aim at more general calibration. Moreover, even for Gaussian noise,
the above form is known to be suboptimal both in the near orthogonal case
(because p could be replaced by p/s, where s is the true sparsity level) and in
the correlated case (where much smaller tuning parameters might be favored),
we refer to [6, 13, 19] and references therein. Finally, even if the above form were
optimal in terms of the standard oracle inequalities for prediction, estimation,
and variable selection, there are no guarantees on their performance in terms of
refitting.
3. Experiments
3.1. General Setup
We measure the numerical performance and the computational speed of AVPr
in two-step prediction. The methods of comparison are cross-validation with 2,
5, 10, and 20 number of folds, which are typically regarded as the standard
calibration schemes.
Variable selection is performed with the lasso. We emphasize that the motivation of this work is not to compare different variable selection methods, but
instead, to compare different calibration schemes in two-step prediction.

Chételat, Lederer, Salmon/Two-Step Prediction

2528

Data: Y, X, δ
Result: σ
b
√
Initialize tuning parameter and variance: λ0 ← 2n log p and σ
b←1
repeat
Save σ
b0 ← σ
b
Update σ
b:
Set λ ← σ
bλ0
Compute β̂ λ as the lasso with tuning parameter λ according to (1.2)
√
Set σ
b ← kY − X β̂ λ k2 / n
0
until |b
σ−σ
b|≤δ

Algorithm 2: Scaled lasso algorithm with early stopping, cf. [35]

The data are generated according to a linear regression model as in (1.1)
with n = p = 100, 200. The first 10 entries of the regression vector β are set
to 1, while all other entries are set to 0. The components of the noise vector
are independently sampled from a univariate standard normal distribution with
mean 0 and variance 1. The rows of the design matrix X are independently
sampled from a multivariate normal distribution with mean 0 and covariance
matrix Σ that is set to Σij = 1 for i = j and to Σij = ρ for i 6= j with
√ ρ = 0.5.
Subsequently, the columns of X are normalized to Euclidean norm n. For all
experiments, we perform 50 repetitions.
In addition to the described parameter settings, we tested various other settings, including different correlation coefficients ρ, regression vectors β, and
tuning parameter grids. As the conclusions were similar across all settings, we
restrict our presentation to the ones described. All computations are conducted
with the standard implementations of the lasso from Python scikit-learn (version
0.16) [29], and our code is available at https://github.com/josephsalmon/AVp.
3.2. Practical Choice of a
We follow the suggestions after Theorem 2.3. Specifically, if σ 2 is known, we
recommend using Algorithm 1 with a = σ 2 as suggested by Theorem 2.3 (regarding the term t log r as a superfluous term coming from our proof technique).
In practice, however, the noise variance σ 2 is often unknown. We then advocate
using a = σ
b2 with a rough estimate σ
b2 of σ 2 . Such a rough estimate can be easily
obtained by using a (very) small number of iterations of the algorithms for the
square-root lasso [9] or the scaled lasso [36]. For our simulations, we have opted
for the latter, which consists of an alternating minimization for estimating both
the regression parameter and the noise level. Algorithm 2 states our concrete
implementation. We set the tolerance to δ = 10−2 , which typically leads to less
than five iterations of the loop in the algorithm and therefore, as illustrated
below, to very low computational costs.

Chételat, Lederer, Salmon/Two-Step Prediction

2529

3.3. Choice of the Tuning Parameter Grids
The tuning parameter grid is chosen as a default grid in the lasso function in
scikit-learn. More precisely, we take a geometric grid of size r = 100 starting
from λmax = kX > Y k∞ , the smallest tuning parameter that leads to a lasso
solution of all zeros, and ending at λmax /1000.
3.4. Computational and Statistical Performance
We first report in Figure 3 the computational times for each of the following:
lasso path:
lslasso path:
lassoAVPr :

lslassoBIC:

lassoCV:

Computation of one tuning parameter path of
lasso.
Computation of one tuning parameter path of
lasso with least-squares refitting.
Least-squares refitted lasso with tuning parameter selected by AVPr with a = σ
b2 as detailed
above.
Least-squares refitted lasso with tuning parameter selected by a BIC-type procedure [2], detailed in Appendix B.
Least-squares refitted lasso with tuning parameter selected by cross-validation on the estimates
of the (one-step) lasso.

lslassoCV:

Least-squares refitted lasso with tuning parameter selected by cross-validation on the estimates
of least-squares refitted lasso.
We then also report in Figure 4 the prediction performances of the last three
methods.
In conclusion, our simulations demonstrate that AVPr is competitive both in
computational speed and in prediction performance.
4. Discussion
The standard scheme for calibrating the lasso is cross-validation. However, crossvalidation entails two main deficiencies: it is computationally inefficient and
lacks finite sample guarantees. In contrast, AVPr is fast and satisfies optimal
bounds for prediction with the refitted lasso. We therefore propose AVPr as an
alternative scheme for prediction with lasso followed by refitting. Moreover, our
work can be readily extended to the square-root lasso and to ridge regression
with thresholding.
A direction for further research could be sharpening the theoretical bounds.
The current result contains a term that grows logarithmically in the number
of tuning parameters under consideration. Such artifacts are common in the

2530

CV
20

as
so

CV
10
lsl

5

as
so

CV
lsl

2
as
so
lsl

lsl

as
so

CV

BI
C
so

lsl

as

AV
p

h

so
las

lsl

as
s

o

pa
t

pa
t
so
las

r

p = 200, n = 200, s = 10, σ = 1, ρ = 0.5, r = 100

90
80
70
60
50
40
30
20
10
0

h

Timing in seconds

Chételat, Lederer, Salmon/Two-Step Prediction

Fig 3. Computation times of the lasso path, lslasso path, lslassoCV, and lassoAVPr (with σ
b).
Cross-validation is performed using a refitting step (lslassoCV) for different numbers of folds.

non-parametric literature. In a sense, one might consider our approach as a
non-parametric version of [11]. Thus, it is not surprising that such a term appears. Empirically, however, there are no indications that this term is needed.
An improved understanding of the bound could especially strengthen the connections between the theoretical and the practical choice of a.
The exact specification of the method is an issue that appears more generally in tuning parameter calibration. In our case, there is flexibility in how to
estimate the noise variance; in cross-validation, one has to specify the number
of folds; in BIC-type approaches, the constant in front of the log penalty needs
to be adjusted; in Q-aggregation, there is a trade-off between the KL regularization and the quadratic term. We believe that although out of the scope of
this contribution, a comprehensive sensitivity analysis studying the selections
in each of the methods would be of interest.
Acknowledgements
We thank Michaël Chichignoud for the insightful remarks and for the inspiring
discussions, and we thank Pierre Bellec for providing us with valuable insights
about numerical aspects of his work. We also thank the editors and reviewers
for the valuable comments that have improved the paper.

Chételat, Lederer, Salmon/Two-Step Prediction

0.8

2531

Lasso for Equi-correlations scenario: p = 100, n = 100, s = 10, σ = 1, ρ = 0.5, r = 100

Prediction Error

0.7
0.6
0.5
0.4
0.3
0.2
0.1

0.6

20
lsl

lsl

as

so

CV
so
as

CV

10

5
CV
lsl

lsl

as

so
as

so

CV

2

10
V
oC
ss

BI
C
so
as

la

la

lsl

ss

oA
V

pr

0.0

Lasso for Equi-correlations scenario: p = 200, n = 200, s = 10, σ = 1, ρ = 0.5, r = 100

Prediction Error

0.5
0.4
0.3
0.2
0.1

20
so
as
lsl

CV
so
as
lsl

CV

10

5
CV
so
as

so
as
lsl

lsl

CV

2

10
V
oC
ss
la

BI
C
so
as
lsl

la

ss

oA
V

pr

0.0

Fig 4. Prediction errors of the lassoCV, lslassoCV, and lassoAVPr (with σ
b). Cross-validation
is performed using (lslassoCV) or not using (lassoCV) a refitting step and using different
numbers of folds.

References
[1] A. Antoniadis. Comments on: `1 -penalization for mixture regression models. TEST, 19(2):257–258, 2010.
[2] P. Bellec. Aggregation of supports along the Lasso path. In COLT, pages
488–529, 2016.
[3] A. Belloni and V. Chernozhukov. Least squares after model selection in
high-dimensional sparse models. Bernoulli, 19(2):521–547, 2013.
[4] A. Belloni, V. Chernozhukov, and L. Wang. Square-root Lasso: Pivotal
recovery of sparse signals via conic programming. Biometrika, 98(4):791–
806, 2011.
[5] P. Bickel, Y. Ritov, and A. Tsybakov. Simultaneous analysis of Lasso and
Dantzig selector. Ann. Statist., 37(4):1705–1732, 2009.

Chételat, Lederer, Salmon/Two-Step Prediction

2532

[6] P. Bühlmann and S. van de Geer. Statistics for high-dimensional data.
Springer Series in Statistics. Springer, Heidelberg, 2011. Methods, theory
and applications.
[7] F. Bunea, A. Tsybakov, and M. Wegkamp. Sparsity oracle inequalities for
the Lasso. Electron. J. Stat., 1:169–194 (electronic), 2007.
[8] F. Bunea, Y. She, H. Ombao, A. Gongvatana, K. Devlin, and R. Cohen. Penalized least squares regression methods and applications to neuroimaging.
Neuroimage, 55, 2011.
[9] F. Bunea, J. Lederer, and Y. She. The group square-root Lasso: Theoretical
properties and fast algorithms. IEEE Trans. Inf. Theory, 60(2):1313–1325,
2014.
[10] S. Chatterjee and J. Jafarov. Prediction error of cross-validated lasso.
arXiv:1502.06291, 2015.
[11] M. Chichignoud and J. Lederer. A robust, adaptive M-estimator for pointwise estimation in heteroscedastic regression. Bernoulli, 20(3):1560–1599,
2014.
[12] M. Chichignoud, J. Lederer, and M. Wainwright. Tuning Lasso for supnorm optimality. J. Mach. Learn. Res., 17, 2016.
[13] A. Dalalyan, M. Hebiri, and J. Lederer. On the prediction performance of
the Lasso. Bernoulli, 23(1):552–581, 2017.
[14] J. Friedman, T. Hastie, and R. Tibshirani. Regularization paths for generalized linear models via coordinate descent. J. Stat. Softw., 33(1):1–22,
2010.
[15] C. Giraud, S. Huet, and N. Verzelen. High-dimensional regression with
unknown variance. Statist. Sci., 27(4):500–518, 2012.
[16] T. Golub, D. Slonim, P. Tamayo, C. Huard, M. Gaasenbeek, J. Mesirov,
H. Coller, M. Loh, J. Downing, M. Caligiuri, C. Bloomfield, and E. Lander.
Molecular classification of cancer: class discovery and class prediction by
gene expression monitoring. Science, 286(5439):531–537, 1999.
[17] B. Grünbaum. Convex Polytopes. Springer-Verlag, New York, second edition, 2003.
[18] N. Harris and A. Sepehri. The accessible lasso models. arXiv:1501.02559,
2015.
[19] M. Hebiri and J. Lederer. How correlations influence Lasso prediction.
IEEE Transactions on Information Theory, 59:1846–1854, 2013.
[20] V. Koltchinskii. Oracle inequalities in empirical risk minimization and
sparse recovery problems, volume 2033 of Lecture Notes in Mathematics.
Springer, Heidelberg, 2011.
[21] J. Lederer. Trust, but verify: benefits and pitfalls of least-squares refitting
in high dimensions. arXiv:1306.0113 [stat.ME], 2013.
[22] J. Lederer and C. Müller. Don’t fall for tuning parameters: Tuning-free
variable selection in high dimensions with the trex. In Proceedings of the
Twenty-Ninth AAAI Conference on Artificial Intelligence, 2015.
[23] D. Lee, J.and Sun and Y. Sun. Exact post-selection inference, with applications to the lasso. Preprint arXiv:1311.6238v5, 2015.
[24] O. Lepski. On a problem of adaptive estimation in gaussian white noise.

Chételat, Lederer, Salmon/Two-Step Prediction

2533

Theory Probab. Appl., 35(3):454–466, 1990.
[25] O. Lepski, E. Mammen, and V. Spokoiny. Optimal spatial adaptation to
inhomogeneous smoothness: an approach based on kernel estimates with
variable bandwidth selectors. Ann. Statist., 25(3):929–947, 1997.
[26] J. Mairal and B. Yu. Complexity analysis of the lasso regularization path.
Proceedings of the 29th International Conference on Machine Learning,
2012.
[27] N. Meinshausen and P. Bühlmann. Stability selection. J. Roy. Statist. Soc.
Ser. B, 72(4):417–473, 2010.
[28] A. Owen. A robust hybrid of lasso and ridge regression. Contemporary
Mathematics, 443:59–72, 2007.
[29] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel,
M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay. Scikit-learn:
Machine learning in Python. J. Mach. Learn. Res., 12:2825–2830, 2011.
[30] J. Sabourin, W. Valdar, and A. Nobel. A permutation approach for selecting the penalty parameter in penalized model selection. Biometrics, 71:
1185–1194, 2015.
[31] R. Schneider. Convex bodies: the Brunn–Minkowski theory, volume 151 of
Encyclopedia of Mathematics and its Applications. Cambridge University
Press, second edition, 2013.
[32] R. Shah and R. Samworth. Variable selection with error control: another
look at stability selection. J. Roy. Statist. Soc. Ser. B, 75(1):55–80, 2013.
[33] J. Shao and X. Deng. Estimation in high-dimensional linear models with
deterministic design matrices. Ann. Statist., 40(2):812–831, 2012.
[34] N. Städler, P. Bühlmann, and Sara s van de Geer. `1 -penalization for
mixture regression models. TEST, 19(2):209–256, 2010.
[35] T. Sun and C.-H. Zhang. Scaled sparse linear regression. Biometrika, 99
(4):879–898, 2012.
[36] T. Sun and C.-H. Zhang. Sparse matrix inversion with scaled lasso. J.
Mach. Learn. Res., 14:3385–3418, 2013.
[37] R. Tibshirani. Regression shrinkage and selection via the lasso. J. Roy.
Statist. Soc. Ser. B, 58(1):267–288, 1996.
[38] R. Tibshirani and J. Taylor. Degrees of freedom in lasso problems. Ann.
Statist., 40(2):1198–1232, 2012.
[39] S. van de Geer and P. Bühlmann. On the conditions used to prove oracle
results for the Lasso. Electron. J. Stat., 3:1360–1392, 2009.
[40] X. Wang, D. Dunson, and C. Leng. No penalty no tears: Least squares in
high-dimensional linear models. arXiv:1506.02222, 2015.
[41] L. Wasserman and K. Roeder. High dimensional variable selection. Ann.
Stat., 37(5A):2178, 2009.
[42] G. Ziegler. Lectures on polytopes, volume 152. Springer, 1995.

Chételat, Lederer, Salmon/Two-Step Prediction

2534

Appendix A: Proofs
A.1. Definitions
A subtlety in the definition of the lasso variable selection scheme is that it is
defined as the solution to the minimization problem (1.2), but the solution is
not necessarily unique. Different lasso algorithms can yield different solutions
to the lasso problem, and all could reasonably called the lasso estimator.
Define the lasso equicorrelation set [38] to be
n
o
E[λ; Y ] := i ∈ [p] : |Xi> (Y − X β̂ λ )| = λ .
This set is unique and contains the support of any lasso solution β̂ λ . In common
cases, there is at least one lasso solution β̂ λ whose support equals the equicorrelation set. This property turns out to be quite valuable in our analysis and
consequently, we will always assume from now on that the support set equals
the equicorrelation set.
Denote the sets of lasso outputs that lead to the same sign vectors η ∈
{1, 0, −1}p by W λ (η) := {Y ∈ Rn : sgn [β̂ λ ] = η}, cf. [23]. The closures of the
sets of equal sign vector will be called regions, and the collection of all regions
will be written V = {cl W 1 (η) : η ⊂ {−1, 0, 1}p }. To simplify notation, we
will write the target as ξ = Xβ. When relevant, this will also be written ξn to
emphasize the dependence on n.
A.2. Lemmas
Recall that in convex geometry, a polyhedron is a finite intersection of closed
half-spaces – more details can be found in Appendix A.
Lemma A.1. The lasso β̂ fulfills the following:
1. it is scale-symmetric, in the sense that supp[β̂ λ (Y )] = supp[β̂ 1 (Y /λ)] for
all λ ∈ Λ, Y ∈ Rn , and X ∈ Rn×p ;
2. for all λ ∈ Λ and η ∈ {1, 0, −1}p , the closure of its regions of equal sign
vector cl [W λ (η)] are polyhedra.
Proof of Lemma A.1. We prove each condition in order.
Part i) The scale-symmetry follows from consideration of the dual problem to
(1.2). Let β̂λ be a lasso solution whose active set equals the equicorrelation set
E[λ; Y ]. Let Cλ stand for the polyhedron {x ∈ Rp : kX > xk∞ ≤ λ}, and let
PCλ denote the Euclidean projection on this set. Notice that for any x ∈ Cλ ,
we have x/λ ∈ C1 and therefore
 
 
Y
Y
Y
x Y
− Y k2 ≤ λkPC1
− k2 ≤ λk − k2 ≤ kx − Y k2 .
kλPC1
λ
λ
λ
λ
λ

Chételat, Lederer, Salmon/Two-Step Prediction

2535


Since this is true for all x ∈ Cλ , and that λPC1 Yλ ∈ Cλ , we conclude that

λPC1 Yλ = PCλ (Y ). As shown in [38], the residual from the lasso satisfies
Y − X β̂λ = PCλ (Y ).
Therefore, for any λ > 0 the active set of β̂λ , which is the equicorrelation set
here, satisfies E[λ; Y ] = E [1; Y /λ] since

 


Y
=1 ,
E[λ; Y ] = i ∈ [p] : Xi> PCλ (Y ) = λ = i ∈ [p] : Xi> PC1
λ
as desired.
Part ii) The polyhedron C1 = {x ∈ Rn : kX > xk∞ ≤ 1} has an irreducible
decomposition into half-spaces
!
!
p
p
\
\ \
n
>
n
>
C1 =
{x ∈ R : Xi x − 1 ≤ 0}
{x ∈ R : −Xi x − 1 ≤ 0} ,
i=1

i=1

so the facets of C1 are Fi± = C1 ∩ {x ∈ Rn : ±Xi> x − 1 ≤ 0} – see [17, Sec.
2.6]. Since by assumption, the active set of the lasso estimate coincides with
the equicorrelation set, by the Karush-Kuhn-Tucker conditions and Y − X β̂1 =
PC1 Y , see [38, Equations (13)-(14) and Lemma 3], we have
Xi> PC1 (Y ) = sgn β̂i ∈ {−1, 1} ⇔ PC1 Y ∈ Fisgn β̂i
for i ∈ E[1; Y ]. Moreover,


|Xi> PC1 (Y )| < 1 ⇔ PC1 Y ∈ C1 \Fi+ ∩ C1 \Fi−
for i 6∈ E[1; Y ]. In light of this, sgn β̂ = η if and only if
\
\
\
+
−
PC1 (Y ) ∈
C1 \Fi+ ∩ C1 \Fi−
Fi ∩
Fi ∩
i∈[p]
ηi =1

i∈[p]
ηi =-1

"
= relint

\

i∈[p]
ηi 6=0

!
Fi+

∩
i∈[p]
ηi =-1

i∈[p]
ηi =1

So, W (η) =

PC−1
(relint Fη )
1

Fi−




for the face Fη =

!#

\

T

+
i∈[p] Fi
ηi =1

.


∩



T

−
i∈[p] Fi
ηi =-1
p

.

Let V ∈ V - then by definition of V, there must be an η ⊂ {−1, 0, 1} such that
V = cl W (η) = cl PC−1
(relint Fη ). According to [31, Equation (2.3) and Page 83],
1
we have N (C1 , Fη ) + relint Fη = PC−1
(relint Fη ) = W (η), where N (C1 , Fη ) is the
1
normal cone of C1 to the face Fη .
Now, since cl A + cl B ⊂ cl (A + B) for any sets A, B,
N (C1 , Fη ) + relint Fη ⊂ N (C1 , Fη ) + Fη ⊂ cl PC−1
(relint Fη ).
1

Chételat, Lederer, Salmon/Two-Step Prediction

2536

But N (C1 , Fη ) + Fη is the sum of two polyhedra, hence a polyhedron, hence
closed. Thus N (C1 , Fη ) + Fη = cl PC−1
(relint Fη ) = V is a polyhedron.
1
Notice that N (C1 , Fη ) + relint Fη is convex as a sum of convex sets, and thus
relint W (η) =relint PC−1 (relint Fη ) = relint (N (C1 , Fη ) + relint Fη )
=relint cl (N (C1 , Fη ) + relint Fη ) = relint V = int V,
since V is n-dimensional. Thus int V ⊂ W (η), and sgn β̂ is constant over int V ,
as desired.
Theorem A.2. For any t ≥ 6, there exists
√ a deterministic integer N such that
with probability at least 1 − t−1 − Rn−1 1 + t log Rn ,
X β̄ i,j − Xβ

2
2

≤ σ 2 (1 + t log r)|Ŝ i,j | + (In −PŜ i,j )Xβ

2
2

(A.1)

for all i, j and all n ≥ N .
Proof of Lemma A.2. By (1.1), the loss of these estimators can be broken down
as
Xβ

i,j

2

− Xβ
2

=
=

= PŜ i,j (Xβ

PŜ i,j (XXŜ+i,j Y
PŜ i,j ε

2
2

i,j

2

− Xβ)
2

+ (In −PŜ i,j )(Xβ

2

− Y + ε)

(In −PŜ i,j )(XXŜ+i,j Y

+
2

+ (In −PŜ i,j )Xβ

2
2

i,j

,

2

− Xβ)
2
2

− Xβ)
2

(A.2)

where XŜ+i,j is the Moore-Penrose pseudo-inverse of the submatrix of X comprising the columns with indexes in Ŝ i,j . Our goal is to control the noise term
kPŜ i,j εk22 in (A.2). We do this in four parts. We first show that on an appropriate scale, the response Y must be close to the target Xβ with high probability.
This is then shown to imply, using the scale-symmetry property of the lasso,
that the ordered active sets must be unique with high probability. This allows us
to control the noise term by showing that the projected noise behaves like a chisquare distribution, construct an appropriate event, and bound its probability.
Finally, on this event the inequality of the theorem is shown to hold.
1. We first use a Gaussian tail bound to show that Y is close to the target
Xβ; precisely, we show that the event
Ω̄1 :=

∞ n
o
\
p
kY − Xβk∞ ≤ 6σ 2 log n
n=n1

fulfills the bound
 
P Ω̄1 ≥ 1 − 1/t,
(A.3)
p
√
where n1 := min{n : 6σ 2 log n > 1/ 2π} ∨ d4te.
For this, write ξn := Xβ (the subscript
the sampleosize depenn n highlightsp
dence) and define the event Ωn := kY − ξn k∞ ≤ 6σ 2 log n for ease of

Chételat, Lederer, Salmon/Two-Step Prediction

2537

notation. Using a union bound, the Gaussian tail bound P[N(0, 1) > t] ≤
√
√
2
2
e−t /2 / 2πt < e−t /2 (note that t > 6 > 1/ 2π), and the definition of
n1 , we find that for the complements ΩC
n of the sets Ωn ,
∞
X

∞
∞
X
X
 
1
P ΩC
2ne−3 log n = 2
n ≤
2
n
n=n1
n=n1
n=n1
Z ∞
2
1
2
2
4
1
≤ 2 +2
dw = 2 +
≤
≤ .
2
n1
w
n
n
n
t
1
1
n1
1

From this and the definition of Ω̄1 , the bound (A.3) follows.
2. We now use Part 1 to deduce that the active sets are deterministic if the
response is close to the target, or more specifically, we derive that on Ω̄1 ,
Ŝ i [Y ] = Ŝ i [ξn ]

(A.4)

for 1 ≤ i ≤ r.
From Part 1, we deduce that on Ω̄1 and for n ≥ n1 ,
r
p
6σ 2 log n/n
6σ 2 log n
kY − ξn k∞
√
≤
=
D(ξn ).
n
D(ξn )
n
Then, by Assumption 2.1, there must be an
√ n2 , without loss of generality
satisfying n2 ≥ n1 , such that kY − ξn k∞ / n < D(ξn ). But by definition
of D(ξn ), this means that for all λ > 0, there must be a λ0 > 0 such
that supp[β̂ 1 (Y /λ)] = supp[β̂ 1 (ξn /λ0 )]. For a given y ∈ Rn , we now define
the collection of active sets by Q̂[y] := {supp[β̂ 1 (y/λ)] : λ > 0}. That is,
Q̂[Y ] is the collection of active sets along the tuning parameter path of
the estimator for given data (Y, X). There are at most 2p different subsets
of [p], so these are always finite sets. We therefore obtain
h
i
P Q̂[Y ] = Q̂[ξn ] ; ∀n ≥ n2 : Ω̄1 = 1.
In particular, the random cardinality r := |Q̂[Y ]| and deterministic cardinality r̄ := |Q̂[ξn ]| coincide almost surely on this event, that is,


P r = r̄ ; ∀n ≥ n2 : Ω̄1 = 1.
Because we follow a fixed ordering rule, we must then have
h
i
P (Ŝ 1 [Y ], ..., Ŝ r [Y ]) = (Ŝ 1 [ξn ], ..., Ŝ r [ξn ]) ∀n ≥ n2 : Ω̄1 = 1.
This finishes the proof of Equation (A.4).
3. Let us define the sets Ŝ i,j := Ŝ i,j [Y ] := Ŝ i [Y ] ∪ Ŝ j [Y ], the random ranks
ri,j := rk XŜ i,j and deterministic ranks r̄i,j := rk XŜ i,j [ξn ] . Our next step
is to show that Part 2 provides a chi-square bound for the noise part
in (A.2), that is, we prove on Ω̄1 the relations
kPŜ i,j εk22 ∼ σ 2 χ2r̄i,j if r̄i,j ≥ 1 and kPŜ i,j εk22 = 0 if r̄i,j = 0.

(A.5)

Chételat, Lederer, Salmon/Two-Step Prediction

2538

To show this, we apply (A.2) to our two-step method β̄ i,j to get
X β̄ i,j − Xβ

2
2

2
2

= PŜ i,j ε

+ (In −PŜ i,j )Xβ

2
2

.

According to Part 2, the sets Ŝ i,j satisfy
h
i
P Ŝ i,j = Ŝ i,j [ξn ] ∀n ≥ n2 : Ω̄1 = 1
for all 1 ≤ i, j ≤ r. This has two consequences. First, the random ranks
ri,j equal the deterministic ranks r̄i,j almost surely on the event Ω̄1 :
h
i
P ri,j = r̄i,j ∀n ≥ n2 : Ω̄1 = 1.
Second, the matrices PŜ i,j are indeed projection matrices on Ω̄1 of rank
ri,j , since on this event the active sets Ŝ i,j (and, therefore, the subspaces
spanned by XŜ i,j ) are constant. Formally,
h
i
P PŜ i,j = PŜ i,j [ξn ] ∀n ≥ n2 : Ω̄1 = 1.
Combining this with ε ∼ N (0, σ 2 In ) yields the results in (A.5). We can
now control the noise part in (A.2) with a chi-square Chernoff bound.
More specifically, we obtain the bound
p
(A.6)
P[Ω̄2 ] ≥ 1 − t−1 − Rn−1 1 + t log Rn

for Ω̄2 := Ω̄1 ∩ kPŜ i,j εk22 ≤ (1 + t log r)r̄i,j for all 1 ≤ i, j ≤ r . To this
end, recall that by definition, the integers r̄i,j = |Ŝ i,j [ξn ]| and r̄ = |Q̂[ξn ]|
are deterministic. According to Part 2, it also holds that on Ω̄1 , we have
r̄i,j = |Ŝ i,j [Y ]| and r̄ = r. We use this, result (A.5), a union bound, and
the bound P[Ω̄1 ] ≥ 1 − 1/t stated in (A.3) to deduce
 max 2


 
χr̄i,j
P Ω̄2 =P Ω̄1 ∩ 1≤i,j≤r̄
≤
1
+
t
log
r̄
r̄ i,j ≥1 r̄ i,j
"
#
X
2
i,j
−1
P χr̄i,j > r̄ (1 + t log r̄) .
≥1 − t −
1≤i,j≤r̄
r̄ i,j ≥1

Now, using the chi-square Chernoff bound P[χ2k > k(1 + a)] < [(1 +
a)e−a ]k/2 , we obtain
X 
 
r̄i,j /2
P Ω̄2 ≥ 1 − t−1 −
(1 + t log r̄)e−t log r̄
.
1≤i,j≤r̄
r̄ i,j ≥1

As r̄ ≥ 1, (1 + t log r̄)e−t log r̄ ≤ 1 and so
X 
 
1/2
P Ω̄2 ≥ 1 − t−1 −
(1 + t log r̄)e−t log r̄
1≤i,j≤r̄
r̄ i,j ≥1

≥ 1 − t−1 − r̄2−t/2

p

1 + t log r̄.

Chételat, Lederer, Salmon/Two-Step Prediction

2539

We now use that that 2 − t/2 ≤ −1 for all t ≥ 6 and the fact that r̄ ≥ Rn
to find that
p
 
P Ω̄2 ≥ 1 − t−1 − Rn−1 1 + t log Rn ,
which concludes Part 3.
4. We finally √
collect the pieces to deduce that with probability at least 1 −
t−1 − Rn−1 1 + t log Rn , the bound
X β̄ i,j − Xβ

2
2

≤ σ 2 (1 + t log r)|Ŝ i,j | + (In −PŜ i,j )Xβ

2
2

(A.7)

holds for all n ≥ n2 .
For this, we assume that indeed n ≥ n2 and then combine the initial
bound (A.2) and the results of Part 3 to find that on Ω̄2 ,
X β̄ i,j − Xβ

2
2

≤ σ 2 (1 + t log r)r̄i,j + (In −PŜ i,j )Xβ

2
2

.

Recalling that r̄i,j = |Ŝ i,j | according to Part 2, the desired bound (A.7)
now follows from Inequality (A.6) derived in Part 3.

A.3. Proofs of Results From Section 2
Recall that the path of active sets is |Ŝ 1 | ≤ |Ŝ 2 | ≤ · · · ≤ |Ŝ r |. The cardinality r
is typically random, but we can always bound it almost surely by some constant
Rn , so that 1 ≤ Rn ≤ r almost surely. This constant should be independent of
the data but can be chosen to vary with n and p. For example, we might have
agreed a priori with considering 50 sets, or our variable selection method might
always yield at least min(n, p) different sets by construction.
Corollary A.3 (Oracle benchmark). Say the oracle set exists, that the design
satisfies Assumption 2.1. Then, for any t ≥ 6, there exists
√ a deterministic integer N such that with probability at least 1 − t−1 − Rn−1 1 + t log Rn , the oracle
estimator satisfies
kXβ ∗ − Xβk22 ≤ σ 2 (1 + t log r)|S ∗ |
for all n ≥ N .
i∗

Proof of Corollary A.3. The oracle estimator is β ∗ = β , the refitted estima∗
tor on the oracle set S ∗ = Ŝ i . The result therefore follows immediately from
Theorem A.2 applied to i = j = i∗ .
Corollary A.4 (Oracle inequality for AVPr ). Say the oracle set exists, that
the design satisfies Assumption 2.1, and that the AVp parameter a is such that
a ≥ 2σ 2 (1 + t log r). Then, for any t ≥ 6, there exists a deterministic integer

Chételat, Lederer, Salmon/Two-Step Prediction

2540

√
N such that with probability at least 1 − t−1 − Rn−1 1 + t log Rn , it holds that
|Ŝ| ≤ |S ∗ | and
h
i
kXβ − Xβk22 ≤ 6a + 4σ 2 (1 + t log r) |S ∗ |
for all n ≥ N .
Proof of Corollary A.4. By Theorem A.2, there exists an N such that the event

2
Ω = X β̄ i,j − Xβ 2 ≤ σ 2 (1 + t log r)|Ŝ i,j |

2
+ (In −PŜ i,j )Xβ 2 , ∀i, j and n ≥ N
√
holds with probability at least 1 − t−1 − Rn−1 1 + t log Rn .

Claim (i): On Ω, it holds that ī ≤ i∗ = min i ∈ [r] Ŝ i ⊃ S .
We prove this claim by contradiction and, therefore, assume that ī > i∗ . Then,
by the definition of our estimator, there must be an integer k ∈ [r] such that
|Ŝ k | ≥ |S ∗ | and
kXβ ∗ − Xβ

k,i∗ 2
k2

> a|S ∗ | + a|Ŝ k ∪ S ∗ | .

(A.8)

The fact that |Ŝ k | ≥ |S ∗ | ≥ |S|, together with the bound (A.1) and S ∗ ⊃ S,
yields
kXβ ∗ − Xβ

k,i∗ 2
k2

≤ 2kXβ ∗ − Xβk22 + 2kXβ − Xβ

k,i∗ 2
k2

≤ 2σ 2 (1 + t log r)|S ∗ | + 2k(In −PS ∗ )Xβk22 + 2σ 2 (1 + t log r)|Ŝ k ∪ S ∗ |
+ 2k(In −PŜ k ∪S ∗ )Xβk22
≤ 2σ 2 (1+t log r)|S ∗ | + 2k(In −PS )Xβk22 + 2c|Ŝ k ∪ S ∗ | + 2k(In −PS )Xβk22
= 2σ 2 (1 + t log r)|S ∗ | + 0 + 2σ 2 (1 + t log r)|Ŝ k ∪ S ∗ | + 0 .
Since a ≥ 2σ 2 (1 + t log r), this contradicts (A.8) and, therefore, concludes the
proof of Claim (i).
Claim (ii): On Ω, it holds that kXβ − Xβk22 ≤ (6a + 4σ 2 (1 + t log r))|S ∗ |.
To prove this claim, we note that by Claim 1, we have ī ≤ i∗ and, therefore,
|S̄| ≤ |S ∗ |. Hence, the definition of the estimator implies for ī = r that
kXβ − Xβ

ī,i∗ 2
k2

ī,ī

= kXβ − Xβ k22 = 0

and otherwise, if ī < r, that (recall that S ∗ ⊃ S)
kXβ − Xβ

ī,i∗ 2
k2

≤ a|S̄| + a|S̄ ∪ S ∗ | ≤ 3a|S ∗ | .

Chételat, Lederer, Salmon/Two-Step Prediction

2541

The bound (A.1), on the other hand, yields
kXβ

ī,i∗

− Xβk22 ≤ σ 2 (1 + t log r)|S̄ ∪ S ∗ | + 0 ≤ 2σ 2 (1 + t log r)|S ∗ | .

Combining these two inequalities, we finally obtain
kXβ − Xβk22 ≤ 2kXβ − Xβ

ī,i∗ 2
k2

+ 2kXβ

ī,i∗

− Xβk22

≤ (6a + 4σ 2 (1 + t log r))|S ∗ | ,
which concludes the proof of Claim (ii).
In particular, this yields the results of Section 2.
Proof of Proposition 2.1. Let t = max( α2 , 6), and let R be large enough so that
r
p
α
2 log R
1
−1
−1
t +R
1 + t log R ≤
+
+
≤ α.
2
R2
α R2
Using Corollary A.3 with t and Rn = R gives the result.
Proof of Theorem 2.3. Let t = max( α2 , 6), and let R be large enough so that
r
p
α
2 log R
1
−1
−1
1 + t log R ≤
+
+
≤ α.
t +R
2
2
R
α R2
Using Corollary A.4 with t and Rn = R gives the result.
A.4. Proof of Theorem 2.2
The proof of Theorem 2.2 rely on various convex geometry notions, so we first
remind the reader of some background on the subject.
The affine hull of a set A, denoted aff A is the intersection of all affine spaces
that contain A, or alternatively, the unique affine set of minimal dimension that
contains A. We write, by extension, dim A = dim aff A. We denote the interior,
closure, and boundary of A by int A, cl A, and ∂A, respectively. The relative
interior and boundary of A, denoted relint A and relbd A, are respectively the
interior and the boundary when A is seen as a subset of its affine hull. We write
A ⊂ B if A is a (not necessarily strict) subset of B.
A half-space H + is a set of the form {x ∈ Rn : α> x ≤ b} for α ∈ Rn ,
b ∈ R. Its boundary H = ∂H + = {x ∈ Rn : α> x = b} is aThyperplane in
Rn . A polyhedron is a finite intersection of half-spaces, X = i∈I Hi+ . Such
decompositions
are usually not unique; we call a decomposition irreducible if
T
+
H
=
6
X
for
all i ∈ I. Given an irreducible decomposition, a facet of X is a
j6=i j
set Fi = X ∩Hi . The faces are the intersections of (potentially many) facets. The
normal cone to a point x0 ∈ X is N (X, x) = {y : y > (x−x0 ) ≤ 0 for all x ∈ X}.
One can show [31, Page 83] that for a given face F , all x0 ∈ F have the same
normal cone; hence we define the normal cone to F to be N (X, F ) = N (X, x0 )
for any x0 ∈ F .

Chételat, Lederer, Salmon/Two-Step Prediction

2542

Next, we make the following remarks. Recall that V is the collection of regions,
namely the closures of sets of points that have the same sign vector under the
lasso.
By Lemma A.1, then the regions must have disjoint interiors. Indeed, for
V 6= V 0 ∈ V we must have V = cl W 1 (η), V 0 = cl W 1 (η 0 ) for some η 6= η 0 ⊂
{−1, 0, 1}p . Since sgn β̂ is constant on int V , we conclude that int V ⊂ W 1 (η),
int V 0 ⊂ W 1 (η 0 ). But W 1 (η) ∩ W (η 0 ) = {Y : η = sgn β̂ 1 = η 0 } = ∅, so
int V ∩ V 0 = ∅.
Take V 6= V 0 ∈ V again: being polyhedra, their boundaries ∂V , ∂V 0 can be
partitioned by the relative interiors of their proper faces [31, Theorem 2.1.2].
Let F(V ) denote the set of proper faces of a polyhedron V and C = {relint F ∩
relint F 0 : F ∈ F(V ), F 0 ∈ F(V 0 ), V 6= V 0 ∈ V} be the collection of “boundary
pieces”. We enumerate, for reference, two properties of C:
i) For two distinct V, V 0 ∈ V and an x ∈ ∂V ∩ ∂V 0 , there is a unique B ∈ C
such that x ∈ B.
ii) Each B ∈ C has dimension at most n − 1.
Indeed, for the first statement we notice that by partitioning, there exists unique
faces F ∈ F(V ) and F 0 ∈ F(V 0 ) such that x ∈ relint F ∩ relint F 0 , hence a
unique B ∈ C such that x ∈ B. For the second, since the interiors of V, V 0 are
disjoint, for B = relint F ∩ relint F 0 6= ∅ we must have F ∈ ∂V ,F 0 ∈ ∂V 0 , hence
dim B ≤ dim relint F ∧ dim relint F 0 ≤ n − 1. In addition to these observations,
we will need the following two lemmas and one supporting proposition.
Lemma A.5. If B ∈ C is of dimension at most n−2, then R+ B is of dimension
at most n − 1.
Proof. Let S = aff B be the affine hull of B, of dimension at most n − 2. Being
an affine subspace, it can be written S = {x : Ax + b = 0} for some matrix
A ∈ Rn×n with rank A ≤ n − 2 and some vector b ∈ Rn . Now, B ⊂ S implies
R+ B ⊂ RS, and RS = {x : ∃t s.t. Atx + b = 0} = [In 0] {(x, t) : Ax + tb = 0} =
[In 0] Ker[A b]. Now, since rank A ≤ n − 2, [A b] has rank at most n − 1, and thus
Ker[A b] is a subspace of Rn+1 of rank at most n − 1. Hence, S = [In 0] Ker[A b]
has dimension at most n − 1. Since R+ B ⊂ RS, and RS is affine (actually, a
subspace), then aff R+ B ⊂ RS. Consequently, dim R+ B ≤ dim RS ≤ n − 1, as
desired.
Proposition A.1. Let K ⊂TRn be a polyhedron of full dimension n with irreducible decomposition K = j Hj+ . Denote the hyperplanes by Hi = ∂Hi+ as
usual. Then for any facet Fi = K ∩ Hi of K and any point x ∈ relint Fi , there
exists an ε > 0 such that int Hi+ ∩ B(x, ε) ⊂ int K, Hi ∩ B(x, ε) ⊂ Fi and
Hi+C ∩ B(x, ε) ⊂ K C .
Proof. With the irreducible decomposition, x ∈ K ∩ Hj+ for all j; say for j 6= i
we had x ∈ Fj = K ∩ Hj . As argued in [17, Page 27], x ∈ Fi ∩ Fj implies
that x is in a facet of Fi , hence in ∂Fi , a contradiction with x ∈ relint Fi . Thus
x ∈ K ∩ int Hj+ for all j 6= i.

Chételat, Lederer, Salmon/Two-Step Prediction

2543

T
T
Since x ∈ j6=i int Hj+ , we can find ε > 0 such that B(x, ε) ⊂ j6=i int Hj+ .
T
T
Then Hi+ ∩ B(x, ε) ⊂ j int Hj+ ⊂ int K, Hi ∩ B(x, ε) ⊂ Hi ∩ j6=i int Hj+ ⊂
S
Hi ∩ K = Fi , and finally Hi+C ∩ B(x, ε) ⊂ j Hj+C = K C , as desired.
Lemma A.6. Let B ∈ C be of dimension n − 1 and L be a ray centered at
the origin such that B ∩ L 6= ∅. Then either R+ B has dimension n − 1, or
L ∩ int V 6= ∅ and L ∩ int V 0 6= ∅.
Proof. Write B = relint F ∩ relint F 0 . Since relint F , relint F 0 both have dimension at most n − 1, and B has dimension n − 1, then they must have exactly
dimension n − 1. They are thus facets of their respective polyhedra V, V 0 ∈ V,
which must have dimension n, and have exactly one supporting hyperplane.
Thus aff B = aff F = aff F 0 , which we might denote S. Let b ∈ L ∩ S, write
S = ha1 , ..., an−1 i + b for some linearly independent vectors a1 , ..., an−1 , and let
an be such that L = R+ an . Consider the affine transformation
φ(x) = [a1 , ..., an−1 , an ] x + b = Ax + b,
which maps vectors (x1 , ..., xn−1 , 0) bijectively to S, and vectors
(0, ..., 0, xn ) with xn ≥ −kbk/kan k bijectively to L. (Recall that b ∈ S ∩ L.)
There are then two possibilities.
Case i) A has rank n − 1. Then an ⊂ ha1 , ...an−1 i and L ⊂ S. But this means
that 0 ∈ L ⊂ S, that is, that S is a subspace. Then R+ B ⊂ RS = S, that is,
R+ B has dimension n − 1.
Case ii) A has rank n. Then φ is bijective and L∩S is theT
singleton {b}.TRecall
that V and V 0 are polyhedra of dimension n, and let V = i Hi+ , V 0 = i Hi0+
be irreducible representations into half-spaces Hi+ , Hj0+ with boundary Hi , Hj0 .
By irreducibility, there are unique indices i, j such that aff F = Hi = aff F 0 =
Hj0 = S. Since b ∈ relint F ∩ relint F 0 , by Proposition A.1 there must be an
ε > 0 small enough so that int Hi+ ∩ B(b, ε) ⊂ int V , Hi+C ∩ B(b, ε) ⊂ V C ,
int Hj0 ∩ B(b, ε) ⊂ int V 0 and Hj0+C ∩ B(b, ε) ⊂ V 0C .
But clearly Hi+ 6= Hj0+ , as otherwise ∅ =
6 B(b, ε)∩int Hi+ = B(b, ε)∩int Hj0+ ⊂
0
int V ∩ int V = ∅, a contradiction since interiors of regions are disjoint. Thus it
must hold that int Hi+ = Hj0+C . In light of this, we may simplify the notation
by writing S + = int Hi+ and S − = Hj0+C . The state of affairs is then that
B(b, ε) ∩ S + ⊂ int V , B(b, ε) ∩ S ⊂ B and B(b, ε) ∩ S − ⊂ int V 0 .
Write R(n−1)+ = {x ∈ Rn : xn > 0} and R(n−1)− = {x ∈ Rn : xn < 0}.
Since φ is open, it must map connected components to connected components,
and being bijective it must hold that φ(R(n−1)+ ) = S + and φ(R(n−1)− ) = S − ,
or vice versa. Fix the former by considering x 7→ −Ax + b instead of x 7→ Ax + b
if necessary. Since φ−1 (B(b, ε)) is an open set, there must be an ε0 such that
B(0, ε0 ) ⊂ φ−1 (B(b, ε)).
Let ε00 = ε0 ∧ kbk2 /kan k2 , so that |t| < ε00 implies tan + b ∈ L. Then, we
have that ∅ 6= φ({(0, ..., 0, t), t ∈ (0, ε00 )}) ⊂ S + ∩ L ∩ B(b, ε) ⊂ L ∩ int V and
∅=
6 φ({(0, ..., 0, t), t ∈ (−ε00 , 0)}) ⊂ S − ∩ L ∩ B(b, ε) ⊂ L ∩ int V 0 . Thus, L ∩ int V
and L ∩ int V 0 are non-empty, as desired.

Chételat, Lederer, Salmon/Two-Step Prediction

2544

We may now turn to the proof of the theorem.
Proof of Theorem 2.2. Every region V ∈ V is a polyhedron, so has a decomposition V = P (V ) + C(V ) into a polytope P (V ) and a cone C(V ) by Minkowski’s
theorem
[42, Theorem 1.2]. Define V0 = {V ∈ V : int V ∩ R+ ξ 6= ∅}, T =
S
V
and
V ∈V0
R=

[

R+ B

B∈C
dim R+ B
≤n−1

[

∪

C(V )

V ∈V
dim C(V )
≤n−1

[

∪

∂C(V ).

V ∈V

The set R is a finite union of closed sets of dimension at most n − 1, so is closed
and has measure zero. We argue that if ξS∈ RC , then R+ ξ ⊂ int T . Indeed, say
that tξ ∈ ∂T for some t > 0. Since ∂T ⊂ V ∈V0 ∂V , there is a V0 ∈ V0 such that
S
S
tξ ∈ ∂V0 . Since Rn = V ∈V V and tξ ∈ ∂T , we have tξ ∈ cl (T C ) = V 6∈V0 V .
Thus there must also be a V1 6= V0 , V1 6∈ V0 such that tξ ∈ V1 . But since the
interiors are disjoint, if tξ ∈ int V1 there would be a contradiction with tξ ∈ ∂V0 ;
hence tξ ∈ ∂V1 . Thus tξ ∈ ∂V0 ∩ ∂V1 and there must be a unique B ∈ C such
that tξ ∈ B. That piece, like all elements of C, must be of dimension n − 1 or
lower.
We argue that B ⊂ R. If it has dimension n − 2, then dim R+ B ≤ n − 1 by
Lemma A.5, so B is indeed a subset of R. Now say instead it has dimension
n − 1 and recall that tξ ∈ R+ ξ ∩ B. Let F ∈ F(V0 ) and F 0 ∈ F(V1 ) be such that
B = relint F ∩ relint F 0 . By Lemma A.6, we must have either dim R+ B = n − 1,
or R+ ξ ∩ int V0 6= ∅ and R+ ξ ∩ int V1 6= ∅. But if the latter was the case, then
V1 ⊂ T by definition, which is impossible; thus we must have dim R+ B = n − 1,
hence B ⊂ R again.
Thus in all cases, B ⊂ R. Let d(x, A) := inf y∈A kx−yk2 denote the Euclidean
distance between a point x and a set A. Since tξ ∈ B, we have d(tξ, R) = 0. But
at the same time, since ξ ∈ RC and R is closed we must have d(ξ, R) > 0, and
since R is invariant under positive multiplication,
d(tξ, R) = inf ktξ −yk2 = t inf kξ −yk2 = t inf kξ −yk2 = td(ξ, R) > 0.
y∈R

y/t∈R

y∈R

This is a contradiction, and we conclude that R+ ξ ⊂ int T .
Now, V0 is finite, since C has only a finite number of faces. Let h be the
continuous map t 7→ tξ, and consider for each V ∈ V0 the closed set h−1 (V ). This
set must be convex, since for s, t ∈ h−1 (V ), h(γs + (1 − γ)t) = [γs + (1 − γ)t]ξ =
γ[sξ] + (1 − γ)[tξ] ∈ V by convexity of V for any γ ∈ [0, 1]. The only closed,
convex sets of R are the closed intervals: thus h−1 (V ) = [α, β] for some α ≤ β.
Enumerate arbitrarily the VS ∈ V0 as V1 , ..., Vm , and
for Vi ∈ V0 let [αi , βi ] =
Sm
m
h−1 (Vi ). Now, R+ ξ ⊂ T = i=1 Vi , so h−1 (T ) = i=1 [αi , βi ] = R+ . Then
some βi must equal ∞, otherwise the union would be bounded. Moreover, since
the interiors of the V ’s are disjoint, (αi , βi ) ∩ (αj , βj ) = ∅ for i 6= j and the
βi = ∞ must be unique, all the others finite. By reordering if necessary, take
0 = α1 < β1 ≤ α2 < β2 ≤ ... ≤ αm < βm = ∞.

Chételat, Lederer, Salmon/Two-Step Prediction

2545

The region Vm is a polyhedron, so has a decomposition as Vm = P (Vm ) +
C(Vm ) for some polytope P (Vm ) and cone C(Vm ) by Minkowski’s theorem [42,
Theorem 1.2]. Fix a point t0 ∈ (αm , ∞); then since t0 ξ + R+ ξ = (t0 , ∞)ξ ⊂
int Vm , by [17, 2.5.1], we conclude that R+ ξ ⊂ C(Vm ), so t0 ξ, ξ ∈ C(Vm ). If
C(Vm ) has dimension at most n − 1, or t0 ξ ∈ ∂C(Vm ) ⇔ ξ ∈ ∂C(Vm ), then
ξ ∈ R, which contradicts ξ ∈ RC – thus C(Vm ) must have dimension n and
t0 ξ, ξ ∈ int C(Vm ). Let ε1 be small enough so that B(ξ, ε1 ) ⊂ int C(Vm ). Then
for any s > 0, B([t0 + s]ξ, sε1 ) = t0 ξ + sB(ξ, ε1 ) ⊂ int Cm ⊂ int Vm . Thus for
all s > 2t0 kξk2 /ε1 , B(sξ, sε1 /2) ⊂ B([t0 + s]ξ, sε1 ) ⊂ int Vm ⊂ int T .
Next, notice that the segment [0, 2t0 kξk2 /ε1 ]ξ is compact and in int T . This
implies that d([0, 2t0 kξk2 /ε1 ]ξ, int T C ) is strictly positive and also that there
must be an ε2 -neighborhood such that B([0, 2t0 kξk2 /ε1 ]ξ, ε2 ) ⊂ int T . Hence,
for all s ∈ [0, 2t0 kξk2 /ε1 ], it holds B(sξ, sε2 ε1 /2t0 kξk2 ) ⊂ B(sξ, ε2 ) ⊂ int T .
Finally, let ε = min(ε1 /2, ε1 ε2 /2t0 kξk2 ). Then for all s ≥ 0, we have B(sξ, sε) ⊂
int T . Let |y − ξ| < ε and define η = sgn β̂(sy).
S Then cl W (η) = V0 for some
V0 ∈ V0 , since otherwise cl W (η) ⊂ cl T C = V 6∈V0 V , which would contradict sy ∈ int T . But there is a t > 0 such that tξ ∈ int V0 , since V0 ∈ V0 ,
and since sgn β̂ is constant over int V0 , sgn β̂(tξ) = sgn β̂(sy). Thus in particular
Ŝ[sy] = Ŝ[tξ]. Since this is true for all |y−ξ| < ε, we conclude that D(ξ) ≥ ε > 0,
as desired.
Appendix B: Description of the lslassoBIC
In this section, we provide details for the lslassoBIC implementation that we
have used. This method is similar to applying a BIC procedure over the refitted
models obtained by a lasso path, and was recently analyzed in [2].
We consider the same collection of tuning parameters Λ = {λ1 , . . . , λr } as
before, and we denote the associated supports by (Ŝ 1 , . . . , Ŝ r ), where Ŝ i :=
supp[β̂ λi ]. Following Equation (1.3), we write (β̄ λ1 , . . . , β̄ λr ) for the estimated
least-squares over these supports. Let us introduce for each j ∈ [r] a prior πj
via
 

−1
p
j
πj = Hp
exp(|
Ŝ
|)
|Ŝ j |
with Hp = (e − e−p )/(e − 1). Then, the lslassoBIC is defined by


1
β̄ lslassoBIC = β̄ λj? with j ? ∈ arg min kY − X β̄ λj k22 + 14σ̂ 2 log
,
πj
j∈[r]
where σ̂ is a standard deviation estimate of the (Gaussian) noise. As the practical
estimation of σ is not discussed further in [2], we have used the same estimator
as for AVPr , namely the one defined in Algorithm 2.
Note that we have adapted the method proposed by [2] to the case of a
predetermined number of lasso parameters. This is because the number of kinks
over the lasso path can be as large as (3p + 1)/2 [26], making an estimator based
on the entire collection of kinks intractable.

