Dual Lasso Selector

The use of microarray technologies have become popular to monitor genome-wide expression changes in health and disease. Typically, a microarray data set is high dimensional in the sense, it usually has tens of thousands of gene expression profile(variables) but only tens or hundreds of subjects(observations). In microarray analysis, a group of genes sharing the same biological pathway tend to have highly correlated expression levels Segal et al. (2003) and the goal is to identify all(rather than a few) of them if they are related to the underlying biological process. This is one example where the need to select groups of correlated variables arises. In many applications it is required to identify all relevant correlated variables. In this paper, we consider the problem of model selection and estimation in sparse high dimensional linear regression models with strongly correlated variables.

We start with the standard linear regression model as

with response vector $\textbf{Y}_{n\times 1}$, design matrix $\textbf{X}_{n\times p}$, true underlying coefficient vector $\beta_{p\times 1}$ and error vector $\epsilon_{n\times 1}\sim N_{n}(0,I_{n})$. In particular, we consider the case of sparse high dimensional linear model $(p\gg n)$ with strong empirical correlation among few variables. The Lasso is a widely used regularized regression method to find sparse solutions, the lasso estimator is defined as

where $\lambda\geq 0$ is the regularization parameter that controls the amount of regularization. It is known that the Lasso tends to select a single variable from a group of strongly correlated variables even if many or all of these variables are important.

In presence of correlated predictors, the concept of clustering or grouping correlated predictors and then pursuing group-wise model fitting was proposed, see Bühlmann et al. (2012) and Gauraha (2016). When the dimension is very high or in case of overlapping clusters, finding an appropriate group structure remains as difficult as the original problem. We note that clustering followed by model fitting is computationally expensive, not reliable and do not scale for large, high-dimensional data sets, so we do not consider it further in this paper. An alternatively approach is simultaneous clustering and model fitting that involves combination of two different penalties. For example, Elastic Net (Zou and Hastie (2005)) is a combination of two regularization techniques, the $\ell_{2}$ regularization provides grouping effects and $\ell_{1}$ regularization produces sparse models. Therefore, the eNet selects or drops highly correlated variables together that depends on the amount of $\ell_{1}$ and $\ell_{2}$ regularization.

The influence of correlations on Lasso prediction has been studied in Hebiri and Lederer (2013) and van de Geer and Lederer (2013), and it is shown that Lasso prediction works well in presence of any degree of correlations with an appropriate amount of regularization. However, studies show that correlations are problematic for parameter estimation and variable selection. It has been proven that the design matrix must satisfy the following two conditions for the Lasso to perform exact variable selection: irrepresentability(IC) condition(Zhao and Yu (2006)) and beta-min condition(Bühlmann and van de Geer (2011)). Having highly correlated variables implies that the design matrix violates the IC, and the Lasso solution is not stable. When active covariates are highly correlated the Lasso solution is not unique and Lasso randomly selects one variable from correlated group. However, even in case of highly correlated variables the corresponding dual Lasso solution is always unique. The dual of the Lasso problem (2), as shown in Wang et al. (2015) is given by

The intuition drawn from the articles OSBORNE et al. (2000) and Wang et al. (2015) further motivates us to consider the Lasso optimal and its dual optimal solution together, that yields in selecting correlated active predictors.

