Merge pull request #56 from michaelweylandt/mw/update_armadillo_call_in_lmBenchmark

Update call to RcppArmadillo_fastLm_impl

Author: eddelbuettel

.travis.yml

@@ -13,9 +13,9 @@ before_install:
   - ./run.sh bootstrap
 
 install:
-  - ./run.sh install_aptget r-cran-rcpp r-cran-matrix r-cran-inline r-cran-runit r-cran-pkgkitten
+  - ./run.sh install_aptget r-cran-rcpp r-cran-matrix r-cran-inline r-cran-runit r-cran-pkgkitten r-cran-microbenchmark
 
-script: 
+script:
   - ./run.sh run_tests
 
 #after_success:
@@ -28,4 +28,3 @@ notifications:
   email:
     on_success: change
     on_failure: change
-

ChangeLog

@@ -1,3 +1,9 @@
+2018-05-30  Michael Weylandt  <michael.weylandt@gmail.com>
+
+	* inst/examples/lmBenchmark.R: Update benchmark script to use
+	microbenchmark and to use exposed fastLm functions from Rcpp
+	packages rather than invoking .Call directly
+
 2018-05-25  Ralf Stubner  <ralf.stubner@daqana.com>
 
 	* inst/include/RcppEigenWrap.h: Use Rf_xlength and R_xlen_t to

DESCRIPTION

@@ -25,6 +25,6 @@ Depends: R (>= 2.15.1)
 LazyLoad: yes
 LinkingTo: Rcpp
 Imports: Matrix (>= 1.1-0), Rcpp (>= 0.11.0), stats, utils
-Suggests: inline, RUnit, pkgKitten
+Suggests: inline, RUnit, pkgKitten, microbenchmark
 URL: http://dirk.eddelbuettel.com/code/rcpp.eigen.html
 BugReports: https://github.com/RcppCore/RcppEigen/issues

inst/examples/lmBenchmark.R

@@ -5,61 +5,78 @@
 ## This file is part of RcppEigen.
 
 require("stats", character=TRUE, quietly=TRUE)
-require("rbenchmark", character=TRUE, quietly=TRUE)
 require("RcppEigen", character=TRUE, quietly=TRUE)
 
