Bolasso: model consistent Lasso estimation through the bootstrap

Introduction

In this paper, we first derive a detailed asymptotic analysis of sparsity pattern selection of the Lasso estimation procedure, that extends previous analysis [Zhaoyu, yuanlin, zou], by focusing on a specific decay of the regularization parameter. We show that when the decay is proportional to $n^{-1/2}$, where $n$ is the number of observations, then the Lasso will select all the variables that should enter the model (the relevant variables) with probability tending to one exponentially fast with $n$, while it selects all other variables (the irrelevant variables) with strictly positive probability. If several datasets generated from the same distribution were available, then the latter property would suggest to consider the intersection of the supports of the Lasso estimates for each dataset: all relevant variables would always be selected for all datasets, while irrelevant variables would enter the models randomly, and intersecting the supports from sufficiently many different datasets would simply eliminate them. However, in practice, only one dataset is given; but resampling methods such as the bootstrap are exactly dedicated to mimic the availability of several datasets by resampling from the same unique dataset [efron]. In this paper, we show that when using the bootstrap and intersecting the supports, we actually get a consistent model estimate, without the consistency condition required by the regular Lasso. We refer to this new procedure as the Bolasso (bootstrap-enhanced least absolute shrinkage operator). Finally, our Bolasso framework could be seen as a voting scheme applied to the supports of the bootstrap Lasso estimates; however, our procedure may rather be considered as a consensus combination scheme, as we keep the (largest) subset of variables on which all regressors agree in terms of variable selection, which is in our case provably consistent and also allows to get rid of a potential additional hyperparameter.

The paper is organized as follows: in Section 2, we present the asymptotic analysis of model selection for the Lasso; in Section 3, we describe the Bolasso algorithm as well as its proof of model consistency, while in Section 4, we illustrate our results on synthetic data, where the true sparse generating model is known, and data from the UCI machine learning repository. Sketches of proofs can be found in Appendix A.

Asymptotic Analysis of Model Selection for the Lasso

In this section, we describe existing and new asymptotic results regarding the model selection capabilities of the Lasso.

We let denote ${\mathbf{J}}=\{j,{\mathbf{w}}_{j}\neq 0\}$ the sparsity pattern of ${\mathbf{w}}$, ${\mathbf{s}}=\mathop{\rm sign}({\mathbf{w}})$ the sign pattern of ${\mathbf{w}}$, and $\varepsilon=Y-X^{\top}{\mathbf{w}}$ the additive noise. Throughout this paper, we use boldface fonts for population quantities. Note that our assumption regarding cumulant generating functions is satisfied when $X$ and $\varepsilon$ have compact support, and also, when the densities of $X$ and $\varepsilon$ have light tails.

Note that the i.i.d. assumption, together with (A(A1)-(A3)), are the simplest assumptions for studying the asymptotic behavior of the Lasso; and it is of course of interest to allow more general assumptions, in particular growing number of variables $p$, more general random variables, etc. (see, e.g., (?)), which are outside the scope of this paper.

2 Lasso Estimation

where $\mu_{n}\geqslant 0$ is the regularization parameter. We denote $\hat{w}$ any global minimum of Eq. (1)—it may not be unique in general, but will with probability tending to one exponentially fast under assumption (A(A1)).

3 Model Consistency - General Results

In this section, we detail the asymptotic behavior of the Lasso estimate $\hat{w}$, both in terms of the difference in norm with the population value ${\mathbf{w}}$ (i.e., regular consistency) and of the sign pattern $\mathop{\rm sign}(\hat{w})$, for all asymptotic behaviors of the regularization parameter $\mu_{n}$. Note that information about the sign pattern includes information about the support, i.e., the indices $i$ for which $\hat{w}_{i}$ is different from zero; moreover, when $\hat{w}$ is consistent, consistency of the sign pattern is in fact equivalent to the consistency of the support.

We now consider five mutually exclusive possible situations which explain various portions of the regularization path (we assume (A(A1)-(A3))); many of these results appear elsewhere [yuanlin, Zhaoyu, fu, zou, grouplasso] but some of the finer results presented below are new (see Section 2.4).

If $\mu_{n}$ tends to infinity, then $\hat{w}=0$ with probability tending to one.

If $\mu_{n}$ tends to a finite strictly positive constant $\mu_{0}$, then $\hat{w}$ converges in probability to the unique global minimum of $\frac{1}{2}(w-{\mathbf{w}})^{\top}\mathbf{Q}(w-{\mathbf{w}})+\mu_{0}\|w\|_{1}$. Thus, the estimate $\hat{w}$ never converges in probability to ${\mathbf{w}}$, while the sign pattern tends to the one of the previous global minimum, which may or may not be the same as the one of ${\mathbf{w}}$.Here and in the third regime, we do not take into account the pathological cases where the sign pattern of the limit in unstable, i.e., the limit is exactly at a hinge point of the regularization path.