Exploiting the fact about uniqueness of the dual Lasso solution, we propose a new variable selection procedure, the Dual Lasso Selector (DLS). For a given $\lambda$ and a Lasso estimator $\hat{\beta}_{Lasso}$, we can compute the corresponding Dual Lasso solution using KKT conditions. Basically, the DLS active set corresponds to the predictors that satisfies dual Lasso feasible boundary conditions (we discuss it in details in later section). We argue that correlations among active predictors are not problematic, and we define a new weaker condition on the design matrix that allows for correlation among active predictors, called Pseudo Irrepresentable Condition (PIC). We prove that the Pseudo Irrepresentable Condition is a necessary and sufficient condition for the proposed dual Lasso selector to select the true active set (under assumption of beta-min condition) with high probability. Moreover, we use the $\ell_{2}$ penalty (the Ridge regression, Hoerl and Kennard (1970)) which is known to perform best in case of correlated variables, to estimate the coefficients of the predictors selected by the dual Lasso selector. We call the combination of the two, the DLSelect+Ridge. Though, DLSelect+Ridge resembles the ”Ridge post Lasso” but it is conceptually different and behaves differently than the Lasso followed by the Ridge, especially in the presence of highly correlated variables. For example DLSelect+Ridge looks like Elastic-net, since both are combination of $\ell_{1}$ and $\ell_{2}$ penalties but Elastic-net is a combination of the Ridge Regression followed by the Lasso. In addition, Enet needs to cross-validate on a two-dimensional surface $O(k^{2})$ to select its the optimal regularization parameters, whereas DLSelect+Ridge needs to cross validated twice on one-dimensional surface $O(k)$, where k is the length of the search space for a regularization parameter.

Our contribution is summarized as follows:

1. 1.

   We briefly review the state-of-the-art methods of simultaneous clustering and model fitting using combination of penalties such as Elastic-net, OSCAR and Fused Lasso etc.

2. 2.

   We study the theoretical properties of the Lasso and its dual optimal solution together and we show that selection of active correlated variables is related to the dual feasible boundary conditions.

3. 3.

   By further exploiting the uniqueness property of the dual Lasso solution, we develop a variable selection algorithm to efficiently select the true active predictors (including correlated active predictors). we call this selection technique as the Dual Lasso Selector.

4. 4.

   We derive the Pseudo Irrepresentable Conditions (PIC) for the design matrix that allow for the correlation between active covariates, and we show that under assumption of PIC the dual Lasso selector is variable selection consistent.

5. 5.

   We propose a new combined approach, the DLSelect+Ridge: Dual Lasso selecting predictors and the Ridge estimating their coefficients.

6. 6.

   We study the theoretical properties of the combination DLSelect+Ridge.

7. 7.

   We implement the DLSelect+Ridge method and empirically compare it with existing methods like Lasso and Enet etc. in terms of variable selection consistency, prediction accuracy, estimation accuracy and time complexity (using various simulations and real data examples).

We have organized the rest of the article in the following manner. We start with background in section 2. In section 3, we present Dual Lasso Selector. We define PIC and discuss variable selection consistency under this assumption on the design matrix, in section 4. Section 5 is concerned with illustration of the proposed method on real and simulated data sets. Section 6 gives computational details. We shall provide some concluding remarks in section 7.

In this section, we state notations and assumptions, used throughout the paper.

We consider usual sparse high-dim linear regression model as given in 1 with $p\gg n$. For the design matrix $\textbf{X}\in\mathcal{R}^{n\times p}$, we represent rows by $x_{i}^{T}\in\mathbb{R}^{p},\ i=1,...,n$, and columns by $X_{j}^{T}\in\mathbb{R}^{n},\ i=1,...,p$. We assume that the design matrix $\textbf{X}_{n\times p}$ is fixed, the data is centred and the predictors are standardized, so that $\sum_{i=1}^{n}\textbf{Y}_{i}=0$, $\sum_{i=1}^{n}({X}_{j}){i}=0$ and $\frac{1}{n}\textbf{X}^{T}_{j}\textbf{X}_{j}=1$ for all $j=1,...,p$. We denote by

the true active set and cardinality of the set $s=|S|$, is called sparsity index. We assume that the true coefficient vector $\beta$ is sparse, that is $s\ll p$. We denote $\textbf{X}_{S}$ as the restriction of X to columns in $S$, and $\beta_{S}$ is the vector $\beta$ restricted to the support $S$, with 0 outside the support $S$. Without loss of generality we can assume that the first $s$ variables are the active variables, and we partition the covariance matrix, $C=\frac{1}{n}\textbf{X}^{T}\textbf{X}$, for the active and the redundant variables as follows.

