Ridge Regression problem with output

I am trying to implement ridge regression in R, but the results are wrong.

$$\hat\beta^{Ridge} = argmin\sum_{i=1}^N(y-\beta_0 - \sum_{j=1}^{p} x_{i,j}\beta_j)^2 + \lambda\sum_{j=1}^{p}\beta_j^2$$

Hence,

$$RSS(\lambda) = (y - X\beta)^T(y-X\beta) + \lambda\beta^T\beta = |y|^2-2|X^Ty\beta|+|X\beta|^2 + \lambda|\beta|^2$$

$$dRSS(\lambda)/d\beta = -2X^Ty + 2X^TX\beta+2\lambda\beta$$, so $$2\beta(X^TX+\lambda I) = 2X^Ty => \hat\beta^{Ridge} = (X^TX+\lambda I)^{-1}X^Ty$$

Now I want to use this in a code to get $$\hat\beta$$'s for Ridge regression

 library(MASS)
y <- longley$GNP.deflator x <- longley[,2:ncol(longley)] X <- model.matrix(~as.matrix(x)) BetaRidge <- function(lam){ solve( crossprod(X) + diag(ncol(X))*lam) %*% crossprod(X,y) } lm.ridge(longley$GNP.deflator~.,
data=longley,lambda = 1)
BetaRidge(1)


The output is considerably different! What is wrong with the logic?

• there is at least one error in your code and calculations: model.matrix will include an intercept, in ridge regression however you don't penalize the intercept. second thing which comes to my mind is that in penalized regression methods one often estimates the coefficients using standardized regressors (implying $\beta_0 = 0$). – BloXX Oct 22 '18 at 10:02
• I had tried without intercept too, but in that way if I remember it right the difference between betas was 10 times bigger. – Sargis Iskandaryan Oct 22 '18 at 11:17
• I will try to standardize then use the model. – Sargis Iskandaryan Oct 22 '18 at 11:19

1.) you don't penalize the intercept in ridge regression and

2.) penalized least square methods like ridge regression or lasso are most meaningful applied to standardized regressors and a centered dependent variable (to get rid of the intercept). So, internally, lm.ridge will perform it's optimization task using $$X_{i,j}^*$$ and $$y_i^*$$, where $$y_i^* = y_i - \overline{y}$$ and $$x_{i,j}^* = (x_{i,j}-\overline{x}_j)/\widehat{\sigma}_j,$$ where (in the case of lm.ridge) $$\widehat{\sigma}_j^2 = \frac{1}{n}\sum_{i=1}^n (x_{i,j}-\overline{x}_j)^2.$$

Indeed: using standardized variables, your R-Code leads to identical results (besides some minor rounding errors, coming from different inversion methods):

 library(MASS)
BetaRidge <- function(lam){solve( crossprod(X) +diag(ncol(X))*lam) %*% crossprod(X,y)}
y <- longley\$GNP.deflator
X <- as.matrix(longley[,2:ncol(longley)])
X<-diag(rep(sqrt(16/15),16))%*%scale(X) #standardization such that \sum_{i}(x_{i,j}^*)^2 = n
y<-y-mean(y)


all(round(coef(lm.ridge(y~X,lambda = 1))[-1],12)== round(BetaRidge(1),12))

true

of course, after estimating the coefficients you need (could!) to transform them back on their original scale.