Réplication des résultats de la régression linéaire glmnet à l'aide d'un optimiseur générique

Comme l'indique le titre, j'essaie de reproduire les résultats de glmnet linear en utilisant l'optimiseur LBFGS de la bibliothèque lbfgs. Cet optimiseur nous permet d'ajouter un terme de régularisateur L1 sans avoir à se soucier de la différentiabilité, tant que notre fonction objectif (sans le terme de régularisateur L1) est convexe.

Le problème de régression linéaire nette élastique dans le papier glmnet est donné par où X \ in \ mathbb {R} ^ {n \ fois p} est la matrice de conception, y \ in \ mathbb {R} ^ p est le vecteur des observations, \ alpha \ in [0,1] est le paramètre net élastique et \ lambda> 0 est le paramètre de régularisation. L'opérateur \ Vert x \ Vert_p désigne la norme Lp habituelle. X∈Rn×py∈Rpα∈[0,1]λ>0‖x‖

Le code ci-dessous définit la fonction, puis inclut un test pour comparer les résultats. Comme vous pouvez le constater, les résultats sont acceptables alpha = 1, mais sont loin des valeurs de alpha < 1.l'erreur se aggrave que nous passons de alpha = 1la alpha = 0, comme le montre le tracé suivant (la « comparaison métrique » est la distance euclidienne moyenne entre les estimations des paramètres de glmnet et lbfgs pour un chemin de régularisation donné).

D'accord, voici donc le code. J'ai ajouté des commentaires chaque fois que possible. Ma question est: pourquoi mes résultats sont-ils différents de ceux des glmnetvaleurs de alpha < 1? Cela a clairement à voir avec le terme de régularisation L2, mais pour autant que je sache, j'ai mis en œuvre ce terme exactement selon le document. Toute aide serait très appréciée!

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 J'ai essayé votre code, voici quelques tests (j'ai fait quelques ajustements mineurs pour correspondre à la structure de glmnet - notez que nous ne régularisons pas le terme d'interception, et les fonctions de perte doivent être mises à l'échelle). C'est pour alpha = 0, mais vous pouvez en essayer alpha- les résultats ne correspondent pas.

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; version dr:

L'objectif contient implicitement un facteur d'échelle , où est l'écart type de l'échantillon.sd(y$\hat{s} = sd(y)$$sd(y)$

Version plus longue

Si vous lisez les petits caractères de la documentation de glmnet, vous verrez:

Notez que la fonction objectif pour «gaussien» est

               1/2  RSS/nobs + lambda*penalty,                  

et pour les autres modèles c'est

               -loglik/nobs + lambda*penalty.                   

Notez également que pour «gaussien», «glmnet» standardise y pour avoir une variance unitaire avant de calculer sa séquence lambda (puis standardise les coefficients résultants); si vous souhaitez reproduire / comparer les résultats avec d'autres logiciels, il est préférable de fournir un y standardisé.

Cela signifie maintenant que l'objectif est en fait et que glmnet signale . ~ β = s

Maintenant, lorsque vous utilisiez un lasso pur ( ), alors la normalisation du de glmnet signifie que les réponses sont équivalentes. D'un autre côté, avec une crête pure, vous devez ensuite mettre la pénalité à l'échelle d'un facteur pour que le chemin soit d'accord, car un facteur supplémentaire de apparaît du carré dans la pénalité . Pour intermédiaire , il n'y a pas de moyen facile de mettre à l'échelle la pénalité des coefficients pour reproduire la sortie.~ ß 1 / s de ℓ 2$\alpha=1$$\tilde{\beta}$$1/\hat{s}$glmnet$\hat{s}$$\ell_2$$\alpha$glmnets

Une fois que j'ai mis à l'échelle le pour avoir la variance unitaire, je trouve $y$

qui ne correspond toujours pas exactement. Cela semble être dû à deux choses:

1. La séquence lambda peut être trop courte pour que l'algorithme de descente de coordonnées cycliques à démarrage à chaud soit entièrement convergent.
2. Il n'y a pas de terme d'erreur dans vos données (le de la régression est 1).$R^2$
3. Notez également qu'il y a un bogue dans le code tel qu'il est fourni dans lequel il faut lambda[2]pour l'ajustement initial, mais cela devrait l'être lambda[1].

Une fois que j'ai rectifié les éléments 1-3, j'obtiens le résultat suivant (bien que YMMV en fonction de la graine aléatoire):