Similarly the coefficient vector $\beta$ can be partitioned as $\left[\begin{array}[]{c}\beta_{1}\\ \beta_{2}\end{array}\right].$
The $\ell_{1}$-norm and $\ell_{2}$-norm (square) are defined as

Throughout the paper, we use the notation $\lambda_{1}>0$ for $\ell_{1}$ penalty and $\lambda_{2}>0$ for other penalty functions. For a vector $a\in\mathbb{R}^{p}$, we denote its sign vector as

We denote sub-gradient of $\ell_{1}$-norm evaluated at $\beta\in\mathbb{R}^{p}$, as $\tau\in\partial\|\beta\|_{1}$, where $\tau$ satisfies the following.

Given the huge literature on the use of Lasso-type penalties for variable selection, we provide only a brief overview here, with focus on previous approaches which are closely related to our work. In particular, we briefly review the Lasso, the Ridge and the state-of-the-art in simultaneous clustering and model fitting using combination of penalties for high-dimensional sparse linear models. In general, we define a penalized least squares method as follows.

where the penalty terms $\mathcal{P}_{method}(\beta,.)$ can be different for different methods depending on the type and number of penalties used. In the following we define various penalized least squares estimators in terms of penalties used by them, and we also mention their computational complexity, variable selection consistency and grouping effects of selecting and dropping highly correlated predictors together.

1. 1.

   Lasso: The Lasso method was proposed by Tibshirani (1996) and the lasso penalty is defined as

   $\displaystyle\mathcal{P}_{Lasso}(\beta,\lambda_{1})=\lambda_{1}\|\beta\|_{1}$

   (16)

   It uses the single $\ell_{1}$ penalty, and due to nature of the $\ell_{1}$ penalty it simultaneously performs variable selection and estimation. The whole regularization path can be computed efficiently with the computational effort of a single OLS fit, by some modification of the LARS algorithm, see Efron et al. (2004). It does not provide grouping effect, in fact the Lasso tends to select a single predictor from a group of highly correlated predictors.

2. 2.

   Ridge Regression (RR): The ridge method was proposed by Hoerl and Kennard (1970) and the ridge penalty is defined as follows.

   $\displaystyle\mathcal{P}_{Ridge}(\beta,\lambda_{2})=\lambda_{2}\|\beta\|_{2}$

   (17)

   It uses the single $L_{2}$ penalty, and it always has a unique solution for a fixed regularization parameter $\lambda_{2}$. Though it is known to correctly detect the variable signs with reduced mean square error with correlated variables, but does not provide variable selection. It provides grouping effect with the highly correlated variables, and the computational complexity of the ridge is same as the computational effort of a single OLS fit.

3. 3.

   Elastic-net (Enet): The Enet method was introduced by Zou and Hastie (2005). The Enet penalty is a combination of $\ell_{1}$ and $L_{2}$ penalties and it is defined as follows.

   $\displaystyle\mathcal{P}_{Enet}(\beta,\lambda_{1},\lambda_{2})=\lambda_{1}\|\beta\|_{1}+\lambda_{2}\|\beta\|_{2}$

   (18)

   Enet addresses both the limitations of the Lasso, that is it can select correlated predictors as well as it can handle the $s>n$ case. It provides grouping effect, but requires to search in two-dimensional space for choosing optimal values of its regularization parameters. Hence its effective time complexity depends on the length of the search space for the regularization parameters $\lambda_{1}$ and $\lambda_{2}$.

