How to find smallest $\lambda$ such all lasso coefficients are set to 0, depending on the intercept?

Question

While bouilding a LASSO-penalized model it is well known that $\lambda =\left\lVert X^ty\right\lVert_\infty$ is the minimum value for which all the $\beta$ coefficients of the model are 0.

Consider for instance the BostonHousing dataset in R

library(mlbench, glmnet)
data(BostonHousing)
x <- data.matrix(BostonHousing[,-14]); y <- BostonHousing[,14]
lm = glmnet(x, y, nlambda=10, alpha=1, standardize=FALSE, intercept=FALSE)
lm$lambda[1]
> 8473.908

We can compute this by our own as

(1/length(y))*max(abs(t(x)%*%y))
> 8473.908

However, I dont know how is glmnet computing the lambda value in the case in which intercept is set to true:

lm2 = glmnet(x, y, nlambda=10, alpha=1, standardize=FALSE, intercept=TRUE)
lm2$lambda[1]
> 724.80000

Answer 1

You can obtain the same value

• when you first center the X and Y.

  > x <- x - matrix(rep(1,506)) %*% colMeans(x)
  > y <- y - mean(y)
  > (1/length(y))*max(abs(t(x)%*%y))
  [1] 724.8204
  

  note that I get a slightly different value than 724.80000 but this is the same when I run glmnet.

• when you manually include an intercept with penalty zero (or at least that is what I had expected, but it is a bit tricky and you need a hack)

  x <- cbind(x,rep(1,506))
  colnames(x)[14] <- 'intercept'
  lm = glmnet(x, y, alpha=1,
      standardize=FALSE, intercept=FALSE,
      penalty.factor = c(rep(1,13),0))
  lm$lambda[1]
  lm[]
  

  , but glmnet does not seem to include this added intercept. Also another strange thing is that when intercept=FALSE only 5 values/steps for beta and lambda are returned (you can manipulate this a bit by manually telling which ranges and numbers of lambda to calculate).

  So somehow this intercept=FALSE makes elnet ignore the intercept also when a vector 1 is inside the matrix x , but to find out how this works you would have to dig into the FORTRAN code, since that is where the work is done and not in the R code.

  You can hack this a bit by not adding a 'pure' intercept term. Such that it will not be ignored:

  x <- cbind(x,rep(1,506))
  x[1,14] <- 1.0000000001     #so we change slighly one value such that the 14-th column won't be recognized as an intercept term
  colnames(x)[14] <- 'intercept'
  lm = glmnet(x, y, alpha=1,
          standardize=FALSE, intercept=FALSE,
          penalty.factor = c(rep(1,13),0))
  lm$lambda[1]*14/13  # note the factor 14/13 which is to correct for a rescaling of the penalty factor based on the number of variables
  

  output is 724.8204