I'm trying to learn some basic Machine Learning and some basic R. I have made a very naive implementation of $L_2$ regularization in R based on the formula:

$\hat w^{ridge} = (X^TX +\lambda I)^{-1} X^T y$

My code looks like this:

fitRidge <- function(X, y, lambda) {
     # Add intercept column to X:
  X <- cbind(1, X)
     # Calculate penalty matrix:
  lambda.diag <- lambda * diag(dim(X)[2])
     # Apply formula for Ridge Regression:
  return(solve(t(X) %*% X + lambda.diag) %*% t(X) %*% y)
}

Note that I'm not yet trying to find an optimal $\lambda$, I'm simply estimating $\hat w^{ridge}$ for a given $\lambda$. However, something seems off. When I enter $\lambda = 0$ I get the expected OLS result. I checked this by applying lm.ridge(lambda = 0) on the same dataset and it gives me the same coefficients. However, when I input any other penalty, like $\lambda=2$ or $\lambda=5$ my coefficients and the coefficients given by lm.ridge disagree wildly. I tried looking at the implementation of lm.ridge but I couldn't work out what it does (and therefore what it does differently).

Could anyone explain why there is a difference between my results and the results from lm.ridge? Am I doing something wrong in my code? I've tried playing around with scale() but couldn't find an answer there.

EDIT:

To see what happens, run the following:

library(car)
X.prestige <- as.matrix.data.frame(Prestige[,c(1,2,3,5)])
y.prestige <- Prestige[,4]

fitRidge(X.prestige, y.prestige, 0)
coef(lm.ridge(formula = prestige~education+income+women+census, data = Prestige, lambda = 0))
fitRidge(X.prestige, y.prestige, 2)
coef(lm.ridge(formula = prestige~education+income+women+census, data = Prestige, lambda = 2))

EDIT2:

Okay, so based on responses below, I've gotten a somewhat clearer understanding of the problem. I've also closely re-read the section about RR in TESL by Hastie, Tibshirani and Friedman, where I discovered that the intercept is often estimated simply as the mean of the response. It seems that many sources on RR online are overly vague. I actually suspect many writers have never implemented RR themselves and might not have realized some important things as many of them leave out 3 important facts:

1. Intercept is not penalized in the normal case, the formula above only applies to the other coefficients.
2. RR is not equivariant under scaling, i.e. different scales gives different results even for the same data.
3. Following from 1, how one actually estimates intercept.

I tried altering my function accordingly:

fitRidge <- function(X, Y, lambda) {
  # Standardize X and Y
  X <- scale(X)
  Y <- scale(Y)
  # Generate penalty matrix
  penalties <- lambda * diag(ncol(X))
  # Estimate intercept
  inter <- mean(Y)
  # Solve ridge system
  coeff <- solve(t(X) %*% X + penalties, t(X) %*% Y)
  # Create standardized weight vector
  wz <- c(inter, coeff )
  return(wz)
}

I still don't get results equivalent to lm.ridge though, but it might just be a question of translating the formula back into the original scales. However, I can't seem to work out how to do this. I thought it would just entail multiplying by the standard deviation of the response and adding the mean, as usual for standard scores, but either my function is still wrong or rescaling is more complex than I realize.

Any advice?