4. 4.

   Correlation Based Penalty (CP): The correlation based penalized least squares method was proposed by Tutz and Ulbricht (2009), which uses the following correlation-based penalty term

   $\displaystyle\mathcal{P}_{CP}(\beta,\lambda_{2})=\lambda_{2}\sum_{i=1}^{p-1}\sum_{j>i}\left\{\frac{(\beta_{i}-\beta_{j})^{2}}{1-\rho_{ij}}+\frac{(\beta_{i}+\beta_{j})^{2}}{1+\rho_{ij}}\right\}$

   It uses the single $CP$ penalty norm, and it always has a unique solution for a fixed regularization parameter and the grouping effect strongly depends on the convexity of the penalty term. It does not provide variable selection. However, a boosted version of the penalized estimator allows to select variables. But the major drawback is that it is not scalable for large high dimensional problems.

5. 5.

   Fused Lasso: The Fused Lasso method was given by Tibshirani et al. (2005), and the Fused Lasso penalty is defined as

   $\displaystyle\mathcal{P}_{Fused}(\beta,\lambda_{1},\lambda_{2})=\lambda_{1}\|\beta\|_{1}+\lambda_{2}\sum_{j=2}^{p}|\beta_{j}-\beta_{j-1}|$

   (19)

   The first constraint encourages sparsity in the coefficients and the second constraint encourages sparsity in their differences. The major drawback of this method is that it requires the covariates to be in some order. It does not perform automated variable clustering to unordered features.

6. 6.

   OSCAR: The OSCAR was invented by Bondell and Reich (2008), and the OSCAR penalty is given as follows.

   $\displaystyle\mathcal{P}_{Oscar}(\beta,\lambda_{1},\lambda_{2})=\lambda_{1}\|\beta\|_{1}+\lambda_{2}\sum_{i<j}max\{|\beta_{i}|,|\beta_{j}|\}$

   with $|\beta_{1}|\leq...\leq|\beta_{p}|$. The first constraint is to encourage sparsity in the coefficients and the second constraint encourages equi-sparsity in $|\beta|$. The time complexity limits its scalability on ultra high-dimensional problems, moreover it requires two-dimensional grid search over the two parameters $(\lambda_{1},\lambda_{2})$

7. 7.

   L1CP : The L1CP penalty term is given by as follows, see Anbari and Mkhadri (2014).

   $\displaystyle\mathcal{P}_{L1CP}(\beta,\lambda_{1},\lambda_{2})=\lambda_{1}\|\beta\|_{1}+\lambda_{2}\sum_{i=1}^{p-1}\sum_{j>i}\left\{\frac{(\beta_{i}-\beta_{j})^{2}}{1-\rho_{ij}}+\frac{(\beta_{i}+\beta_{j})^{2}}{1+\rho_{ij}}\right\}.$

   (20)

   It performs variable selection with grouping effect and estimation together but is not scalable to the large scale problems due to expensive computation time, and it also requires two-dimensional grid search over the two parameters $(\lambda_{1},\lambda_{2})$.

8. 8.

   Clustered Lasso: The Clustered Lasso penalty is defined as

   $\displaystyle\mathcal{P}_{CL}(\beta,\lambda_{1},\lambda_{2})=\lambda_{1}\|\beta\|_{1}+\lambda_{2}\sum_{i<j}|\beta_{i}-\beta_{j}|$

   The first constraint encourages sparsity in the coefficients and the second constraint encourages equi-sparsity in $|\beta|$, It is similar as the Fused lasso but does not require ordering of variables, see She (2010). It provides grouping effect, but it requires two-dimensional grid search over the two parameters $(\lambda_{1},\lambda_{2})$. It is computationally expensive since it has to check equi-sparsity pattern for each pair of variables.

In the following table we summarize the properties discussed above for various regularization methods.

In this section, we present the dual Lasso selector, a new variable selection method for sparse high-dim regression models with correlated variables. First, we study the theoretical properties of the Lasso and dual Lasso solutions. Then, we show that the magnitude of correlations between the predictors and the dual vector determines the set of active predictors. This is the basis for our correlated variable selection.

The dual problem of the Lasso problem (2) can be given as follows (we provide the detailed derivation of the Lasso’s dual in the appendix A.1):

