b_tpsob()

Author: CoryMcCartan

.air.toml

@@ -1,4 +1,4 @@
 [format]
-line-width = 85
+line-width = 96
 indent-width = 4
 indent-style = "space"

DESCRIPTION

@@ -1,13 +1,13 @@
 Package: bases
 Title: Basis Expansions for Regression Modeling
 Version: 0.1.2.9999
-Authors@R: 
+Authors@R:
     person("Cory", "McCartan", , "mccartan@psu.edu", role = c("aut", "cre", "cph"),
            comment = c(ORCID = "0000-0002-6251-669X"))
 Description: Provides various basis expansions for flexible regression modeling,
     including random Fourier features (Rahimi & Recht, 2007)
     <https://proceedings.neurips.cc/paper_files/paper/2007/file/013a006f03dbc5392effeb8f18fda755-Paper.pdf>,
-    exact kernel / Gaussian process feature maps, prior features for Bayesian 
+    exact kernel / Gaussian process feature maps, prior features for Bayesian
     Additive Regression Trees (BART) (Chipman et al., 2010) <doi:10.1214/09-AOAS285>,
     and a helpful interface for n-way interactions. The provided functions may
     be used within any modeling formula, allowing the use of kernel methods and
@@ -20,13 +20,14 @@ Depends:
 Imports:
     rlang,
     stats
-Suggests: 
+Suggests:
     recipes,
+    Sieve,
     tibble,
     testthat (>= 3.0.0),
     knitr,
     rmarkdown
-LinkingTo: 
+LinkingTo:
     cpp11
 License: MIT + file LICENSE
 Encoding: UTF-8

NAMESPACE

@@ -9,11 +9,13 @@ S3method(makepredictcall,b_inter)
 S3method(makepredictcall,b_ker)
 S3method(makepredictcall,b_nn)
 S3method(makepredictcall,b_rff)
+S3method(makepredictcall,b_tpsob)
 S3method(predict,b_bart)
 S3method(predict,b_inter)
 S3method(predict,b_ker)
 S3method(predict,b_nn)
 S3method(predict,b_rff)
+S3method(predict,b_tpsob)
 S3method(predict,ridge)
 S3method(print,ridge)
 S3method(print,step_basis)
@@ -25,6 +27,7 @@ export(b_inter)
 export(b_ker)
 export(b_nn)
 export(b_rff)
+export(b_tpsob)
 export(bart_depth_prior)
 export(k_lapl)
 export(k_matern)

NEWS.md

@@ -2,6 +2,8 @@
 
 * `mgcv` smooth interface via `s()` for more flexible penalization TODO
 * New `b_nn()` for neural network basis expansion
+* New `b_tpsob()` for tensor product Sobolev space basis expansion (Zhang and
+     Simon 2023). Requires the `Sieve` package to be installed.
 * 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

R/b_conv.R

@@ -31,4 +31,6 @@ b_conv <- function(
     stride = size,
     activation = c("max", "mean", "ppv"),
     kernels = NULL
-) {}
+) {
+    abort("Not yet implemented.")
+}

R/b_tpsob.R

@@ -8,11 +8,79 @@
 #'      \phi_1(x) = 1 \quad\text{and}\quad
 #'      \phi_j(x) = \sqrt{2}\cos(\pi (j-1) x).
 #' }
+#' The multi-indices \eqn{\mathbf{j}} are generated in a specific order to
+#' maximize statistical efficiency.
+#' All inputs are standardized to lie in the unit hypercube \eqn{[0, 1]^d}.
 #'
