Extended Comparisons of Best Subset Selection, Forward Stepwise Selection, and the Lasso
Trevor Hastie Robert Tibshirani Ryan J. Tibshirani
stat.ME stat.CO

A Jointly Learned Deep Architecture for Facial Attribute Analysis and Face Detection in the Wild
    Keke He, Yanwei Fu, Xiangyang Xue

Fudan University 

{kkhe15, yanweifu, xyxue}@fudan.edu.cn

    December 30, 2023
=================================================================================================




In exciting new work, <cit.> showed that the
classical best subset selection problem in regression modeling can be
formulated as a mixed integer    
optimization (MIO) problem. Using recent
advances in MIO algorithms, they demonstrated that  best subset selection can
now be solved at much larger problem sizes that what was thought possible 
in the statistics community.  They presented
empirical comparisons of best subset selection with other popular  
variable selection procedures, in particular, the lasso and forward
stepwise selection. Surprisingly (to us), their simulations suggested
that best subset selection consistently outperformed both methods in
terms of prediction accuracy.  Here we present an expanded set of
simulations to shed more light on these comparisons.
The summary is roughly as follows: 


  * neither best subset selection
nor the lasso uniformly dominate the other, with best subset
selection generally performing better in high signal-to-noise (SNR)
ratio regimes, and the lasso better in low SNR regimes; 

  * best subset selection and forward stepwise perform quite
  similarly throughout;

  * the relaxed lasso (actually, a simplified
version of the original relaxed estimator defined in
) 
 is the overall winner, performing just about as well
as the lasso in low SNR scenarios, and as well as best subset selection
in high SNR scenarios. 














.tocmain


§ INTRODUCTION


Best subset selection, forward stepwise selection, and the lasso are
popular methods for selection and estimation of the parameters in a
linear model.  The first two are classical methods in statistics,
dating back to at least
<cit.> for  
best subset selection and <cit.>
for forward selection; the lasso is (relatively speaking) more recent,
due to <cit.>.

Given a response vector Y ∈^n, predictor matrix 
X ∈^n× p, and a subset size k between 0 and min{n,p}, best     
subset selection finds the subset of k predictors that produces the best fit 
in terms of squared error, solving the nonconvex problem

    _β∈^p Y-Xβ_2^2     β_0 ≤ k,

where β_0=∑_i=1^p 1{β_i ≠ 0} is
the ℓ_0 norm of β. (Here and throughout, for notational simplicity, 
we omit the intercept term from the regression model.) 

Forward stepwise selection is less ambitious: starting with the empty model,  
it iteratively adds the variable that best improves the fit.[Other ways
  of defining the variable j_k that “best improves the fit” are possible,
  but the entry criterion is (<ref>) is the standard one in
  statistics.]  It hence yields a subset of each size k=0,…,min{n,p}, 
but none of these are generally globally optimal in the sense of
(<ref>). Formally, the procedure starts with an empty active 
set A_0={0}, and for k=1,…,min{n,p}, selects the variable
indexed by  

    j_k = _j ∉ A_k-1 Y - P_A_k-1∪{j_k} Y
    _2^2 = _j ∉ A_k-1 X_j^T P_A_k-1^⊥
      Y/P_A_k-1^⊥ X_j_2

that leads to the lowest squared error when added to A_k-1, or equivalently,
such that X_j_k, achieves the maximum absolute correlation with Y, 
after we project out the contributions from X_A_k-1. A note on
notation: here we write X_S ∈^n× |S| for the submatrix of X
whose columns are indexed by a set S (and when S={j}, we 
simply use X_j). We also write P_S for the projection matrix onto the column   
span of X_S, and P_S^⊥=I-P_S for the projection matrix onto the 
orthocomplement.  At the end of step 
k of the procedure, the active set is updated, A_k=A_k-1∪{j_k}, and the forward stepwise estimator of the 
regression coefficients is defined by the least squares fit onto
X_A_k.      














The lasso solves a convex relaxation of (<ref>) where we
replace the ℓ_0 norm by the ℓ_1 norm, namely 

    _β∈^p Y-Xβ_2^2      β_1 ≤ t,

where β_1=∑_i=1^p |β_i|, and t ≥ 0 
is a tuning parameter. By convex duality, the above problem
is equivalent to the more common (and more easily solveable) penalized
form   

    _β∈^p Y-Xβ_2^2  + λβ_1

where now λ≥ 0 is a tuning parameter. This is the form 
that we focus on in this paper.

The lasso problem (<ref>) is convex (and highly
structured) and there is by now a sizeable literature in statistics,
machine learning, and optimization dedicated to efficient algorithms
for this problem.  On the other hand, the best subset selection
problem (<ref>) is nonconvex and is known to be NP-hard
<cit.>.  The accepted view in statistics for many
years has been that this problem is not solveable
beyond (say) p in the 30s, this view being shaped by the
available software for best subset selection (e.g., in the R language,
the leaps package implements a branch-and-bound algorithm for
best subset selection of ).  

For a much more detailed introduction to best subset selection, forward stepwise 
selection, and the lasso, see, e.g., Chapter 3 of <cit.>. 



 §.§ An exciting new development