-## define different versions of lm
-exprs <- list()
-
-## These versions use rank-revealing decompositions and thus can
-## handle rank-deficient cases.
-
-                                        # default version used in lm()
-exprs$lm.fit <- expression(stats::lm.fit(mm, y))
-                                        # versions from RcppEigen
-## column-pivoted QR decomposition - similar to lm.fit
-exprs$PivQR <- expression(.Call("RcppEigen_fastLm_Impl", mm, y, 0L, PACKAGE="RcppEigen"))
-## LDLt Cholesky decomposition with rank detection
-exprs$LDLt <- expression(.Call("RcppEigen_fastLm_Impl", mm, y, 2L, PACKAGE="RcppEigen"))
-## SVD using the Lapack subroutine dgesdd and Eigen support
-exprs$GESDD <- expression(.Call("RcppEigen_fastLm_Impl", mm, y, 6L, PACKAGE="RcppEigen"))
-## SVD (the JacobiSVD class from Eigen)
-exprs$SVD <- expression(.Call("RcppEigen_fastLm_Impl", mm, y, 4L, PACKAGE="RcppEigen"))
-## eigenvalues and eigenvectors of X'X
-exprs$SymmEig <- expression(.Call("RcppEigen_fastLm_Impl", mm, y, 5L, PACKAGE="RcppEigen"))
-
-## Non-rank-revealing decompositions.  These work fine except when
-## they don't.
-
-## Unpivoted  QR decomposition
-exprs$QR <- expression(.Call("RcppEigen_fastLm_Impl", mm, y, 1L, PACKAGE="RcppEigen"))
-## LLt Cholesky decomposition
-exprs$LLt <- expression(.Call("RcppEigen_fastLm_Impl", mm, y, 3L, PACKAGE="RcppEigen"))
-
-if (suppressMessages(require("RcppArmadillo", character=TRUE, quietly=TRUE))) {
-    exprs$arma <- expression(.Call("RcppArmadillo_fastLm", mm, y, PACKAGE="RcppArmadillo"))
-}
+if(require("microbenchmark", character=TRUE, quietly=TRUE)){
 
-if (suppressMessages(require("RcppGSL", character=TRUE, quietly=TRUE))) {
-    exprs$GSL <- expression(.Call("RcppGSL_fastLm", mm, y, PACKAGE="RcppGSL"))
-}
+    ## define different versions of lm
+    exprs <- list()
 
-do_bench <- function(n=100000L, p=40L, nrep=20L, suppressSVD=(n > 100000L)) {
-    mm <- cbind(1, matrix(rnorm(n * (p - 1L)), nc=p-1L))
-    y <- rnorm(n)
-    if (suppressSVD) exprs <- exprs[!names(exprs) %in% c("SVD", "GSL")]
-    cat("lm benchmark for n = ", n, " and p = ", p, ": nrep = ", nrep, "\n", sep='')
-    do.call(benchmark, c(exprs,
-                         list(order="relative",
-                              columns = c("test", "relative",
-                              "elapsed", "user.self", "sys.self"),
-                              replications = nrep)))
-}
+    ## These versions use rank-revealing decompositions and thus can
+    ## handle rank-deficient cases.
+
+    # default version used in lm()
+    exprs["lm.fit"] <- alist(stats::lm.fit(mm, y))
+
+    # versions from RcppEigen
+    ## column-pivoted QR decomposition - similar to lm.fit
+    exprs["PivQR"] <- alist(RcppEigen::fastLmPure(mm, y, 0L))
+    ## LDLt Cholesky decomposition with rank detection
+    exprs["LDLt"] <- alist(RcppEigen::fastLmPure(mm, y, 2L))
+    ## SVD using the Lapack subroutine dgesdd and Eigen support
+    exprs["GESDD"] <- alist(RcppEigen::fastLmPure(mm, y, 6L))
+    ## SVD (the JacobiSVD class from Eigen)
+    exprs["SVD"] <- alist(RcppEigen::fastLmPure(mm, y, 4L))
+    ## eigenvalues and eigenvectors of X'X
+    exprs["SymmEig"] <- alist(RcppEigen::fastLmPure(mm, y, 5L))
+
+    ## Non-rank-revealing decompositions.  These work fine except when
+    ## they don't.
+
+    ## Unpivoted  QR decomposition
+    exprs["QR"] <- alist(RcppEigen::fastLmPure(mm, y, 1L))
+    ## LLt Cholesky decomposition
+    exprs["LLt"] <- alist(RcppEigen::fastLmPure(mm, y, 3L))
+
+    if (suppressMessages(require("RcppArmadillo", character=TRUE, quietly=TRUE))) {
+        exprs["arma"] <- alist(RcppArmadillo::fastLmPure(mm, y))
+    }
 
-print(do_bench())
+    if (suppressMessages(require("RcppGSL", character=TRUE, quietly=TRUE))) {
+        exprs["GSL"] <- alist(RcppGSL::fastLmPure(mm, y))
+    }
 
-sessionInfo()
+    do_bench <- function(n=100000L, p=40L, nrep=20L, suppressSVD=(n > 100000L)) {
+        mm <- cbind(1, matrix(rnorm(n * (p - 1L)), nc=p-1L))
+        y <- rnorm(n)
+        if (suppressSVD) exprs <- exprs[!names(exprs) %in% c("SVD", "GSL")]
 
-.Call("RcppEigen_eigen_version", FALSE, PACKAGE="RcppEigen")
+        cat("lm benchmark for n = ", n, " and p = ", p, ": nrep = ", nrep, "\n", sep='')
+        cat("RcppEigen: Included Eigen version", paste(RcppEigen:::eigen_version(FALSE), collapse="."), "\n")
+        cat("RcppEigen: Eigen SSE support", RcppEigen:::Eigen_SSE(), "\n")
 
-.Call("RcppEigen_Eigen_SSE", PACKAGE="RcppEigen")
+        mb <- microbenchmark(list=exprs, times = nrep)
+
+        op <- options(microbenchmark.unit="relative")
+        on.exit(options(op))
+
+        mb_relative <- summary(mb)
+        levels(mb_relative$expr) <- names(exprs)
+
+        options(microbenchmark.unit=NULL)
+        mb_absolute <- summary(mb)
+        levels(mb_absolute$expr) <- names(exprs)
+
+        mb_combined <- merge(mb_relative[, c("expr", "median")],
+                             mb_absolute[, c("expr", "median")],
+                             by="expr")
+
+        colnames(mb_combined) <- c("Method",
+                                   "Relative",
+                                   paste0("Elapsed (", attr(mb_absolute, "unit"), ")"))
+
+        mb_combined[order(mb_combined$Relative),]
+    }
+
+    print(do_bench())
+}