where $\theta$ is the dual vector, as defined in equation (44). For a fixed $\lambda\geq 0$, let $\hat{\beta}_{lasso}(\lambda)$ and $\hat{\theta}(\lambda)$ denote the optimal solutions of the Lasso and its dual problem respectively. Since it is implicit that the Lasso and its dual optimal depends on the $\lambda$, we drop the term $\lambda$ from the expression for notational simplicity. From KKT conditions (derivation of the KKT conditions is given in the appendix A.2) we get the following primal dual relationship:

It is worth mentioning the basic properties of the Lasso and its dual, which has already been derived and studied by various authors (see Tibshirani and Taylor (2011) and Wang et al. (2015) for more insights).

1. 1.

   Uniqueness of the Lasso-fit: There may not be a unique solution for the Lasso problem because for the criterion (42) is not strictly convex in $\beta$. But the least square loss is strictly convex in $\textbf{X}\beta$, hence there is always a unique fitted value $\textbf{X}\hat{\beta}$.

2. 2.

   Uniqueness of the dual vector: The dual problem is strictly convex in $\theta$, therefore the dual optimal $\hat{\theta}$ is unique. Another argument for the uniqueness of $\hat{\theta}$ is that it is a function of $\textbf{X}\hat{\beta}$ (23) which itself is unique. The fact that the DLS can achieve consistent variable selection for situations (with correlated active predictors) when the Lasso is unstable for estimation of the true active set is related to the uniqueness of the dual Lasso solution.

3. 3.

   Uniqueness of the Sub-gradient: Sub-gradient of $\ell_{1}$ norm of any Lasso solution $\hat{\beta}$ is unique because it is a function of $\textbf{X}\hat{\beta}$ (see Appendix A.2). More specifically, suppose $\hat{\beta}$ and $\tilde{\beta}$ are two lasso solutions for a fixed $\lambda$ value, then they must have the same signs $sign(\hat{\beta})=sign(\tilde{\beta})$, it is not possible that $\hat{\beta}_{j}>0$ and $\hat{\beta}_{j}<0$ for some $j$.

Let $\hat{S}_{lasso}$ denote the support set or active set of the Lasso estimator $\hat{\beta}$ which is given as

Similarly, we define the active set of the dual Lasso vector that corresponds to the active constraints of the dual optimization problem. We note that constraints are said to be active at a feasible point if that point lies on a boundary formed by the constraint.

Now, we define the following lemmas that will be used later for our mathematical derivations.

It is known that Irrepresentable condition (assuming beta-min conditions holds) is necessary and sufficient condition for the Lasso to select true model, see Zhao and Yu (2006) (for completeness we have proved it in appendix A.4).

The proof is worked out in Appendix A.3.

The IC may fail to hold due to violation of any one of (or both) the following two conditions:

1. 1.

   When $C_{11}$ is not invertible, that implies there is strong correlation among variables of the true active set.

2. 2.

   The active predictors are correlated with the noise features (this situation is better explained in terms of irrepresentable condition).

When there is strong correlation among variables of the active set, then $C_{11}$ is not invertible and the IC does not hold, and the Lasso fails to do variable selection. But we argue that the dual Lasso can still perform variable selection consistently even when $C_{11}$ not invertible, when we impose some milder condition on the design matrix, we call it Pseudo Irrepresentable Condition (PIC). The Pseudo Irrepresentable Condition is defined as follows.

The following lemma gives the sufficient condition for the dual Lasso for support recovery. This lemma is similar in spirit of the Lemma 2 define in Omidiran and Wainwright (2010). Here, we do not assume that $\Sigma_{11}$ is invertible.

We have shown that the dual Lasso optimal $\hat{\theta}$ is always unique, and it remains to show that the $\hat{S}_{dual}$ recovers the true active set. Under the assumption (32), we can derive that $|X_{j}^{T}\hat{\theta}|<\lambda$ for all $j\in S^{c}$. Therefore $\hat{S}_{dual}=S$.

