Hasil replikasi untuk regresi linier glmnet menggunakan pengoptimal generik

Seperti yang dinyatakan judul, saya mencoba mereplikasi hasil dari glmnet linear menggunakan pengoptimal LBFGS dari perpustakaan lbfgs. Pengoptimal ini memungkinkan kita untuk menambahkan istilah regularizer L1 tanpa harus khawatir tentang diferensiabilitas, selama fungsi objektif kami (tanpa istilah regularizer L1) adalah cembung.

min_{β \in R^{p}} \frac{1}{2 n} ‖ β_{0} + X β - y ‖_{2}^{2} + α λ ‖ β ‖_{1} + \frac{1}{2} (1 - α) λ ‖ β ‖_{2}^{2}

$\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$

X \in R^{n \times p}

$X \in \mathbb{R}^{n \times p}$

y \in R^{p}

$y \in \mathbb{R}^p$

α \in [0, 1]

$\alpha \in [0,1]$

λ > 0

$\lambda > 0$

‖ x ‖_{p}

$\Vert x \Vert_p$

Kode di bawah ini mendefinisikan fungsi, dan kemudian menyertakan tes untuk membandingkan hasilnya. Seperti yang Anda lihat, hasilnya dapat diterima kapan alpha = 1, tetapi jauh untuk nilai-nilai alpha < 1.Kesalahan semakin memburuk saat kita beralih alpha = 1ke alpha = 0, seperti yang ditunjukkan plot berikut ("metrik perbandingan" adalah jarak Euclidean rata-rata antara estimasi parameter glmnet dan lbfgs untuk jalur regularisasi yang diberikan).

Oke, jadi ini kodenya. Saya telah menambahkan komentar sedapat mungkin. Pertanyaan saya adalah: Mengapa hasil saya berbeda glmnetdengan nilai untuk alpha < 1? Ini jelas ada hubungannya dengan istilah regularisasi L2, tapi sejauh yang saya tahu, saya sudah menerapkan istilah ini persis seperti yang tertulis di kertas. Bantuan apa pun akan sangat dihargai!

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 Saya mencoba kode Anda, berikut adalah beberapa tes (saya membuat beberapa perubahan kecil agar sesuai dengan struktur glmnet - perhatikan kami tidak mengatur istilah intersep, dan fungsi kerugian harus ditingkatkan). Ini untuk alpha = 0, tetapi Anda dapat mencoba alpha- hasilnya tidak cocok.

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)))

r regression regularization glmnet elastic-net pengguna3294195
sumber

Saya tidak yakin pertanyaan Anda tentang topik (saya pikir itu mungkin, karena ini tentang teknik optimasi yang mendasarinya), dan saya tidak dapat benar-benar memeriksa kode Anda sekarang, tetapi lbfgsmemunculkan poin tentang orthantwise_cargumen tentang glmnetkesetaraan.

Firebug

Masalahnya tidak benar-benar dengan lbfgsdan orthantwise_c, seperti kapan alpha = 1, solusinya hampir sama persis dengan glmnet. Ini ada hubungannya dengan sisi regularisasi L2 hal yaitu kapan alpha < 1. Saya pikir membuat semacam modifikasi untuk definisi SSEdan SSE_grharus memperbaikinya, tapi saya tidak yakin apa yang harus modifikasi - sejauh yang saya tahu, fungsi-fungsi tersebut didefinisikan persis seperti yang dijelaskan dalam kertas glmnet.

user3294195

Ini mungkin lebih dari stackoverflow, pertanyaan pemrograman.

Matthew Gunn

Saya pikir ini lebih berkaitan dengan optimasi & regularisasi daripada kode itu sendiri, itulah sebabnya saya diposting di sini.

user3294195

Untuk pertanyaan optimasi murni, scicomp.stackexchange.com juga merupakan opsi. Saya tidak yakin apakah pertanyaan khusus bahasa ada di topik di sana? (mis. "lakukan ini dalam R")

GeoMatt22

Hasil replikasi untuk regresi linier glmnet menggunakan pengoptimal generik

Jawaban:

tl; versi dr:

Versi yang lebih panjang