CVXR与glmnet (弹性网)的比较

Question

我在模拟数据集上拟合一个模型，比较glmnet和CVXR的结果。

如果我没有错误的代码，结果是非常不同的。

显式的glmnet产生的结果非常接近真实的参数。

为什么是这种情况？

library(CVXR)
library(glmnet)

set.seed(571)
n = 500
p = 9
x = matrix(rnorm(n*p), ncol=p)
b = c(0.5, 0, 25, -25, 125, -125, rep(0, 3))
y = x %*% b + rnorm(n, sd=.05)

n = nrow(x); p = ncol(x)

lam = 0.4
al = 0.3

# glmnet

glmnet_res = coef(glmnet(x,y,alpha=al,standardize=F,intercept=F),s=lam)[-1]

# CVXR

elastic_reg = function(beta, lambda = 0, alpha = 0) {
  ridge = 0.5*(1 - alpha) * sum(beta^2)
  lasso = alpha * p_norm(beta, 1)
  lambda * (lasso + ridge)
}

beta = Variable(p)  
loss = sum((y - x %*% beta)^2)/(2*n)

## Elastic-net regression
obj = loss + elastic_reg(beta, lam, al)
prob = Problem(Minimize(obj))
result = solve(prob)
beta_vals = result$getValue(beta)

cvxr_res = round(beta_vals,7)

cbind(glmnet_res,cvxr_res)

结果

      glmnet_res    cvxr_res         
 [1,]    0.00000   0.2417734
 [2,]    0.00000   0.0000475
 [3,]   23.39102  19.0372445
 [4,]  -23.26282 -18.6020795
 [5,]  121.59156  96.7286536
 [6,] -121.17658 -95.0466518
 [7,]    0.00000  -1.8589296
 [8,]    0.00000   0.2651426
 [9,]    0.00000   1.0167725

Answer 1

对于连续的结果，glmnet用其标准差来衡量结果(y)。将glmnet中的解决方案与其他软件进行比较的最简单方法是显式缩放y。此外，您还需要根据标准偏差来缩放在CVXR中使用的相应的惩罚值(lam)，因为您提供给coef()的惩罚值也是由y的标准差自动缩放的。从CVXR的估计，然后可以不标准化后，拟合。我还对您的代码做了另外两个小更改：

• 改变glmnet和CVXR的收敛阈值以提高精度
• 增加惩罚值(lam)，因为在CVXR中，对于较大的值，该解决方案更稳定(我发现对于较小的值，它没有达到最优解)

修改代码

library(CVXR)
library(glmnet)

# simulate data
set.seed(571)
n <- 500
p <- 9
x <- matrix(rnorm(n*p), ncol=p)
b <- c(0.5, 0, 25, -25, 125, -125, rep(0, 3))
y <- x %*% b + rnorm(n, sd = .5)
sd_y <- drop(sqrt(var(y) * (n - 1) / n))
y_stnd <- y / sd_y

# fix penalty value and EN parameter
lam <- 20
al <- 0.3

# fit EN in glmnet
fit_glmnet <- glmnet(x = x,
                     y = y,
                     alpha = al,
                     standardize = FALSE, 
                     intercept = FALSE,
                     thresh = 1e-20)

betas_glmnet <- as.vector(coef(fit_glmnet, 
                               s = lam, 
                               exact = TRUE, 
                               x = x, 
                               y = y)[-1])

# fit EN in CVXR (using standardized y and rescaled penalty, lambda / sd_y)
beta <- Variable(p)  
obj <- Minimize(sum((y_stnd - x %*% beta)^2) / (2 * n) + 
                (lam / sd_y) * ((1 - al) * sum_squares(beta) / 2 + al * p_norm(beta, 1)))
prob <- Problem(obj)
result <- solve(prob, solver = "ECOS", verbose = TRUE, ABSTOL = 1e-12, RELTOL = 1e-10)
betas_cvxr <- drop(result$getValue(beta))

# Compare results (unstandardize estimates for CVXR)
round(cbind(betas_glmnet, sd_y * betas_cvxr), 6)

结果

 [1,]      0.00000    0.00000
 [2,]      0.00000    0.00000
 [3,]     17.84706   17.84706
 [4,]    -17.28221  -17.28221
 [5,]    109.82539  109.82539
 [6,]   -108.07262 -108.07262
 [7,]      0.00000    0.00000
 [8,]      0.00000    0.00000
 [9,]      0.00000    0.00000