R
使用通用優化器複製 glmnet 線性回歸的結果
正如標題所述，我正在嘗試使用 library 中的 LBFGS 優化器從 glmnet linear 複製結果
lbfgs。只要我們的目標函數（沒有 L1 正則項）是凸的，這個優化器就允許我們添加 L1 正則項而不必擔心可微性。glmnet論文中的彈性網絡線性回歸問題由下式給出
在哪裡是設計矩陣，是觀察向量，是彈性網絡參數，並且是正則化參數。運營商表示通常的 Lp 範數。 下面的代碼定義了函數，然後包含一個比較結果的測試。如您所見，結果在 時是可以接受的，
alpha = 1但對於和給定正則化路徑的 lbfgs）。alpha < 1.``alpha = 1``alpha = 0
好的，這裡是代碼。我盡可能地添加了評論。我的問題是：為什麼我的結果與 的
glmnet值不同alpha < 1？它顯然與 L2 正則化術語有關，但據我所知，我已經完全按照論文實現了這個術語。任何幫助將非常感激！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 我嘗試了你的代碼，這裡有一些測試（我做了一些小的調整以匹配 glmnet 的結構 - 注意我們沒有規範截距項，並且必須縮放損失函數）。這是為
alpha = 0，但您可以嘗試任何alpha- 結果不匹配。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;博士版本：
目標隱含地包含一個比例因子， 在哪裡是樣本標準差。
更長的版本
如果您閱讀 glmnet 文檔的細則，您將看到：
請注意，“高斯”的目標函數是
1/2 RSS/nobs + lambda*penalty,對於其他型號，它是
-loglik/nobs + lambda*penalty.另請注意，對於 ‘“gaussian”'，‘glmnet’ 在計算其 lambda 序列之前將 y 標準化為具有單位方差（然後對結果係數進行非標準化）；如果您希望使用其他軟件重現/比較結果，最好提供標準化的 y。
現在這意味著目標實際上是
那glmnet報告. 現在，當您使用純套索時（)，然後是 glmnet 的非標準化表示答案是等價的。另一方面，如果是純嶺，那麼您需要將懲罰擴大一個因子為了讓
glmnet路徑一致，因為一個額外的因素從廣場中彈出懲罰。對於中級，沒有一種簡單的方法來衡量係數的懲罰來重現glmnets輸出。一旦我縮放有單位方差，我發現
這仍然不完全匹配。這似乎是由於兩件事：
- 對於熱啟動循環坐標下降算法而言，lambda 序列可能太短而無法完全收斂。
- 您的數據中沒有錯誤術語（的回歸是 1)。
- 另請注意，提供的代碼中存在一個錯誤，它需要
lambda[2]初始擬合，但應該是lambda[1].一旦我糾正了項目 1-3，我得到以下結果（儘管 YMMV 取決於隨機種子）：
