Why do Lars and Glmnet give different solutions for the Lasso problem?

Question

I want to better understand the R packages Lars and Glmnet, which are used to solve the Lasso problem: $$min_{(\beta_0 \beta) \in R^{p+1}} \left[\frac{1}{2N}\sum_{i=1}^{N}(y_i-\beta_0-x_i^T\beta)^2 + \lambda||\beta ||_{l_{1}} \right]$$ (for $p$ Variables and $N$ samples, see www.stanford.edu/~hastie/Papers/glmnet.pdf on page 3)

Therefore, I applied them both on the same toy dataset. Unfortunately, the two methods do not give the same solutions for the same data input. Does anybody have an idea where the difference comes from?

I obtained the results as follows: After generating some data (8 samples, 12 features, Toeplitz design, everything centered), I computed the whole Lasso path using Lars. Then, I ran Glmnet using the sequence of lambdas computed by Lars (multiplied by 0.5) and hoped to obtain the same solution, but I did not.

[edit:Please ignore the following paragraph, the diffierent intercepts were caused by a bug in the code, which I fixed now, sorry for that. The rest of the question stays relevant, as fixing the bug did not solve the problem adressed.]
When I saw the (non-identical) outputs, I recognized that glmnet computed an intercept for each lambda, while lars did not. So I subtracted by hand the intercept of some glmnet solution from Y and recentered Y. Then I computed the Lasso path again using Lars, and then with glmnet using the new lambdas from Lars.

One can see that the solutions are similar. But how can I explain the differences? Please find my code below. There is a related question here: GLMNET or LARS for computing LASSO solutions? , but it does not contain the answer to my question.

### Load packages
library(lars)
library(glmnet)
library(MASS)

###Set parameters 
nbFeatures<-12
nbSamples<-8
nbRelevantIndices<-3
snr<-1
nbLambdas<-10
nlambda<-nbLambdas


### Create data, not really important 
sigma<-matrix(0,nbFeatures,nbFeatures)
for (i in (1:nbFeatures))
  for (j in (1:nbFeatures))
    sigma[i,j]<-0.99^(abs(i-j))
X<-mvrnorm(n = nbSamples, rep(0, nbFeatures), sigma, tol = 1e-6, empirical = FALSE)
relevantIndices <- sample(1:nbFeatures, nbRelevantIndices, replace=F)
X<-scale(X) #bug fixed, replaced 'scale(X)' by 'X<-scale(X)'
beta <- rep(0, times=nbFeatures)
beta[relevantIndices] <- runif(nbRelevantIndices,0,1)
epsilon<-matrix(rnorm(nbSamples),nbSamples,1)
simulatedSnr<-(norm(X %*% beta,type="F"))/(norm(epsilon,type="F"))
epsilon<-epsilon*(simulatedSnr/snr)
Y <- X %*% beta + epsilon;
Y<-scale(Y) #bug fixed, replaced 'scale(Y)' by 'Y<-scale(Y)'

### launch lars 
la <- lars(X,Y,intercept=TRUE, max.steps=1000, use.Gram=FALSE)
coLars <- as.matrix(coef(la,,mode="lambda"))
print(round(coLars,4))



### Launch glmnet with lambda=1/2*lambda_lars 
glm2 <- glmnet(X,Y,family="gaussian",lambda=0.5*la$lambda,thresh=1e-16)
coGlm2 <- as.matrix(t(coef(glm2,mode="lambda")))
print(round(coGlm2,4))

##########################
#edit: The rest of the program is not relevant any more, as the issue of the intercept was due to a bug.    
##########################

### Remove intercept computed by glmnet for some lambda (adjusted by hand, depending on output) 
Y<-Y--0.0863
Y<-Y/norm(Y,"F")

### Launch Lars 
la <- lars(X,Y,intercept=TRUE, max.steps=1000, use.Gram=FALSE)
coLars <- as.matrix(coef(la,,mode="lambda"))
print(round(coLars,4))

### Launch Glmnet with lambda=1/2*lambda_lars 
glm2 <- glmnet(X,Y,family="gaussian",lambda=0.5*la$lambda,thresh=1e-16)
coGlm2 <- as.matrix(t(coef(glm2,mode="lambda")))
print(round(coGlm2,4))

Answer 1

Finally we were able to produce the same solution with both methods! First issue is that glmnet solves the lasso problem as stated in the question, but lars has a slightly different normalization in the objective function, it replaces $\frac{1}{2N}$by $\frac{1}{2}$. Second, both methods normalize the data differently, so the normalization must be swiched off when calling the methods.

To reproduce that, and see that the same solutions for the lasso problem can be computed using lars and glmnet, the following lines in the code above must be changed:

la <- lars(X,Y,intercept=TRUE, max.steps=1000, use.Gram=FALSE)

to

la <- lars(X,Y,intercept=TRUE, normalize=FALSE, max.steps=1000, use.Gram=FALSE)

and

glm2 <- glmnet(X,Y,family="gaussian",lambda=0.5*la$lambda,thresh=1e-16)

to

glm2 <- glmnet(X,Y,family="gaussian",lambda=1/nbSamples*la$lambda,standardize=FALSE,thresh=1e-16)

Answer 2

Obviously if the methods use different models you will get different answers. Subtracting off the intercept terms does not lead to the model without the intercept because the best fitting coefficients will change and you do not change them the way you are approaching it. You need to fit the same model with both methods if you want the same or nearly the same answers.