Massive difference between R's glmnet and Python's sklearn regarding Lasso regression

Question

I have a burning question. First, in Python:

import os
import time
import numpy as np
from sklearn.linear_model import LassoCV
from sklearn.datasets import make_regression


X, y = make_regression(1000, 5000, noise = 4, random_state = 123)
if not os.path.exists("tmp"): os.mkdir("tmp")


# Save on disk to be read in R session.
np.save("tmp/X.npy", X)
np.save("tmp/y.npy", y)


# Let first 70% rows be training examples.
Nrow = int(round(len(X) * 0.7))


tik = time.time()
reg = LassoCV(cv = 5).fit(X[:Nrow, :], y[:Nrow])
time.time() - tik 
# 6.59s


# Compute mse on the validation examples:
np.mean((reg.predict(X[Nrow:, :]) - y[Nrow:]) ** 2)
# 19.036

Then in R:

X = RcppCNPy::npyLoad("tmp/X.npy") # Load the data saved from Python session.
y = RcppCNPy::npyLoad("tmp/y.npy")


# First 70% rows are training data.
trainInd = 1:as.integer(round(nrow(X) * 0.7))


system.time({
reg = glmnet::cv.glmnet(
  X[trainInd, ], y[trainInd], type.measure = "mse", nfolds = 5)
})
# 1.047s


mean((predict(reg, X[-trainInd, ], s = "lambda.min") - y[-trainInd]) ^ 2)
# 32.01

Why does R glmnet result in an MSE 68% higher than that from Python sklearn? I anticipated glmnet would be faster, but a 6.6x speedup over sklearn seems to suggest glmnet skipped something important, which might contribute to its significantly worse error metric? I couldn't figure out where I did wrong..

Answer 1

As always, it helps to make plots. Here I plot the cv.glmnet predictions on the left and the LassoCV predictions on the right, against the true values in the hold out set.

We see immediately that both implementations do a great job: they explain most of the variability in $y$ and the mean squared error (MSE) is small compared to the total variability in $y$.

We can quantify the difference between the implementations in terms of the proportion of variance explained (PVE) = 1 - MSE / Var$(y_{\text{holdout}})$.

• PVE_cv.glmnet = .9989
• PVE_LassoCV = .9993

Note: PVE is not the same as the R-squared, which is computed on the training data. And it's not equivalent to comparing the Lasso to the "naive" model that predicts the average $y$ observed in the training data.

So there is a difference but it is small relative to the variability of the outcome variable. I wouldn't use this example to conclude that LassoCV is a better implementation than cv.glmnet in general.

PS. An increase of .0004 could be a meaningful improvement and worthwhile to put in practice. However, to conclude that a new model is better than the production model by such a small increment, we would want to compare the models on more than 300 examples.