4
$\begingroup$

I'm trying to understand why ridge regression coefficient estimates (through the glmnet package in R) do not match the ordinary least squares (OLS) estimates in the case that $\lambda$ = 0. I have seen a couple of other posts on this topic, but none of them answered my question.

Here is a minimal reprex:

library(glmnet)
set.seed(1)

X <- matrix(rnorm(90), ncol = 9, nrow = 10, byrow = TRUE)
y <- matrix(rnorm(10), nrow = 10, ncol = 1)
X_scaled <- scale(X)

ridge1 <- glmnet(X_scaled, y, alpha = 0, lambda = 0)
lm1 <- lm(y~X_scaled)

This results in:

> coef(lm1)
(Intercept)   X_scaled1   X_scaled2   X_scaled3   X_scaled4   X_scaled5   X_scaled6   X_scaled7   X_scaled8   X_scaled9 
  0.1123413   4.4105824  -4.1680260   4.9959933   2.2281174   3.0542372   3.8673192  -2.5323069   0.4444550   5.0073531
 
> coef(ridge1)
10 x 1 sparse Matrix of class "dgCMatrix"
                    s0
(Intercept)  0.1123413
V1           4.1667913
V2          -3.9353740
V3           4.7692778
V4           2.1239412
V5           2.8683159
V6           3.6622262
V7          -2.3987696
V8           0.4305574
V9           4.7282300

The coefficient estimates from ridge regression should match the OLS coefficients when $\lambda$=0, however, these do not match (except for the intercept). What is going on here?

Improve this question
$\endgroup$
4
  • $\begingroup$ There ia a warning in the manual for glmnet(): "Avoid supplying a single value for lambda .... Supply instead a decreasing sequence of lambda values. glmnet relies on its warms starts for speed, and its often faster to fit a whole path than compute a single fit." See what happens if you instead supply a sequence of $\lambda$ values that ends with 0. If that removes the problem, write that up as an answer here; it's OK to answer your own question. If not, I'm perplexed. Clearly all of the coefficients are slightly penalized. $\endgroup$
    – EdM
    Commented Aug 4, 2020 at 20:59
  • $\begingroup$ Thank you, @EdM. Tried that, but to no avail. I updated the ridge regression line to ridge1 <- glmnet(X_scaled, y, alpha = 0, lambda = c(50, 40, 30, 20, 10, 9, 8, 7, 6, 5, 4, 3, 2.5, 2, 1.5, 1, .9, .8, .7, .6, .4, .3, .2, .1, .05, .01, .001, .0001, 0)), yet coef(ridge1, s = 0) results in the same coefficient values as before (to the fourth decimal). $\endgroup$
    – bob
    Commented Aug 4, 2020 at 21:19
  • $\begingroup$ The ratio of the square of the sums of the non-intercept coefficients between ridge1 and lm1 is suspiciously close to 0.9: sum(coef(ridge1)[2:10])^2/sum(coef(lm1)[2:10])^2 = 0.8995228. I wonder if the function somehow imposes a minimum 10% penalty. One is always welcome to inspect the source code, but it might be simpler just to ask the package authors about that choice. $\endgroup$
    – EdM
    Commented Aug 4, 2020 at 23:51
  • $\begingroup$ Your original example has no degrees of freedom and lm1 no residuals: you would get closer results if you tried something like X <- matrix(rnorm(50), ncol = 5, nrow = 10, byrow = TRUE) $\endgroup$
    – Henry
    Commented Aug 5, 2020 at 7:38

1 Answer 1

Reset to default
5
$\begingroup$

glmnet finds an approximate solution using coordinate descent. You can get the solution closer by changing the thresh parameter which is the threshold for the algorithm to stop converging to the solution.

library(glmnet)
set.seed(1)

### data
X <- matrix(rnorm(90), ncol = 9, nrow = 10, byrow = TRUE)
y <- matrix(rnorm(10), nrow = 10, ncol = 1)
X_scaled <- scale(X)

### perform fitting
ridge1 <- glmnet(X_scaled, y, alpha = 0, lambda = 0, thresh = 10^-7)   # this is the default
ridge2 <- glmnet(X_scaled, y, alpha = 0, lambda = 0, thresh = 10^-14)
lm1 <- lm(y~X_scaled)

### output
out <- cbind(lm1$coefficients,
             coefficients(ridge1),
             coefficients(ridge2))
colnames(out) <- c("lm", "ridge1", "ridge2")
out

### the difference is in the duration of the coordinate descent algorithm
### it is 2466 passes versus 9238 passes
ridge1$npasses
ridge2$npasses

which gives:

> out
10 x 3 sparse Matrix of class "dgCMatrix"
                    lm     ridge1     ridge2
(Intercept)  0.1123413  0.1123413  0.1123413
V1           4.4105824  4.1667913  4.4105053
V2          -4.1680260 -3.9353740 -4.1679524
V3           4.9959933  4.7692778  4.9959216
V4           2.2281174  2.1239412  2.2280845
V5           3.0542372  2.8683159  3.0541784
V6           3.8673192  3.6622262  3.8672543
V7          -2.5323069 -2.3987696 -2.5322647
V8           0.4444550  0.4305574  0.4444506
V9           5.0073531  4.7282300  5.0072648

You can get the result closer by decreasing the parameter thresh even further. An exact result could be computed if you use alpha = 0 (Tikhonov regularization can be computed directly using a matrix equation).

I am not sure why glmnet doesn't do that direct computation, but it is a much more general function. Probably, when your only interest is plain ridge regression (and not the generalizations) then you might find some other packages that use the direct (and exact) computation. Although maybe there is some speed advantage when you need to compute a path of solutions for many different $\lambda$, which EdM notices in the comments has a speed advantage over separate calculations.

Improve this answer
$\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.