If $\mu_{n}$ tends to zero slower than $n^{-1/2}$, then $\hat{w}$ converges in probability to ${\mathbf{w}}$ (regular consistency) and the sign pattern converges to the sign pattern of the global minimum of $\frac{1}{2}v^{\top}\mathbf{Q}v+v_{\mathbf{J}}^{\top}\mathop{\rm sign}({\mathbf{w}}_{\mathbf{J}})+\|v_{{\mathbf{J}}^{c}}\|_{1}$. This sign pattern is equal to the population sign vector ${\mathbf{s}}=\mathop{\rm sign}({\mathbf{w}})$ if and only if the following consistency condition is satisfied:

Thus, if Eq. (2) is satisfied, the probability of correct sign estimation is tending to one, and to zero otherwise [yuanlin].

If $\mu_{n}=\mu_{0}n^{-1/2}$ for $\mu_{0}\in(0,\infty)$, then the sign pattern of $\hat{w}$ agrees on ${\mathbf{J}}$ with the one of ${\mathbf{w}}$ with probability tending to one, while for all sign patterns consistent on ${\mathbf{J}}$ with the one of ${\mathbf{w}}$, the probability of obtaining this pattern is tending to a limit in $(0,1)$ (in particular strictly positive); that is, all patterns consistent on ${\mathbf{J}}$ are possible with positive probability. See Section 2.4 for more details.

Among the five previous regimes, the only ones with consistent estimates (in norm) and a sparsity-inducing effect are $\mu_{n}$ tending to zero and $\mu_{n}n^{1/2}$ tending to a limit $\mu_{0}\in(0,\infty]$ (i.e., potentially infinite). When $\mu_{0}=+\infty$, then we can only hope for model consistent estimates if the consistency condition in Eq. (2) is satisfied. This somewhat disappointing result for the Lasso has led to various improvements on the Lasso to ensure model consistency even when Eq. (2) is not satisfied [yuanlin, zou]. Those are based on adaptive weights based on the non regularized least-square estimate. We propose in Section 3 an alternative way which is based on resampling.

In this paper, we now consider the specific case where $\mu_{n}=\mu_{0}n^{-1/2}$ for $\mu_{0}\in(0,\infty)$, where we derive new asymptotic results. Indeed, in this situation, we get the correct signs of the relevant variables (those in ${\mathbf{J}}$) with probability tending to one, but we also get all possible sign patterns consistent with this, i.e., all other variables (those not in ${\mathbf{J}}$) may be non zero with asymptotically strictly positive probability. However, if we were to repeat the Lasso estimation for many datasets obtained from the same distribution, we would obtain for each $\mu_{0}$, a set of active variables, all of which include ${\mathbf{J}}$ with probability tending to one, but potentially containing all other subsets. By intersecting those, we would get exactly ${\mathbf{J}}$.

However, this requires multiple copies of the samples, which are not usually available. Instead, we consider bootstrapped samples which exactly mimic the behavior of having multiple copies. See Section 3 for more details.

4 Model Consistency with Exact Root-n𝑛n Regularization Decay

In this section we present detailed new results regarding the pattern consistency for $\mu_{n}$ tending to zero exactly at rate $n^{-1/2}$ (see proofs in Appendix A):

Assume (A(A1)-(A3)) and $\mu_{n}=\mu_{0}n^{-1/2}$, $\mu_{0}>0$. Then, for any pattern $s\in\{-1,0,1\}^{p}$ such that $s_{\mathbf{J}}\neq\mathop{\rm sign}({\mathbf{w}}_{\mathbf{J}})$, there exist a constant $A(\mu_{0})>0$ such that

The last two propositions state that we get all relevant variables with probability tending to one exponentially fast, while we get exactly get all other patterns with probability tending to a limit strictly between zero and one. Note that the results that we give in this paper are valid for finite $n$, i.e., we could derive actual bounds on probability of sign pattern selections with known constants that explictly depend on ${\mathbf{w}}$, $\mathbf{Q}$ and $P_{XY}$.

Bolasso: Bootstrapped Lasso

The asymptotic analysis from Section 2 suggests to estimate the supports $J_{k}=\{j,\ \hat{w}^{k}_{j}\neq 0\}$ of the Lasso estimates $\hat{w}^{k}$ for the bootstrap samples, $k=1,\dots,m$, and to intersect them to define the Bolasso model estimate of the support: $J=\bigcap_{k=1}^{m}J_{k}$. Once $J$ is selected, we estimate $w$ by the unregularized least-square fit restricted to variables in $J$. The detailed algorithm is given in Algorithm 1. The algorithm has only one extra parameter (the number of bootstrap samples $m$). Following Proposition 3, $\log(m)$ should be chosen growing with $n$ asymptotically slower than $n$. In simulations, we always use $m=128$ (except in Figure 3, where we exactly study the influence of $m$).

