A comparison between CVXR and glmnet (elastic net)

Question

I fit a model on a simulated data set to compare glmnet and CVXR results.

If I don't have a code mistake, the results are very different.

Explicitly glmnet yields results which are very close to the true parameters.

Why is this the case?

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)

Results

      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

skijunkie · Answer 1 · 2018-11-16T21:12:17.333

For continuous outcomes, glmnet scales the outcome (y) by its standard deviation. The easiest way to compare solutions in glmnet to other software is to explicitly scale y. Additionally, you need to scale the corresponding penalty value (lam) you use in CVXR by the standard deviation, because the penalty value that you provide to coef() is automatically scaled by the standard deviation of y as well. The estimates from CVXR can then be unstandardized after fitting. I also made two other small changes to your code:

Changed the convergence threshold for both glmnet and CVXR to increase accuracy
Increased the penalty value (lam) as the solution is more stable in CVXR for larger values (I found that it was not reaching an optimal solution for small values)

Modified Code

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)

Results

 [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

A comparison between CVXR and glmnet (elastic net)

1 Answers1