Ridge

ridge <- function(w, X, y, lambda = .1) {
  # X: model matrix; 
  # y: target; 
  # lambda: penalty parameter; 
  # w: the weights/coefficients
  
  crossprod(y - X %*% w) + lambda * length(y) * crossprod(w)
}


set.seed(8675309)
N = 500
p = 10
X = scale(matrix(rnorm(N * p), ncol = p))
b = c(.5, -.5, .25, -.25, .125, -.125, rep(0, 4))
y = scale(X %*% b + rnorm(N, sd = .5))

#' Note, if `lambda=0`, result is the same as  `lm.fit`.
#' 
#' 
result_ridge = optim(
  rep(0, ncol(X)),
  ridge,
  X = X,
  y = y,
  lambda = .1,
  method = 'BFGS'
)

#' Analytical result.
#' 
result_ridge2 =  solve(crossprod(X) + diag(length(y)*.1, ncol(X))) %*% crossprod(X, y)

#' Alternative with augmented data (note sigma ignored as it equals 1, but otherwise
#' X/sigma and y/sigma).
#' 
X2 = rbind(X, diag(sqrt(length(y)*.1), ncol(X)))
y2 = c(y, rep(0, ncol(X)))
result_ridge3 = solve(crossprod(X2)) %*% crossprod(X2, y2)




#' `glmnet` is by default a mixture of ridge and lasso penalties, setting alpha
#' = 1 reduces to lasso, while alpha=0 would be ridge.


library(glmnet)

Le chargement a nécessité le package : Matrix

Loaded glmnet 4.1-6

glmnet_res = coef(
  glmnet(
    X,
    y,
    alpha = 0,
    lambda = c(10, 1, .1),
    thresh = 1e-12,
    intercept = F
  ), 
  s = .1
)

#' # Comparison

data.frame(
  lm     = coef(lm(y ~ . - 1, data.frame(X))),
  ridge  = result_ridge$par,
  ridge2 = result_ridge2,
  ridge3 = result_ridge3,
  glmnet = glmnet_res[-1, 1],
  truth  = b
)

              lm        ridge       ridge2       ridge3       glmnet  truth
X1   0.534988063  0.485323748  0.485323748  0.485323748  0.485368766  0.500
X2  -0.529993422 -0.480742032 -0.480742032 -0.480742032 -0.480786661 -0.500
X3   0.234376590  0.209412833  0.209412833  0.209412833  0.209435147  0.250
X4  -0.294350608 -0.268814168 -0.268814168 -0.268814168 -0.268837476 -0.250
X5   0.126037566  0.114963716  0.114963716  0.114963716  0.114973801  0.125
X6  -0.159386728 -0.145880488 -0.145880488 -0.145880488 -0.145892837 -0.125
X7  -0.016718534 -0.021658889 -0.021658889 -0.021658889 -0.021655033  0.000
X8   0.009894575  0.006956965  0.006956965  0.006956965  0.006959470  0.000
X9  -0.005441959  0.001392244  0.001392244  0.001392244  0.001386661  0.000
X10  0.010561128  0.010985385  0.010985385  0.010985385  0.010985102  0.000