Note that in practice, the Bolasso estimate can be computed simultaneously for a large number of regularization parameters because of the efficiency of the Lars algorithm (which we use in simulations), that allows to find the entire regularization path for the Lasso at the (empirical) cost of a single matrix inversion [lars]. Thus computational complexity of the Bolasso is $O(m(p^{3}+p^{2}n))$.

The following proposition (proved in Appendix A) shows that the previous algorithm leads to consistent model selection.

where $A_{1},A_{2},A_{3},A_{4}$ are strictly positive constants.

Therefore, if $\log(m)$ tends to infinity slower than $n$ when $n$ tends to infinity, the Bolasso asymptotically selects with overwhelming probability the correct active variable, and by regular consistency of the restricted least-square estimate, the correct sign pattern as well. Note that the previous bound is true whether the condition in Eq. (2) is satisfied or not, but could be improved on if we suppose that Eq. (2) is satisfied. See Section 4.1 for a detailed comparison with the Lasso on synthetic examples.

Simulations

In this section, we illustrate the consistency results obtained in this paper with a few simple simulations on synthetic examples similar to the ones used by (?) and some medium scale datasets from the UCI machine learning repository [UCI].

In Figure 1, we sampled two distributions $P_{XY}$ with $p=16$ and $r=8$ relevant variables, one for which the consistency condition in Eq. (2) is satisfied (left), one for which it was not satisfied (right). For a fixed number of sample $n=1000$, we generated 256 replications and computed the empirical frequencies of selecting any given variable for the Lasso as the regularization parameter $\mu$ varies. Those plots show the various asymptotic regimes of the Lasso detailed in Section 2. In particular, on the right plot, although no $\mu$ leads to perfect selection (i.e., exactly variables with indices less than $r=8$ are selected), there is a range where all relevant variables are always selected, while all others are selected with probability within $(0,1)$.

In Figure 2, we plot the results under the same conditions for the Bolasso (with a fixed number of bootstrap replications $m=128$). We can see that in the Lasso-consistent case (left), the Bolasso widens the consistency region, while in the Lasso-inconsistent case (right), the Bolasso “creates” a consistency region.

In Figure 3, we selected the same two distributions and compared the probability of exactly selecting the correct support pattern, for the Lasso, and for the Bolasso with varying numbers of bootstrap replications (those probabilities are computed by averaging over 256 experiments with the same distribution). In Figure 3, we can see that in the Lasso-inconsistent case (right), the Bolasso indeed allows to fix the unability of the Lasso to find the correct pattern. Moreover, increasing $m$ looks always beneficial; note that although it seems to contradict the asymptotic analysis in Section 3 (which imposes an upper bound for consistency), this is due to the fact that not selecting (at least) the relevant variables has very low probability and is not observed with only 256 replications.

Finally, in Figure 4, we compare various variable selection procedures for linear regression, to the Bolasso, with two distributions where $p=64$, $r=8$ and varying $n$. For all the methods we consider, there is a natural way to select exactly $r$ variables with no free parameters (for the Bolasso, we select the most stable pattern with $r$ elements, i.e., the pattern which corresponds to most values of $\mu$). We can see that the Bolasso outperforms all other variable selection methods, even in settings where the number of samples becomes of the order of the number of variables, which requires additional theoretical analysis, subject of ongoing research. Note in particular that we compare with bagging of least-square regression [bagging] followed by a thresholding of the loading vector, which is another simple way of using bootstrap samples: the Bolasso provides a more efficient way to use the extra information, not for usual stabilization purposes [stabilizing], but directly for model selection. Note finally, that the bagging of Lasso estimates requires an additional parameter and is thus not tested.

2 UCI datasets

The previous simulations have shown that the Bolasso is succesful at performing model selection in synthetic examples. We now apply it to several linear regression problems and compare it to alternative methods for linear regression, namely, ridge regression, Lasso, bagging of Lasso estimates [bagging], and a soft version of the Bolasso (referred to as Bolasso-S), where instead of intersecting the supports for each bootstrap replications, we select those which are present in at least $90\%$ of the bootstrap replications. In Table 1, we consider data randomly generated as in Section 4.1 (with $p=32$, $r=8$, $n=64$), where the true model is known to be composed of a sparse loading vector, while in Table 2, we consider regression datasets from the UCI machine learning repository. For all of those, we perform 10 replications of 10-fold cross validation and for all methods (which all have one free regularization parameter), we select the best regularization parameter on the 100 folds and plot the mean square prediction error and its standard deviation.

