r - R 中的套索回归实现有什么问题？

Question

我正在尝试从头开始编写一个函数，该函数使用具有协调下降(Gauss-Seidel)和软阈值的LASSO来估计回归系数。

我的代码如下:

library(MASS)
set.seed(1)
n = 200
p = 200
V = matrix(0.2, p, p)
diag(V) = 1
X = as.matrix(mvrnorm(n, mu = rep(0, p), Sigma = V))
y = X[, 1] + 0.5*X[, 2] + 0.25*X[, 3] + rnorm(n)
X = scale(X)
y = scale(y)

soft_th <- function(b, lambda){
  if (b > lambda){
    return (b - lambda)
  }
  else if (b < -lambda){
    return(b + lambda)
  }
  else {
    return (0)
  }
}

myLasso <- function(X,y, lambda=0.3,tol=1e-5,maxitr=100){
  beta_old <- rep(0,p)
  beta_new <- rep(0,p)
  for(i in 1:maxitr){
    beta_old <- beta_new
    for (j in (1:ncol(X)))
      {
      X_j <- X[,j]
      y_pred <- t(X)%*%beta_old
      rho <- t(X_j)%*%(y - y_pred + beta_old[j]*X_j)
      beta_new[j] <- soft_th(rho,0.7)
    }

    l1 <- sum(abs(beta_old-beta_new))
    print(l1)
    r <-  y - t(X)%*%beta_old

    if (l1<tol){
      print('Convergence reached')
      break
  }
  }

}
myLasso(X,y)

我遇到的问题是 beta_old 和 beta_new 之间的 L1 范数在每次迭代之间增加(很多!)。我遵循这篇文章中的内容:

https://stats.stackexchange.com/questions/123672/coordinate-descent-soft-thresholding-update-operator-for-lasso/351134#351134

我认为我在某种程度上没有正确实现下降更新规则。

如有任何帮助，我们将不胜感激。提前致谢。

Answer 1

我做了更多研究，看来我没有对 X 矩阵进行归一化。定义X后添加X <- X/norm(X,type='2')后，问题解决。

现在我遇到的一个新问题是这个新函数没有复制 LASSO 回归的 glmnet 实现的结果。会是什么呢？我使用 glmnet 得到的 RMSE 为 0.6，使用我的实现得到的 RMSE 为 0.997。如果有人能指导我如何改进我的功能，我会很高兴。

提前致谢。