Faster GLM submodel fits for non-Gaussian families

NEWS.md

@@ -7,6 +7,10 @@ If you read this from a place other than <https://mc-stan.org/projpred/news/inde
 
 ## Major changes
 
+* Faster projections onto submodels without group-level or additive terms: the projected draws (or clusters) share one model matrix, and for the Gaussian family with the identity link, one linear solve gives the submodel fits for all draws. For example, `varsel()` with forward search for a Gaussian model with 50 observations and 200 predictor terms is about 7 times faster. (GitHub: #574)
+
+* Faster projections onto binomial and Poisson GLM submodels (and Gaussian GLM submodels with a non-identity link): the IRLS fits compute their pseudo-observations in C++, and one C++ call fits all projected draws (or clusters). For example, `cv_varsel()` with `validate_search = TRUE` for binomial and Poisson models is about 4 to 6 times faster. (GitHub: #581)
+
 ## Minor changes
 
 ## Bug fixes and maintenance

R/RcppExports.R

@@ -9,3 +9,11 @@ glm_ridge_c <- function(x, pseudo_obs, lambda, intercept, penalty, beta_init, w_
     .Call('_projpred_glm_ridge_c', PACKAGE = 'projpred', x, pseudo_obs, lambda, intercept, penalty, beta_init, w_init, thresh, qa_updates_max, ls_iter_max, debug)
 }
 
+glm_ridge_native_c <- function(x, y, weights, offset, family_code, link_code, lambda, intercept, penalty, beta_init, w_init, thresh, qa_updates_max, ls_iter_max = 100L, debug = FALSE) {
+    .Call('_projpred_glm_ridge_native_c', PACKAGE = 'projpred', x, y, weights, offset, family_code, link_code, lambda, intercept, penalty, beta_init, w_init, thresh, qa_updates_max, ls_iter_max, debug)
+}
+
+glm_ridge_native_multi_c <- function(x, y, weights, offset, family_code, link_code, lambda, intercept, penalty, beta_init, w_init, thresh, qa_updates_max, ls_iter_max = 100L, debug = FALSE) {
+    .Call('_projpred_glm_ridge_native_multi_c', PACKAGE = 'projpred', x, y, weights, offset, family_code, link_code, lambda, intercept, penalty, beta_init, w_init, thresh, qa_updates_max, ls_iter_max, debug)
+}
+

R/glmfun.R

@@ -207,10 +207,12 @@ glm_ridge <- function(x, y, family = gaussian(), lambda = 0, thresh = 1e-7,
                       qa_updates_max = NULL, weights = NULL, offset = NULL,
                       obsvar = 0, intercept = TRUE, penalty = NULL,
                       normalize = TRUE, beta_init = NULL, beta0_init = NULL,
-                      ls_iter_max = 30) {
+                      ls_iter_max = 30, print_msgs = TRUE) {
   #
   # Fits GLM with ridge penalty on the regression coefficients.
   # Does not handle any dispersion parameters.
+  # The warning messages of the C++ fitter are returned in element `msgs`. With
+  # `print_msgs = TRUE`, they are also printed to the standard output.
   #
   if (is.null(x)) {
     x <- matrix(ncol = 0, nrow = length(y))
@@ -251,22 +253,21 @@ glm_ridge <- function(x, y, family = gaussian(), lambda = 0, thresh = 1e-7,
       # constant predictor)
       x <- matrix(rep(1, length(y)), ncol = 1)
       w0 <- weights
-      pseudo_obs <- function(f, wprev)
-        pseudo_data(f, y, family, offset = offset, weights = weights,
-                    obsvar = obsvar, wprev = wprev)
-      out <- glm_ridge_c(x, pseudo_obs, lambda, FALSE, 0, beta_start, w0,
-                         thresh, qa_updates_max, ls_iter_max)
+      out <- glm_ridge_fit_c(x, y, family, weights, offset, obsvar, lambda,
+                             FALSE, 0, beta_start, w0, thresh, qa_updates_max,
+                             ls_iter_max)
+      print_ridge_msgs(out[[6]], print_msgs)
       return(list(beta = matrix(integer(length = 0)),
                   beta0 = as.vector(out[[1]]), w = out[[3]], loss = out[[4]],
-                  qa_updates = out[[5]]))
+                  qa_updates = out[[5]], msgs = out[[6]]))
     } else {
       # null model with no predictors and no intercept
       pseudo_obs <- function(f, wprev)
         pseudo_data(f, y, family, offset = offset, weights = weights,
                     obsvar = obsvar, wprev = wprev)
       pobs <- pseudo_obs(rep(0, length(y)), weights)
       return(list(beta = matrix(integer(length = 0)), beta0 = 0, w = pobs$wobs,
-                  qa_updates = 0))
+                  qa_updates = 0, msgs = character()))
     }
   }
 
