Merge pull request #123 from ModelOriented/performance

Performance improvements

Author: mayer79

NEWS.md

@@ -1,9 +1,10 @@
 # kernelshap 0.4.1
 
-## Other changes
+## Performance improvements
 
-- Slight speed-up of `kernelshap()` and `permshap()` for single-output predictions.
-- Slight speed-up of `kernelshap()` and `permshap()` for factor-valued predictions.
+- Significant speed-up for data objects with *single class "data.frame"*, i.e., no data.tables or tibbles or grouped data etc. This change makes it almost as fast to work with data.frames as with matrices.
+- Slight speed-up for single-output predictions.
+- Slight speed-up for factor-valued predictions.
 - Slight speed-up of `permshap()` by caching calculations for the two special permutations of all 0 and all 1. Consequently, the `m_exact` component in the output is reduced by 2.
 
 # kernelshap 0.4.0

R/kernelshap.R

@@ -233,14 +233,14 @@ kernelshap.default <- function(object, X, bg_X, pred_fun = stats::predict,
     }
     m_exact <- nrow(precalc[["Z"]])
     prop_exact <- sum(precalc[["w"]])
-    precalc[["bg_X_exact"]] <- bg_X[rep(seq_len(bg_n), times = m_exact), , drop = FALSE]
+    precalc[["bg_X_exact"]] <- rep_rows(bg_X, rep.int(seq_len(bg_n), m_exact))
   } else {
     precalc <- list()
     m_exact <- 0L
     prop_exact <- 0
   }
   if (!exact) {
-    precalc[["bg_X_m"]] <- bg_X[rep(seq_len(bg_n), times = m), , drop = FALSE]  
+    precalc[["bg_X_m"]] <- rep_rows(bg_X, rep.int(seq_len(bg_n), m))
   }
 
   # Some infos

R/permshap.R