Recently, <cit.> presented a mixed integer
optimization (MIO) formulation for the best subset selection problem
(<ref>).  This allows one to use highly optimized MIO  
solvers, like Gurobi (based on branch-and-cut methods, 
hybrids of branch-and-bound and cutting plane algorithms), to solve 
(<ref>). Using these MIO solvers, problems with p in the
hundreds and even thousands are not out of reach, and this presents us
with exciting new ground on which to perform empirical 
comparisons. Simulation studies in 
<cit.> demonstrated that best subset selection
generally gives superior prediction accuracy compared to forward
stepwise selection and the lasso, over a variety of problem setups.  

In what follows, we replicate and expand these simulations to shed 
more light on such comparisons.
For convenience, we made an R package bestsubset for   
optimizing the best subset selection problem using the Gurobi MIO
solver (after this problem has been translated into a mixed 
integer quadratic program as in ). This
package, as well as R code for reproducing all of the results in this 
paper, are available at  
<https://github.com/ryantibs/best-subset/>. 



§ PRELIMINARY DISCUSSION




 §.§ Is best subset selection the holy grail?


Various researchers throughout the years have viewed best subset 
selection as the “holy grail” of estimators 
for sparse modeling in regression, suggesting (perhaps implicitly)
that it should be used whenever possible, and that other methods for
sparse regression—such as forward stepwise selection and the lasso—should  
be seen as approximations or heuristics, used only out of necessity when best  
subset selection is not computable.       
However, as we will demonstrate in the simulations that
follow, this is not the case.  Different procedures have different 
operating characteristics, i.e., give rise to different bias-variance 
tradeoffs as we vary their respective tuning parameters.  In fact, depending on 
the problem setting, the bias-variance tradeoff provided by best
subset selection may be more or less useful than the tradeoff
provided by the lasso. 

As a brief interlude, let us inspect the “noiseless” versions of the
best subset and lasso optimization problems, namely  

    _β∈^p β_0      Xβ=Y, 
    _β∈^p β_1      Xβ=Y,

respectively.  Suppose that our goal is to find the sparsest
solution to the linear system Xβ=Y. Problem (<ref>), 
by definition of the ℓ_0 norm, produces
it. Problem (<ref>), in which the criterion has been  
convexified, does not generally give the sparsest solution and so 
in this sense we may rightly view it as a heuristic for the 
nonconvex problem (<ref>).  Indeed, much of the 
literature on compressed sensing (in which
(<ref>), (<ref>) have
been intensely studied) uses this language.
However, in the noiseless setting, there is no bias-variance tradeoff,
because (trivially) there is no bias and no variance; both of
the estimators defined by (<ref>),
(<ref>) have zero mean squared error owing 
to the linear constraint Xβ=Y (and the fact that Y=(Y|X),
as there is no noise).  

The noisy setting—which is the traditional and
most practical setting for statistical estimation, and that studied in 
this paper—is truly different. Here, it is no longer appropriate
to view the estimator defined by the ℓ_1-regularized problem   
(<ref>) 
as a heuristic for that defined by the ℓ_0-regularized problem 
(<ref>) (or (<ref>) as a heuristic for
(<ref>)). 
Generally speaking, the lasso and best subset selection 
differ in terms of their “aggressiveness” in selecting and estimating the 
coefficients in a linear model, with the lasso being less aggressive
than best subset selection; meanwhile, forward stepwise lands somewhere in the
middle, in terms of its aggressiveness. There are various ways to make this
vague but intuitive comparison more explicit.  For example:


  * forward stepwise can be seen as a “locally optimal”
  version of best subset selection, updating the active
  set by one variable at each step, instead of re-optimizing
  over all possible subsets of a given size; in turn, the lasso can
  be seen as a more “democratic” version of forward stepwise,
  updating the coefficients so as 
  maintain equal absolute correlation of all active variables with
  the residual <cit.>;

  * the lasso applies shrinkage to its nonzero estimated
  coefficients (e.g., see (<ref>) with
  γ=1) but forward stepwise and best subset selection do
  not, and simply perform least squares on their respective active
  sets;  

  * thanks to such shrinkage, the fitted values from the lasso (for
  any fixed λ≥ 0) are continuous functions of y
  <cit.>, whereas the fitted
  values from forward stepwise and best subset selection
  (for fixed k ≥ 1) jump discontinuously as y moves across
  a decision boundary for the active set; 

  * again thanks to shrinkage, the effective degrees of freedom of
  the lasso (at any fixed λ≥ 0) is equal to the expected
  number of selected variables
  <cit.>,  
  whereas the degrees of freedom of both forward stepwise and
  best subset selection can greatly exceed k at any given step k
  ≥ 1 <cit.>.  

Figure <ref> uses the latter perspective of effective degrees 
of freedom to contrast the aggressiveness of the three methods.  



When the signal-to-noise ratio (SNR) is low, and also depending  
on other factors like the correlations between predictor variables,  
the more aggressive best subset and  
forward stepwise methods can already have quite high
variance at the start of their model paths (i.e., for small
step numbers k). Even after optimizing over the tuning 
parameter k (using say, an external validation set or an
oracle which reveals the true risk), we can arrive at an
estimator with unwanted variance and worse accuracy than a 
properly tuned lasso estimator. On the other hand, for high SNR
values, and other configurations for the correlations between
predictors, etc., the story can be completely flipped and the
shrinkage applied by the lasso estimator can result in unwanted bias  
and worse accuracy than best subset selection and forward stepwise
selection. See Figure <ref> for empirical evidence.   



