Merge pull request #450 from fweber144/merge_getters

Reduce peak memory usage during performance evaluation

NEWS.md

@@ -34,6 +34,7 @@ If you read this from a place other than <https://mc-stan.org/projpred/news/inde
     + In case of an **rstanarm** reference model, the default for arguments `weightsnew` and `offsetnew` (see `proj_linpred()`, `proj_predict()`, and `predict.refmodel()`) now causes the original observation weights and offsets to be used (instead of ones and zeros, respectively) if possible. For **brms** reference models, this behavior had already been implemented before.
     + An error is now thrown if a length-zero element `weights` or `offset` is returned by the function supplied to argument `extract_model_data` of `init_refmodel()` (before, a vector of ones or zeros was used silently for the observation weights and offsets, respectively).
 * Added the helper function `force_search_terms()` which allows to construct `search_terms` where certain predictor terms are forced to be included (i.e., they are forced to be selected first) whereas other predictor terms are optional (i.e., they are subject to the variable selection, but only after the inclusion of the "forced" terms). (GitHub: #346)
+* Reduced peak memory usage during performance evaluation (more precisely, during the re-projections done for the performance evaluation). This reduction is considerable especially for multilevel submodels, but possibly also for additive submodels. (GitHub: #440, #450)
 
 ## Bug fixes
 

R/cv_varsel.R

@@ -518,24 +518,23 @@ loo_varsel <- function(refmodel, method, nterms_max, ndraws,
     )
     verb_out("-----", verbose = verbose)
 
-    verb_out("-----\nFor performance evaluation: Re-projecting (using the ",
+    verb_out("-----\nPerformance evaluation, step 1: Re-projecting (using the ",
              "full dataset) onto the submodels along the full-data solution ",
-             "path ...", verbose = verbose && refit_prj)
-    submodls <- get_submodls(
-      search_path = search_path,
-      nterms = c(0, seq_along(search_path$solution_terms)),
-      refmodel = refmodel, regul = opt$regul, refit_prj = refit_prj,
-      ndraws = ndraws_pred, nclusters = nclusters_pred, return_p_ref = TRUE, ...
+             "path and evaluating their predictive performance ...",
+             verbose = verbose && refit_prj)
+    perf_eval_out <- perf_eval(
+      search_path = search_path, refmodel = refmodel, regul = opt$regul,
+      refit_prj = refit_prj, ndraws = ndraws_pred, nclusters = nclusters_pred,
+      return_p_ref = TRUE, return_preds = TRUE, indices_test = inds, ...
     )
-    refdist_eval <- submodls$p_ref
-    submodls <- submodls$submodls
-    clust_used_eval <- element_unq(submodls, nm = "clust_used")
-    nprjdraws_eval <- element_unq(submodls, nm = "nprjdraws")
+    clust_used_eval <- perf_eval_out[["clust_used"]]
+    nprjdraws_eval <- perf_eval_out[["nprjdraws"]]
+    refdist_eval <- perf_eval_out[["p_ref"]]
     verb_out("-----", verbose = verbose && refit_prj)
 
-    verb_out("-----\nCalculating the full-data performance evaluation ",
-             "results and weighting those results according to the PSIS-LOO ",
-             "CV weights ...", verbose = verbose)
+    verb_out("-----\nPerformance evaluation, step 2: Weighting the full-data ",
+             "performance evaluation results according to the PSIS-LOO CV ",
+             "weights ...", verbose = verbose)
     if (refmodel$family$for_latent) {
       refdist_eval_mu_offs_oscale <- refmodel$family$latent_ilink(
         t(refdist_eval$mu_offs), cl_ref = refdist_eval$cl,
@@ -665,13 +664,10 @@ loo_varsel <- function(refmodel, method, nterms_max, ndraws,
     # This re-weighting requires a re-normalization (as.array() is applied to
     # have stricter consistency checks, see `?sweep`):
     lw_sub <- sweep(lw_sub, 2, as.array(apply(lw_sub, 2, log_sum_exp)))
-    for (k in seq_along(submodls)) {
-      mu_k <- refmodel$family$mu_fun(submodls[[k]]$outdmin,
-                                     obs = inds,
-                                     offset = refmodel$offset[inds])
-      log_lik_sub <- t(refmodel$family$ll_fun(
-        mu_k, submodls[[k]]$dis, refmodel$y[inds], refmodel$wobs[inds]
-      ))
+    for (k in seq_len(1 + length(search_path$solution_terms))) {
+      # TODO: For consistency, replace `k` in this `for` loop by `j`.
+      mu_k <- perf_eval_out[["mu_by_size"]][[k]]
+      log_lik_sub <- perf_eval_out[["lppd_by_size"]][[k]]
       loo_sub[[k]][inds] <- apply(log_lik_sub + lw_sub, 2, log_sum_exp)
       if (refmodel$family$for_latent) {
         mu_k_oscale <- refmodel$family$latent_ilink(
@@ -750,26 +746,19 @@ loo_varsel <- function(refmodel, method, nterms_max, ndraws,
         verbose = verbose_search, opt = opt, search_terms = search_terms, ...
       )
 
-      # Re-project along the solution path (or fetch the projections from the
-      # search results) of the current fold:
-      submodls <- get_submodls(
-        search_path = search_path,
-        nterms = c(0, seq_along(search_path$solution_terms)),
-        refmodel = refmodel, regul = opt$regul, refit_prj = refit_prj,
-        ndraws = ndraws_pred, nclusters = nclusters_pred,
+      # Run the performance evaluation for the submodels along the predictor
+      # ranking:
+      perf_eval_out <- perf_eval(
+        search_path = search_path, refmodel = refmodel, regul = opt$regul,
+        refit_prj = refit_prj, ndraws = ndraws_pred, nclusters = nclusters_pred,
         reweighting_args = list(cl_ref = cl_pred, wdraws_ref = exp(lw[, i])),
-        ...
+        indices_test = i, ...
       )
-      clust_used_eval <- element_unq(submodls, nm = "clust_used")
-      nprjdraws_eval <- element_unq(submodls, nm = "nprjdraws")
-
-      # Predictive performance at the omitted observation:
-      summaries_sub <- get_sub_summaries(submodls = submodls,
-                                         refmodel = refmodel,
-                                         test_points = i)
 
       return(nlist(predictor_ranking = search_path[["solution_terms"]],
-                   summaries_sub, clust_used_eval, nprjdraws_eval))
+                   summaries_sub = perf_eval_out[["sub_summaries"]],
+                   clust_used_eval = perf_eval_out[["clust_used"]],
+                   nprjdraws_eval = perf_eval_out[["nprjdraws"]]))
     }
     if (!parallel) {
       # Sequential case. Actually, we could simply use ``%do_projpred%` <-
@@ -981,7 +970,7 @@ loo_varsel <- function(refmodel, method, nterms_max, ndraws,
   summaries <- list(sub = summ_sub, ref = summ_ref)
 
   if (!validate_search) {
-    out_list <- nlist(search_path, ce = sapply(submodls, "[[", "ce"))
+    out_list <- nlist(search_path, ce = perf_eval_out[["ce"]])
   } else {
     out_list <- nlist(solution_terms_cv = solution_terms_mat[
       , seq_len(prv_len_soltrms), drop = FALSE
@@ -1034,23 +1023,13 @@ kfold_varsel <- function(refmodel, method, nterms_max, ndraws,
       verbose = verbose_search, opt = opt, search_terms = search_terms, ...
     )
 
-    # For performance evaluation: Re-project (using the training data of the
-    # current fold) along the predictor ranking (or fetch the projections from
-    # the search output) of the current fold:
-    submodls <- get_submodls(
-      search_path = search_path,
-      nterms = c(0, seq_along(search_path$solution_terms)),
-      refmodel = fold$refmodel, regul = opt$regul, refit_prj = refit_prj,
-      ndraws = ndraws_pred, nclusters = nclusters_pred, ...
+    # Run the performance evaluation for the submodels along the predictor
+    # ranking:
+    perf_eval_out <- perf_eval(
+      search_path = search_path, refmodel = fold$refmodel, regul = opt$regul,
+      refit_prj = refit_prj, ndraws = ndraws_pred, nclusters = nclusters_pred,
+      refmodel_fulldata = refmodel, indices_test = fold$omitted, ...
     )
-    clust_used_eval <- element_unq(submodls, nm = "clust_used")
-    nprjdraws_eval <- element_unq(submodls, nm = "nprjdraws")
-
-    # Performance evaluation for the re-projected or fetched submodels of the
-    # current fold:
-    summaries_sub <- get_sub_summaries(submodls = submodls,
-                                       refmodel = refmodel,
-                                       test_points = fold$omitted)
 
     # Performance evaluation for the reference model of the current fold:
     eta_test <- fold$refmodel$ref_predfun(
@@ -1069,7 +1048,9 @@ kfold_varsel <- function(refmodel, method, nterms_max, ndraws,
     )
 
     return(nlist(predictor_ranking = search_path[["solution_terms"]],
-                 summaries_sub, summaries_ref, clust_used_eval, nprjdraws_eval))
+                 summaries_sub = perf_eval_out[["sub_summaries"]],
+                 summaries_ref, clust_used_eval = perf_eval_out[["clust_used"]],
+                 nprjdraws_eval = perf_eval_out[["nprjdraws"]]))
   }
   if (!parallel) {
     # Sequential case. Actually, we could simply use ``%do_projpred%` <-

R/project.R

@@ -262,14 +262,15 @@ project <- function(object, nterms = NULL, solution_terms = NULL,
   }
 
   ## project onto the submodels
-  submodls <- get_submodls(
+  submodls <- perf_eval(
     search_path = nlist(
       solution_terms,
       p_sel = object$search_path$p_sel,
       outdmins = object$search_path$outdmins
     ),
     nterms = nterms, refmodel = refmodel, regul = regul, refit_prj = refit_prj,
-    ndraws = ndraws, nclusters = nclusters, projpred_verbose = verbose, ...
+    ndraws = ndraws, nclusters = nclusters, return_submodls = TRUE,
+    projpred_verbose = verbose, ...
   )
 
   # Output:

R/projfun.R

@@ -45,23 +45,33 @@ get_submodl_prj <- function(solution_terms, p_ref, refmodel, regul = 1e-4,
   ))
 }
 
-# Function to fetch init_submodl() output (of class `submodl`) for each of given
-# model sizes `nterms`, so this gives a list of objects, each of class
-# `submodl`.
-get_submodls <- function(search_path, nterms, refmodel, regul,
-                         refit_prj = FALSE, ndraws, nclusters,
-                         reweighting_args = NULL, return_p_ref = FALSE, ...) {
+# Function to prepare the performance evaluation. For each submodel size along
+# the predictor ranking, it first fetches the init_submodl() output (of class
+# `submodl`) and then calculates the submodel "summary" (precursor quantities)
+# for the actual performance evaluation performed by summary.vsel() later.
+perf_eval <- function(search_path,
+                      nterms = c(0, seq_along(search_path$solution_terms)),
+                      refmodel, regul, refit_prj = FALSE, ndraws, nclusters,
+                      reweighting_args = NULL, return_submodls = FALSE,
+                      return_preds = FALSE, return_p_ref = FALSE,
+                      refmodel_fulldata = refmodel,
+                      indices_test, newdata_test = NULL,
+                      offset_test = refmodel_fulldata$offset[indices_test],
+                      wobs_test = refmodel_fulldata$wobs[indices_test],
+                      y_test = refmodel_fulldata$y[indices_test],
+                      y_oscale_test = refmodel_fulldata$y_oscale[indices_test],
+                      ...) {
   if (!refit_prj) {
     p_ref <- search_path$p_sel
     # In this case, simply fetch the already computed projections, so don't
     # project again.
-    fetch_submodl <- function(nterms, ...) {
+    fetch_submodl <- function(size_j, ...) {
       return(init_submodl(
         # Re-use the submodel fits from the search:
-        outdmin = search_path$outdmins[[nterms + 1]],
+        outdmin = search_path$outdmins[[size_j + 1]],
         p_ref = p_ref,
         refmodel = refmodel,
-        solution_terms = utils::head(search_path$solution_terms, nterms),
+        solution_terms = utils::head(search_path$solution_terms, size_j),
         wobs = refmodel$wobs
       ))
     }
@@ -77,16 +87,64 @@ get_submodls <- function(search_path, nterms, refmodel, regul,
         wdraws = reweighting_args$wdraws_ref, cl = reweighting_args$cl_ref
       )
     }
-    fetch_submodl <- function(nterms, ...) {
+    fetch_submodl <- function(size_j, ...) {
       return(get_submodl_prj(
-        solution_terms = utils::head(search_path$solution_terms, nterms),
+        solution_terms = utils::head(search_path$solution_terms, size_j),
         p_ref = p_ref, refmodel = refmodel, regul = regul, ...
       ))
     }
   }
-  out <- lapply(nterms, fetch_submodl, ...)
+  out_by_size <- lapply(nterms, function(size_j) {
+    # Fetch the init_submodl() output (of class `submodl`) for the submodel at
+    # position `size_j + 1` of the predictor ranking:
+    submodl <- fetch_submodl(size_j, ...)
+    if (return_submodls) {
+      # Currently only called in project().
+      return(submodl)
+    }
+    if (return_preds) {
+      # Currently only called in loo_varsel()'s `validate_search = FALSE` case.
+      mu_j <- refmodel$family$mu_fun(submodl$outdmin, obs = indices_test,
+                                     offset = refmodel$offset[indices_test])
+      lppd_j <- t(refmodel$family$ll_fun(
+        mu_j, submodl$dis, refmodel$y[indices_test], refmodel$wobs[indices_test]
+      ))
+      out_j <- nlist(mu_j, lppd_j)
+    } else {
+      # Calculate precursor quantities for predictive performance statistic(s)
+      # of the submodel at position `size_j + 1` of the predictor ranking:
+      sub_summary <- weighted_summary_means(
+        y_wobs_test = data.frame(y = y_test, y_oscale = y_oscale_test,
+                                 wobs = wobs_test),
+        family = refmodel_fulldata$family,
+        wdraws = submodl$wdraws_prj,
+        mu = refmodel_fulldata$family$mu_fun(submodl$outdmin,
+                                             obs = indices_test,
+                                             newdata = newdata_test,
+                                             offset = offset_test),
+        dis = submodl$dis,
+        cl_ref = submodl$cl_ref,
+        wdraws_ref = submodl$wdraws_ref
+      )
+      out_j <- nlist(sub_summary)
+    }
+    return(c(out_j, list(ce = submodl[["ce"]])))
+  })
+  if (return_submodls) {
+    return(out_by_size)
+  }
+  if (return_preds) {
+    out <- list(mu_by_size = lapply(out_by_size, "[[", "mu_j"),
+                lppd_by_size = lapply(out_by_size, "[[", "lppd_j"))
+  } else {
+    out <- list(sub_summaries = lapply(out_by_size, "[[", "sub_summary"))
+  }
+  out <- c(out, list(ce = sapply(out_by_size, "[[", "ce"),
+                     clust_used = p_ref$clust_used,
+                     nprjdraws = NCOL(p_ref$mu)))
   if (return_p_ref) {
-    out <- list(submodls = out, p_ref = p_ref)
+    # Currently only called in loo_varsel()'s `validate_search = FALSE` case.
+    out <- c(out, nlist(p_ref))
   }
   return(out)
 }

R/search.R

@@ -49,7 +49,7 @@ search_forward <- function(p_ref, refmodel, nterms_max, verbose = TRUE, opt,
   # `chosen`. However, `reduce_models(chosen)` and `chosen` should be identical
   # at this place because select_possible_terms_size() already avoids redundant
   # models. Thus, use `chosen` here because it matches `outdmins` (this
-  # matching is necessary because later in get_submodls()'s `!refit_prj` case,
+  # matching is necessary because later in perf_eval()'s `!refit_prj` case,
   # `outdmins` is indexed with integers which are based on `solution_terms`):
   return(nlist(solution_terms = setdiff(chosen, "1"), outdmins))
 }

R/summary_funs.R

@@ -1,22 +1,3 @@
-get_sub_summaries <- function(submodls, refmodel, test_points, newdata = NULL,
-                              offset = refmodel$offset[test_points],
-                              wobs = refmodel$wobs[test_points],
-                              y = refmodel$y[test_points],
-                              y_oscale = refmodel$y_oscale[test_points]) {
-  lapply(submodls, function(submodl) {
-    weighted_summary_means(
-      y_wobs_test = data.frame(y, y_oscale, wobs),
-      family = refmodel$family,
-      wdraws = submodl$wdraws_prj,
-      mu = refmodel$family$mu_fun(submodl$outdmin, obs = test_points,
-                                  newdata = newdata, offset = offset),
-      dis = submodl$dis,
-      cl_ref = submodl$cl_ref,
-      wdraws_ref = submodl$wdraws_ref
-    )
-  })
-}
-
 # Calculate log posterior(-projection) predictive density values and average
 # them across parameter draws (together with the corresponding expected response
 # values).

R/varsel.R

@@ -286,30 +286,19 @@ varsel.refmodel <- function(object, d_test = NULL, method = NULL,
   )
   verb_out("-----", verbose = verbose)
 
-  # For the performance evaluation: Re-project along the solution path (or fetch
-  # the projections from the search results):
-  verb_out("-----\nFor performance evaluation: Re-projecting onto the ",
-           "submodels along the solution path ...",
+  # "Run" the performance evaluation for the submodels along the predictor
+  # ranking (in fact, we only prepare the performance evaluation by computing
+  # precursor quantities, but for users, this difference is not perceivable):
+  verb_out("-----\nRunning the performance evaluation ...",
            verbose = verbose && refit_prj)
-  submodls <- get_submodls(
-    search_path = search_path,
-    nterms = c(0, seq_along(search_path$solution_terms)),
-    refmodel = refmodel, regul = regul, refit_prj = refit_prj,
-    ndraws = ndraws_pred, nclusters = nclusters_pred, ...
+  perf_eval_out <- perf_eval(
+    search_path = search_path, refmodel = refmodel, regul = regul,
+    refit_prj = refit_prj, ndraws = ndraws_pred, nclusters = nclusters_pred,
+    indices_test = NULL, newdata_test = d_test$data,
+    offset_test = d_test$offset, wobs_test = d_test$weights, y_test = d_test$y,
+    y_oscale_test = d_test$y_oscale, ...
   )
-  clust_used_eval <- element_unq(submodls, nm = "clust_used")
-  nprjdraws_eval <- element_unq(submodls, nm = "nprjdraws")
   verb_out("-----", verbose = verbose && refit_prj)
-  # The performance evaluation itself, i.e., the calculation of the predictive
-  # performance statistic(s) for the submodels along the solution path:
-  sub <- get_sub_summaries(submodls = submodls,
-                           refmodel = refmodel,
-                           test_points = NULL,
-                           newdata = d_test$data,
-                           offset = d_test$offset,
-                           wobs = d_test$weights,
-                           y = d_test$y,
-                           y_oscale = d_test$y_oscale)
 
   # Predictive performance of the reference model:
   if (inherits(refmodel, "datafit")) {
@@ -369,21 +358,21 @@ varsel.refmodel <- function(object, d_test = NULL, method = NULL,
               search_path,
               solution_terms = search_path$solution_terms,
               solution_terms_cv = NULL,
-              ce = sapply(submodls, "[[", "ce"),
+              ce = perf_eval_out[["ce"]],
               type_test = d_test$type,
               y_wobs_test,
               nobs_test,
-              summaries = nlist(sub, ref),
+              summaries = nlist(sub = perf_eval_out[["sub_summaries"]], ref),
               nterms_all,
               nterms_max,
               method,
               cv_method = NULL,
               K = NULL,
               validate_search = NULL,
               clust_used_search = search_path$p_sel$clust_used,
-              clust_used_eval,
+              clust_used_eval = perf_eval_out[["clust_used"]],
               nprjdraws_search = NCOL(search_path$p_sel$mu),
-              nprjdraws_eval,
+              nprjdraws_eval = perf_eval_out[["nprjdraws"]],
               projpred_version = utils::packageVersion("projpred"))
   class(vs) <- "vsel"
 

tests/testthat/test_varsel.R

@@ -622,9 +622,9 @@ if (run_more) {
 
 # In fact, `regul` is already checked in `test_project.R`, so the `regul` tests
 # could be omitted here since varsel() and cv_varsel() also pass `regul` to
-# get_submodl_prj() (usually via get_submodls(), just like project()). This
-# doesn't hold for L1 search, though. So for L1 search, the `regul` tests are
-# still needed.
+# get_submodl_prj() (usually via perf_eval(), just like project()). This doesn't
+# hold for L1 search, though. So for L1 search, the `regul` tests are still
+# needed.
 
 test_that(paste(
   "for GLMs with L1 search, `regul` only has an effect on prediction, not on",