Note that when the generating model is actually sparse (Table 1), the Bolasso outperforms all other models, while in other cases (Table 2) the Bolasso is sometimes too strict in intersecting models, i.e., the softened version works better and is competitive with other methods. Studying the effects of this softened scheme (which is more similar to usual voting schemes), in particular in terms of the potential trade-off between good model selection and low prediction error, and under conditions where $p$ is large, is the subject of ongoing work.

Conclusion

We have presented a detailed analysis of variable selection properties of a boostrapped version of the Lasso. The model estimation procedure, referred to as the Bolasso, is provably consistent under general assumptions. This work brings to light that poor variable selection results of the Lasso may be easily enhanced thanks to a simple parameter-free resampling procedure. Our contribution also suggests that the use of bootstrap samples by L. Breiman in Bagging/Arcing/Random Forests [arcing] may have been so far slightly overlooked and considered a minor feature, while using boostrap samples may actually be a key computational feature in such algorithms for good model selection performances, and eventually good prediction performances on real datasets.

Appendix A Proof of Model Consistency Results

In this appendix, we give sketches of proofs for the asymptotic results presented in Section 2 and Section 3. The proofs rely on the well-known property of the Lasso optimization problems, namely that if the sign pattern of the solution is known, then we can get the solution in closed form.

The optimality conditions for Eq. (3) can be written in terms of the sign pattern $s=s(w)=\mathop{\rm sign}(w)$ and the sparsity pattern $J=J(w)=\{j,\ w_{j}\neq 0\}$ [yuanlin]:

In this paper, we focus on regularization parameters $\mu_{n}$ of the form $\mu_{n}=\mu_{0}n^{-1/2}$. The main idea behind the results is to consider that $(Q,q)$ are distributed according to their limiting distributions, obtained from the law of large numbers and the central limit theorem, i.e., $Q$ converges to $\mathbf{Q}$ a.s. and $n^{1/2}q$ is asymptotically normally distributed with mean zero and covariance matrix $\sigma^{2}\mathbf{Q}$. When assuming this, Propositions 1 and 2 are straightforward. The main effort is to make sure that we can safely replace $(Q,q)$ by their limiting distributions. The following lemmas give sufficient conditions for correct estimation of the signs of variables in ${\mathbf{J}}$ and for selecting a given pattern $s$ (note that all constants could be expressed in terms of $\mathbf{Q}$ and ${\mathbf{w}}$, details are omitted here):

Assume (A(A2)) and $\|Q-\mathbf{Q}\|_{2}\leqslant\lambda_{\min}(\mathbf{Q})/2$. Then $\mathop{\rm sign}(\hat{w}_{\mathbf{J}})\neq\mathop{\rm sign}({\mathbf{w}}_{\mathbf{J}})$ implies $\|\mathbf{Q}^{-1/2}q\|_{2}\geqslant C_{1}-\mu_{n}C_{2}$, where $C_{1},C_{2}>0$.

Assume (A(A2)) and let $s\in\{-1,0,1\}^{p}$ such that $s_{\mathbf{J}}=\mathop{\rm sign}({\mathbf{w}}_{\mathbf{J}})$. Let $J=\{j,s_{j}\neq 0\}\supset{\mathbf{J}}$. Assume

with $C_{4},C_{5},C_{6},C_{7},C_{8}$ are positive constants. Then $\mathop{\rm sign}(\hat{w})=\mathop{\rm sign}({\mathbf{w}})$.

Those two lemmas are interesting because they relate optimality of certain sign patterns to quantities from which we can derive concentration inequalities.

A.2 Concentration Inequalities

A.3 Proof of Proposition 1

By Lemma 2, for any given $A$, and $n$ large enough, the probability that the sign is different from $s$ is upperbounded by

A.4 Proof of Proposition 2

A.5 Proof of Proposition 3

In order to simplify the proof, we made the simplifying assumption that the random variables $X$ and $\varepsilon$ have compact supports. Extending the proofs to take into account the looser condition that $\|X\|^{2}$ and $\varepsilon^{2}$ have non uniformly infinite cumulant generating functions (i.e., assumption (A(A1))) can be done with minor changes. The probability that $\bigcap_{k=1}^{m}J_{k}$ is different from ${\mathbf{J}}$ is upper bounded by the sum of the following probabilities:

(b) Selecting at most variables in 𝐉𝐉{\mathbf{J}}:

the probability that for all replications, the set ${\mathbf{J}}$ is not exactly selected (note that this is not tight at all since on top of the relevant variables which are selected with overwhelming probability, different additional variables may be selected for different replications and cancel out when intersecting).

References