This is a simple point, but is worth emphasizing.  To convey
the idea once more:
Different procedures bring us from the high bias to the high variance
ends of the tradeoff along different model paths; and these paths are 
affected by aspects of the problem setting, like the SNR and predictor
correlations, in different ways. For some classes of problems, some 
procedures admit more fruitful paths, and for other classes, other
procedures admit more fruitful paths.  For example, neither best
subset selection nor the lasso dominates the other, across
all problem settings.   




 §.§ What is a realistic signal-to-noise ratio?



In their simulation studies, <cit.> considered 
SNRs in the range of about 2 to 8 in their low-dimensional cases, and
about 3 to 10 in their high-dimensional cases. Is this a realistic
range that one encounters in practice? 


In our view, inspecting the proportion of variance explained (PVE) 
can help to answer this question.

Let (x_0,y_0) ∈^p × be a pair of predictor and response    
variables, and define f(x_0)=(y_0|x_0) and ϵ_0=y_0-f(x_0), so that  
we may express the relationship between x_0,y_0 as:

    y_0 = f(x_0) + ϵ_0.

The signal-to-noise ratio (SNR) in this model is defined as 

    = (f(x_0))/(ϵ_0).

For a given prediction fuction g—e.g., one trained on n samples
(x_i,y_i), i=1,…,n that are i.i.d. to (x_0,y_0)—its associated
proportion of variance explained (PVE) is defined as 

    (g) = 1 - (y_0-g(x_0))^2/(y_0).

Of course, this is maximized when we take g to be the mean function
f itself, in which case

    (f) = 1 - (ϵ_0)/(y_0) = /1+.

In the second equality we have assumed independence of x_0,
ϵ_0, so (y_0)=(f(x_0))+(ϵ_0).
As the optimal prediction function is f, it sets the gold-standard of
/(1+) for the PVE, so we should always expect to see the attained PVE
be less than /(1+) and greater than 0 (otherwise we could simply
replace our prediction function by g=0.) 

We illustrate using a simulation with n=200 and p=100. The 
predictor autocorrelation was zero and the coefficients followed the beta-type 2 
pattern with s=5; see Section <ref> for details.  We varied the SNR 
in the simulation from 0.05 to 6 in 20 equally spaced values. We 
computed the lasso over 50 values of the tuning parameter λ, and
selected the tuning 
parameter by optimizing prediction error on a separate validation set of size 
n. Figure <ref> shows the PVE of the tuned lasso estimator, 
averaged over 20 repetitions from this simulation setup.  Also shown is the
population PVE, i.e., the maximum possible PVE at any given SNR level, of 
/(1+).   
We see that an SNR of 1.0 corresponds to a PVE of about 0.45 (with a maximum of 
0.5), while an SNR as low as 0.25 yields a PVE of 0.1 (with a maximum of 0.2).   
In our experience, a
PVE of 0.5 is rare for noisy observational data, and 0.2 may be more
typical. A PVE of 0.86, corresponding to an SNR of 6, is unheard of!  With
financial returns data, explaining even 2% of the variance (PVE 
of 0.02) would be considered huge, and the corresponding prediction function
could lead to considerable profits if used in a trading scheme.  Therefore,
based on these observations, we examine a wider lower range of SNRs in our 
simulations, compared to the SNRs studied in <cit.>.   





 §.§ A (simplified) relaxed lasso


In addition to the lasso estimator, we consider a simplified version
of the relaxed lasso estimator as originally defined by
<cit.>.
Let ^lasso(λ) denote the solution in
problem (<ref>), i.e., the lasso estimator at the tuning 
parameter value λ≥ 0. Let A_λ denote its active set, and let 
^LS_A_λ denote
the least squares coefficients from regressing of Y on X_A_λ, the 
submatrix of active predictors. Finally, let
^LS(λ) be the full-sized (p-dimensional)
version of the least squares coefficients, padded with zeros in the
appropriately.  We consider the estimator
^relax(λ,γ) defined by 

    ^relax(λ,γ) =
    γ^lasso(λ) + (1-γ)
    ^LS(λ)

with respect to the pair of tuning parameter values λ≥ 0 
and γ∈ [0,1]. Recall 
<cit.> that when the columns of X are in
general position (a weak condition occurring almost surely for
continuously distributed pedictors, regardless of n,p), it holds
that:  


  * the lasso solution is unique;

  * the submatrix X_A_λ of active predictors has full column rank,
thus _A_λ^LS=(X_A_λ^T
  X_A_λ)^-1 X_A_λ^T Y
is well-defined; 

  * the lasso solution can be written (over its active set) as
^lasso_A_λ(λ)=(X_A_λ^T
  X_A_λ)^-1(X_A_λ^T Y -  λ s), where 
s ∈{-1,1}^|A_λ| contains the signs of the active lasso
coefficients.  

