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$
  • 3
    $\begingroup$ This will be hard to reproduce without code and data. $\endgroup$ – Sycorax May 18 '18 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 May 20 '18 at 19:10
  • $\begingroup$ @Chris why set intercepts to false? $\endgroup$ – shadowtalker May 21 '18 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$ – aginensky Jun 23 '18 at 12:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.