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 budget
of 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