New kernelshap solver

Author: mayer79

DESCRIPTION

@@ -28,7 +28,6 @@ RoxygenNote: 7.3.2
 Imports: 
     doFuture,
     foreach,
-    MASS,
     stats,
     utils
 Suggests: 

NEWS.md

@@ -31,17 +31,23 @@ The list is passed to `[foreach::foreach(.options.future = ...)]`.
 ### Internal changes
 
 - Matrices holding on-off vectors are now consistently of type logical ([#167](https://github.com/ModelOriented/kernelshap/pull/167)).
+- `kernelshap()` solver: Replacing the Moore-Penrose pseudo-inverse by two direct solves, a trick of [Ian Covert](https://github.com/iancovert/shapley-regression/blob/master/shapreg/shapley.py),
+  and ported to R in ([#171](https://github.com/ModelOriented/kernelshap/pull/171)).
 
 ### Changes in parallelism
 
 We have switched from `%dopar%` to `doFuture` ([#170](https://github.com/ModelOriented/kernelshap/pull/170)) with the following impact:
 
 - No need for calling `registerDoFuture()` anymore.
 - Random seeding is properly handled, and respects `seed`, thanks [#163](https://github.com/ModelOriented/kernelshap/issues/163) for reporting.
-- {doFuture} is listed under "imports", not as "suggested".
 - If missing packages or globals have to be specified, this now has to be done through `parallel_args = list(packages = ..., globals = ...)` 
 instead of `parallel_args = list(.packages = ..., .globals = ...)`. The list is passed to `[foreach::foreach(.options.future = ...)]`.
 
+### Dependencies
+
+- {MASS}: Dropped from imports
+- {doFuture}: suggests -> imports
+
 # kernelshap 0.8.0
 
 ### Major improvement

R/utils_kernelshap.R

@@ -95,14 +95,15 @@ kernelshap_one <- function(
 
 # Regression coefficients given sum(beta) = constraint
 # A: (p x p), b: (p x k), constraint: (1 x K)
+# Full credits: https://github.com/iancovert/shapley-regression/blob/master/shapreg/shapley.py
 solver <- function(A, b, constraint) {
-  p <- ncol(A)
-  Ainv <- MASS::ginv(A)
-  dimnames(Ainv) <- dimnames(A)
-  s <- (matrix(colSums(Ainv %*% b), nrow = 1L) - constraint) / sum(Ainv) #  (1 x K)
-  Ainv %*% (b - s[rep.int(1L, p), , drop = FALSE]) #  (p x K)
+  Ainv1 <- solve(A, matrix(1, nrow = nrow(A)))
+  Ainvb <- solve(A, b)
+  num <- rbind(colSums(Ainvb)) - constraint
+  return(Ainvb - Ainv1 %*% num / sum(Ainv1))
 }
 
+
 # Draw m binary vectors z of length p with sum(z) distributed according
 # to Kernel SHAP weights -> (m x p) matrix.
 # The argument S can be used to restrict the range of sum(z).

packaging.R

@@ -39,7 +39,6 @@ use_description(
 
 use_package("doFuture", "Imports")
 use_package("foreach", "Imports")
-use_package("MASS", "Imports")
 use_package("stats", "Imports")
 use_package("utils", "Imports")
 

tests/testthat/test-kernelshap-utils.R

@@ -128,3 +128,23 @@ test_that("input_partly_exact(p, deg) fails for bad p or deg", {
   expect_error(input_partly_exact(2L, deg = 0L, feature_names = LETTERS[1:p]))
   expect_error(input_partly_exact(5L, deg = 3L, feature_names = LETTERS[1:p]))
 })
+
+test_that("new solver gives same results as original one", {
+  solver_old <- function(A, b, constraint) {
+    p <- ncol(A)
+    Ainv <- solve(A) # was actually: Ainv <- MASS::ginv(A)
+    dimnames(Ainv) <- dimnames(A)
+    s <- (matrix(colSums(Ainv %*% b), nrow = 1L) - constraint) / sum(Ainv) #  (1 x K)
+    Ainv %*% (b - s[rep.int(1L, p), , drop = FALSE]) #  (p x K)
+  }
+
+  A <- matrix(seq(0.1, 0.20, length.out = 25), ncol = 5)
+  diag(A) <- 0.5
+  b <- cbind(1:5)
+  constraint <- rbind(8)
+  expect_equal(solver_old(A, b, constraint), solver(A, b, constraint))
+
+  b <- cbind(1:5, seq(2, 10, by = 2))
+  constraint <- rbind(1:2)
+  expect_equal(solver_old(A, b, constraint), solver(A, b, constraint))
+})