@@ -284,21 +285,91 @@ glm_ridge <- function(x, y, family = gaussian(), lambda = 0, thresh = 1e-7,
 
   # compute the solution
   w0 <- weights
-  pseudo_obs <- function(f, wprev)
-    pseudo_data(f, y, family, offset = offset, weights = weights,
-                obsvar = obsvar, wprev = wprev)
-  out <- glm_ridge_c(x, pseudo_obs, lambda, intercept, penalty, beta_start, w0,
-                     thresh, qa_updates_max, ls_iter_max)
+  out <- glm_ridge_fit_c(x, y, family, weights, offset, obsvar, lambda,
+                         intercept, penalty, beta_start, w0, thresh,
+                         qa_updates_max, ls_iter_max)
   beta <- out[[1]]
   beta0 <- as.vector(out[[2]])
   w <- out[[3]]
   loss <- out[[4]]
+  print_ridge_msgs(out[[6]], print_msgs)
 
   # return the intercept and the coefficients on the original scale
   beta_orig <- beta / transf$scale
   beta0_orig <- beta0 - sum(transf$shift * beta_orig)
 
   out <- nlist(beta = beta_orig, beta0 = beta0_orig, w,
-               qa_updates = out[[5]])
+               qa_updates = out[[5]], msgs = out[[6]])
   return(out)
 }
+
+# Calls the C++ fitter: glm_ridge_native_c() if it supports the family and the
+# link (see native_pseudo_codes()), otherwise glm_ridge_c() with pseudo_data()
+# as the R callback. Both give the same results.
+glm_ridge_fit_c <- function(x, y, family, weights, offset, obsvar, lambda,
+                            intercept, penalty, beta_start, w0, thresh,
+                            qa_updates_max, ls_iter_max) {
+  codes <- native_pseudo_codes(family)
+  if (!is.null(codes)) {
+    return(glm_ridge_native_c(
+      x, as.numeric(y), as.numeric(weights), as.numeric(offset), codes[1],
+      codes[2], lambda, intercept, penalty, beta_start, w0, thresh,
+      qa_updates_max, ls_iter_max
+    ))
+  }
+  pseudo_obs <- function(f, wprev) {
+    pseudo_data(f, y, family, offset = offset, weights = weights,
+                obsvar = obsvar, wprev = wprev)
+  }
+  return(glm_ridge_c(x, pseudo_obs, lambda, intercept, penalty, beta_start,
+                     w0, thresh, qa_updates_max, ls_iter_max))
+}
+
+# Family and link codes for glm_ridge_native_c(), or `NULL` if the
+# pseudo-observations have to be computed by pseudo_data() in R. The native
+# code reproduces the family objects of the stats package, so other family
+# objects with the same names are not covered. Global option
+# `projpred.glm_native` (default `TRUE`) switches the native code off.
+native_pseudo_codes <- function(family) {
+  if (!getOption("projpred.glm_native", TRUE)) {
+    return(NULL)
+  }
+  fam_code <- match(family$family, c("gaussian", "binomial", "poisson"))
+  link_code <- match(family$link, c("identity", "log", "logit", "probit",
+                                    "cloglog", "cauchit"))
+  if (is.na(fam_code) || is.na(link_code)) {
+    return(NULL)
+  }
+  # The function bodies of the stats family, cached per family and link:
+  key <- paste(family$family, family$link)
+  ref_bodies <- native_pseudo_cache[[key]]
+  if (is.null(ref_bodies)) {
+    stats_fam <- get(family$family, envir = asNamespace("stats"),
+                     mode = "function")
+    stats_fam <- tryCatch(stats_fam(link = family$link),
+                          error = function(e) NULL)
+    if (is.null(stats_fam)) {
+      ref_bodies <- FALSE
+    } else {
+      ref_bodies <- lapply(stats_fam[c("linkinv", "mu.eta", "variance")], body)
+    }
+    assign(key, ref_bodies, envir = native_pseudo_cache)
+  }
+  if (isFALSE(ref_bodies) ||
+      !identical(body(family$linkinv), ref_bodies$linkinv) ||
+      !identical(body(family$mu.eta), ref_bodies$mu.eta) ||
+      !identical(body(family$variance), ref_bodies$variance)) {
+    return(NULL)
+  }
+  return(c(fam_code, link_code))
+}
+native_pseudo_cache <- new.env(parent = emptyenv())
+
+# Prints the warning messages of glm_ridge_c() to the standard output, one line
+# per message:
+print_ridge_msgs <- function(msgs, print_msgs) {
+  if (print_msgs && length(msgs) > 0) {
+    cat(paste0(msgs, "\n"), sep = "")
+  }
+  return(invisible())
+}

