使用通用优化器复制 glmnet 线性回归的结果

机器算法验证 r 回归 正则化 网络 弹性网
2022-03-16 21:47:36

正如标题所述,我正在尝试使用 library 中的 LBFGS 优化器复制 glmnet linear 的结果lbfgs只要我们的目标函数(没有 L1 正则项)是凸的,这个优化器就允许我们添加 L1 正则项而不必担心可微性。

glmnet论文中的弹性网络线性回归问题由下式给出

minβRp12nβ0+Xβy22+αλβ1+12(1α)λβ22
其中XRn×p是设计矩阵,yRp是观测向量,α[0,1]是弹性网络参数,λ>0是正则化参数。运算符xp表示通常的 Lp 范数。

下面的代码定义了函数,然后包含一个比较结果的测试。如您所见,结果在 时是可以接受alpha = 1对于和给定正则化路径的 lbfgs)。alpha < 1.alpha = 1alpha = 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)))
1个回答

tl;博士版本:

目标隐含地包含一个比例因子,其中是样本标准差。s^=sd(y)sd(y)

更长的版本

如果您阅读 glmnet 文档的细则,您将看到:

请注意,“高斯”的目标函数是

               1/2  RSS/nobs + lambda*penalty,                  

对于其他型号,它是

               -loglik/nobs + lambda*penalty.                   

另请注意,对于 '"gaussian"','glmnet' 在计算其 lambda 序列之前将 y 标准化为具有单位方差(然后对结果系数进行非标准化);如果您希望使用其他软件重现/比较结果,最好提供标准化的 y。

现在这意味着目标实际上是 并且 glmnet 报告

12ny/s^Xβ22+λαβ1+λ(1α)β22,
β~=s^β

现在,当您使用纯套索()时,glmnet 的的非标准化意味着答案是等价的。另一方面,对于纯脊,则需要将惩罚缩放系数以使路径一致,因为从正方形中弹出了一个额外的系数在惩罚中。对于中间,没有一种简单的方法来缩放系数的惩罚以重现输出。α=1β~1/s^glmnets^2αglmnets

一旦我缩放以获得单位方差,我发现 y在此处输入图像描述

这仍然不完全匹配。这似乎是由于两件事:

  1. 对于热启动循环坐标下降算法而言,lambda 序列可能太短而无法完全收敛。
  2. 您的数据中没有错误项(回归的为 1)。R2
  3. 另请注意,提供的代码中存在一个错误,它需要lambda[2]初始拟合,但应该是lambda[1].

一旦我纠正了项目 1-3,我得到以下结果(尽管 YMMV 取决于随机种子):

在此处输入图像描述