LASSO with two predictors

Question

I have a question regarding LASSO with two predictors, somewhat related to another one of mine posted here.

I am trying to illustrate equation (6) of the original paper by Tibshirani, JRSSB 1996, which says that the LASSO estimates $(\hat\beta_1,\hat\beta_2)$ in the case of two regressors will be

$$\hat\beta_1=[s/2+(\hat\beta^{(ols)}_1-\hat\beta^{(ols)}_2)/2]^+\qquad\qquad(1)$$

and

$$\hat\beta_2=[s/2-(\hat\beta^{(ols)}_1-\hat\beta^{(ols)}_2)/2]^+\qquad\qquad(2)$$

where $(x)^+$ means "positive part of $x$" (i.e. $(x)^+=\max(x,0)$) and $s$ is our "budget"

$$|\hat\beta_1|+|\hat\beta_2| \leq s$$

The OLS estimates $\hat\beta^{(ols)}_1$, $\hat\beta^{(ols)}_2$ need to be positive and

$$\hat\beta^{(ols)}_1+\hat\beta^{(ols)}_2 \geq s.$$

The code below attempts to numerically replicate this finding, with only very partial success.

More specifically, in the setting in which both LASSO coefficients are positive (which you can generate by commenting out the second time y is generated in line 15 of the code), the LASSO coefficients indeed sum to the budget $s$ of 0.5. In the more interesting case in which one of the coefficients is shrunk to zero, however (obtained here by using the second y sample), the other coefficient is larger than the allocated budgetof 0.5, obviously leading to a sum of the coefficients that exceeds the budget $s$.

Where's my mistake?

rm(list=ls())
library(glmnet)
library(mvtnorm)
set.seed(4) 

N = 50
K = 2
X = rmvnorm(N, mean=rep(0,K))
u = rnorm(N,sd=.1)  
X = scale(X)
vX = var(X)
X = sqrt(1/(N-1))*X%*%solve(chol(vX)) # even generates orthonormal X, which should not even be necessary according to Tibshirani 

y = .9*X[,1]+.6*X[,2]+u # yields a sample in which both LASSO coefficients are nonzero - works fine
y = .9*X[,1]+.1*X[,2]+u # yields a sample in which only one LASSO coefficient is  nonzero - but larger than budget!

ytilde = y-mean(y)

reg = lm(ytilde~X-1)

budget = .5           # how large the sum of the absolute LASSO coefficients is allowed to be (since they are not negative here, abs plays no role)
reg$coefficients      # must both be positive from JRSSB, which they are
sum(reg$coefficients) # must be bigger than budget from JRSSB, which it is

# what LASSO estimates should be according to Tibshirani (JRSSB 1996), eq (6) 
beta_1_LASSO = max(c(0,budget/2+(reg$coefficients[1]-reg$coefficients[2])/2))   # corresponds to equation (1) in the post
beta_2_LASSO = max(c(0,budget/2-(reg$coefficients[1]-reg$coefficients[2])/2))   # corresponds to equation (2) in the post
beta_1_LASSO+beta_2_LASSO # only sums to budget if both coefficients are nonzero

Answer 1

The two formulas

$$\hat\beta_1=[s/2+(\hat\beta^{(ols)}_1-\hat\beta^{(ols)}_2)/2]^+\qquad\qquad(1)$$

and

$$\hat\beta_2=[s/2-(\hat\beta^{(ols)}_1-\hat\beta^{(ols)}_2)/2]^+\qquad\qquad(2)$$

are inconsistent with the budget equation

$$\hat\beta_1+\hat\beta_2 = s$$

in the case where only one predictor is in the model. Note that I switched out your inequality with an equality because it is always better to use up your full budget to drive down some squared error. In the case where both predictors are in the model, this is consistent

$$\hat\beta_1 + \hat\beta_2 = s/2+(\hat\beta^{(ols)}_1-\hat\beta^{(ols)}_2)/2 + s/2-(\hat\beta^{(ols)}_1-\hat\beta^{(ols)}_2)/2 = s$$

but for only one predictor it is not.

In the one predictor case, you just observe that since you have to use up your entire budget, it must be entirely spent on the one predictor in the model, so your working estimates are

$$\hat\beta_1 = s $$

and

$$\hat\beta_2 = 0 $$

This happens until the continuity equation is satisfied

$$ s = s/2 + (\hat\beta^{(ols)}_1-\hat\beta^{(ols)}_2)/2 $$

after which the second parameter enters the model. There's a continuity equation for the second parameter here, but it is automatically satisfied because

$$ s/2 - (\hat\beta^{(ols)}_1-\hat\beta^{(ols)}_2)/2 = (\hat\beta^{(ols)}_1-\hat\beta^{(ols)}_2)/2 - (\hat\beta^{(ols)}_1-\hat\beta^{(ols)}_2)/2 = 0$$

So after this point, we can move on with your original set of equations and both predictors in the model.

It's been a while since I've revisited this paper, but it does seem to me that there are some errors of omission in the discussion of the two predictor case. I'll try to work through the whole discussion in detail and make sure I agree with myself.

I went through it and I do agree with myself. Here's the argument.

Tibshirani's three working equations in this section are:

$$ \hat\beta_j = (\hat\beta^o_j - \lambda)^+ \quad \text{for} \ j=1,2 $$

which I'll call the lambda equations and

$$ \hat\beta_1 + \hat\beta_2 = s$$

which I'll call the budget equation. He also makes the assumption that $\hat\beta^o_1 \geq \hat\beta^o_2 \geq 0$. We can analyze this case by case.

Case 1: $\hat\beta^o_2 \leq \hat\beta^o_1 \leq \lambda$

Here both lambda equations reduce to $\hat\beta_j = 0$. So $s = 0$ as well, and we have no budget.

Case 2: $\hat\beta^o_2 < \lambda \leq \hat\beta^o_1$

Here, one lambda equation is $\hat\beta_2 = 0$, so the beget equation reads $\hat\beta_1 = s$.

Case 3: $\lambda < \hat\beta^o_2 \leq \hat\beta^o_1$

Here both lambda equations are non-trivial, so they can be subtracted, canceling the $\lambda$s and resulting in

$$ \hat\beta_1 - \hat\beta_2 = \hat\beta^o_1 - \hat\beta^o_2 $$

This can be added to the budget equation to clear out $\hat\beta^o_2$ and get

$$ 2 \hat\beta_1 = s + \hat\beta^o_1 - \hat\beta^o_2 $$

Solving this gets equation (1), and then using the budget equation again gives equation (2).