28
$\begingroup$

I'm using glmnet to calculate ridge regression estimates. I got some results that made me suspicious in that glmnet is really doing what I think it does. To check this I wrote a simple R script where I compare the result of ridge regression done by solve and the one in glmnet, the difference is significant:

n    <- 1000
p.   <-  100
X.   <- matrix(rnorm(n*p,0,1),n,p)
beta <- rnorm(p,0,1)
Y    <- X%*%beta+rnorm(n,0,0.5)

beta1 <- solve(t(X)%*%X+5*diag(p),t(X)%*%Y)
beta2 <- glmnet(X,Y, alpha=0, lambda=10, intercept=FALSE, standardize=FALSE, 
                family="gaussian")$beta@x
beta1-beta2

The norm of the difference is usually around 20 which cannot be due to numerically different algorithms, I must be doing something wrong. What are the settings I have to set in glmnet in order to obtain the same result as with ridge?

$\endgroup$
  • 1
    $\begingroup$ Have you seen this question? $\endgroup$ – cdeterman Dec 15 '14 at 16:00
  • 1
    $\begingroup$ Yes, but I still don't get the same result using the normalization. $\endgroup$ – John Dec 15 '14 at 16:46
  • $\begingroup$ Could you post your code then? $\endgroup$ – shadowtalker Dec 16 '14 at 14:30
  • $\begingroup$ I've just had the same problem! a = data.frame(a=jitter(1:10), b=jitter(1:10), c=jitter(1:10), d=jitter(1:10), e=jitter(1:10), f=jitter(1:10), g=sample(jitter(1:10)), y=seq(10,100,10)); coef(lm.ridge(y~a+b+c+d+e+f+g, a, lambda=2.57)); coef(glmnet(as.matrix(a[,1:7]),a$y,family = "gaussian", alpha=0, lambda=2.57/10)) The results differ quite a bit and become much more similar when I use much higher lambdas for glmnet. $\endgroup$ – a11msp Oct 5 '16 at 21:13
  • $\begingroup$ Intriguing. The coefficients seem to differ roughly by the factor of 10. $\endgroup$ – tomka Dec 14 '16 at 0:06
27
+50
$\begingroup$

The difference you are observing is due to the additional division by the number of observations, N, that GLMNET uses in their objective function and implicit standardization of Y by its sample standard deviation as shown below.

$$ \frac{1}{2N}\left\|\frac{y}{s_y}-X\beta\right\|^2_{2}+\lambda\|\beta\|^2_{2}/2 $$

where we use $1/n$ in place of $1/(n-1)$ for $s_y$, $$ s_y=\frac{\sum_i(y_i-\bar{y})^2}{n} $$

By differentiating with respect to beta, setting the equation to zero,

$$ X^TX\beta-\frac{X^Ty}{s_y}+N\lambda\beta =0 $$

And solving for beta, we obtain the estimate,

$$ \tilde{\beta}_{GLMNET}= (X^TX+N\lambda I_p)^{-1}\frac{X^Ty}{s_y} $$

To recover the estimates (and their corresponding penalties) on the original metric of Y, GLMNET multiplies both the estimates and the lambdas by $s_y$ and returns these results to the user,

$$ \hat{\beta}_{GLMNET}=s_y\tilde{\beta}_{GLMNET}= (X^TX+N\lambda I_p)^{-1}X^Ty $$ $$ \lambda_{unstd.}=s_y\lambda $$

Compare this solution with the standard derivation of ridge regression.

$$ \hat{\beta}= (X^TX+\lambda I_p)^{-1}X^Ty $$

Notice that $\lambda$ is scaled by an extra factor of N. Additionally, when we use the predict() or coef() function, the penalty is going to be implicitly scaled by $1/s_y$. That is to say, when we use these functions to obtain the coefficient estimates for some $\lambda^*$, we are effectively obtaining estimates for$\lambda=\lambda^*/s_y$.

Based on these observations, the penalty used in GLMNET needs to be scaled by a factor of $s_y/N$.

set.seed(123)

n    <- 1000
p   <-  100
X   <- matrix(rnorm(n*p,0,1),n,p)
beta <- rnorm(p,0,1)
Y    <- X%*%beta+rnorm(n,0,0.5)

sd_y <- sqrt(var(Y)*(n-1)/n)[1,1]

beta1 <- solve(t(X)%*%X+10*diag(p),t(X)%*%(Y))[,1]