Thus, under the general position assumption on X, the simplified
relaxed lasso can be rewritten as 

    ^relax_A_λ(λ,γ)    = (X_A_λ^T X_A_λ)^-1 X_A_λ^T
    Y - γλ (X_A_λ^T X_A_λ)^-1 s 
    ^relax_-A_λ (λ,γ)   = 0,

so we see that γ∈ [0,1] acts as a multiplicative factor
applied directly to the “extra” shrinkage term apparent in the lasso
coefficients.  Henceforth, we will drop the word “simplified” and 
will just refer to this estimator as the relaxed lasso.

The relaxed lasso tries to undo the shrinkage inherent in the lasso estimator,
to a varying degree, depending on γ.  In this sense, we would expect it
to be more aggressive than the lasso, and have a larger effective degrees of
freedom.  However, even in its most aggressive mode, γ=0, it 
is typically less aggressive than both forward stepwise selection
and best subset selection, in that it often has a smaller degrees of
freedom than these two.  See Figure <ref> for an example.  





 §.§ Other estimators


Many other sparse estimators for regression could be
considered, for example, ℓ_1-penalized alternatives to 
the lasso, like the Dantzig selector
<cit.> and square-root lasso
<cit.>; greedy alternatives to forward stepwise
algorithm, like matching pursuit <cit.> and
orthogonal matching pursuit <cit.>;
nonconvex-penalized methods, such as SCAD <cit.>,
MC+ <cit.>, and SparseNet <cit.>;  
hybrid lasso/stepwise approaches like FLASH <cit.>; 
and many others.   

It would be interesting to include all of these estimators in our
comparisons, though that would make for a huge simulation suite and
would dilute 
the comparisons between best subset selection, forward stepwise, and  
the lasso that we would like to highlight. Roughly speaking, we would
expect the Dantzig selector and square-root lasso to perform similarly 
to the lasso; the matching pursuit variants to perform similarly
to forward stepwise; and the nonconvex-penalized methods to perform
somewhere in between the lasso and best subset selection.  (It is
worth noting that our R package is structured in such a way to
make further simulations and comparisons straightforward. We
invite interested readers to use it to perform comparisons to other
methods.) 



 §.§ Brief discussion of computational costs



Computation of the lasso solution in 
(<ref>) has been a
popular topic of research, and there are by now many efficient lasso  
algorithms.  In our simulations, we use coordinate descent with warm starts
over a sequence of tuning parameter values λ_1 > ⋯  
> λ_m > 0, as implemented in the glmnet R package 
<cit.>.  The base
code for this is 
written in Fortran, and warm starts—plus additional tricks like 
active set optimization and screening rules
<cit.>—make this implementation highly
efficient.  For example, for a problem with n=500 observations and
p=100 variables, glmnet delivers the lasso solutions across 
100 values of λ in less than 0.01 seconds, on a standard
laptop computer.  The relaxed lasso in (<ref>) comes
at only a slight increase in computational cost, seeing as we must
only additionally compute the least squares coefficients on each
active set.  We provide an implementation in the bestsubset R
package accompanying this paper, which
just uses an R wrapper around glmnet. 
For the same example with n=500 and p=100, computing the relaxed
lasso path over 100 values of λ and 10 values of γ 
again took less than 0.01 seconds.

For forward stepwise selection, we implemented our own version 
in the bestsubset R package.  The core matrix
manipulations for this method are written in C, and the rest 
is written in R.  The forward stepwise path is highly
structured and this greatly aids its computation: at step k, we have 
k-1 active variables included in the model, and we seek the variable
among the remaining p-k+1 that—once orthogonalized with respect
to the current active set of variables—achieves the greatest
absolute correlation with Y, as in (<ref>). Suppose that we 
have maintained a QR decomposition of the active submatrix
X_A_k-1 of predictors, as well as the orthogonalization of the
remaining p-k-1 predictors with respect to X_A_k-1.
We can compute the necessary
correlations in O(n(p-k+1)) operations, update the  QR  
factorization of X_A_k in constant time, and orthogonalize
the remaining predictors with respect to the one just
included in O(n(p-k))  
operations (refer to the modified Gram-Schmidt algorithm in
).  Hence, the forward stepwise 
selection path can be seen as a certain guided QR decomposition for 
computing the least squares coefficients on all p variables (or, on some 
subset of n variables when p>n).  For the same
example with n=500 and p=100, our implementation computes the
forward stepwise path in less than 0.5 seconds.

Best subset selection (<ref>) is the most
computationally challenging, by a large margin.
<cit.> 
describe two reformulations of (<ref>) as a mixed integer
quadratic program, one that is preferred when n ≥ p, and the
other when p>n, and recommend using the Gurobi commercial MIO
solver (which is free for academic use).  They also describe a
proximal gradient descent method for computing approximate solutions
in (<ref>), and recommend using the best output from this
algorithm over many randomly-initialized runs to warm start the
Gurobi solver.  See <cit.> for details.
We have implemented the method of these authors[We thank  
  the third author Rahul Mazumder for his help and guidance.]—which 
transforms the best subset selection problem into one of two MIO formulations
depending on the relative sizes of n and p, uses proximal gradient
to compute a warm start, and then calls Gurobi through its R
interface—in our accompanying R package bestsubset.  

