Elastic Net

Introduction

Often in applications, we encounter problems that require regularization to prevent overfitting, introduce sparsity, facilitate variable selection, or impose prior distributions on parameters. Two of the most common regularization functions are the \(l_1\)-norm and squared \(l_2\)-norm, combined in the elastic net regression model (H. Zou 2005, Friedman, Hastie, and Tibshirani (2010)),

\[ \begin{array}{ll} \underset{\beta}{\mbox{minimize}} & \frac{1}{2m}\|y - X\beta\|_2^2 + \lambda(\frac{1-\alpha}{2}\|\beta\|_2^2 + \alpha\|\beta\|_1). \end{array} \]

Here \(\lambda \geq 0\) is the overall regularization weight and \(\alpha \in [0,1]\) controls the relative \(l_1\) versus squared \(l_2\) penalty. Thus, this model encompasses both ridge (\(\alpha = 0\)) and lasso (\(\alpha = 1\)) regression.

Example

To solve this problem in CVXR, we first define a function that calculates the regularization term given the variable and penalty weights.

suppressMessages(suppressWarnings(library(CVXR)))
elastic_reg <- function(beta, lambda = 0, alpha = 0) {
    ridge <- (1 - alpha) * sum(beta^2)
    lasso <- alpha * p_norm(beta, 1)
    lambda * (lasso + ridge)
}

Later, we will add it to the scaled least squares loss as shown below.

loss <- sum((y - X %*% beta)^2)/(2*m)
obj <- loss + elastic_reg(beta, lambda, alpha)

The advantage of this modular approach is that we can easily incorporate elastic net regularization into other regression models. For instance, if we wanted to run regularized Huber regression, CVXR allows us to reuse the above code with just a single changed line.

loss <- huber(y - X %*% beta, M)

We generate some synthetic sparse data for this example.

set.seed(1)

# Problem data
n <- 20
m <- 1000
DENSITY <- 0.25    # Fraction of non-zero beta
beta_true <- matrix(rnorm(n), ncol = 1)
idxs <- sample.int(n, size = floor((1 - DENSITY) * n), replace = FALSE)
beta_true[idxs] <- 0

sigma <- 45
X <- matrix(rnorm(m * n, sd = 5), nrow = m, ncol = n)
eps <- matrix(rnorm(m, sd = sigma), ncol = 1)
y <- X %*% beta_true + eps

We fit the elastic net model for several values of \(\lambda\) .

TRIALS <- 10
beta_vals <- matrix(0, nrow = n, ncol = TRIALS)
lambda_vals <- 10^seq(-2, log10(50), length.out = TRIALS)
beta <- Variable(n)  
loss <- sum((y - X %*% beta)^2)/(2*m)

## Elastic-net regression
alpha <- 0.75
for(i in 1:TRIALS) {
    lambda <- lambda_vals[i]
    obj <- loss + elastic_reg(beta, lambda, alpha)
    prob <- Problem(Minimize(obj))
    result <- solve(prob)
    beta_vals[,i] <- result$getValue(beta)
}

We can now get a table of the coefficients.

d <- data.frame(lambda_vals, t(beta_vals))
names(d) <- c("$\\lambda$", sprintf("$\\beta_{%d}$", seq_len(n)))
rownames(d) <- paste0("s", 0:9)
library(kableExtra)
knitr::kable(t(d), format = "html", caption = "Elastic net fits from `CVXR`", digits = 3) %>%
    kable_styling("striped") %>%
    column_spec(1:11, background = "#ececec")

We plot the coefficients against the regularization.

plot(0, 0, type = "n", main = "Regularization Path for Elastic-net Regression",
     xlab = expression(lambda), ylab = expression(beta),
     ylim = c(-0.75, 1.25), xlim = c(0, 50))
matlines(lambda_vals, t(beta_vals))

We can also compare with the glmnet results.

library(glmnet)

## Loading required package: Matrix

## Loading required package: foreach

## Loaded glmnet 2.0-13

model_net <- glmnet(X, y, family = "gaussian", alpha = alpha,
                    lambda = lambda_vals, standardize = FALSE,
                    intercept = FALSE)
## Reverse order to match beta_vals
coef_net <- data.frame(lambda_vals, t(as.matrix(coef(model_net)[-1, seq(TRIALS, 1, by = -1)])))
names(coef_net) <- c("$\\lambda$", sprintf("$\\beta_{%d}$", seq_len(n)))
knitr::kable(t(coef_net), format = "html", digits = 3, caption = "Coefficients from `glmnet`") %>%
    kable_styling("striped") %>%
    column_spec(1:11, background = "#ececec")

Session Info

sessionInfo()

## R version 3.4.3 (2017-11-30)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS High Sierra 10.13.3
## 
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] methods   stats     graphics  grDevices datasets  utils     base     
## 
## other attached packages:
## [1] glmnet_2.0-13    foreach_1.4.4    Matrix_1.2-12    kableExtra_0.7.0
## [5] CVXR_0.95       
## 
## loaded via a namespace (and not attached):
##  [1] gmp_0.5-13.1      Rcpp_0.12.15      highr_0.6        
##  [4] plyr_1.8.4        compiler_3.4.3    pillar_1.1.0     
##  [7] iterators_1.0.9   R.methodsS3_1.7.1 R.utils_2.6.0    
## [10] tools_3.4.3       digest_0.6.15     bit_1.1-12       
## [13] viridisLite_0.3.0 evaluate_0.10.1   tibble_1.4.2     
## [16] lattice_0.20-35   pkgconfig_2.0.1   rlang_0.2.0      
## [19] yaml_2.1.16       blogdown_0.5.4    xfun_0.1         
## [22] xml2_1.2.0        Rmpfr_0.7-0       ECOSolveR_0.4    
## [25] stringr_1.3.0     httr_1.3.1        knitr_1.20       
## [28] hms_0.4.1         rprojroot_1.3-2   bit64_0.9-7      
## [31] grid_3.4.3        R6_2.2.2          rmarkdown_1.8.10 
## [34] bookdown_0.7      readr_1.1.1       magrittr_1.5     
## [37] codetools_0.2-15  scales_0.5.0      backports_1.1.2  
## [40] htmltools_0.3.6   scs_1.1-1         rvest_0.3.2      
## [43] colorspace_1.3-2  stringi_1.1.6     munsell_0.4.3    
## [46] R.oo_1.21.0

Source

R Markdown