[avehtari]

This PR was assisted by Claude.

Summary

This PR builds on #581. After that PR, the IRLS fits of a non-Gaussian GLM projection take about 35 % of the time of cv_varsel(validate_search = TRUE). About 40 % goes to R code in proj_to_submodl() and divmin() that builds, for each projection, the submodel formula (update(), flatten_formula()), a data frame with one response column per draw, the draw-wise formulas (validate_response_formula()) and the model matrix (model.frame(), model.matrix()). For one candidate with 5 terms and 20 clusters, these steps took 1.4 ms of 5.3 ms. They depend only on the predictor terms, the number of draws and the reference data, not on the projected draws. A validated LOO search projects onto almost the same term sets in each fold, and the evaluation projects onto the same submodels.

This PR caches these objects during a call of varsel(), cv_varsel() or project():

• proj_to_submodl() keeps a store per reference model (identified by the environment of refmodel$fetch_data) with the key (predictor terms, number of draws). On a miss, it computes the entry with the existing code. On a hit, it calls the multi-draw fitters of divmin() (new divmin_multi(), which divmin() now also uses) directly with the cached model matrix and the projected draws p_ref$mu.
• proj_to_submodl() uses the fast path only for the traditional projection with refmodel$div_minimizer = divmin(), the default projpred.glm_fitter, no NA in p_ref$mu, and no arguments of divmin() itself in .... In all other cases, it runs the existing code.
• The three functions activate the cache for the duration of their call. Nested calls share it (a depth counter), and the outermost call empties it at its exit. No state remains between calls, and the returned objects do not contain the cache. Parallel workers (foreach) have no active cache and use the existing code. The global option projpred.proj_cache = FALSE switches the cache off.
• After a multi-draw fit, divmin_multi() calls check_conv() only if a draw has a message, because check_conv_s() returns TRUE for subfit objects.

Results

The results are the same as without the cache: in the bench below, the rankings, the suggested sizes and the CV ranking proportions are identical, and the largest relative difference in the ELPD estimates is 1.9e-16.

Bench as in #581. Times in seconds. The first speedup compares this PR with #581. The second compares this PR together with #581 and #574 (on which #581 builds) with master.

model	call	master	#581	this PR	speedup vs #581	speedup vs master
collinear (Poisson)	cvt	596.4	117.7	74.7	1.6	8.0
collinear (Poisson)	cvt50	170.6	28.5	22.8	1.3	7.5
sim. Bernoulli, logit	cvt50	728.9	169.2	99.7	1.7	7.3
sim. Bernoulli, probit	cvt50	853.2	219.0	161.5	1.4	5.3
student, binomial	cvt50	716.0	154.3	100.3	1.5	7.1
student, Poisson	cvt50	802.4	159.6	104.2	1.5	7.7
sim. Bernoulli, logit	vs	25.7	6.3	5.5	1.1	4.7
student, Poisson	vs	28.7	5.9	5.7	1.0	5.0

A single varsel() search repeats few term sets, so the cache gives little there: in this bench, varsel() with this PR ran 0.9 to 1.2 times as fast as with #581 alone.

Tests

The test suite passes with NOT_CRAN=true: 263 tests, 55 934 expectations, 0 failures, 0 errors, the same 4 warnings as master. Tested with this change applied to master 73cf0ee (projpred 2.11.0.9000) together with #574 and #581, R 4.6.1, rstanarm 2.32.2, brms 2.23.0, on CSC Roihu (Linux, AMD EPYC).

New test in tests/testthat/test_varsel.R: varsel() with the cache switched off gives the same results as the setup's varsel() calls (cache on) for the GLM setups of the gaussian, Bernoulli, binomial and poisson families; cv_varsel(method = "forward", validate_search = TRUE, nloo = 5) gives the same results with and without the cache for two Bernoulli/binomial/poisson GLM setups; the cache is empty after the calls. The test fails on master (option and functions do not exist; the depth check errors).

Overlap with other PRs

• Leave the left-out observation out of the fold projections in validated LOO CV #580: this PR and Leave the left-out observation out of the fold projections in validated LOO CV #580 change the same lines in two places. The second PR to merge needs a rebase.
  • In the draw-wise loop of divmin(), the resolution is the same as in Faster GLM submodel fits for non-Gaussian families #581: capt_fun(soutdmin <- do.call(sdivmin, c(list(...), dot_args))).
  • In proj_to_submodl(), keep the cached fast path of this PR. In the fallback, multiply projpred_ws_aug by the observation weights as Leave the left-out observation out of the fold projections in validated LOO CV #580 does, and pass projpred_ws_aug instead of p_ref$mu. The fast path does not apply to the augmented-data projection, so the fast path needs no change.
  • The cache does not depend on the observation weights. With Leave the left-out observation out of the fold projections in validated LOO CV #580, the folds of the validated LOO CV have different observation weights and use the same cache.