Сравнение между 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. Кроме того, вам нужно масштабировать соответствующее значение штрафа (lam) вы используете в CVXR стандартное отклонение, потому что значение штрафа, которое вы предоставляете 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