Appendix A — Comparing glmnet and CVX

CVX is a software package for convex optimization, primarily used in Matlab and python (via the cvxpy package). R users can access CVX functionality through the CVXR package.

Since the optimization problems defining the lasso and elastic net models are convex optimization problems, it is natural to try to solve them in CVXR and compare the output with that from glmnet. This note explains what you need to do to get the outputs to match.

Let’s start with a dataset:

data(QuickStartExample)
x <- QuickStartExample$x
y <- QuickStartExample$y

Next, let’s say we want to solve the following problem:

\[ \min_{\mathbf{\beta}} \| \mathbf{y} - \mathbf{X\beta} \|_2^2 + \lambda_0 \|\beta \|_1, \tag{A.1}\]

with \(\lambda_0 = 8\).

A.1 Solving with glmnet

First, recall that in the minimization problem that glmnet solves (Equation 2.4), there is a \(1/(2n)\) factor in front of the first term (the RSS term). The problem we want to solve is equivalent to

\[ \min_{\mathbf{\beta}} \frac{1}{2n} \| \mathbf{y} - \mathbf{X\beta} \|_2^2 + \frac{\lambda_0}{2n} \|\beta \|_1, \]

so we need to modify the \(\lambda\) hyperparameter accordingly: lambda = lambda_0 / (2 * n).

There are two other important points to note:

There is no intercept term in Equation A.1 but there is one in glmnet’s minimization problem (@eq-lasso). To account for this, we need to set intercept = FALSE.
Equation A.1 uses the \(\mathbf{X}\) matrix as-is without any standardization of columns. Thus, we must set standardize = FALSE.

This is the glmnet call we would use to solve Equation A.1:

n <- nrow(x)
fit <- glmnet(x, y, intercept = FALSE, standardize = FALSE, lambda = 8 / (2 * n),
              thresh = 1e-20)
#> Warning in .resolve_control(control = control, nvars = nvars, deprecated =
#> list(thresh = if (!missing(thresh)) thresh, : Passing 'thresh' to glmnet()
#> is deprecated. Use control = list(thresh = ...) instead.

In the above, we set thresh = 1e-20 to tighten the convergence criteria and get a more accurate solution for comparison. The default value of thresh = 1e-7 is sufficient for most purposes.

The line of code below extracts the coefficients, dropping the intercept term:

beta_glmnet <- coef(fit)[-1,]

A.2 Solving with CVXR

The code below solves Equation A.1 with \(\lambda_0 = 8\):

library(CVXR)
#> 
#> Attaching package: 'CVXR'
#> The following objects are masked from 'package:Matrix':
#> 
#>     diag, norm
#> The following objects are masked from 'package:stats':
#> 
#>     convolve, power, sd, var
#> The following objects are masked from 'package:base':
#> 
#>     diag, norm, outer
n <- nrow(x)
p <- ncol(x)
lambda_0 <- 8

# create the optimization variables
beta <- Variable(p)

# create the loss function
rss <- sum((y-x %*% beta)^2)
lassoPenalty <- function(beta, lambda) lambda * p_norm(beta, 1)
obj <- rss + lassoPenalty(beta, lambda = lambda_0)

# define and solve the problem
prob <- Problem(Minimize(obj))
obj_result <- psolve(prob)

# extract results
beta_CVXR <- value(beta)

A.3 Comparing the results

The results are listed below up to 6 decimal digits (due to convergence thresholds). We see that all three packages give the same result.

cmp <- round(cbind(beta_glmnet, beta_CVXR), digits = 6)
colnames(cmp) <- c("beta_glmnet", "beta_CVXR")
cmp
#>     beta_glmnet beta_CVXR
#> V1     1.389118  1.389118
#> V2     0.007991  0.007991
#> V3     0.731234  0.731234
#> V4     0.031119  0.031119
#> V5    -0.866793 -0.866793
#> V6     0.564867  0.564867
#> V7     0.069678  0.069678
#> V8     0.358346  0.358346
#> V9     0.000000  0.000000
#> V10    0.070565  0.070565
#> V11    0.173464  0.173464
#> V12   -0.027472 -0.027472
#> V13   -0.017960 -0.017960
#> V14   -1.138053 -1.138053
#> V15   -0.070990 -0.070990
#> V16    0.000000  0.000000
#> V17    0.000000  0.000000
#> V18    0.000000  0.000000
#> V19    0.000000  0.000000
#> V20   -1.097528 -1.097528