Gurobi uses branch-and-cut techniques 
(a combination of branch-and-bound and cutting plane methods),
along with many other sophisticated optimization tools, for MIO
problems.  Compared to the pure branch-and-bound method from 
the leaps R package, 
its speed can be impressive: for example, in one run with n=500 and 
p=100, it returned the best subset selection solution of size k=8
in about 3 minutes (brute-force search for this problem would need to
have looked at about 186 billion candidates!).  But for most problems
of this size (n=500 and p=100) it has been our experience that
Gurobi typically requires 1 hour or longer to complete its optimization.    
The third author Rahul Mazumder of
<cit.> suggested to us that for these problem
sizes, it is often the case that Gurobi has found the solution in less
than 3 minutes, though it takes much longer to certify its optimality. For
our simulations in the next section, we used a time limit of 3 minutes for
Gurobi to optimize the best subset selection problem (<ref>) at any
particular value of the subset size k (once the time limit has been reached,
the solver returns its best iterate).  For more discussion on this
choice and its implications, see Section <ref>.  We note that this is
already quite a steep computational cost for “regular” practical usage: at 3
minutes per value of k, if we wanted to use 10-fold cross-validation to choose
between the subset sizes k=0,…,50, then we are already facing 25 hours of   
computation time.  



§ SIMULATIONS




 §.§ Setup



We present simulations, basically following the simulation setup of  
<cit.>, except that we consider a wider range of SNR 
values. Given n,p (problem dimensions), s
(sparsity level), beta-type (pattern of sparsity), ρ 
(predictor autocorrelation level), and ν (SNR level), our process can be
described as follows:    


  i. we define coefficients β_0 ∈^p according to s
  and the beta-type, as described below;

  ii. we draw the rows of the predictor matrix X ∈^n × p i.i.d. from N_p(0,Σ), where Σ∈^p× p has entry (i,j) equal to ρ^|i-j|; 

  iii. we draw the response vector Y ∈^n from 
N_n(Xβ_0, σ^2 I), with σ^2 defined to meet the   
  desired SNR level, i.e., σ^2=β_0^T Σβ_0 /ν;  

  iv. we run the lasso, relaxed lasso, forward stepwise
  selection, and best subset selection on the data X,Y, each over a 
  wide range of tuning parameter values; for each method we choose
  the tuning parameter by minimizing prediction error on  
  validation set X∈^n× p, 
    Y∈^n that is generated independently of and  
  identically to X,Y, as in steps ii–iii above;

  v. we record several metrics of interest, as specified below; 

  vi. we repeat steps ii-v a total of 10 times, and average the
  results. 

Below we describe some aspects of the simulation process in more
detail. 



  
Coefficients.  We considered four settings for the
coefficients β_0 ∈^p: 


  * beta-type 1: β_0 has s components equal to 1, 
  occurring at (roughly) equally-spaced indices between 1 and p,
  and the rest equal to 0; 

  * beta-type 2: β_0 has its first s components equal
  to 1, and the rest equal to 0;

  * beta-type 3: β_0 has its first s components
  taking nonzero values equally-spaced between 10 and 0.5, and the  
  rest equal to 0; 

  * beta-type 5: β_0 has its first s components equal
  to 1, and the rest decaying exponentially to 0, 
  specifically, β_0i=0.5^i-s, for i=s+1,…,p. 

The first three types were studied in
<cit.>.  They also defined a fourth type 
that we did not include here, as we found it yielded basically the 
same results as beta-type 3.  The last type above is new: we
included it to investigate the effects of weak sparsity and call it 
beta-type 5, to avoid confusion.



  
Evaluation metrics. Let x_0 ∈^p denote test  
predictor values drawn from N_p(0,Σ) (as in the rows of the 
training predictor matrix X) and let y_0 ∈ denote its
associated response value drawn from N(x_0^T β_0,
σ^2). Also let  denote estimated coefficients
from one of the regression procedures.  We considered the following
evaluation metrics:  


  * Relative risk: this is the accuracy metric studied in
  <cit.>[Actually, these authors used
    an “in-sample” version of this metric defined as
    X-Xβ_0_2^2/X_2^2, whereas our 
    definition is “out-of-sample”, with an expectation over the new
    test predictor value x_0 taking the place of the sample average 
    over the training values x_i, i=1,…,n.], defined as

    RR() = 
    (x_0^T-x_0^Tβ_0)^2/(x_0^Tβ_0)^2 = 
    (-β_0)^TΣ(-β_0)/β_0^TΣβ_0.

The expectations here and below are taken over the test point (x_0,y_0), 
with all training data and validation data (and thus ) held
fixed. A perfect score is 0 (if =β_0)
and the null score is 1 (if =0).






























  * Relative test error: this measures the expected test error  
  relative to the Bayes error rate, 

    RTE() =
    (y_0-x_0^T)^2/σ^2 =
    (-β_0)^TΣ(-β_0) 
      +σ^2/σ^2.

A perfect score is 1 and the null score is
(β_0^TΣβ_0+σ^2)/σ^2=+1.  


  * Proportion of variance explained: 
