Introduction to the parglm package
Benjamin Christoffersen
2019-01-19
The motivation for the parglm package is a parallel version of the glm function. It solves the iteratively re-weighted least squares using a QR decomposition with column pivoting with DGEQP3 function from LAPACK. The computation is done in parallel as in the bam function in the mgcv package. The cost is an additional \(O(Mp^2 + p^3)\) where \(p\) is the number of coefficients and \(M\) is the number chunks to be computed in parallel. The advantage is that you do not need to compile the package with an optimized BLAS or LAPACK which supports multithreading. The package also includes a method that computes the Fisher information and then solves the normal equation as in the speedglm package. This is faster but less numerically stable.
Example of computation time
Below, we perform estimate a logistic regression with 1000000 observations and 50 covariates. We vary the number of cores being used with the nthreads argument to parglm.control. The method arguments sets which method is used. LINPACK yields the QR method and FAST yields a similar method as in the speedglm package.
#####
# simulate
n # number of observations
#> [1] 1000000
p # number of covariates
#> [1] 50
set.seed(68024947)
X <- matrix(rnorm(n * p, 1/p, 1/sqrt(p)), n, ncol = p)
df <- data.frame(y = 1/(1 + exp(-(rowSums(X) - 1))) > runif(n), X)#####
# compute and measure time. Setup call to make
library(microbenchmark)
library(speedglm)
library(parglm)
cl <- list(
quote(microbenchmark),
glm = quote(glm (y ~ ., binomial(), df)),
speedglm = quote(speedglm(y ~ ., family = binomial(), data = df)),
times = 11L)
tfunc <- function(method = "LINPACK", nthreads)
parglm(y ~ ., binomial(), df, control = parglm.control(method = method,
nthreads = nthreads))
cl <- c(
cl, lapply(1:n_threads, function(i) bquote(tfunc(nthreads = .(i)))))
names(cl)[5:(5L + n_threads - 1L)] <- paste0("parglm-LINPACK-", 1:n_threads)
cl <- c(
cl, lapply(1:n_threads, function(i) bquote(tfunc(
nthreads = .(i), method = "FAST"))))
names(cl)[(5L + n_threads):(5L + 2L * n_threads - 1L)] <-
paste0("parglm-FAST-", 1:n_threads)
cl <- as.call(cl)
cl # the call we make
#> microbenchmark(glm = glm(y ~ ., binomial(), df), speedglm = speedglm(y ~
#> ., family = binomial(), data = df), times = 11L, `parglm-LINPACK-1` = tfunc(nthreads = 1L),
#> `parglm-LINPACK-2` = tfunc(nthreads = 2L), `parglm-LINPACK-3` = tfunc(nthreads = 3L),
#> `parglm-LINPACK-4` = tfunc(nthreads = 4L), `parglm-LINPACK-5` = tfunc(nthreads = 5L),
#> `parglm-LINPACK-6` = tfunc(nthreads = 6L), `parglm-LINPACK-7` = tfunc(nthreads = 7L),
#> `parglm-LINPACK-8` = tfunc(nthreads = 8L), `parglm-LINPACK-9` = tfunc(nthreads = 9L),
#> `parglm-LINPACK-10` = tfunc(nthreads = 10L), `parglm-LINPACK-11` = tfunc(nthreads = 11L),
#> `parglm-LINPACK-12` = tfunc(nthreads = 12L), `parglm-LINPACK-13` = tfunc(nthreads = 13L),
#> `parglm-LINPACK-14` = tfunc(nthreads = 14L), `parglm-LINPACK-15` = tfunc(nthreads = 15L),
#> `parglm-LINPACK-16` = tfunc(nthreads = 16L), `parglm-LINPACK-17` = tfunc(nthreads = 17L),
#> `parglm-LINPACK-18` = tfunc(nthreads = 18L), `parglm-FAST-1` = tfunc(nthreads = 1L,
#> method = "FAST"), `parglm-FAST-2` = tfunc(nthreads = 2L,
#> method = "FAST"), `parglm-FAST-3` = tfunc(nthreads = 3L,
#> method = "FAST"), `parglm-FAST-4` = tfunc(nthreads = 4L,
#> method = "FAST"), `parglm-FAST-5` = tfunc(nthreads = 5L,
#> method = "FAST"), `parglm-FAST-6` = tfunc(nthreads = 6L,
#> method = "FAST"), `parglm-FAST-7` = tfunc(nthreads = 7L,
#> method = "FAST"), `parglm-FAST-8` = tfunc(nthreads = 8L,
#> method = "FAST"), `parglm-FAST-9` = tfunc(nthreads = 9L,
#> method = "FAST"), `parglm-FAST-10` = tfunc(nthreads = 10L,
#> method = "FAST"), `parglm-FAST-11` = tfunc(nthreads = 11L,
#> method = "FAST"), `parglm-FAST-12` = tfunc(nthreads = 12L,
#> method = "FAST"), `parglm-FAST-13` = tfunc(nthreads = 13L,
#> method = "FAST"), `parglm-FAST-14` = tfunc(nthreads = 14L,
#> method = "FAST"), `parglm-FAST-15` = tfunc(nthreads = 15L,
#> method = "FAST"), `parglm-FAST-16` = tfunc(nthreads = 16L,
#> method = "FAST"), `parglm-FAST-17` = tfunc(nthreads = 17L,
#> method = "FAST"), `parglm-FAST-18` = tfunc(nthreads = 18L,
#> method = "FAST"))
out <- eval(cl)s <- summary(out) # result from `microbenchmark`
print(s[, c("expr", "min", "mean", "median", "max")], digits = 3,
row.names = FALSE)
#> expr min mean median max
#> glm 14.87 15.35 15.28 16.20
#> speedglm 4.13 4.57 4.62 4.75
#> parglm-LINPACK-1 13.23 13.42 13.39 13.67
#> parglm-LINPACK-2 7.78 7.93 7.97 8.03
#> parglm-LINPACK-3 5.97 6.12 6.09 6.35
#> parglm-LINPACK-4 5.13 5.29 5.28 5.51
#> parglm-LINPACK-5 4.61 4.80 4.84 4.90
#> parglm-LINPACK-6 4.26 4.44 4.43 4.60
#> parglm-LINPACK-7 4.02 4.26 4.32 4.41
#> parglm-LINPACK-8 3.89 4.04 4.01 4.26
#> parglm-LINPACK-9 3.71 3.92 4.00 4.11
#> parglm-LINPACK-10 3.58 3.77 3.72 3.99
#> parglm-LINPACK-11 3.60 3.70 3.68 3.96
#> parglm-LINPACK-12 3.54 3.68 3.66 3.83
#> parglm-LINPACK-13 3.57 3.73 3.72 3.87
#> parglm-LINPACK-14 3.47 3.57 3.57 3.67
#> parglm-LINPACK-15 3.31 3.64 3.60 4.38
#> parglm-LINPACK-16 3.14 3.57 3.58 3.85
#> parglm-LINPACK-17 3.38 3.60 3.56 3.80
#> parglm-LINPACK-18 3.19 3.48 3.50 3.70
#> parglm-FAST-1 4.03 4.16 4.15 4.30
#> parglm-FAST-2 2.88 3.08 2.97 3.35
#> parglm-FAST-3 2.33 2.62 2.65 2.84
#> parglm-FAST-4 2.32 2.51 2.49 2.67
#> parglm-FAST-5 2.36 2.50 2.51 2.58
#> parglm-FAST-6 2.00 2.32 2.31 2.56
#> parglm-FAST-7 2.19 2.36 2.40 2.44
#> parglm-FAST-8 1.97 2.23 2.26 2.36
#> parglm-FAST-9 1.96 2.06 2.04 2.19
#> parglm-FAST-10 1.94 2.14 2.18 2.34
#> parglm-FAST-11 1.90 2.12 2.08 2.30
#> parglm-FAST-12 1.99 2.13 2.07 2.32
#> parglm-FAST-13 1.92 2.12 2.16 2.24
#> parglm-FAST-14 1.92 2.08 2.04 2.24
#> parglm-FAST-15 1.88 2.07 2.03 2.26
#> parglm-FAST-16 1.88 2.04 2.04 2.19
#> parglm-FAST-17 1.91 2.08 2.08 2.30
#> parglm-FAST-18 1.97 2.15 2.19 2.30The plot below shows median run times versus the number of cores. The dashed line is the median run time of glm and the dotted line is the median run time of speedglm. We could have used glm.fit and parglm.fit. This would make the relative difference bigger as both call e.g., model.matrix and model.frame which do take some time. To show this point, we first compute how much times this takes and then we make the plot. The continuous line is the computation time of model.matrix and model.frame.
modmat_time <- microbenchmark(
modmat_time = {
mf <- model.frame(y ~ ., df); model.matrix(terms(mf), mf)
}, times = 10)
modmat_time # time taken by `model.matrix` and `model.frame`
#> Unit: milliseconds
#> expr min lq mean median uq
#> modmat_time 976.246588 1177.963766 1183.842476 1182.462718 1196.209505
#> max neval
#> 1356.438417 10par(mar = c(4.5, 4.5, .5, .5))
o <- aggregate(time ~ expr, out, median)[, 2] / 10^9
ylim <- range(o, 0); ylim[2] <- ylim[2] + .04 * diff(ylim)
o_linpack <- o[-c(1:2, (n_threads + 3L):length(o))]
o_fast <- o[-(1:(n_threads + 2L))]
plot(1:n_threads, o_linpack, xlab = "Number of cores", yaxs = "i",
ylim = ylim, ylab = "Run time", pch = 16)
points(1:n_threads, o_fast, pch = 1)
abline(h = o[1], lty = 2)
abline(h = o[2], lty = 3)
abline(h = median(modmat_time$time) / 10^9, lty = 1)
The open circles are the FAST method and the other circles are the LINPACK method. Again, the FAST method and the speedglm package compute the Fisher information and then solves the normal equation. This is advantages in terms of computation cost but may lead to unstable solutions. You can alter the number of observations in each parallel chunk with the block_size argument of parglm.control.
The single threaded performance of parglm may be slower when there are more coefficients. The cause seems to be the difference between the LAPACK and LINPACK implementation. This presumably due to either the QR decomposition method and/or the qr.qty method. On Windows, the parglm do seems slower when build with Rtools and the reason seems so be the qr.qty method in LAPACK, dormqr, which is slower then the LINPACK method, dqrsl. Below is an illustration of the computation times on this machine.
qr1 <- qr(X)
qr2 <- qr(X, LAPACK = TRUE)
microbenchmark::microbenchmark(
`qr LINPACK` = qr(X),
`qr LAPACK` = qr(X, LAPACK = TRUE),
`qr.qty LINPACK` = qr.qty(qr1, df$y),
`qr.qty LAPACK` = qr.qty(qr2, df$y),
times = 11)
#> Unit: milliseconds
#> expr min lq mean median
#> qr LINPACK 3110.641142 3115.4115140 3139.2611351 3133.540966
#> qr LAPACK 2493.754534 2497.5766770 2516.8324145 2498.076688
#> qr.qty LINPACK 474.270743 478.8815830 494.5624567 492.309339
#> qr.qty LAPACK 529.304587 530.3935825 530.7427405 530.972688
#> uq max neval
#> 3138.1808365 3263.890541 11
#> 2512.5776900 2644.825526 11
#> 507.2853105 526.827080 11
#> 531.3701745 531.790984 11Session info
sessionInfo()
#> R version 3.4.1 (2017-06-30)
#> Platform: x86_64-redhat-linux-gnu (64-bit)
#> Running under: Amazon Linux AMI 2018.03
#>
#> Matrix products: default
#> BLAS/LAPACK: /usr/lib64/R/lib/libRblas.so
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
#> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
#> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] speedglm_0.3-2 MASS_7.3-47 Matrix_1.2-10
#> [4] microbenchmark_1.4-6 parglm_0.1.0
#>
#> loaded via a namespace (and not attached):
#> [1] Rcpp_0.12.18 codetools_0.2-15 lattice_0.20-35 digest_0.6.15
#> [5] rprojroot_1.3-2 grid_3.4.1 backports_1.1.2 magrittr_1.5
#> [9] evaluate_0.11 stringi_1.2.4 rmarkdown_1.10 tools_3.4.1
#> [13] stringr_1.3.1 yaml_2.2.0 compiler_3.4.1 htmltools_0.3.6
#> [17] knitr_1.20