fit_glmnet <- glmnet(X,Y, alpha=0, standardize = F, intercept = FALSE, thresh = 1e-20)
beta2 <- as.vector(coef(fit_glmnet, s = sd_y*10/n, exact = TRUE))[-1]
cbind(beta1[1:10], beta2[1:10])

           [,1]        [,2]
[1,]  0.23793862  0.23793862
[2,]  1.81859695  1.81859695
[3,] -0.06000195 -0.06000195
[4,] -0.04958695 -0.04958695
[5,]  0.41870613  0.41870613
[6,]  1.30244151  1.30244151
[7,]  0.06566168  0.06566168
[8,]  0.44634038  0.44634038
[9,]  0.86477108  0.86477108
[10,] -2.47535340 -2.47535340

The results generalize to the inclusion of an intercept and standardized X variables. We modify a standardized X matrix to include a column of ones and the diagonal matrix to have an additional zero entry in the [1,1] position (i.e. do not penalize the intercept). You can then unstandardize the estimates by their respective sample standard deviations (again ensure you are using 1/n when computing standard deviation).

$$ \hat\beta_{j}=\frac{\tilde{\beta_j}}{s_{x_j}} $$

$$ \hat\beta_{0}=\tilde{\beta_0}-\bar{x}^T\hat{\beta} $$

mean_x <- colMeans(X)
sd_x <- sqrt(apply(X,2,var)*(n-1)/n)
X_scaled <- matrix(NA, nrow = n, ncol = p)
for(i in 1:p){
    X_scaled[,i] <- (X[,i] - mean_x[i])/sd_x[i] 
}
X_scaled_ones <- cbind(rep(1,n), X_scaled)

beta3 <- solve(t(X_scaled_ones)%*%X_scaled_ones+1000*diag(x = c(0, rep(1,p))),t(X_scaled_ones)%*%(Y))[,1]
beta3 <- c(beta3[1] - crossprod(mean_x,beta3[-1]/sd_x), beta3[-1]/sd_x)

fit_glmnet2 <- glmnet(X,Y, alpha=0, thresh = 1e-20)
beta4 <- as.vector(coef(fit_glmnet2, s = sd_y*1000/n, exact = TRUE))

cbind(beta3[1:10], beta4[1:10])
             [,1]        [,2]
 [1,]  0.24534485  0.24534485
 [2,]  0.17661130  0.17661130
 [3,]  0.86993230  0.86993230
 [4,] -0.12449217 -0.12449217
 [5,] -0.06410361 -0.06410361
 [6,]  0.17568987  0.17568987
 [7,]  0.59773230  0.59773230
 [8,]  0.06594704  0.06594704
 [9,]  0.22860655  0.22860655
[10,]  0.33254206  0.33254206

Added code to show standardized X with no intercept:

set.seed(123)

n <- 1000
p <-  100
X <- matrix(rnorm(n*p,0,1),n,p)
beta <- rnorm(p,0,1)
Y <- X%*%beta+rnorm(n,0,0.5)

sd_y <- sqrt(var(Y)*(n-1)/n)[1,1]

mean_x <- colMeans(X)
sd_x <- sqrt(apply(X,2,var)*(n-1)/n)

X_scaled <- matrix(NA, nrow = n, ncol = p)
for(i in 1:p){
    X_scaled[,i] <- (X[,i] - mean_x[i])/sd_x[i] 
}

beta1 <- solve(t(X_scaled)%*%X_scaled+10*diag(p),t(X_scaled)%*%(Y))[,1]

fit_glmnet <- glmnet(X_scaled,Y, alpha=0, standardize = F, intercept = 
FALSE, thresh = 1e-20)
beta2 <- as.vector(coef(fit_glmnet, s = sd_y*10/n, exact = TRUE))[-1]
cbind(beta1[1:10], beta2[1:10])

             [,1]        [,2]
 [1,]  0.23560948  0.23560948
 [2,]  1.83469846  1.83469846
 [3,] -0.05827086 -0.05827086
 [4,] -0.04927314 -0.04927314
 [5,]  0.41871870  0.41871870
 [6,]  1.28969361  1.28969361
 [7,]  0.06552927  0.06552927
 [8,]  0.44576008  0.44576008
 [9,]  0.90156795  0.90156795
