复现glmnet的LASSO回归结果与官方输出存在差异求助

差异原因排查及修正方案 #

核心问题点 #

• 初始值设置错误：LASSOstandardized函数中beta_start默认值为NULL，未显式初始化为全零向量，导致初始残差计算逻辑错误，需要将默认值改为rep(0, ncol(Xtilde))。
• 特征标准化逻辑不匹配：glmnet默认对输入特征做自动标准化（standardize=TRUE），先将每列特征缩放到标准差为1，拟合完成后再把系数转换回原始特征尺度，你的代码直接使用原始尺度的特征计算，两者的lambda惩罚量级完全不对等。
• 截距拟合逻辑不匹配：glmnet默认自动拟合截距项（intercept=TRUE），你构造的输入x已经通过model.matrix(~.-1)去掉了截距，且自定义函数没有处理截距拟合逻辑，两者的拟合目标不一致。
• lambda计算精度问题：glmnet默认如果传入的lambda不在自动生成的lambda序列范围内，会通过插值返回近似结果，而非精确拟合指定的lambda值，需要开启exact=TRUE参数获得精确结果。

对齐后的测试代码 #

首先修正自定义函数的初始值问题：

soft <- function(a, lambda){
  # calculate new beta
  new_beta = sign(a) * max(abs(a) - lambda, 0)
  return(new_beta)
}

lasso <- function(Xtilde, Ytilde, beta, lambda){
  return((crossprod(Ytilde - Xtilde %*% beta)) /(2 * length(Ytilde)) + lambda * sum(abs(beta)))
}

LASSOstandardized <- function(Xtilde, Ytilde, lambda, eps, beta_start = NULL){
  # 补充默认初始值为全零向量
  if(is.null(beta_start)) beta_start = rep(0, ncol(Xtilde))
  beta_last = beta_start
  beta_new = beta_start
  r = Ytilde - Xtilde %*% beta_last
  
  repeat {
    beta_last = beta_new
    for (i in 1:ncol(Xtilde)) {
      beta_new[i] = soft((beta_last[i] + crossprod(r, Xtilde[, i])/nrow(Xtilde)), lambda )
      r = r + Xtilde[, i] * (beta_last[i] - beta_new[i])
    }
    
    loss_old = lasso(Xtilde, Ytilde, beta_last, lambda)
    loss_new = lasso(Xtilde, Ytilde, beta_new, lambda)
    loss_diff = loss_old - loss_new
    
    if (abs(loss_diff) < eps){
      break
    }
  }
  return(beta_new)
}

调用时对齐glmnet的参数：

library(ISLR)
library(glmnet)
Hitters=na.omit(Hitters)[1:259,c(1:13,19)]
x=model.matrix(Salary~.-1,data=Hitters) 
y=Hitters$Salary
# 自定义函数计算结果
custom_res <- LASSOstandardized(x,y,2,0.0000001)
# glmnet关闭截距拟合、关闭标准化、开启精确拟合
glmnet_res <- glmnet(x,y, alpha=1,lambda=2, intercept=FALSE, standardize=FALSE, exact=TRUE, x=x, y=y)$beta
# 验证结果一致
all.equal(as.vector(glmnet_res), custom_res)

内容的提问来源于stack exchange，提问作者mbeauch