How does glmnet handle larger datasets?

Question

I'm looking to fit a model with about 1k-40k variables and up to a few million observations.

Can anyone with a bit more experience speak to its performance for larger datasets?

It looks like I can use R's sparse matrix support to represent the data in a much smaller form (>90% of variables will have zero values for any given observation). However, the few examples I've found that use glmnet with a sparse matrix appear to also have very few observations.

I'm left a bit unsure whether to prep my data for glmnet if it won't be able to converge to a solution in a reasonable amount of time, so any advice would be appreciated.

Answer 1

glmnet is very strong in this respect. Of course it depends on the situation but to give you an impression of its performance, I made an example with 3 mio lines and a sparse model matrix of dimension 40'001 (one factor with 40k levels and one normally distributed variable). It takes about half a minute to run on my normal laptop to find the regularization path of length 44.

library(glmnet)
library(Matrix)

n <- 3e6 # 3 mio

set.seed(483487)
y <- sample(0:1, n, replace = TRUE) 
x1 <- rnorm(n)
x2 <- sample(1:40000, n, replace = TRUE)

X <- sparse.model.matrix(~ x1 + factor(x2))
dim(X)  # 3000000   40001

# Runs 25.32 seconds on normal laptop
system.time(fit <- glmnet(y = y, x = X, family = "binomial", alpha = 0.5))
fit 

Unfortunately, glmnet does not return standard errors etc. for its coefficients, even in the unregularized mode. (Using regularization invalidates quantities like this.)

An alternative is the h2o R package that connects to h2o.ai backend, a small JAVA program designed for high performance machine learning. Their elastic net GLM implementation is able to provide p-values, standard errors etc. for your coefficients as long as you use the slow "IRWLS" optimizer and no regularization. If you are not interested in p values etc., it is as fast as glmnet. You can fit data sets of any size as long as it fits into memory.

library(h2o)
h2o.init(max_mem_size = "4G", nthreads = -1)

# 20 seconds
train <- as.h2o(data.frame(y = y, x1 = x1, x2 = factor(x2)))

# still at 2% done after 1h
system.time(fit_h2o <- h2o.glm(x = c("x1", "x2"), y = "y", 
                               training_frame = train, 
                               family = "binomial",
                               compute_p_values = TRUE,
                               lambda = 0,
                               alpha = 0,
                               solver = "IRLSM"))