5
$\begingroup$

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_)
$\endgroup$
4
  • 3
    $\begingroup$ This will be hard to reproduce without code and data. $\endgroup$
    – Sycorax
    Commented May 18, 2018 at 20:48
  • $\begingroup$ 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) $\endgroup$
    – Chris
    Commented May 20, 2018 at 19:10
  • $\begingroup$ @Chris why set intercepts to false? $\endgroup$ Commented May 21, 2018 at 4:16
  • 1
    $\begingroup$ 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 ? $\endgroup$
    – meh
    Commented Jun 23, 2018 at 12:59

1 Answer 1

1
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.