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

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],,[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_)

• This will be hard to reproduce without code and data.
– Sycorax
May 18, 2018 at 20:48
• I updated the post with code. I am newer to this so the Python code is particularly messy. The toy data chosen was super simple. To note I want intercepts set to False. (In my own data I normalized the data prior to running in R and Python so that I have standardization/Normalization also set to False) May 20, 2018 at 19:10
• @Chris why set intercepts to false? May 21, 2018 at 4:16
• Let me address your statement "At first I thought it would be a simple factor of 2 between the two λ's, however, this does not seem to be the case." . You are missing parenthesis. The correct expression is $\frac{1}{ 2\cdot n_{samples}} (||y - Xw||^2_2 + \lambda ||w||_1)$. The normalizing constant out front is simply to make the result after differentiation look nice. That is why the coefficients should be the same. Python has a funny habit of using $C = \frac{1}{\lambda}$ for regularizing constants. Any chance that explains it ?
– meh
Jun 23, 2018 at 12:59

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],,[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.