If multi-collinearity is high, would LASSO coefficients shrink to 0?

Question

Given $x_2 = 2 x_1$, what's the theoretical behavior of LASSO coefficients and why?

Would one of $x_1$ or $x_2$ shrink to $0$ or both of them?

require(glmnet)
x1 = runif(100, 1, 2)
x2 = 2*x1
x_train = cbind(x1, x2)
y = 100*x1 + 100 + runif(1)
ridge.mod = cv.glmnet(x_train, y, alpha = 1)
coef(ridge.mod)

#3 x 1 sparse Matrix of class "dgCMatrix"
#                       1
#(Intercept) 1.057426e+02
#x1          9.680073e+01
#x2          3.122502e-15

Answer 1

Notice that \begin{align*} \|y-X\beta\|_2^2 + \lambda \|\beta\|_1 & = \|y - \beta_1 x_1 - \beta_2 x_2 \|_2^2 + \lambda \left( |\beta_1| + |\beta_2| \right) \\ & = \|y - (\beta_1 + 2 \beta_2) x_1 \|_2^2 + \lambda \left( |\beta_1| + |\beta_2| \right). \end{align*}

For any fixed value of the coefficient $\beta_1 + 2\beta_2$, the penalty $|\beta_1| + |\beta_2|$ is minimized when $\beta_1 = 0$. This is because the penalty on $\beta_1$ is twice as weighted! To put this in notation, $$\tilde\beta = \arg\min_{\beta \, : \, \beta_1 + 2\beta_2 = K}|\beta_1| + |\beta_2|$$ satisfies $\tilde\beta_1 = 0$ for any $K$. Therefore, the lasso estimator \begin{align*} \hat\beta & = \arg\min_{\beta \in \mathbb{R}^p} \|y - X \beta\|_2^2 + \lambda \|\beta\|_1 \\ & = \arg\min_{\beta \in \mathbb{R}^p} \|y - (\beta_1 + 2 \beta_2) x_1 \|_2^2 + \lambda \left( |\beta_1| + |\beta_2| \right) \\ & = \arg_\beta \min_{K \in \mathbb{R}} \, \min_{\beta \in \mathbb{R}^p \, : \, \beta_1 + 2 \beta_2 = K} \, \|y - K x_1 \|_2^2 + \lambda \left( |\beta_1| + |\beta_2| \right) \\ & = \arg_\beta \min_{K \in \mathbb{R}} \, \left\{ \|y - K x_1 \|_2^2 + \lambda \min_{\beta \in \mathbb{R}^p \, : \, \beta_1 + 2 \beta_2 = K} \, \left\{ \left( |\beta_1| + |\beta_2| \right) \right\} \right\} \end{align*} satisfies $\hat\beta_1 = 0$. The reason why the comments to OP's question are misleading is because there's a penalty on the model: those $(0, 50)$ and $(100,0)$ coefficients give the same error, but different $\ell_1$ norm! Further, it's not necessary to look at anything like LARs: this result follows immediately from the first principles.

As pointed out by Firebug, the reason why your simulation shows a contradictory result is that glmnet automatically scales to unit variance the features. That is, due to the use of glmnet, we're effectively in the case that $x_1 = x_2$. There, the estimator is no longer unique: $(100,0)$ and $(0,100)$ are both in the arg min. Indeed, $(a,b)$ is in the $\arg\min$ for any $a,b \geq 0$ such that $a+b = 100$.

In this case of equal features, glmnet will converge in exactly one iteration: it soft-thresholds the first coefficient, and then the second coefficient is soft-thresholded to zero.

This explains why the the simulation found $\hat\beta_2 = 0$ in particular. Indeed, the second coefficient will always be zero, regardless of the ordering of the features.

Proof: Assume WLOG that the feature $x \in \mathbb{R}^n$ satisfies $\|x\|_2 = 1$. Coordinate descent (the algorithm used by glmnet) computes for it's first iteration: $$\hat\beta_1^{(1)} = S_\lambda(x^T y)$$ followed by \begin{align*} \hat\beta_2^{(1)} & = S_\lambda \left[ x^T \left( y - x S_\lambda (x^T y) \right) \right] \\ & = S_\lambda \left[ x^T y - x^T x \left( x^T y + T \right) \right] \\ & = S_\lambda \left[ - T \right] \\ & = 0, \end{align*} where $T = \begin{cases} - \lambda & \textrm{ if } x^T y > \lambda \\ \lambda & \textrm{ if } x^T y < -\lambda \\ 0 & \textrm{ otherwise} \end{cases}$. Then, since $\hat\beta_2^{(1)}= 0$, the second iteration of coordinate descent will repeat the computations above. Inductively, we see that $\hat\beta_j^{(i)} = \hat\beta_j^{(i)}$ for all iterations $i$ and $j \in \{1,2\}$. Therefore glmnet will report $\hat\beta_1 = \hat\beta_1^{(1)}$ and $\hat\beta_2 = \hat\beta_2^{(1)}$ since the stopping criterion is immediately reached.

Answer 2

When I re-run your code, I get that the coefficient of $x_2$ is numerically indistinguishable from zero.

To understand better why LASSO sets that coefficient to zero, you should look at the relationship between LASSO and Least Angle Regression (LAR). LASSO can be seen as a LAR with a special modification.

The algorithm of LAR is roughly like this: Start with an empty model (except for an intercept). Then add the predictor variable that is the most correlated with $y$, say $x_j$. Incrementally change that predictor's coefficient $\beta_j$, until the residual $y - c - x_j\beta_j$ is equally correlated with $x_j$ and another predictor variable $x_k$. Then change the coefficients of both $x_j$ and $x_k$ until a third predictor $x_l$ is equally correlated with the residual $y - c - x_j\beta_j -x_k\beta_k$ and so on.

LASSO can be seen as LAR with the following twist: as soon as the coefficient of a predictor in your model (an "active" predictor) hits zero, drop that predictor from the model. This is what happens when you regress $y$ on the collinear predictors: both will get added to the model at the same time and, as their coefficients are changed, their respective correlation with the residuals will change proportionately, but one of the predictors will get dropped from the active set first because it hits zero first. As for which of the two collinear predictors it will be, I don't know. [EDIT: When you reverse the order of $x_1$ and $x_2$, you can see that the coefficient of $x_1$ is set to zero. So the glmnet algorithm simply seems to set those coefficients to zero first that are ordered later in the design matrix.]

A source that explains these things more in detail is Chapter 3 in "The Elements of Statistical Learning" by Friedman, Hastie and Tibshirani.