kernel reprediction speedup option

Author: CoryMcCartan

NEWS.md

@@ -5,7 +5,7 @@
 * New `b_conv()` for random convolutional features for regression on images TODO im2col()
 * New `b_rocket()` for random convolutional features for regression on time series TODO https://arxiv.org/pdf/2012.08791
 * New `b_echo()` for echo state network features for time series forecasting TODO https://www.ai.rug.nl/minds/uploads/PracticalESN.pdf
-* More efficient `b_ker()` option for many predictions TODO
+* More efficient `b_ker()` option for many predictions
 * New vignette on other packages that help produce basis expansions or embeddings.
 
 # bases 0.1.2

R/b_ker.R

@@ -13,6 +13,10 @@
 #' @inheritParams b_rff
 #' @param x The (training) data points at which to evaluate the kernel. If
 #'   provided, overrides `...`.
+#' @param L_inv The inverse of the Cholesky factor of the kernel matrix at the
+#'   training points. Will be automatically computed if not provided, but in
+#'   order to avoid recomputing it for new predictions, pass `L_inv = TRUE`,
+#'   which will save and re-use this matrix for future calls.
 #'
 #' @returns A matrix of kernel features.
 #'
@@ -39,19 +43,26 @@
 #' @export
 b_ker <- function(..., kernel = k_rbf(),
                   stdize = c("scale", "box", "symbox", "none"),
-                  x = NULL, shift = NULL, scale = NULL) {
+                  x = NULL, shift = NULL, scale = NULL, L_inv = NULL) {
     y = as.matrix(cbind(...))
     std = do_std(y, stdize, shift, scale)
     y = std$x
     if (is.null(x)) x = y
 
-    II = diag(nrow(x))
-    K = kernel(x, x) + 1e-9 * II
-    m = kernel(y, x) %*% backsolve(chol(K), II)
+    save = isTRUE(L_inv)
+    if (is.null(L_inv) || save) {
+        II = diag(nrow(x))
+        K = kernel(x, x) + 1e-9 * II
+        L_inv = backsolve(chol(K), II)
+    }
+    m = kernel(y, x) %*% L_inv
 
     attr(m, "x") = x
     attr(m, "shift") = std$shift
     attr(m, "scale") = std$scale
+    if (save) {
+        attr(m, "L_inv") = L_inv
+    }
     attr(m, "call") = rlang::current_call()
     class(m) = c("b_ker", "matrix", "array")
 
@@ -70,7 +81,7 @@ predict.b_ker <- function (object, newdata, ...)  {
 makepredictcall.b_ker <- function(var, call) {
     if (as.character(call)[1L] == "b_ker" ||
         (is.call(call) && identical(eval(call[[1L]]), b_ker))) {
-        at = attributes(var)[c("x", "shift",  "scale")]
+        at = attributes(var)[c("x", "shift",  "scale", "L_inv")]
         call[names(at)] = at
     }
     call

man/b_ker.Rd

@@ -23,3 +23,21 @@ test_that("predict() method works correctly", {
     pred_m = suppressWarnings(predict(m, newdata=list(x=xn)))
     expect_equal(unname(pred_m[2:21]), unname(fitted(m)[1:20]), tolerance=1e-4)
 })
+
+test_that("predict() works with saved L_inv", {
+    y = c(BJsales)
+    x = seq_along(y)
+    xn = c(0:20, 150:200)
+
+    m0 = lm(y ~ b_ker(x))
+    m1 = lm(y ~ b_ker(x, L_inv=TRUE))
+
+    expect_equal(predict(m1), predict(m0))
+    expect_equal(predict(m1), fitted(m1), tolerance=1e-4)
+    expect_equal(predict(m1), predict(m1, list(x=x)))
+
+    skip_on_cran() # don't run timing tests on CRAN
+    t0 = system.time(for (i in 1:50) predict(m0, list(x = x)))["elapsed"]
+    t1 = system.time(for (i in 1:50) predict(m1, list(x = x)))["elapsed"]
+    expect_lt(t1, t0)
+})