[10,] -2.43163420 -2.43163420
$\endgroup$
  • 3
    $\begingroup$ +6. Welcome to CV and thanks for answering this old question in such a clear way. $\endgroup$ – amoeba says Reinstate Monica Feb 8 '17 at 9:55
  • 1
    $\begingroup$ It should be the identity matrix instead of $\beta$ in the solution of $\tilde{\beta}$, correct? $\endgroup$ – user1769197 Jul 4 '17 at 16:35
  • $\begingroup$ I also notice that for the second part where you said "The results generalize to the inclusion of an intercept and standardized X variables"; for this part, if you exclude the intercept, then following the same calculations, the results of glmnet become different from manual computation. $\endgroup$ – user1769197 Jul 4 '17 at 17:38
  • $\begingroup$ Correct, I have updated the solution with the identity matrix in place of $\beta$ as needed. I checked the solution for standardized X with no intercept and still obtain identical results (see additional code above). $\endgroup$ – skijunkie Jul 14 '17 at 20:25
3
$\begingroup$

According to https://web.stanford.edu/~hastie/glmnet/glmnet_alpha.html, when the family is gaussian, glmnet() should minimize $$\frac{1}{2n} \sum_{i=1}^n (y_i-\beta_0-x_i^T\beta)^2 +\lambda\sum_{j=1}^p(\alpha|\beta_j| +(1-\alpha)\beta_j^2/2). \tag{1}$$

When using glmnet(x, y, alpha=1) to fit the lasso with the columns in $x$ standardized, the solution for the reported penalty $\lambda$ is the solution for minimizing $$\frac{1}{2n} \sum_{i=1}^n (y_i-\beta_0-x_i^T\beta)^2 +\lambda \sum_{j=1}^p |\beta_j|.$$ However, at least in glmnet_2.0-13, when using glmnet(x, y, alpha=0) to fit ridge regression, the solution for a reported penalty $\lambda$ is the solution for minimizing $$\frac{1}{2n} \sum_{i=1}^n (y_i-\beta_0-x_i^T\beta)^2 +\lambda \frac{1}{2s_y} \sum_{j=1}^p \beta_j^2.$$ where $s_y$ is the standard deviation of $y$. Here, the penalty should have been reported as $\lambda/s_y$.

What might happen is that the function first standardizes $y$ to $y_0$ and then minimizes $$\frac{1}{2n} \sum_{i=1}^n (y_{0i}-x_i^T\gamma)^2 +\eta \sum_{j=1}^p(\alpha|\gamma_j| +(1-\alpha)\gamma_j^2/2), \tag{2}$$ which effectively is to minimize $$\frac{1}{2n s_y^2} \sum_{i=1}^n (y_i-\beta_0-x_i^T\beta)^2 +\eta \frac{\alpha}{s_y} \sum_{j=1}^p |\beta_j| +\eta \frac{1-\alpha}{2s_y^2} \sum_{j=1}^p \beta_j^2,$$ or equivalently, to minimize $$\frac{1}{2n} \sum_{i=1}^n (y_i-\beta_0-x_i^T\beta)^2 +\eta s_y \alpha \sum_{j=1}^p |\beta_j| +\eta (1-\alpha) \sum_{j=1}^p \beta_j^2/2.$$

For the lasso ($\alpha=1$), scaling $\eta$ back to report the penalty as $\eta s_y$ makes sense. Then for all $\alpha$, $\eta s_y$ has to be reported as the penalty to maintain continuity of the results across $\alpha$. This probably is the cause of the problem above. This is partly due to using (2) to solve (1). Only when $\alpha=0$ or $\alpha=1$ there is some equivalence between problems (1) and (2) (i.e., a correspondence between the $\lambda$ in (1) and the $\eta$ in (2)). For any other $\alpha\in(0,1)$, problems (1) and (2) are two different optimization problems, and there is no one-to-one correspondence between the $\lambda$ in (1) and the $\eta$ in (2).

$\endgroup$
  • 1
    $\begingroup$ I can't see where does your answer differ from the previous one. Could you explain, please? $\endgroup$ – Firebug Feb 23 '18 at 13:50
  • 1
    $\begingroup$ @Firebug I wanted to shed light on why the function reports the lambda this way, which appears unnatural when viewed solely from the perspective of ridge regression, but makes sense (or has to be this way) when viewed from the perspective of the whole spectrum including both ridge and the lasso. $\endgroup$ – Chun Li Feb 24 '18 at 18:18

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.