When $C_{11}$ is invertible the PIC coincides with the IC, and under assumption of the IC we have already shown that $\hat{S}_{dual}=\hat{S}_{lasso}$. In Appendix A.4, we prove that the PIC is necessary and sufficient condition (beta-min condition is implicit) for the dual Lasso to consistently select the true active set. The PIC may hold even when $C_{11}$ is not invertible, which implies that the PIC is weaker than the IC. We illustrate it with the following examples:

Let $S=\{1,2,3,4\}$ be the active set, the covariance matrix $C=\frac{\textbf{X}^{T}\textbf{X}}{n}$ and is given as

where the active variables are uncorrelated and the noise variable is equally correlated with all active covariates. First of all, it is easy to check that only for $|\rho|\leq\frac{1}{2}$, $C$ is positive semi definite, and for $|\rho|<\frac{1}{4}$, $C$ satisfies the IC.

Now, we augment this matrix with two additional columns, one copy of the first and second active variables, and we rearrange the columns such that we get the following covariance matrix, and we redefine the set of active variables as $S=\{1,2,3,4,5,6\}$ and we assume that $|\rho|<\frac{1}{4}$.

We partition $\Sigma$ as (7), and it is clear that $C_{11}$ is not invertible and IC does not hold, hence the Lasso does not perform variable selection. The rank of the $C_{11}$ is 4. Let us consider any $(4\times 4)$ sub matrix of $C_{11}$ such that its rank is four ($S_{1}\subset S,andrank(S_{1})=4,S_{1}=\{\{1,3,5,6\},\{1,4,5,6\},\{2,3,5,6\},\{2,4,5,6\}$). Further, we partition $C_{11}$ as

where $rank(C_{rr})$ has full column rank and $rank(C_{rr})=rank(C_{11})$, let $C_{R}$ be set of are four such possible invertible sub matrices of $C_{11}$. Then considering the generalized inverses corresponding to them as

where $C_{A}\in C_{R}$ is invertible. With the above inverse $C_{11}^{+}$ the PIC holds for the design matrix X. It can be also viewed as the IC is satisfied for each reduced active set $S^{\prime}\in S_{1}$ and the corresponding reduced design matrix $\textbf{X}_{S^{\prime}}$, and hence the Lasso picks randomly one element from the set $S_{1}$ and sets the coefficient of the noise variable to zero (with high probability). Also, since PIC holds, the dual Lasso will select the true active set $S$ with high probability and will set zero for the coefficient of noise feature.

In this section, we apply the DLSelect+Ridge for variable selection and estimation on simulations and real data and compare the results with that of the Lasso, Ridge and Elastic-net. We consider the True Positive Rate (power) and False Discovery Rate (FDR) as the measure of performances for variable selection, which are defined as follows.

For prediction performance we consider the Mean Squared Prediction Error, which is defined as

where $\hat{\textbf{Y}}$ is the predicted response vector or an estimate $\textbf{X}\hat{\beta}$ based on an estimator $\hat{\beta}$.

Since our aim is to avoid false negatives, we do not report false positives, and ridge does not perform variable selection therefore TPR is not reported for the Ridge. The Ridge is considered as a competitor because its prediction performance is better than the Lasso for correlated designs.

Statistical analysis was performed in R 3.2.2. We used the package “glmnet” for penalized regression method(the Lasso).

The main achievements of this work are summarized as follows: We argued that the correlations among active predictors is not problematic, as long as the PIC is satisfied by the design matrix. In particular, we proved that the dual Lasso performs consistent variable selection under assumption of PIC. Exploiting this result we proposed the dual Lasso+Ridge method. We illustrated DLSelection+Ridge method on simulated and real high dimensional data sets. The numerical studies based on the simulations and real examples show clearly that the proposed method is very competitive in terms of variable selection, prediction accuracy, estimation accuracy and computation speed.