What is the difference between Lasso regression in glmnet (in R) and Sklearn lasso (in Python)?

Question

A similar post was discussed here regarding Ridge Regression: What are the differences between Ridge regression using R's glmnet and Python's scikit-learn?

My question is what is this difference for Lasso? In R my data gives a corresponding lambda value of $0.02$ from the cross-validation plot. This is also almost the same value that is returned from Python as well. However, the coefficients returned from R return several variable coefficients while Python returns $0$'s for the coefficients.

The expression which Python lasso optimizes is:

$\frac{1}{ 2\cdot n_{samples}} ||y - Xw||^2_2 + \lambda ||w||_1$.

While R minimizes:

$\frac{1}{ n_{samples}} ||y - Xw||^2_2 + \lambda ||w||_1$.

At first I thought it would be a simple factor of $2$ between the two $\lambda$'s, however, this does not seem to be the case.

I am new to coding in general, so any ideas would be appreciated.

Ok, so as an edit I am adding code with just some simple toy data. Both R and Python still give different values for the Lambda chosen. Although if you notice the coefficient plot it looks to be shifted by some factor... this is also the case in my full dataset.

Code for R:

#initialize libraries
library(xlsx)
library(glmnet)
library(plotmo)

#Data Input
x<-matrix(c(-1,-4,0,0,1,16),nrow=3,ncol=2,byrow=TRUE)
y<-matrix(c(-2.2,0,3.8),nrow=3,ncol=1,byrow=TRUE)

#lasso fit and plots
yfit=glmnet(x,y,family="gaussian",standardize.response=FALSE,intercept = FALSE,lambda=(exp(seq(log(0.0001), log(100), length.out=100))),upper.limits = 1000,lower.limits = -1000)
ycvfit=cv.glmnet(x,y,family="gaussian",nfolds=3,standardize.response=FALSE,intercept = FALSE,grouped = FALSE,lambda=(exp(seq(log(0.0001), log(100), length.out=100))),upper.limits = 1000,lower.limits = -1000)
plot.cv.glmnet(ycvfit)
plot.glmnet(yfit,xvar="lambda")


#writing results 
lassocoeff=as.matrix(coef(yfit,s=0.1))
print(lassocoeff)

And the Python Code:

import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import Lasso, LassoCV
from __future__ import division
import time


x=np.array([[-1,-4],[0,0],[1,16]])
y=np.array([[-2.2],[0],[3.8]]).reshape(3,)

#set alphas to be tested, note this is in logspace
alphas = np.logspace(-7,5,num=1000,base=np.e)

#plot coefficent size vs log(alpha) graph. 
lasso = Lasso(max_iter=10000, normalize=False, fit_intercept=False)
coefs = []
for a in alphas:
    lasso.set_params(alpha=a)
    lasso.fit((x), (y))
    coefs.append(lasso.coef_)
ax = plt.gca()
ax.plot(np.log(alphas), coefs)
ax.set_xscale('linear')
plt.axis('tight')
plt.xlabel('log(alpha)')
plt.ylabel('weights')


print("Computing regularization path using the coordinate descent lasso...")
t1 = time.time()
model = LassoCV(alphas=None,cv=3, normalize=False, fit_intercept=False).fit((x), (y))
t_lasso_cv = time.time() - t1
# Display results
m_log_alphas = np.log(model.alphas_)
plt.figure()
#To properly show, the ymin and ymax may need adjusted
ymin, ymax = 0, 25
plt.plot(m_log_alphas, model.mse_path_.mean(axis=-1), 'k',
         label='Average across the folds', linewidth=2)
plt.axvline(np.log(model.alpha_), linestyle='--', color='k',
            label='alpha: CV estimate')
plt.legend()
plt.xlabel('log(alpha)')
plt.ylabel('Mean square error')
plt.title('Mean square error on each fold: coordinate descent '
          '(train time: %.2fs)' % t_lasso_cv)
plt.axis('tight')
plt.ylim(ymin, ymax)
plt.xlim(-10,10)

#Calculate Lasso coefficents from proper alpha. Skipping for now, as we are just comparing 1 alpha value.
#lassocv = LassoCV(eps=.000000001,max_n_alphas=1000, cv=N, max_iter=100000, normalize=False, fit_intercept=False)
#lassocv.fit(x, y)
#minimumalpha = lassocv.alpha_
lasso=Lasso(alpha=.1,normalize=False,fit_intercept=False,max_iter=100000)
lasso.fit(x, y)
print(lasso.coef_)

Answer 1

Instead of using standardize.response=FALSE you should be doing standardize=FALSE.

Python

import numpy as np
import sklearn.linear_model
X=np.array([[-1,-4],[0,0],[1,16]])
y=np.array([[-2.2],[0],[3.8]]).reshape(3,)
lasso=sklearn.linear_model.Lasso(alpha=0.1, fit_intercept=False)
lasso.fit(X, y)
lasso.coef_

R

X<-matrix(c(-1,-4,0,0,1,16),nrow=3,ncol=2,byrow=TRUE)
y<-matrix(c(-2.2,0,3.8),nrow=3,ncol=1,byrow=TRUE)
fit = glmnet::glmnet(X, y, standardize=FALSE, intercept=FALSE, lambda=0.1)
fit$beta

The results are: $$ \begin{aligned} \beta_\mathrm{py} &= \begin{pmatrix} 1.14197818 \\ 0.17081043 \end{pmatrix} \\ \beta_\mathrm{R} &= \begin{pmatrix} 1.1436378 \\ 0.1706884 \end{pmatrix} \end{aligned} $$

I calculate $X'(X\beta - y)$ which should be $\pm \lambda n$ (see Ex 3.27 in the Elements of Statistical Learning 2nd edition).

$$ \begin{aligned} \nabla_\mathrm{py} &= \begin{pmatrix} -0.29983508 \\ -0.3 \end{pmatrix} \\ \nabla_\mathrm{R} &= \begin{pmatrix} -0.2989565 \\ -0.3 \end{pmatrix} \end{aligned} $$

So Python seems to be a bit more accurate on this "dataset". You can tweak the convergence settings for both algorithms if you like.