as defined in Section <ref>, this is

    PVE(β̂)
    =1-(y_0-x_0^T )^2/(y_0) 
    =1-(-β_0)^TΣ(-β_0)+σ^2/β_0^TΣβ_0+σ^2.

A perfect score is /(1+) and the null score is 0. 


  * Number of nonzeros: unlike the last three metrics
  which measure predictive accuracy, this metric simply records the
  number of nonzero estimated coefficients, 
  _0 = ∑_i=1^p 1{_i≠0}.


It is worth noting that, in addition to metrics based on
predictive accuracy, it would be useful to consider a metric
that measures proper variable recovery, i.e., the extent to which the 
sparsity pattern in the estimated  matches that in
β_0. We briefly touch on this in the discussion. Here we mention
one advantage to studying predictive accuracy: any of the metrics
defined above are still relevant when (y|x) is
no longer assumed to be linear, making the predictive angle more
broadly practically relevant than a study of proper variable recovery 
(which necessarily requires linearity of the true mean). 



  
Configurations.  We considered the following
four problem settings:
 

  * low: n=100, p=10, s=5;

  * medium: n=500, p=100, s=5; 

  * high-5: n=50, p=1000, s=5;

  * high-10: n=100, p=1000, s=10.

In each setting, we considered ten values for the SNR ranging
from 0.05 to 6 on a log scale, namely

    [      0.05 0.09 0.14 0.25 0.42 0.71 1.22 2.07 3.52 6.00;     ;      0.05 0.08 0.12 0.20 0.30 0.42 0.55 0.67 0.78 0.86 ]

(For convenience we provide the corresponding population PVE as well.)
In each setting, we also considered three values for the
predictor autocorrelation ρ, namely 0, 0.35, and 0.7.



  
Tuning of procedures.
In the low setting, the lasso was tuned over 50       
values of λ ranging from λ_max=X^T Y_∞
to a small fraction of λ_max on a log scale, as per the 
default in glmnet, and the relaxed lasso was tuned over the same 50 
values of λ, and 10 values of γ equally spaced from 1 to
0 (hence a total of 500 tuning parameter values).  Also in the low setting,  
forward stepwise and best subset selection were tuned over steps
k=0,…,10.  In all other problem settings (medium, high-5, and high-10),
the lasso was tuned over 100 values of λ, the relaxed lasso was tuned
over the same 100 values of λ and 10 values of γ (hence 1000
tuning parameter values total), and forward stepwise and best subset selection
were tuned over steps k=0,…,50. In all cases, tuning was performed by  
by minimizing prediction error on an external validation set of size n, which
we note approximately matches the precision of leave-one-out cross-validation.  



 §.§ Time budget for Gurobi



As mentioned in Section <ref>, for each problem instance
and subset size k, we used a time limit of 3
minutes for Gurobi to optimize the best subset selection problem.  In
comparison, <cit.> 
used much larger time budgets: 15 minutes (per problem per k) for
problems with p=100 as in our medium setup, and 66 minutes (per
problem per k) for problems with p ≥ 2000 as in our high-5 and high-10
setups.  Their simulations however were not as extensive, as they
looked at fewer combinations of beta-types, SNR levels, and correlation levels.  
Another important difference worth mentioning: the Gurobi optimizer,
when run through its Python/Matlab interface,
 as in <cit.>, 
automatically takes advantage of multithreading capabilities; this does not
appear to be the case when run through its R interface, as in our simulations. 

The third author of <cit.>, Rahul Mazumder, suggested 
in personal communication that the MIO solver in the medium setting will often
arrive at the best subset selection solution in less than 3 minutes, but it can
take much longer to certify its optimality[Gurobi
  constructs a sequence of lower and upper bounds on the criterion in  
  (<ref>); typically the lower bounds come from convex 
  relaxations and the upper bounds from the current iterates, and it
  is the lower bounds that take so long to converge.] (usually over 1 hour, in
absence of extra speedup tricks as described in ).  
Meanwhile, in the high-5 and high-10 settings, this author also pointed out that
3 minutes may no longer be enough.  For practical reasons, we have kept 
the 3 minute budget per problem instance per subset size.  
Note that this amounts to 150 minutes per path of 50 solutions, 1500 minutes or
25 hours per set of 10 repetitions, and in total 750 hours or 31.25 days for
any given setting, once we go through the 10 SNR levels and 3 correlation levels.    



 §.§ Results: computation time


In Table <ref>, we report the time in seconds taken by each 
method to compute one path of solutions, averaged over 10 repetitions and all
SNR and predictor correlation levels in the given setting.  All
timings were recorded on a Linux cluster.
As explained above, the lasso path consisted of 50 tuning parameter values in
the low setting and 100 in all other settings, the relaxed lasso path consisted
of 500 tuning parameter values in the low setting and 1000 in all other
settings, and the forward stepwise and best subset selection paths each
consisted of min{p,50} tuning parameter values.  



We can see that the lasso and relaxed lasso are very fast,
requiring less than 25 milliseconds in every case.
Forward stepwise selection is also fast, though not quite as fast as the lasso
(some of the differences here might be due to the fact that 
our forward stepwise  
algorithm is implemented partly in R).  Moreover, it should be noted that 
when n and p is large, and one wants to explore models with a sizeable
number of variables (we limited our search to models of size 50), forward
stepwise has to plod through its path one variable at a time, but the lasso 
can make jumps over subset sizes bigger than one by varying λ and
leveraging warm starts.