R/misc.R

@@ -690,6 +690,26 @@ capt_mssgs_warns <- function(expr) {
   return(mssgs_warns_capt)
 }
 
+# Like capt_mssgs_warns(), but collects the messages and warnings by condition
+# handlers. Faster, but does not capture text that is written directly to the
+# message stream, so use it only for fitters that throw R conditions.
+capt_conds <- function(expr) {
+  capt <- character()
+  withCallingHandlers(
+    expr,
+    warning = function(w) {
+      capt <<- c(capt, conditionMessage(w))
+      invokeRestart("muffleWarning")
+    },
+    message = function(m) {
+      capt <<- c(capt, sub("\n$", "", conditionMessage(m)))
+      invokeRestart("muffleMessage")
+    }
+  )
+  capt <- as.character(unlist(strsplit(capt, "\n", fixed = TRUE)))
+  return(setdiff(capt, ""))
+}
+
 # Parse the argument containing the observation weights (`wobs` or `weights`)
 # for the <family_object>$ppd() functions used by proj_predict():
 parse_wobs_ppd <- function(wobs, n_obs) {

src/RcppExports.cpp

@@ -56,10 +56,62 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// glm_ridge_native_c
+List glm_ridge_native_c(arma::mat x, arma::vec y, arma::vec weights, arma::vec offset, int family_code, int link_code, double lambda, bool intercept, arma::vec penalty, arma::vec beta_init, arma::vec w_init, double thresh, int qa_updates_max, int ls_iter_max, bool debug);
+RcppExport SEXP _projpred_glm_ridge_native_c(SEXP xSEXP, SEXP ySEXP, SEXP weightsSEXP, SEXP offsetSEXP, SEXP family_codeSEXP, SEXP link_codeSEXP, SEXP lambdaSEXP, SEXP interceptSEXP, SEXP penaltySEXP, SEXP beta_initSEXP, SEXP w_initSEXP, SEXP threshSEXP, SEXP qa_updates_maxSEXP, SEXP ls_iter_maxSEXP, SEXP debugSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< arma::mat >::type x(xSEXP);
+    Rcpp::traits::input_parameter< arma::vec >::type y(ySEXP);
+    Rcpp::traits::input_parameter< arma::vec >::type weights(weightsSEXP);
+    Rcpp::traits::input_parameter< arma::vec >::type offset(offsetSEXP);
+    Rcpp::traits::input_parameter< int >::type family_code(family_codeSEXP);
+    Rcpp::traits::input_parameter< int >::type link_code(link_codeSEXP);
+    Rcpp::traits::input_parameter< double >::type lambda(lambdaSEXP);
+    Rcpp::traits::input_parameter< bool >::type intercept(interceptSEXP);
+    Rcpp::traits::input_parameter< arma::vec >::type penalty(penaltySEXP);
+    Rcpp::traits::input_parameter< arma::vec >::type beta_init(beta_initSEXP);
+    Rcpp::traits::input_parameter< arma::vec >::type w_init(w_initSEXP);
+    Rcpp::traits::input_parameter< double >::type thresh(threshSEXP);
+    Rcpp::traits::input_parameter< int >::type qa_updates_max(qa_updates_maxSEXP);
+    Rcpp::traits::input_parameter< int >::type ls_iter_max(ls_iter_maxSEXP);
+    Rcpp::traits::input_parameter< bool >::type debug(debugSEXP);
+    rcpp_result_gen = Rcpp::wrap(glm_ridge_native_c(x, y, weights, offset, family_code, link_code, lambda, intercept, penalty, beta_init, w_init, thresh, qa_updates_max, ls_iter_max, debug));
+    return rcpp_result_gen;
+END_RCPP
+}
+// glm_ridge_native_multi_c
+List glm_ridge_native_multi_c(arma::mat x, arma::mat y, arma::vec weights, arma::vec offset, int family_code, int link_code, double lambda, bool intercept, arma::vec penalty, arma::vec beta_init, arma::vec w_init, double thresh, int qa_updates_max, int ls_iter_max, bool debug);
+RcppExport SEXP _projpred_glm_ridge_native_multi_c(SEXP xSEXP, SEXP ySEXP, SEXP weightsSEXP, SEXP offsetSEXP, SEXP family_codeSEXP, SEXP link_codeSEXP, SEXP lambdaSEXP, SEXP interceptSEXP, SEXP penaltySEXP, SEXP beta_initSEXP, SEXP w_initSEXP, SEXP threshSEXP, SEXP qa_updates_maxSEXP, SEXP ls_iter_maxSEXP, SEXP debugSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< arma::mat >::type x(xSEXP);
+    Rcpp::traits::input_parameter< arma::mat >::type y(ySEXP);
+    Rcpp::traits::input_parameter< arma::vec >::type weights(weightsSEXP);
+    Rcpp::traits::input_parameter< arma::vec >::type offset(offsetSEXP);
+    Rcpp::traits::input_parameter< int >::type family_code(family_codeSEXP);
+    Rcpp::traits::input_parameter< int >::type link_code(link_codeSEXP);
+    Rcpp::traits::input_parameter< double >::type lambda(lambdaSEXP);
+    Rcpp::traits::input_parameter< bool >::type intercept(interceptSEXP);
+    Rcpp::traits::input_parameter< arma::vec >::type penalty(penaltySEXP);
+    Rcpp::traits::input_parameter< arma::vec >::type beta_init(beta_initSEXP);
+    Rcpp::traits::input_parameter< arma::vec >::type w_init(w_initSEXP);
+    Rcpp::traits::input_parameter< double >::type thresh(threshSEXP);
+    Rcpp::traits::input_parameter< int >::type qa_updates_max(qa_updates_maxSEXP);
+    Rcpp::traits::input_parameter< int >::type ls_iter_max(ls_iter_maxSEXP);
+    Rcpp::traits::input_parameter< bool >::type debug(debugSEXP);
+    rcpp_result_gen = Rcpp::wrap(glm_ridge_native_multi_c(x, y, weights, offset, family_code, link_code, lambda, intercept, penalty, beta_init, w_init, thresh, qa_updates_max, ls_iter_max, debug));
+    return rcpp_result_gen;
+END_RCPP
+}
 
 static const R_CallMethodDef CallEntries[] = {
     {"_projpred_glm_elnet_c", (DL_FUNC) &_projpred_glm_elnet_c, 14},
     {"_projpred_glm_ridge_c", (DL_FUNC) &_projpred_glm_ridge_c, 11},
+    {"_projpred_glm_ridge_native_c", (DL_FUNC) &_projpred_glm_ridge_native_c, 15},
+    {"_projpred_glm_ridge_native_multi_c", (DL_FUNC) &_projpred_glm_ridge_native_multi_c, 15},
     {NULL, NULL, 0}
 };