11
$\begingroup$

I am going through the LAB section §6.6 on Ridge Regression/Lasso in the book 'An Introduction to Statistical Learning with Applications in R' by James, Witten, Hastie, Tibshirani (2013).

More specifically, I am trying to do apply the scikit-learn Ridge model to the 'Hitters' dataset from the R package 'ISLR'. I have created the same set of features as shown in the R code. However, I cannot get close to the results from the glmnet() model. I have selected one L2 tuning parameter to compare. ('alpha' argument in scikit-learn).

Python:

regr = Ridge(alpha=11498)
regr.fit(X, y)

http://nbviewer.ipython.org/github/JWarmenhoven/ISL-python/blob/master/Notebooks/Chapter%206.ipynb

R:

Note that the argument alpha=0 in glmnet() means that a L2 penalty should be applied (Ridge regression). The documentation warns not to enter a single value for lambda, but the result is the same as in ISL, where a vector is used.

ridge.mod <- glmnet(x,y,alpha=0,lambda=11498)

What causes the differences?

Edit:
When using penalized() from the penalized package in R, the coefficients are the same as with scikit-learn.

ridge.mod2 <- penalized(y,x,lambda2=11498)

Maybe the question could then also be: 'What is the difference between glmnet() and penalized() when doing Ridge regression?

New python wrapper for actual Fortran code used in R package glmnet
https://github.com/civisanalytics/python-glmnet

$\endgroup$

migrated from stackoverflow.com Jul 6 '15 at 9:15

This question came from our site for professional and enthusiast programmers.

  • 5
    $\begingroup$ Totally unfamiliar with glmnet ridge regression. But by default, sklearn.linear_model.Ridge does unpenalized intercept estimation (standard) and the penalty is such that ||Xb - y - intercept||^2 + alpha ||b||^2 is minimized for b. There can be factors 1/2 or 1/n_samples or both in front of the penalty, making results different immediately. To factor out the penalty scaling problem, set the penalty to 0 in both cases, resolve any discrepancies there and then check what adding back the penalty does. And btw IMHO here IS the right place to ask this question. $\endgroup$ – eickenberg Jul 4 '15 at 12:51
9
$\begingroup$

My answer is missing a factor of $\frac{1}{N}$, please see @visitors answer below for the correct comparison.


Here are two references that should clarify the relationship.

The sklearn documentation says that linear_model.Ridge optimizes the following objective function

$$ \left| X \beta - y \right|_2^2 + \alpha \left| \beta \right|_2^2 $$

The glmnet paper says that the elastic net optimizes the following objective function

$$ \left| X \beta - y \right|_2^2 + \lambda \left( \frac{1}{2} (1 - \alpha) \left| \beta \right|_2^2 + \alpha \left| \beta \right|_1 \right) $$

Notice that the two implementations use $\alpha$ in totally different ways, sklearn uses $\alpha$ for the overall level of regularization while glmnet uses $\lambda$ for that purpose, reserving $\alpha$ for trading between ridge and lasso regularization.

Comparing the formulas, it look like setting $\alpha = 0$ and $\lambda = 2 \alpha_{\text{sklearn}}$ in glmnet should recover the solution from linear_model.Ridge.

$\endgroup$
  • $\begingroup$ And I totally missed that in @eickenberg 's comment as well. I do have to use standardize = FALSE in glmnet() to get the identical results. $\endgroup$ – Jordi Jul 7 '15 at 8:30
  • $\begingroup$ @Jordi You should definitely standardized if using linear_model.Ridge for any real world analysis. $\endgroup$ – Matthew Drury Jul 7 '15 at 11:56
  • $\begingroup$ I understand that sklearn linear_model.Ridge model standardizes the features automatically. Normalization is optional. I wonder why I then need to deactivate standardization in glmnet() to get the models to produce identical results. $\endgroup$ – Jordi Jul 7 '15 at 13:26
10
$\begingroup$

Matthew Drury's answer should have a factor of 1/N. More precisely...

The glmnet documentation states that the elastic net minimizes the loss function

$$ \frac{1}{N} \| X\beta - y \|_2^2 + \lambda \left( \frac{1}{2} (1 - \alpha) \, \| \beta \|_2^2 + \alpha \| \beta \|_1 \right) $$

The sklearn documentation says that linear_model.Ridge minimizes the loss function

$$ \| X\beta - y \|_2^2 + \alpha \| \beta \|_2^2 $$

which is equivalent to minimizing

$$ \frac{1}{N} \| X\beta - y \|_2^2 + \frac{\alpha}{N} \| \beta \|_2^2 $$

To obtain the same solution from glmnet and sklearn, both of their loss functions must be equal. This means setting $\alpha = 0$ and $\displaystyle{\lambda = \frac{2}{N} \alpha_{\text{sklearn}}}$ in glmnet.

library(glmnet)
X = matrix(c(1, 1, 2, 3, 4, 2, 6, 5, 2, 5, 5, 3), byrow = TRUE, ncol = 3)
y = c(1, 0, 0, 1)
reg = glmnet(X, y, alpha = 0, lambda = 2 / nrow(X))
coef(reg)

glmnet output: –0.03862100, –0.03997036, –0.07276511, 0.42727955

import numpy as np
from sklearn.linear_model import Ridge
X = np.array([[1, 1, 2], [3, 4, 2], [6, 5, 2], [5, 5, 3]])
y = np.array([1, 0, 0, 1])
reg = Ridge(alpha = 1, fit_intercept = True, normalize = True)
reg.fit(X, y)
np.hstack((reg.intercept_, reg.coef_))

sklearn output: –0.03862178, –0.0399697, –0.07276535, 0.42727921

$\endgroup$
  • 4
    $\begingroup$ The different definitions of parameters and their scaling used in different libraries is a common source of confusion. $\endgroup$ – AaronDefazio Mar 30 '17 at 3:31
  • 1
    $\begingroup$ I wouldn't expect that both Gung and I would get this wrong. $\endgroup$ – Michael Chernick Mar 30 '17 at 6:10
  • 2
    $\begingroup$ Yes, both of you got it wrong. Your reasons for rejecting my edit make it clear that both of you didn't see my comment "Missing factor of 1/N" at stats.stackexchange.com/review/suggested-edits/139985 $\endgroup$ – visitor Mar 30 '17 at 8:04
  • $\begingroup$ Your edit was probably rejected because it changed much more than only what you claim. If youd like to edit my post and only change the missing factor, please do, but changing my links and wording and code as well is overkill. The comments about your unfair treatment in your answer are innapropriate, and unrelated to the content of the question, please remove them. Your wording also plagarized my answer, this is not the right way to respond to a rejected edit. We would love your valuable contributions to our community, but please aquaint yourself with our norms before eviscerating us. $\endgroup$ – Matthew Drury Mar 31 '17 at 16:53
  • 1
    $\begingroup$ @visitor Sorry if I came off a bit gruff. I really should just be trying to communicate that you seem like a good potential contributor to the site, and I want you to have a good experience. We have some social norms, just like any other group, and you will have a better experience if you stay aware of them. I still think "Matthew Drury's answer is wrong" is quite harsh, there are surely better ways to communicate that my answer is erroneously missing a factor of $\frac{1}{N}$. "X's answer is wrong" reads as a personal attack. $\endgroup$ – Matthew Drury Mar 31 '17 at 20:41

Your Answer

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

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