+#' @inheritParams b_rff
+#' @param p The number of basis functions to generate.
+#'
+#' @returns A matrix of tensor-product Sobolev space basis features.
+#'
+#' @references
+#' Zhang, T., & Simon, N. (2023). Regression in tensor product spaces by the
+#' method of sieves. _Electronic journal of statistics_, 17(2), 3660.
+#'
+#' @examples
+#' data(quakes)
+#'
+#' m = ridge(depth ~ b_tpsob(lat, long, p = 100), quakes)
+#' plot(fitted(m), quakes$depth)
+#'
+#' x = 1:150
+#' y = as.numeric(BJsales)
+#' m = lm(y ~ b_tpsob(x, p = 10))
+#' plot(x, y)
+#' lines(x, fitted(m), col="blue")
+#' @export
 b_tpsob <- function(
     ...,
     p = 100,
-    stdize = c("scale", "box", "symbox", "none"),
     shift = NULL,
     scale = NULL
-) {}
+) {
+    x = as.matrix(cbind(...))
+    n = nrow(x)
+    d = ncol(x)
+
+    std = do_std(x, "box", shift, scale)
+    x = std$x
+
+    rlang::check_installed("Sieve", "for this basis function.")
+    idx = Sieve::create_index_matrix(d, p + 1L, interaction_order = p)
+    idx = idx[1L + seq_len(p), -1, drop = FALSE]
+
+    m = matrix(nrow = n, ncol = p)
+    for (j in seq_len(p)) {
+        resc = 2^(sum(idx[j, ] > 1L) / 2)
+        m[, j] = resc * row_prod(cos(pi * x * rep(idx[j, ] - 1L, each = n)))
+    }
+
+    attr(m, "shift") = std$shift
+    attr(m, "scale") = std$scale
+    attr(m, "call") = rlang::current_call()
+    class(m) = c("b_tpsob", "matrix", "array")
+
+    m
+}
+
+
+#' @export
+predict.b_tpsob <- function(object, newdata, ...) {
+    if (missing(newdata)) {
+        return(object)
+    }
+    rlang::eval_tidy(makepredictcall(object, attr(object, "call")), newdata)
+}
+
+#' @export
+makepredictcall.b_tpsob <- function(var, call) {
+    if (
+        as.character(call)[1L] == "b_tpsob" ||
+            (is.call(call) && identical(eval(call[[1L]]), b_tpsob))
+    ) {
+        at = attributes(var)[c("shift", "scale")]
+        call[names(at)] = at
+    }
+    call
+}

R/cpp11.R

@@ -12,6 +12,10 @@ dist_l1 <- function(x, y) {
   .Call(`_bases_dist_l1`, x, y)
 }
 
+row_prod <- function(x) {
+  .Call(`_bases_row_prod`, x)
+}
+
 im2col <- function(x, h, w, c, size, stride) {
   .Call(`_bases_im2col`, x, h, w, c, size, stride)
 }

man/b_tpsob.Rd

@@ -26,6 +26,13 @@ extern "C" SEXP _bases_dist_l1(SEXP x, SEXP y) {
     return cpp11::as_sexp(dist_l1(cpp11::as_cpp<cpp11::decay_t<const doubles_matrix<>>>(x), cpp11::as_cpp<cpp11::decay_t<const doubles_matrix<>>>(y)));
   END_CPP11
 }
+// dist.cpp
+doubles row_prod(const doubles_matrix<> x);
+extern "C" SEXP _bases_row_prod(SEXP x) {
+  BEGIN_CPP11
+    return cpp11::as_sexp(row_prod(cpp11::as_cpp<cpp11::decay_t<const doubles_matrix<>>>(x)));
+  END_CPP11
+}
 // im2col.cpp
 doubles_matrix<> im2col(const doubles& x, int h, int w, int c, int size, int stride);
 extern "C" SEXP _bases_im2col(SEXP x, SEXP h, SEXP w, SEXP c, SEXP size, SEXP stride) {
@@ -40,6 +47,7 @@ static const R_CallMethodDef CallEntries[] = {
     {"_bases_dist_l2",    (DL_FUNC) &_bases_dist_l2,    2},
     {"_bases_forest_mat", (DL_FUNC) &_bases_forest_mat, 4},
     {"_bases_im2col",     (DL_FUNC) &_bases_im2col,     6},
+    {"_bases_row_prod",   (DL_FUNC) &_bases_row_prod,   1},
     {NULL, NULL, 0}
 };
 }

src/cpp11.cpp

@@ -51,3 +51,27 @@ doubles_matrix<> dist_l1(const doubles_matrix<> x, const doubles_matrix<> y) {
 
     return out;
 }
+
+
+/*
+ * L1 distance between `x` and `y`, as a matrix with a row for each row of `x`
+ * and a column for each row of `y`.
+ */
+[[cpp11::register]]
+doubles row_prod(const doubles_matrix<> x) {
+    int n = x.nrow();
+    int p = x.ncol();
+    writable::doubles out(n);
+
+    for (int i = 0; i < n; i++) {
+        out[i] = x(i, 0);
+    }
+    for (int j = 1; j < p; j++) {
+        for (int i = 0; i < n; i++) {
+            out[i] *= x(i, j);
+        }
+    }
+
+    return out;
+}
+