@@ -90,7 +90,7 @@ permshap.default <- function(object, X, bg_X, pred_fun = stats::predict,
   precalc <- list(
     Z = Z,
     Z_code = rowpaste(Z), 
-    bg_X_rep = bg_X[rep(seq_len(bg_n), times = m_exact), , drop = FALSE]
+    bg_X_rep = rep_rows(bg_X, rep.int(seq_len(bg_n), m_exact))
   )
 
   if (m_exact * bg_n > 2e5) {

R/utils.R

@@ -1,3 +1,26 @@
+#' Fast Row Subsetting
+#' 
+#' Internal function used to row-subset data.frames.
+#' Brings a massive speed-up for data.frames. All other classes (tibble, data.table,
+#' matrix) are subsetted in the usual way.
+#' 
+#' @noRd
+#' @keywords internal
+#' 
+#' @param x A matrix-like object.
+#' @param i Logical or integer vector of rows to pick.
+#' @returns Subsetted version of `x`.
+rep_rows <- function(x, i) {
+  if (!(all(class(x) == "data.frame"))) {
+    return(x[i, , drop = FALSE])  # matrix, tibble, data.table, ...
+  }
+  # data.frame
+  out <- lapply(x, function(z) if (length(dim(z)) != 2L) z[i] else z[i, , drop = FALSE])
+  attr(out, "row.names") <- .set_row_names(length(i))
+  class(out) <- "data.frame"
+  out
+}
+
 #' Weighted Version of colMeans()
 #' 
 #' Internal function used to calculate column-wise weighted means.

R/utils_kernelshap.R

@@ -11,11 +11,11 @@ kernelshap_one <- function(x, v1, object, pred_fun, feature_names, bg_w, exact,
     bg_X_exact <- precalc[["bg_X_exact"]]                         #  (m_ex*n_bg x p)
     Z <- precalc[["Z"]]                                           #  (m_ex x p)
     m_exact <- nrow(Z)
-    v0_m_exact <- v0[rep(1L, m_exact), , drop = FALSE]            #  (m_ex x K)
-    
+    v0_m_exact <- v0[rep.int(1L, m_exact), , drop = FALSE]        #  (m_ex x K)
+
     # Most expensive part
     vz <- get_vz(                                                 #  (m_ex x K)
-      X = x[rep.int(1L, times = nrow(bg_X_exact)), , drop = FALSE],# (m_ex*n_bg x p)
+      X = rep_rows(x, rep.int(1L, nrow(bg_X_exact))),             #  (m_ex*n_bg x p)
       bg = bg_X_exact,                                            #  (m_ex*n_bg x p)
       Z = Z,                                                      #  (m_ex x p)
       object = object, 
@@ -35,8 +35,8 @@ kernelshap_one <- function(x, v1, object, pred_fun, feature_names, bg_w, exact,
 
   # Iterative sampling part, always using A_exact and b_exact to fill up the weights
   bg_X_m <- precalc[["bg_X_m"]]                                   #  (m*n_bg x p)
-  X <- x[rep(1L, times = nrow(bg_X_m)), , drop = FALSE]           #  (m*n_bg x p)
-  v0_m <- v0[rep(1L, m), , drop = FALSE]                          #  (m x K)
+  X <- rep_rows(x, rep.int(1L, nrow(bg_X_m)))                     #  (m*n_bg x p)
+  v0_m <- v0[rep.int(1L, m), , drop = FALSE]                      #  (m x K)
 
   est_m = list()
   converged <- FALSE
@@ -91,7 +91,7 @@ solver <- function(A, b, constraint) {
   Ainv <- ginv(A)
   dimnames(Ainv) <- dimnames(A)
   s <- (matrix(colSums(Ainv %*% b), nrow = 1L) - constraint) / sum(Ainv)     #  (1 x K)
-  Ainv %*% (b - s[rep(1L, p), , drop = FALSE])                               #  (p x K)
+  Ainv %*% (b - s[rep.int(1L, p), , drop = FALSE])                           #  (p x K)
 }
 
 ginv <- function (X, tol = sqrt(.Machine$double.eps)) {
@@ -138,7 +138,7 @@ sample_Z <- function(p, m, feature_names, S = 1:(p - 1L)) {
   # t(out)
 
   # Vectorized by Mathias Ambuehl
-  out <- rep(rep(0:1, m), as.vector(rbind(p - N, N)))
+  out <- rep(rep.int(0:1, m), as.vector(rbind(p - N, N)))
   dim(out) <- c(p, m)
   ord <- order(col(out), sample.int(m * p))
   out[] <- out[ord]
@@ -180,7 +180,7 @@ input_sampling <- function(p, m, deg, paired, feature_names) {
   }
   w_total <- if (deg == 0L) 1 else 1 - 2 * sum(kernel_weights(p)[seq_len(deg)])
   w <- w_total / m
-  list(Z = Z, w = rep(w, m), A = crossprod(Z) * w)
+  list(Z = Z, w = rep.int(w, m), A = crossprod(Z) * w)
 }
 
 # Functions required only for handling (partly) exact cases
@@ -262,7 +262,7 @@ input_partly_exact <- function(p, deg, feature_names) {
     Z[[k]] <- partly_exact_Z(p, k = k, feature_names = feature_names)
     n <- nrow(Z[[k]])
     w_tot <- kw[k] * (2 - (p == 2L * k))
-    w[[k]] <- rep(w_tot / n, n)
+    w[[k]] <- rep.int(w_tot / n, n)
   }
   w <- unlist(w, recursive = FALSE, use.names = FALSE)
   Z <- do.call(rbind, Z)

R/utils_permshap.R

@@ -28,7 +28,7 @@ shapley_weights <- function(p, ell) {
 permshap_one <- function(x, v1, object, pred_fun, bg_w, v0, precalc, ...) {
   Z <- precalc[["Z"]]                                                    # ((m_ex+2) x K)
   vz <- get_vz(                                                          # (m_ex x K)
-    X = x[rep.int(1L, times = nrow(precalc[["bg_X_rep"]])), , drop = FALSE], # (m_ex*n_bg x p)
+    X = rep_rows(x, rep.int(1L, times = nrow(precalc[["bg_X_rep"]]))),   # (m_ex*n_bg x p)
     bg = precalc[["bg_X_rep"]],                                          # (m_ex*n_bg x p)
     Z = Z[2:(nrow(Z) - 1L), , drop = FALSE],                             # (m_ex x p)
     object = object,

README.md

@@ -206,7 +206,7 @@ sv_dependence(sv_nn, "clarity")
 
 ## Parallel computing
 
-Parallel computing is supported via `foreach`, at the price of losing the progress bar. Note that this does not work with Keras models (and some others).
+Parallel computing is supported via {foreach}, at the price of losing the progress bar. Note that this does not work with Keras models (and some others).
 
 ### Example: Linear regression continued
 
@@ -226,7 +226,7 @@ system.time(
 
 ### Example: Parallel GAM
 
-On Windows, sometimes not all packages or global objects are passed to the parallel sessions. In this case, the necessary instructions to `foreach` can be specified through a named list via `parallel_args`, see the following example:
+On Windows, sometimes not all packages or global objects are passed to the parallel sessions. In this case, the necessary instructions to {foreach} can be specified through a named list via `parallel_args`, see the following example:
 
 ```r
 library(mgcv)

backlog/2023-11-11 Permutation-SHAP.R

@@ -3,19 +3,21 @@ library(ranger)
 
 differences <- numeric(4)
 
+system.time({
 set.seed(1)
-
 for (depth in 1:4) {
+  print(depth)
   fit <- ranger(
     Sepal.Length ~ ., 
     mtry = 3,
     data = iris, 
     max.depth = depth
   )
-  ps <- permshap(fit, iris[2:5], bg_X = iris)
-  ks <- kernelshap(fit, iris[2:5], bg_X = iris)
+  ps <- permshap(fit, iris[2:5], bg_X = iris, verbose = FALSE)
+  ks <- kernelshap(fit, iris[2:5], bg_X = iris, verbose = FALSE)
   differences[depth] <- mean(abs(ks$S - ps$S))
 }
+})
 
 differences  # for tree depth 1, 2, 3, 4
 # 5.053249e-17 9.046443e-17 2.387905e-04 4.403375e-04

tests/testthat/test-utils.R

@@ -32,6 +32,26 @@ test_that("exact_Z() works for both kernel- and permshap", {
   expect_equal(dim(r1), c(2^p, p))
 })
 
+test_that("rep_rows() gives the same as usual subsetting (except rownames)", {
+  setrn <- function(x) {rownames(x) <- 1:nrow(x); x}
+  
+  expect_equal(rep_rows(iris, 1), iris[1, ])
+  expect_equal(rep_rows(iris, 2:1), setrn(iris[2:1, ]))
+  expect_equal(rep_rows(iris, c(1, 1, 1)), setrn(iris[c(1, 1, 1), ]))
+  
+  ir <- iris[1, ]
+  ir$y <- list(list(a = 1, b = 2))
+  expect_equal(rep_rows(ir, c(1, 1)), setrn(ir[c(1, 1), ]))
+})
+
+test_that("rep_rows() gives the same as usual subsetting for matrices", {
+  ir <- data.matrix(iris[1:4])
+  
+  expect_equal(rep_rows(ir, c(1, 1, 2)), ir[c(1, 1, 2), ])
+  expect_equal(rep_rows(ir, 1), ir[1, , drop = FALSE])
+})
+
+
 # Unit tests copied from {hstats}
 
 test_that("rep_each() works", {