The oem package provides estimation for various penalized linear models using the Orthogonalizing EM algorithm. Documentation for the package can be found here: oem site.
Install using the devtools package (RcppEigen must be installed first as well):
devtools::install_github("jaredhuling/oem")
or by cloning and building using R CMD INSTALL
To cite oem please use:
Xiong, S., Dai, B., Huling, J., Qian, P. Z. G. (2016) Orthogonalizing
EM: A design-based least squares algorithm, Technometrics, Volume 58,
Pages 285-293,
http://dx.doi.org/10.1080/00401706.2015.1054436.
Huling, J.D. and Chien, P. (2018), Fast Penalized Regression and Cross Validation for Tall Data with the OEM Package, Journal of Statistical Software, to appear, URL: https://arxiv.org/abs/1801.09661.
library(microbenchmark) library(glmnet) library(oem) # compute the full solution path, n > p set.seed(123) n <- 1000000 p <- 100 m <- 25 b <- matrix(c(runif(m), rep(0, p - m))) x <- matrix(rnorm(n * p, sd = 3), n, p) y <- drop(x %*% b) + rnorm(n) lambdas = oem(x, y, intercept = TRUE, standardize = FALSE, penalty = "elastic.net")$lambda[[1]] microbenchmark( "glmnet[lasso]" = {res1 <- glmnet(x, y, thresh = 1e-10, standardize = FALSE, intercept = TRUE, lambda = lambdas)}, "oem[lasso]" = {res2 <- oem(x, y, penalty = "elastic.net", intercept = TRUE, standardize = FALSE, lambda = lambdas, tol = 1e-10)}, times = 5 )
## Unit: seconds
## expr min lq mean median uq max neval cld
## glmnet[lasso] 5.325385 5.374823 5.859432 6.000302 6.292411 6.304239 5 a
## oem[lasso] 1.539320 1.573569 1.600241 1.617136 1.631450 1.639730 5 b
# difference of results max(abs(coef(res1) - res2$beta[[1]]))
## [1] 1.048243e-07
res1 <- glmnet(x, y, thresh = 1e-12, standardize = FALSE, intercept = TRUE, lambda = lambdas) # answers are now more close once we require more precise glmnet solutions max(abs(coef(res1) - res2$beta[[1]]))
## [1] 3.763507e-09
library(sparsenet) library(ncvreg) # compute the full solution path, n > p set.seed(123) n <- 5000 p <- 200 m <- 25 b <- matrix(c(runif(m, -0.5, 0.5), rep(0, p - m))) x <- matrix(rnorm(n * p, sd = 3), n, p) y <- drop(x %*% b) + rnorm(n) mcp.lam <- oem(x, y, penalty = "mcp", gamma = 2, intercept = TRUE, standardize = TRUE, nlambda = 200, tol = 1e-10)$lambda[[1]] scad.lam <- oem(x, y, penalty = "scad", gamma = 4, intercept = TRUE, standardize = TRUE, nlambda = 200, tol = 1e-10)$lambda[[1]] microbenchmark( "sparsenet[mcp]" = {res1 <- sparsenet(x, y, thresh = 1e-10, gamma = c(2,3), #sparsenet throws an error #if you only fit 1 value of gamma nlambda = 200)}, "oem[mcp]" = {res2 <- oem(x, y, penalty = "mcp", gamma = 2, intercept = TRUE, standardize = TRUE, nlambda = 200, tol = 1e-10)}, "ncvreg[mcp]" = {res3 <- ncvreg(x, y, penalty = "MCP", gamma = 2, lambda = mcp.lam, eps = 1e-7)}, "oem[scad]" = {res4 <- oem(x, y, penalty = "scad", gamma = 4, intercept = TRUE, standardize = TRUE, nlambda = 200, tol = 1e-10)}, "ncvreg[scad]" = {res5 <- ncvreg(x, y, penalty = "SCAD", gamma = 4, lambda = scad.lam, eps = 1e-7)}, times = 5 )
## Unit: milliseconds
## expr min lq mean median uq
## sparsenet[mcp] 1466.54465 1492.72548 1527.32113 1517.19926 1579.70827
## oem[mcp] 95.71381 98.09740 105.90083 105.76415 110.31668
## ncvreg[mcp] 5196.48035 5541.69429 5669.10010 5611.31491 5865.06723
## oem[scad] 70.74110 71.46554 80.21926 78.76494 84.25458
## ncvreg[scad] 5289.59790 5810.69254 5801.60997 5950.84377 5964.01276
## max neval cld
## 1580.42800 5 a
## 119.61209 5 b
## 6130.94372 5 c
## 95.87013 5 b
## 5992.90288 5 c
diffs <- array(NA, dim = c(2, 1)) colnames(diffs) <- "abs diff" rownames(diffs) <- c("MCP: oem and ncvreg", "SCAD: oem and ncvreg") diffs[,1] <- c(max(abs(res2$beta[[1]] - res3$beta)), max(abs(res4$beta[[1]] - res5$beta))) diffs
## abs diff
## MCP: oem and ncvreg 1.725859e-07
## SCAD: oem and ncvreg 5.094648e-08
In addition to the group lasso, the oem package offers computation for the group MCP, group SCAD, and group sparse lasso penalties. All aforementioned penalties can also be augmented with a ridge penalty.
library(gglasso) library(grpreg) library(grplasso) # compute the full solution path, n > p set.seed(123) n <- 10000 p <- 200 m <- 25 b <- matrix(c(runif(m, -0.5, 0.5), rep(0, p - m))) x <- matrix(rnorm(n * p, sd = 3), n, p) y <- drop(x %*% b) + rnorm(n) groups <- rep(1:floor(p/10), each = 10) grp.lam <- oem(x, y, penalty = "grp.lasso", groups = groups, nlambda = 100, tol = 1e-10)$lambda[[1]] microbenchmark( "gglasso[grp.lasso]" = {res1 <- gglasso(x, y, group = groups, lambda = grp.lam, intercept = FALSE, eps = 1e-8)}, "oem[grp.lasso]" = {res2 <- oem(x, y, groups = groups, intercept = FALSE, standardize = FALSE, penalty = "grp.lasso", lambda = grp.lam, tol = 1e-10)}, "grplasso[grp.lasso]" = {res3 <- grplasso(x=x, y=y, index = groups, standardize = FALSE, center = FALSE, model = LinReg(), lambda = grp.lam * n * 2, control = grpl.control(trace = 0, tol = 1e-10))}, "grpreg[grp.lasso]" = {res4 <- grpreg(x, y, group = groups, eps = 1e-10, lambda = grp.lam)}, times = 5 )
## Unit: milliseconds
## expr min lq mean median uq
## gglasso[grp.lasso] 679.59049 724.16350 858.99280 801.79179 865.83580
## oem[grp.lasso] 59.84769 62.23879 64.11779 63.36026 64.30146
## grplasso[grp.lasso] 3714.92601 3753.18663 4322.32431 4537.50185 4786.80867
## grpreg[grp.lasso] 1216.21794 1248.84647 1270.46132 1269.71047 1287.75969
## max neval cld
## 1223.58241 5 a
## 70.84075 5 b
## 4819.19839 5 c
## 1329.77201 5 a
diffs <- array(NA, dim = c(2, 1)) colnames(diffs) <- "abs diff" rownames(diffs) <- c("oem and gglasso", "oem and grplasso") diffs[,1] <- c( max(abs(res2$beta[[1]][-1,] - res1$beta)), max(abs(res2$beta[[1]][-1,] - res3$coefficients)) ) diffs
## abs diff
## oem and gglasso 1.303906e-06
## oem and grplasso 1.645871e-08
set.seed(123) n <- 500000 p <- 200 m <- 25 b <- matrix(c(runif(m, -0.5, 0.5), rep(0, p - m))) x <- matrix(rnorm(n * p, sd = 3), n, p) y <- drop(x %*% b) + rnorm(n) groups <- rep(1:floor(p/10), each = 10) # fit all group penalties at once grp.penalties <- c("grp.lasso", "grp.mcp", "grp.scad", "grp.mcp.net", "grp.scad.net", "sparse.group.lasso") system.time(res <- oem(x, y, penalty = grp.penalties, groups = groups, alpha = 0.25, # mixing param for l2 penalties tau = 0.5)) # mixing param for sparse grp lasso
## user system elapsed
## 2.043 0.222 2.267
The oem algorithm is quite efficient at fitting multiple penalties simultaneously when p is not too big.
set.seed(123) n <- 100000 p <- 100 m <- 15 b <- matrix(c(runif(m, -0.25, 0.25), rep(0, p - m))) x <- matrix(rnorm(n * p, sd = 3), n, p) y <- drop(x %*% b) + rnorm(n) microbenchmark( "oem[lasso]" = {res1 <- oem(x, y, penalty = "elastic.net", intercept = TRUE, standardize = TRUE, tol = 1e-10)}, "oem[lasso/mcp/scad/ols]" = {res2 <- oem(x, y, penalty = c("elastic.net", "mcp", "scad", "grp.lasso", "grp.mcp", "sparse.grp.lasso", "grp.mcp.net", "mcp.net"), gamma = 4, tau = 0.5, alpha = 0.25, groups = rep(1:10, each = 10), intercept = TRUE, standardize = TRUE, tol = 1e-10)}, times = 5 )
## Unit: milliseconds
## expr min lq mean median uq max
## oem[lasso] 125.6408 126.7870 130.3534 127.3374 133.4962 138.5055
## oem[lasso/mcp/scad/ols] 148.3162 152.1743 153.0176 152.4529 154.4300 157.7144
## neval cld
## 5 a
## 5 b
#png("../mcp_path.png", width = 3000, height = 3000, res = 400);par(mar=c(5.1,5.1,4.1,2.1));plot(res2, which.model = 2, main = "mcp",lwd = 3,cex.axis=2.0, cex.lab=2.0, cex.main=2.0, cex.sub=2.0);dev.off() # layout(matrix(1:4, ncol=2, byrow = TRUE)) plot(res2, which.model = 1, lwd = 2, xvar = "lambda") plot(res2, which.model = 2, lwd = 2, xvar = "lambda") plot(res2, which.model = 4, lwd = 2, xvar = "lambda") plot(res2, which.model = 7, lwd = 2, xvar = "lambda")
