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?

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)

#3 x 1 sparse Matrix of class "dgCMatrix"
#                       1
#(Intercept) 1.057426e+02
#x1          9.680073e+01
#x2          3.122502e-15
  • 2
    $\begingroup$ I'm not sure if this is a good simulation because both coefficients are in fact zero. It's a bit more interesting to look at the behavior of the coefficient estimates when there's a real relationship. $\endgroup$ – dsaxton Oct 20 '16 at 19:57
  • 1
    $\begingroup$ Simulation improved. I provide the simulation because I want to explain what my question is. I only interested in theoretical results of this question. $\endgroup$ – Jill Clover Oct 20 '16 at 20:39
  • 1
    $\begingroup$ I think the behavior will be unpredictable because the model is not identifiable. That is, how can the model fitting procedure possibly know for instance that $\beta_1 = 100$ and $\beta_2 = 0$ rather than $\beta_1 = 0$ and $\beta_2 = 50$? It can't, because either is "correct." $\endgroup$ – dsaxton Oct 20 '16 at 20:51
  • $\begingroup$ I agree with your reasoning. Is there a mathematical way to describe it? $\endgroup$ – Jill Clover Oct 20 '16 at 20:54
  • 1
    $\begingroup$ I think you meant y = 100*x1 + 100 + runif(100), otherwise you get a single random number that is recycled and added uniformly to all other entries. $\endgroup$ – Firebug Sep 11 '17 at 16:53

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.

| cite | improve this answer | |
  • 2
    $\begingroup$ glmnet has feature scaling turned on by default, I'm pretty sure. So $x_1$ and $x_2$ become the same in the model. $\endgroup$ – Firebug Sep 11 '17 at 16:22
  • 2
    $\begingroup$ Try this instead: ridge.mod=cv.glmnet(x_train,y,alpha=1, standardize = FALSE); coef(ridge.mod) $\endgroup$ – Firebug Sep 11 '17 at 16:25
  • 2
    $\begingroup$ That did it! Great thinking, @Firebug! Now the coefficient of $x_1$ does indeed become estimated as zero. Thank you for sharing your insight! $\endgroup$ – user795305 Sep 11 '17 at 16:31

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.

| cite | improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.