Replicar resultados para regresión lineal de glmnet utilizando un optimizador genérico

Como dice el título, estoy tratando de replicar los resultados de glmnet linear usando el optimizador LBFGS de la biblioteca lbfgs. Este optimizador nos permite agregar un término regularizador L1 sin tener que preocuparnos por la diferenciabilidad, siempre que nuestra función objetivo (sin el término regularizador L1) sea convexa.

El problema de la regresión lineal neta elástica en el documento glmnet viene dado por donde X \ in \ mathbb {R} ^ {n \ times p} es la matriz de diseño, y \ in \ mathbb {R} ^ p es el vector de observaciones, \ alpha \ in [0,1] es el parámetro neto elástico y \ lambda> 0 es el parámetro de regularización. El operador \ Vert x \ Vert_p denota la norma Lp habitual.

El siguiente código define la función y luego incluye una prueba para comparar los resultados. Como se puede ver, los resultados son aceptables cuando alpha = 1, pero están muy lejos de los valores de alpha < 1.El error empeora a medida que se pasa de alpha = 1que alpha = 0, como muestra la siguiente trama (la "comparación métrica" es la distancia euclídea medio entre las estimaciones de los parámetros de glmnet y lbfgs para una ruta de regularización dada).

Bien, entonces aquí está el código. He agregado comentarios siempre que sea posible. Mi pregunta es: ¿Por qué mis resultados son diferentes de los de los glmnetvalores de alpha < 1? Claramente tiene que ver con el término de regularización L2, pero por lo que puedo decir, he implementado este término exactamente según el documento. Cualquier ayuda sería muy apreciada!

library(lbfgs)
linreg_lbfgs <- function(X, y, alpha = 1, scale = TRUE, lambda) {
  p <- ncol(X) + 1; n <- nrow(X); nlambda <- length(lambda)

  # Scale design matrix
  if (scale) {
    means <- colMeans(X)
    sds <- apply(X, 2, sd)
    sX <- (X - tcrossprod(rep(1,n), means) ) / tcrossprod(rep(1,n), sds)
  } else {
    means <- rep(0,p-1)
    sds <- rep(1,p-1)
    sX <- X
  }
  X_ <- cbind(1, sX)

  # loss function for ridge regression (Sum of squared errors plus l2 penalty)
  SSE <- function(Beta, X, y, lambda0, alpha) {
    1/2 * (sum((X%*%Beta - y)^2) / length(y)) +
      1/2 * (1 - alpha) * lambda0 * sum(Beta[2:length(Beta)]^2) 
                    # l2 regularization (note intercept is excluded)
  }

  # loss function gradient
  SSE_gr <- function(Beta, X, y, lambda0, alpha) {
    colSums(tcrossprod(X%*%Beta - y, rep(1,ncol(X))) *X) / length(y) + # SSE grad
  (1-alpha) * lambda0 * c(0, Beta[2:length(Beta)]) # l2 reg grad
  }

  # matrix of parameters
  Betamat_scaled <- matrix(nrow=p, ncol = nlambda)

  # initial value for Beta
  Beta_init <- c(mean(y), rep(0,p-1)) 

  # parameter estimate for max lambda
  Betamat_scaled[,1] <- lbfgs(call_eval = SSE, call_grad = SSE_gr, vars = Beta_init, 
                              X = X_, y = y, lambda0 = lambda[2], alpha = alpha,
                              orthantwise_c = alpha*lambda[2], orthantwise_start = 1, 
                              invisible = TRUE)$par

  # parameter estimates for rest of lambdas (using warm starts)
  if (nlambda > 1) {
    for (j in 2:nlambda) {
      Betamat_scaled[,j] <- lbfgs(call_eval = SSE, call_grad = SSE_gr, vars = Betamat_scaled[,j-1], 
                                  X = X_, y = y, lambda0 = lambda[j], alpha = alpha,
                                  orthantwise_c = alpha*lambda[j], orthantwise_start = 1, 
                                  invisible = TRUE)$par
    }
  }

  # rescale Betas if required
  if (scale) {
    Betamat <- rbind(Betamat_scaled[1,] -
colSums(Betamat_scaled[-1,]*tcrossprod(means, rep(1,nlambda)) / tcrossprod(sds, rep(1,nlambda)) ), Betamat_scaled[-1,] / tcrossprod(sds, rep(1,nlambda)) )
  } else {
    Betamat <- Betamat_scaled
  }
  colnames(Betamat) <- lambda
  return (Betamat)
}

# CODE FOR TESTING
# simulate some linear regression data
n <- 100
p <- 5
X <- matrix(rnorm(n*p),n,p)
true_Beta <- sample(seq(0,9),p+1,replace = TRUE)
y <- drop(cbind(1,X) %*% true_Beta)