Recall, the MIO solver for best subset selection was allowed 3 minutes per
subset size k, or 150 minutes for a path of 50 subset sizes.
As the times in Table <ref> suggest, the maximum allotted time was
not reached in all instances, and the MIO solver managed to verify optimality
of some solutions along the path. In the medium setting, on average 
17.55 of the 50 solutions were verified as being optimal. In the
high-5 and high-10 settings, only 1.61 of the 50 were verified on average (note
this count includes the subset of size 1, which is trivial).  These measures may 
be pessismistic, as Gurobi may have found high-quality approximate solutions or 
even exact solutions but was just not able to verify them in time, see the
discussion in the above subsection.



 §.§ Results: accuracy metrics


Here we display a slice of the accuracy results, focusing for concreteness on
the case in which the predictor autocorrelation is ρ=0.35, and the 
population coefficients follow the beta-type 2 pattern.  In a supplementary
document, we display the full set of results, over the whole simulation
design.  

Figure <ref> plots the relative risk, relative test error, PVE, and
number of nonzero coefficients as functions of 
the SNR level, for the low setting.  Figures <ref>, <ref>,
and <ref> show the same for the medium, high-5, and high-10 settings,  
respectively.  Each panel in the figures displays the average of
a given metric over 10 repetitions, for the four methods in question, and
vertical bars denote one standard 
error.  In the relative test error plots, the dotted curve denotes
the performance of the null model (null score); in the PVE plots, it denotes 
the performance of the true model (perfect score); in the number of 
nonzero plots, it marks the true support size s. 

The low and medium settings, Figures <ref> and <ref>, yield
somewhat similar results.  In the relative 
risk and PVE plots (top left and bottom left panels), we see that best subset
and forward stepwise selection lag behind the lasso and relaxed lasso in terms
of accuracy for low SNR levels, and as the SNR increases, we see that all
four methods converge to nearly perfect accuracy.  The relative test error plot
(top right panel) magnifies the differences between the methods.
For low SNR levels, 
we see that the lasso outperforms the more aggressive best subset and forward
stepwise methods, but for high SNR levels, it is outperformed by the
latter two methods.  The critical transition point—the SNR value at which  
their relative test error curves cross—is different for the low and medium 
settings: for the low setting, it is around 1.22, and for the medium setting,
it is earlier, around 0.42.  The relaxed lasso, meanwhile, is competitive across
all SNR levels: at low SNR levels it matches the performance of the lasso at low
SNR levels, and at high SNR levels it matches that of best subset and forward
stepwise selection. It is able to do so by properly tuning the amount of
shrinkage (via its parameter γ) on the validation set.  Lastly, the
number of nonzero estimated coefficients from the four methods (bottom right
panel) is also revealing.  The lasso consistently delivers much denser
solutions; essentially, to optimize prediction error on the validation set, it
is forced to do so, as the sparser solutions along its path feature too
much shrinkage.  The relaxed lasso does not suffer from this issue, again thanks
to its ability to unshrink (move γ away from 1); it delivers
solutions that are just as sparse as those from best subset and forward stepwise
selection, except at the low SNR range.

The high-5 and high-10 settings, Figures <ref> and <ref>,
behave quite differently.  The high-5 setting (smaller n and 
smaller s) is more dire: the PVEs delivered by all methods—especially best  
subset selection—are negative for low SNR values, due to poor tuning 
on the validation set (had we chosen the null model for each 
method, the PVE would have been zero).  In both the high-5 and high-10
settings, we see that there is generally no reason, based on relative risk, 
relative test error, or PVE, to favor best subset selection or forward stepwise 
selection over the lasso.  At low SNR levels, best subset and forward stepwise
selection often have 
worse accuracy metrics (and certainly more erratic metrics); at high SNR
levels, these procedures do not show much of an advantage.  For best subset 
selection, it is quite possible that its performance at the high SNR range
would improve if we gave Gurobi a greater budget (than 3 minutes per problem
instance per subset size).  The relaxed lasso again performs the best overall,
with a noticeable gap in performance at the high SNR levels.  As is confirmed by
the number of nonzero  
coefficients plots, the lasso and best subset/forward stepwise selection achieve
similar accuracy in the high SNR range using two opposite strategies:
the former uses high-bias and low-variance estimates, and the latter uses
low-bias and high-variance estimates.  The relaxed lasso 
is most accurate by striking a favorable balance between these two polar
regimes.  











 §.§ Summary of results


