Replicación de los resultados de la regresión lineal glmnet mediante un optimizador genérico

Question

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 añadir 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 regresión lineal de la red elástica en el papel glmnet viene dada por $$\min_{\beta \in \mathbb{R}^p} \frac{1}{2n}\Vert \beta_0 + X\beta - y \Vert_2^2 + \alpha \lambda \Vert \beta\Vert_1 + \frac{1}{2}(1-\alpha)\lambda\Vert\beta\Vert^2_2$$ 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 de la red elástica, 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 puede ver, los resultados son aceptables cuando alpha = 1 pero están muy lejos para valores de alpha < 1. El error empeora a medida que pasamos de alpha = 1 a alpha = 0 como muestra el siguiente gráfico (la "métrica de comparación" es la distancia euclidiana media entre las estimaciones de los parámetros de glmnet y lbfgs para una ruta de regularización determinada).

Bien, aquí está el código. He añadido comentarios siempre que ha sido posible. Mi pregunta es: ¿Por qué mis resultados son diferentes a los de glmnet para los valores de alpha < 1 ? Está claro que tiene que ver con el término de regularización L2, pero por lo que veo, he implementado este término exactamente como en el artículo. Cualquier ayuda será 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é tu código, aquí hay algunas pruebas (hice algunos ajustes menores para que coincida con la estructura de glmnet - nota que no regularizamos el término de intercepción, y las funciones de pérdida deben ser escaladas). Esto es para alpha = 0 , pero puedes probar con cualquier 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)))

Answer 1

Versión tl;dr:

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

Versión más larga

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

Obsérvese que la función objetivo para el "gaussiano" es

               1/2  RSS/nobs + lambda*penalty,                  

y para los otros modelos es

               -loglik/nobs + lambda*penalty.                   

Obsérvese también que para '"gaussiano", 'glmnet' estandariza y para que tenga una varianza unitaria varianza unitaria antes de calcular su secuencia lambda (y luego los coeficientes resultantes); si desea reproducir/comparar los resultados Si desea reproducir/comparar los resultados con otro software, lo mejor es suministrar un estandarizado.

Ahora bien, esto significa que el objetivo es realmente $$ \frac{1}{2n} \lVert y/\hat{s}-X\beta \rVert_2^2 + \lambda\alpha \lVert \beta \rVert_1 + \lambda(1-\alpha)\lVert \beta \rVert_2^2, $$ y que glmnet informa $\tilde{\beta} = \hat{s} \beta$ .

Ahora bien, cuando se utiliza un lazo puro ( $\alpha=1$ ), entonces la no estandarización de glmnet's $\tilde{\beta}$ significa que las respuestas son equivalentes. Por otro lado, con una cresta pura, entonces hay que escalar la penalización por un factor $1/\hat{s}$ para que el glmnet camino para estar de acuerdo, porque un factor extra de $\hat{s}$ sale de la plaza en la $\ell_2$ pena. Para los intermedios $\alpha$ no hay una manera fácil de escalar la penalización de los coeficientes para reproducir glmnets de salida.

Una vez que escale el $y$ para tener una varianza unitaria, encuentro

que sigue sin coincidir exactamente. Esto parece ser debido a dos cosas:

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

Una vez que rectifico los puntos 1 a 3, obtengo el siguiente resultado (aunque depende de la semilla aleatoria):

Answer 2

Puede que esto no responda a su pregunta, pero debería ser útil: Yo eliminaría el efecto por lbfgs . Y utilizar una cantidad mínima de código para comprobar.

Aquí hay una cantidad mínima de código para hacer la regresión lineal glmnet. Estoy usando optimx con comprobación de gradiente. Incluso puede utilizar su propio código para optimizar el objetivo.

    set.seed(0)
    n_data=1000
    n_features=5
    x=matrix(runif(n_data*n_features),ncol=n_features)
    y=runif(n_data)

    # linear regression loss and gradient
    lr_loss<-function(w,lambda1,lambda2){
      e=x %*% w -y
      v=t(e) %*% e + lambda1 * sum(abs(w)) + lambda2 * crossprod(w)
      return(as.numeric(v))
    }

    lr_loss_gr<-function(w,lambda1,lambda2){
      e=x %*% w -y
      v=2 * t(x) %*% e + lambda1*sign(w)+2*lambda2*w
      return(as.numeric(v))
    }

    optimx::optimx(runif(n_features),lr_loss,lr_loss_gr,
                   lambda1=1,lambda2=1,method="L-BFGS")