[avehtari]

This PR speeds up cv_varsel() with validate_search = TRUE for Gaussian submodels with the identity link, also with the latent projection. It changes the evaluation and search code that #579 added.

Evaluation

For each left-out observation, perf_eval() calls get_p_clust(). get_p_clust() builds observations x clusters matrices of means and predictive variances, but perf_eval_gauss() uses only weighted column sums of them. The evaluation uses ndraws_pred draws of the reference model (default 400), and each of them is its own cluster. (The search uses the nclusters clusters, which are not related.) This PR uses that case:

• The new get_p_clust(as_moments = TRUE) returns the draw index, the weight and the predictive variance of each cluster, and no matrices.
• The new C++ function wcol_moments_c() computes the weighted means and sums of squares from the draws.
• The observation weights now multiply the model matrix before crossprod(), not the observations x draws matrix.
• If perf_eval_gauss() does not apply or get_p_clust(as_moments = TRUE) returns NULL, perf_eval() computes the full get_p_clust() result. NULL occurs, for example, if a cluster has more than one draw. If a projected draw is not finite, perf_eval_gauss() gets the full get_p_clust() result.

Without the observations x clusters and observations x draws matrices, garbage collection takes less time (see "Timings").

Search

Each step got its candidates from select_possible_terms_size() and their new labels from setdiff() on each candidate string. If each search term is "1" or has one term label, the new single_label_info() (once per search) and cands_single_label() (once per step) give the same candidates and labels with vector operations. Other search terms use the old code. labels_gauss() now uses the session cache of #599.

Timings

The timings are for full validated LOO. "Before" is master with #601 merged. For each case, both versions ran in one SLURM job, one run each. Case 1 ran on CSC Roihu with R 4.6.1, case 2 on Triton with R 4.5.3.

case	model	search clusters	evaluation draws	before (s)	after (s)	speedup
1	Gaussian, N 50, 200 predictors	20	400	32.3	6.2	5.2
2	Gaussian, N 382, 26 predictors	400	4000	54.8	32.2	1.7

The search clusters are nclusters, the evaluation draws are ndraws_pred. Garbage collection took 1.7 s before and 0.5 s after in case 1, and 10.2 s before and 4.7 s after in case 2. The predictor ranking is identical in both cases. In case 2, the other elements of the vsel object are identical or differ by at most 1.1e-15 (relative).

Tests

• test_wcol_moments.R:
  • It compares wcol_moments_c() with R code, including NaN and infinite values and a weight of 0.
  • It compares get_p_clust(as_moments = TRUE) with the full get_p_clust() for all draws, permuted draws and thinned draws, and for Gaussian and Student-t families.
  • It checks the cases that return NULL: clusters with more than one draw, a weight of 0, and a user-defined predvar. It checks that wcol_moments_c() gives a non-finite sum of squares for a non-finite projected draw.
  • test_search_gauss.R:
  • It uses trace() to check that validated LOO uses as_moments = TRUE. The ranking and the submodel summaries are equal to those from the full get_p_clust() (tolerance 1e-12).
  • In one reference model, the ranking contains the factor interaction x1:f, so perf_eval_gauss() does not apply and perf_eval() uses the full result.
  • A latent-projection case with thinned draws compares the results with the default options(projpred.search_gauss_crossprod = TRUE) and with FALSE. This option of Gaussian forward search and validated LOO evaluation from cross-products of the model matrix #579 switches the Gaussian search and evaluation code on and off.
  • Other tests compare cands_single_label() with select_possible_terms_size(), and the search with and without the new candidate code (single_label_info() mocked to NULL). The search terms include sets without "1", with a duplicated term, and with one interaction in two orders.
• The full test suite passes on the workstation (R 4.6.1, testthat 3.3.2, 58,119 expectations).

This PR was assisted by Claude.

[avehtari]

I pushed a1d321c. get_p_clust(as_moments = TRUE) no longer checks with range() that the projected draws are finite. range() copied the observations x draws matrix for each left-out observation. perf_eval_gauss() now detects a non-finite draw from the weighted sums of squares that it already computes. For such a draw, it uses the full get_p_clust() result, as before.

In case 2 of the description, full validated LOO took 34.9 s before this commit and 29.4 s with it (Triton, R 4.5.3, one SLURM job, one run each). The vsel objects of the two versions are identical. The new and changed tests fail without the commit, and the full test suite passes (58,146 expectations). I updated the description.

[avehtari]

I pushed c7c6166 and 10b714e. In each search step, search_forward_gauss() looked up the model-matrix columns of each candidate by name. Now it looks up the columns of each search term once per search, if each search term is "1" or has one term label.

In case 1 of the description (200 predictors), full validated LOO took 8.3 s before these commits and 6.7 s with them. A second job gave 9.3 s and 7.5 s.