As mentioned above, the results from our entire simulation suite can be found in  
a supplementary document.  Here is a high-level summary.  



  * An important caveat to emphasize upfront is that the Gurobi MIO algorithm
  for best subset selection was given 3 minutes per problem
  instance per subset size.  This practical restriction may have caused
  best subset selection to underperform, in particular, at the high SNR levels
  in the high-5 and high-10 settings. 


  * Forward stepwise selection and best subset selection perform quite
  similarly throughout (with the former being much faster). This does not agree
  with the results for forward stepwise in  <cit.>, where
  it performed quite poorly in comparison.  In talking with the third author,
  Rahul Mazumder, we have learned that this was due to the fact that forward
  stepwise in their study was tuned using AIC, rather than a separate
  validation set. So, when put on equal footing and allowed to select its
  tuning  parameter using validation data just as the other methods, we see that
  it performs quite comparably. 


  * The lasso gives better accuracy results than best subset selection in the
  low SNR range and worse accuracy than best subset in the high SNR range.  The
  transition point—the SNR level past which best subset outperforms the
  lasso—varies depending on the problem dimensions (n,p) predictor
  autocorrelation (ρ), and beta-type (1 through 5).  For the medium 
  setting, the transition point comes earlier than in the low setting.  For the
  high-5 and high-10 settings, the transition point often does not come at all
  (before an SNR of 6, which is the maximum value we considered).  As the
  predictor autocorrelation level increases, the transition point typically
  appears later (again, in some cases it does not come at all, e.g., for beta-type
  5 and autocorrelation ρ=0.7).


  * The relaxed lasso provides altogether the top accuracy results.  In 
  nearly all cases (across all SNR levels, and in all problem configurations) we
  considered, it performs as well as or better than all other methods.  We 
  conclude that it is able to use its auxiliary shrinkage parameter (γ)
  to get the “best of both worlds”: it accepts the heavy shrinkage from the
  lasso when such shrinkage is helpful, and reverses it when it is not. 


  * The proportion of variance explained plots remind us that, despite what
  may seem like large relative differences, the four methods under consideration
  do not have very different absolute performances in this intuitive
  and important metric.  It thus makes sense overall to favor the methods that
  are easy to compute. 




§ DISCUSSION


The recent work of <cit.>  has enabled the first large-scale
empirical examinations of best subset selection.  In this paper, we have
expanded and refined the simulations in their work, comparing best subset 
selection to forward stepwise selection, the lasso, and the relaxed lasso.  We
have found: (a) forward stepwise selection and best subset selection perform
similarly throughout; (b) best subset selection often loses to the lasso except
in the high SNR range; (c) the relaxed lasso achieves “the best of both
worlds” and performs on par with the best method in each scenario. 
We note that these comparisons are based on (various measures of)
out-of-sample prediction accuracy.  A different target, e.g., a measure of
support recovery, may yield different results.

Our R package bestsubset, designed to easily replicate all of   
the simulations in this work, or forge new comparisons, is available at 
<https://github.com/ryantibs/best-subset/>. 

agsm


160
140
120




Supplement to “Extended Comparisons of Best Subset Selection, 
  Forward Stepwise Selection, and the Lasso”



This supplementary document contains plots from the simulation suite
described in the paper “Extended Comparisons of Best Subset Selection,  Forward 
Stepwise Selection, and the Lasso”. 
The plots in Section 1 precisely follow the simulation format described in the
paper. Those in Section 2 follow an analogous format, except that the
tuning has been done using an “oracle”, rather than a validation set as in
Section 1. Specifically, the tuning parameter for each method in each scenario
is chosen to minimize the average risk over all of the repetitions. 

.tocapp
mainnone
appsubsubsection
‎‎




§ VALIDATION TUNING




 §.§ Low setting: n=100, p=10, s=5
 


  §.§.§ Relative risk (to null model)
 





  §.§.§ Relative test error (to Bayes)






  §.§.§ Proportion of variance explained






  §.§.§ Number of nonzero coefficients






 §.§ Medium setting: n=500, p=100, s=5
 


  §.§.§ Relative risk (to null model)
 





  §.§.§ Relative test error (to Bayes)






  §.§.§ Proportion of variance explained






  §.§.§ Number of nonzero coefficients






 §.§ High-5 setting: n=50, p=1000, s=5
 


  §.§.§ Relative risk (to null model)
 





  §.§.§ Relative test error (to Bayes)






  §.§.§ Proportion of variance explained






  §.§.§ Number of nonzero coefficients






 §.§ High-10 setting: n=100, p=1000, s=10
 


  §.§.§ Relative risk (to null model)
 





  §.§.§ Relative test error (to Bayes)






  §.§.§ Proportion of variance explained






  §.§.§ Number of nonzero coefficients






§ ORACLE TUNING




 §.§ Low setting: n=100, p=10, s=5
 


  §.§.§ Relative risk (to null model)
 





  §.§.§ Relative test error (to Bayes)






  §.§.§ Proportion of variance explained






  §.§.§ Number of nonzero coefficients






 §.§ Medium setting: n=500, p=100, s=5
 


  §.§.§ Relative risk (to null model)
 





  §.§.§ Relative test error (to Bayes)






  §.§.§ Proportion of variance explained






  §.§.§ Number of nonzero coefficients






 §.§ High-5 setting: n=50, p=1000, s=5
 


  §.§.§ Relative risk (to null model)
 





  §.§.§ Relative test error (to Bayes)






  §.§.§ Proportion of variance explained






  §.§.§ Number of nonzero coefficients






 §.§ High-10 setting: n=100, p=1000, s=10
 


  §.§.§ Relative risk (to null model)
 





  §.§.§ Relative test error (to Bayes)






  §.§.§ Proportion of variance explained






  §.§.§ Number of nonzero coefficients