library(glmnet)

# function to compare glmnet vs lbfgs for a given alpha
glmnet_compare <- function(X, y, alpha) {
  m_glmnet <- glmnet(X, y, nlambda = 5, lambda.min.ratio = 1e-4, alpha = alpha)
  Beta1 <- coef(m_glmnet)
  Beta2 <- linreg_lbfgs(X, y, alpha = alpha, scale = TRUE, lambda = m_glmnet$lambda)
  # mean Euclidean distance between glmnet and lbfgs results
  mean(apply (Beta1 - Beta2, 2, function(x) sqrt(sum(x^2))) ) 
}

# compare results
alpha_seq <- seq(0,1,0.2)
plot(alpha_seq, sapply(alpha_seq, function(alpha) glmnet_compare(X,y,alpha)), type = "l", ylab = "Comparison metric")

@ hxd1011 Probé su código, aquí hay algunas pruebas (realicé algunos ajustes menores para que coincidan con la estructura de glmnet; tenga en cuenta que no regularizamos el término de intercepción, y las funciones de pérdida deben ajustarse). Esto es para alpha = 0, pero puede intentarlo alpha: los resultados no coinciden.

rm(list=ls())
set.seed(0)
# simulate some linear regression data
n <- 1e3
p <- 20
x <- matrix(rnorm(n*p),n,p)
true_Beta <- sample(seq(0,9),p+1,replace = TRUE)
y <- drop(cbind(1,x) %*% true_Beta)

library(glmnet)
alpha = 0

m_glmnet = glmnet(x, y, alpha = alpha, nlambda = 5)

# linear regression loss and gradient
lr_loss<-function(w,lambda1,lambda2){
  e=cbind(1,x) %*% w -y
  v= 1/(2*n) * (t(e) %*% e) + lambda1 * sum(abs(w[2:(p+1)])) + lambda2/2 * crossprod(w[2:(p+1)])
  return(as.numeric(v))
}

lr_loss_gr<-function(w,lambda1,lambda2){
  e=cbind(1,x) %*% w -y
  v= 1/n * (t(cbind(1,x)) %*% e) + c(0, lambda1*sign(w[2:(p+1)]) + lambda2*w[2:(p+1)])
  return(as.numeric(v))
}

outmat <- do.call(cbind, lapply(m_glmnet$lambda, function(lambda) 
  optim(rnorm(p+1),lr_loss,lr_loss_gr,lambda1=alpha*lambda,lambda2=(1-alpha)*lambda,method="L-BFGS")$par
))

glmnet_coef <- coef(m_glmnet)
apply(outmat - glmnet_coef, 2, function(x) sqrt(sum(x^2)))

tl; versión dr:

El objetivo contiene implícitamente un factor de escala , donde es la desviación estándar de la muestra.$\hat{s} = sd(y)$$sd(y)$

Versión más larga

Si lee la letra pequeña de la documentación de glmnet, verá:

Tenga en cuenta que la función objetivo para '"gaussiano" es

               1/2  RSS/nobs + lambda*penalty,                  

y para los otros modelos es

               -loglik/nobs + lambda*penalty.                   

Tenga en cuenta también que para '' gaussiano '', 'glmnet' estandariza y para tener una varianza unitaria antes de calcular su secuencia lambda (y luego no estandariza los coeficientes resultantes); Si desea reproducir / comparar resultados con otro software, lo mejor es proporcionar una y estandarizada.

Ahora esto significa que el objetivo es en realidad y ese glmnet informa . ~ β = s

Ahora, cuando estaba usando un lazo puro ( ), entonces la no estandarización de de glmnet significa que las respuestas son equivalentes. Por otro lado, con una cresta pura, entonces debes escalar la penalización por un factor para que la ruta esté de acuerdo, porque un factor adicional de emerge del cuadrado en la penalización . Para intermedio , no hay una manera fácil de escalar la penalización de los coeficientes para reproducir la salida.~ beta 1 / s s ℓ 2$\alpha=1$$\tilde{\beta}$$1/\hat{s}$glmnet$\hat{s}$$\ell_2$$\alpha$glmnets

Una vez que escalo la para tener varianza unitaria, encuentro $y$

que todavía no coincide exactamente Esto parece deberse a dos cosas:

1. La secuencia lambda puede ser demasiado corta para que el algoritmo de descenso de coordenadas cíclicas de inicio en caliente sea totalmente convergente.
2. No hay término de error en sus datos (el de la regresión es 1).$R^2$
3. Tenga en cuenta también que hay un error en el código, según lo previsto, que se necesita lambda[2]para el ajuste inicial, pero eso debería ser lambda[1].

Una vez que rectifico los ítems 1-3, obtengo el siguiente resultado (aunque YMMV depende de la semilla aleatoria):