From af78fc0e0e07e4759d3322193ab4a17c022ff979 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Sun, 25 Oct 2026 15:19:45 +0200 Subject: [PATCH 1/8] Rank the terms of a Gaussian forward search from cross-products of the model matrix --- R/projpred-package.R | 4 + R/search.R | 213 +++++++++++++++++++++++++++++ man/internal-global-options.Rd | 4 + tests/testthat/test_search_gauss.R | 76 ++++++++++ 4 files changed, 297 insertions(+) create mode 100644 tests/testthat/test_search_gauss.R diff --git a/R/projpred-package.R b/R/projpred-package.R index fb1c0508e..f571a7927 100644 --- a/R/projpred-package.R +++ b/R/projpred-package.R @@ -266,6 +266,10 @@ #' [#207](https://github.com/stan-dev/projpred/issues/207) and GitHub pull #' request [#353](https://github.com/stan-dev/projpred/pull/353)). Currently, #' this is an experimental feature. +#' * `projpred.search_gauss_crossprod`: A single logical value indicating +#' whether to compute the forward search for Gaussian submodels with the +#' identity link from cross-products of the model matrix (default `TRUE`). With +#' `FALSE`, \pkg{projpred} projects onto each candidate submodel draw by draw. #' #' @name internal-global-options #' @keywords internal diff --git a/R/search.R b/R/search.R index 1384d2bf2..542eabb5d 100644 --- a/R/search.R +++ b/R/search.R @@ -7,6 +7,18 @@ search_forward <- function(p_ref, refmodel, nterms_max, verbose = TRUE, stop("Did not expect `search_terms` to be `NULL`. Please report this.") } + # Gaussian submodels: rank the terms from cross-products (see + # search_forward_gauss()); `NULL` if this case does not apply: + if (getOption("projpred.search_gauss_crossprod", TRUE)) { + out_gauss <- search_forward_gauss( + p_ref = p_ref, refmodel = refmodel, nterms_max = nterms_max, + verbose = verbose, search_terms = search_terms, ... + ) + if (!is.null(out_gauss)) { + return(out_gauss) + } + } + chosen <- character() outdmins <- c() @@ -142,6 +154,207 @@ search_forward <- function(p_ref, refmodel, nterms_max, verbose = TRUE, return(nlist(predictor_ranking = setdiff(chosen, "1"), outdmins)) } +# Forward search for Gaussian submodels with the identity link, from the +# cross-products of the model matrix of all search terms. Let J be the +# model-matrix columns of a candidate. For draw s, the projection solves the +# weighted ridge problem of glm_ridge(): +# (G[J, J] + regul * diag(sd[J]^2)) beta_s = C[J, s]. +# Here, G = Xc' W Xc and C = Xc' W M. Xc is the model matrix, centered at the +# weighted column means. W holds the observation weights, M the fitted values +# of the projected reference draws, and sd the weighted column standard +# deviations from .weighted_sd(). glm_ridge() centers and scales the columns +# in the same way and does not penalize the intercept. The part of ce_gauss() +# that differs between the candidates is +# -beta_s' (G[J, J] + 2 * regul * diag(sd[J]^2)) beta_s / sum(W). +# Thus, the cost per candidate does not depend on the number of observations. +# This function gives only the ranking. The submodel fits along the ranking come +# from proj_to_submodl(), as in search_forward(). The function returns `NULL` +# if one of its conditions does not hold; then search_forward() continues with +# its own search. +search_forward_gauss <- function(p_ref, refmodel, nterms_max, verbose, + search_terms, ...) { + fam <- refmodel$family + glm_fitter <- getOption("projpred.glm_fitter", "fit_glm_ridge_callback") + # Custom search terms are possible, because the candidates come from + # select_possible_terms_size(), as in search_forward(). + if (fam$family != "gaussian" || + fam$link != "identity" || isTRUE(fam$for_latent) || + isTRUE(fam$for_augdat) || !identical(refmodel$div_minimizer, divmin) || + !identical(glm_fitter, "fit_glm_ridge_callback") || + any(refmodel$offset != 0)) { + return(NULL) + } + # The arguments for the divergence minimizer, as in proj_to_submodl(): + dot_args <- list(...) + divmin_args <- dot_args$search_control %||% + dot_args[setdiff(names(dot_args), "search_control")] + # fit_glm_ridge_callback() uses `regul` and passes only the formal arguments + # of glm_ridge() on. Of these, only `thresh` (not used for this family) is + # allowed here: + args_glm_ridge <- intersect(names(divmin_args), methods::formalArgs(glm_ridge)) + if (!all(args_glm_ridge %in% "thresh")) { + return(NULL) + } + regul <- divmin_args$regul %||% 1e-4 + + # Model matrix of all search terms (a search term can consist of several + # terms, for example "x1 + x2 + x1:x2"): + terms_all <- setdiff(unique(trimws(unlist( + strsplit(search_terms, "+", fixed = TRUE) + ))), "1") + if (length(terms_all) == 0) { + return(NULL) + } + trms <- extract_terms_response(make_formula(terms_all)) + if (length(trms$group_terms) > 0 || length(trms$additive_terms) > 0) { + return(NULL) + } + tt <- terms(as.formula(paste("~", paste(terms_all, collapse = " + ")))) + data <- refmodel$fetch_data() + mf <- model.frame(tt, data = data, drop.unused.levels = TRUE) + if (nrow(mf) != nrow(data)) { + return(NULL) + } + # The dummy coding of a factor in an interaction depends on the other terms + # of the submodel, so these cases are not covered: + fac <- attr(tt, "factors") + is_categ <- vapply(rownames(fac), function(v) { + is.factor(mf[[v]]) || is.character(mf[[v]]) || is.logical(mf[[v]]) + }, logical(1)) + if (any(attr(tt, "order") > 1 & + colSums(fac[is_categ, , drop = FALSE]) > 0)) { + return(NULL) + } + x <- model.matrix(tt, data = mf) + assign_x <- attr(x, "assign") + x <- x[, assign_x > 0, drop = FALSE] + assign_x <- assign_x[assign_x > 0] + norm_label <- function(lbl) { + vapply(strsplit(lbl, ":", fixed = TRUE), function(parts) { + paste(sort(trimws(parts)), collapse = ":") + }, character(1)) + } + cols_by_term <- split(seq_along(assign_x), + norm_label(attr(tt, "term.labels"))[assign_x]) + + # Cross-products: + wobs <- refmodel$wobs + sum_w <- sum(wobs) + n <- nrow(x) + x_c <- sweep(x, 2, colSums(x * wobs) / sum_w) + g_mat <- crossprod(x_c, x_c * wobs) + c_mat <- crossprod(x_c, p_ref$mu * wobs) + sd_x <- sqrt(n / (n - 1) * diag(g_mat) / sum_w) + if (n < 2 || any(!is.finite(sd_x)) || any(sd_x == 0)) { + return(NULL) + } + wdraws <- p_ref$wdraws_prj + lambda_x <- regul * sd_x^2 + a_full <- g_mat + diag(a_full) <- diag(a_full) + lambda_x + + # The term labels of a candidate string. The same strings come back in every + # step, so the labels are computed once: + label_cache <- new.env(hash = TRUE, parent = emptyenv()) + labels_of <- function(cand) { + lbls <- label_cache[[cand]] + if (is.null(lbls)) { + lbls <- character(0) + if (!identical(cand, "1")) { + lbls <- unique(norm_label(trimws(strsplit(cand, "+", fixed = TRUE)[[1]]))) + } + assign(cand, lbls, envir = label_cache) + } + return(lbls) + } + + # The criterion is -sum_s wdraws_s * beta_s' (G[J, J] + 2 * Lambda) beta_s + # (smaller is better, as for the cross-entropy). With + # G[J, J] + Lambda = R' R and y_s = R^-T C[J, s], it equals + # -sum_s wdraws_s * (||y_s||^2 + beta_s' Lambda beta_s). In one step, all + # candidates share the columns of the chosen terms, so the Cholesky factor + # `r_chosen` and `y_chosen` of these columns are computed once per step. Each + # candidate extends them by its new columns `cols_new` (block Cholesky). + crit_cand <- function(cols_new, cols_chosen, r_chosen, y_chosen) { + quad <- 0 + r_cn <- NULL + r_n <- NULL + y_n <- NULL + beta_n <- NULL + if (length(cols_new) > 0) { + a_nn <- a_full[cols_new, cols_new, drop = FALSE] + rhs_n <- c_mat[cols_new, , drop = FALSE] + if (length(cols_chosen) > 0) { + r_cn <- forwardsolve(t(r_chosen), + a_full[cols_chosen, cols_new, drop = FALSE]) + a_nn <- a_nn - crossprod(r_cn) + rhs_n <- rhs_n - crossprod(r_cn, y_chosen) + } + r_n <- chol(a_nn) + y_n <- forwardsolve(t(r_n), rhs_n) + beta_n <- backsolve(r_n, y_n) + quad <- colSums(y_n^2) + colSums(lambda_x[cols_new] * beta_n^2) + } + if (length(cols_chosen) > 0) { + y_c <- y_chosen + if (!is.null(beta_n)) { + y_c <- y_c - r_cn %*% beta_n + } + beta_c <- backsolve(r_chosen, y_c) + quad <- quad + colSums(y_chosen^2) + + colSums(lambda_x[cols_chosen] * beta_c^2) + } + return(list(crit = -sum(wdraws * quad), r_cn = r_cn, r_n = r_n, y_n = y_n)) + } + + chosen <- character() + labels_chosen <- character() + cols_chosen <- integer(0) + r_chosen <- NULL + y_chosen <- NULL + for (size in seq_len(nterms_max + 1L)) { + cands <- select_possible_terms_size(chosen, search_terms, size = size) + if (is.null(cands)) { + next + } + labels_new <- lapply(cands, function(cand) { + setdiff(labels_of(cand), labels_chosen) + }) + if (!all(unlist(labels_new) %in% names(cols_by_term))) { + return(NULL) + } + res_cands <- lapply(labels_new, function(lbls) { + crit_cand(unlist(cols_by_term[lbls], use.names = FALSE), cols_chosen, + r_chosen, y_chosen) + }) + imin <- which.min(vapply(res_cands, `[[`, numeric(1), "crit")) + chosen <- c(chosen, cands[imin]) + cols_new <- unlist(cols_by_term[labels_new[[imin]]], use.names = FALSE) + if (length(cols_new) > 0) { + best <- res_cands[[imin]] + if (length(cols_chosen) == 0) { + r_chosen <- best$r_n + y_chosen <- best$y_n + } else { + r_chosen <- rbind( + cbind(r_chosen, best$r_cn), + cbind(matrix(0, length(cols_new), length(cols_chosen)), best$r_n) + ) + y_chosen <- rbind(y_chosen, best$y_n) + } + cols_chosen <- c(cols_chosen, cols_new) + labels_chosen <- c(labels_chosen, labels_new[[imin]]) + } + } + + # The submodel fits along the ranking (as stored by search_forward()): + outdmins <- lapply(seq_along(chosen), function(j) { + proj_to_submodl(predictor_terms = chosen[seq_len(j)], p_ref = p_ref, + refmodel = refmodel, verbose = verbose, ...)$outdmin + }) + return(nlist(predictor_ranking = setdiff(chosen, "1"), outdmins)) +} + #' Force search terms #' #' A helper function to construct the input for argument `search_terms` of diff --git a/man/internal-global-options.Rd b/man/internal-global-options.Rd index e3b9d2560..001878a63 100644 --- a/man/internal-global-options.Rd +++ b/man/internal-global-options.Rd @@ -40,6 +40,10 @@ projections (see GitHub issue \href{https://github.com/stan-dev/projpred/issues/207}{#207} and GitHub pull request \href{https://github.com/stan-dev/projpred/pull/353}{#353}). Currently, this is an experimental feature. +\item \code{projpred.search_gauss_crossprod}: A single logical value indicating +whether to compute the forward search for Gaussian submodels with the +identity link from cross-products of the model matrix (default \code{TRUE}). With +\code{FALSE}, \pkg{projpred} projects onto each candidate submodel draw by draw. } } \keyword{internal} diff --git a/tests/testthat/test_search_gauss.R b/tests/testthat/test_search_gauss.R new file mode 100644 index 000000000..78c886338 --- /dev/null +++ b/tests/testthat/test_search_gauss.R @@ -0,0 +1,76 @@ +context("forward search from cross-products for Gaussian submodels") + +# A Gaussian reference model with simulated draws of the coefficients (no model +# fit is necessary). `fml` gives the predictor terms. +refm_search_gauss <- function(fml, n = 80, ndraws = 200, seed = 6283) { + set.seed(seed) + dat <- data.frame(x1 = rnorm(n), x2 = rnorm(n), x3 = rnorm(n), + x4 = rnorm(n), + f = factor(sample(c("a", "b", "c"), n, replace = TRUE))) + dat$y <- 0.8 * dat$x1 - 0.4 * dat$x2 + 0.3 * dat$x1 * dat$x2 + + c(a = 0, b = 0.5, c = -0.5)[as.character(dat$f)] + rnorm(n) + mm <- model.matrix(fml, dat) + beta_hat <- drop(solve(crossprod(mm) + diag(1e-6, ncol(mm)), + crossprod(mm, dat$y))) + beta <- beta_hat + matrix(rnorm(ncol(mm) * ndraws, sd = 0.05), ncol(mm)) + fit <- list(beta = beta, data = dat, fml = fml) + ref_predfun <- function(fit, newdata = NULL) { + nd <- if (is.null(newdata)) fit$data else newdata + model.matrix(delete.response(terms(fit$fml)), nd) %*% fit$beta + } + suppressWarnings(init_refmodel(fit, data = dat, formula = fml, + family = gaussian(), ref_predfun = ref_predfun, + dis = rep(1, ndraws))) +} + +fml_search_gauss <- y ~ x1 + x2 + x3 + x4 + f + x1:x2 + poly(x3, 2) + +test_that(paste( + "the cross-product search gives the same results as the search with", + "draw-wise projections" +), { + refm <- refm_search_gauss(fml_search_gauss) + for (nclusters_crr in c(1, 5)) { + vs_crossprod <- suppressWarnings(varsel( + refm, nclusters = nclusters_crr, nclusters_pred = 7, nterms_max = 6, + seed = 3317, verbose = FALSE + )) + opts_old <- options(projpred.search_gauss_crossprod = FALSE) + vs_drawwise <- suppressWarnings(varsel( + refm, nclusters = nclusters_crr, nclusters_pred = 7, nterms_max = 6, + seed = 3317, verbose = FALSE + )) + options(opts_old) + expect_identical(ranking(vs_crossprod)$fulldata, + ranking(vs_drawwise)$fulldata, info = nclusters_crr) + expect_equal(performances(vs_crossprod, deltas = FALSE), + performances(vs_drawwise, deltas = FALSE), + info = nclusters_crr) + expect_equal(vs_crossprod$search_path$outdmins, + vs_drawwise$search_path$outdmins, info = nclusters_crr) + } +}) + +test_that("search_forward_gauss() applies to Gaussian GLM-type submodels", { + refm <- refm_search_gauss(fml_search_gauss) + p_ref <- get_refdist(refm, nclusters = 3) + search_terms <- union("1", split_formula(refm$formula, + data = refm$fetch_data())) + out <- search_forward_gauss(p_ref = p_ref, refmodel = refm, nterms_max = 4, + verbose = FALSE, search_terms = search_terms) + expect_false(is.null(out)) + expect_length(out$predictor_ranking, 4) + expect_length(out$outdmins, 5) +}) + +test_that(paste( + "search_forward_gauss() returns `NULL` for a factor in an interaction term" +), { + refm <- refm_search_gauss(y ~ x1 + x2 + f + x1:f) + p_ref <- get_refdist(refm, nclusters = 3) + search_terms <- union("1", split_formula(refm$formula, + data = refm$fetch_data())) + out <- search_forward_gauss(p_ref = p_ref, refmodel = refm, nterms_max = 3, + verbose = FALSE, search_terms = search_terms) + expect_null(out) +}) From 34b142ecc5a9cfb7ae40e0e91606eee3e34d2968 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Fri, 23 Oct 2026 13:16:29 +0300 Subject: [PATCH 2/8] Evaluate Gaussian submodels in the LOO CV with validated search from cross-products --- R/gauss_crossprod.R | 234 +++++++++++++++++++++++++++++ R/projfun.R | 15 +- R/projpred-package.R | 7 +- R/search.R | 81 ++-------- man/internal-global-options.Rd | 7 +- tests/testthat/test_search_gauss.R | 37 +++++ 6 files changed, 304 insertions(+), 77 deletions(-) create mode 100644 R/gauss_crossprod.R diff --git a/R/gauss_crossprod.R b/R/gauss_crossprod.R new file mode 100644 index 000000000..994ccd2cf --- /dev/null +++ b/R/gauss_crossprod.R @@ -0,0 +1,234 @@ +# Gaussian submodels from cross-products -------------------------------------- +# +# For a Gaussian submodel with the identity link, glm_ridge() solves a weighted +# ridge least-squares problem: the columns of the model matrix are centered and +# scaled with the observation weights, and the intercept is not penalized. Let +# Xc be the model matrix of all terms, centered at the weighted column means, +# W the observation weights, M the fitted values of the projected reference +# draws, sd the weighted column standard deviations (as in .weighted_sd()), and +# Lambda = regul * diag(sd^2). For the columns J of a submodel and draw s, the +# coefficients solve +# (G[J, J] + Lambda[J, J]) beta_s = C[J, s], G = Xc' W Xc, C = Xc' W M, +# and the intercept is ybar_s - m[J]' beta_s, with ybar_s the weighted mean of +# M[, s] and m the weighted column means. With G[J, J] + Lambda[J, J] = R' R and +# y_s = R^-T C[J, s], let q_s = ||y_s||^2 + beta_s' Lambda[J, J] beta_s. Then +# sum(W * (M[, s] - fit_s)^2) = sum(W * (M[, s] - ybar_s)^2) - q_s, +# sum(W * (-2 * M[, s] * fit_s + fit_s^2)) = -ybar_s^2 * sum(W) - q_s. +# These give the projected dispersion (dis_gauss()) and the cross-entropy +# (ce_gauss()) without the n-dimensional fitted values. + +# The condition for the cross-product computations: Gaussian family with the +# identity link (traditional projection), the default divergence minimizer and +# submodel fitter, zero offsets, and no argument for glm_ridge() other than +# `thresh` (not used for this family) in `divmin_args`. Returns `regul` or +# `NULL`. +gauss_crossprod_regul <- function(refmodel, divmin_args) { + fam <- refmodel$family + glm_fitter <- getOption("projpred.glm_fitter", "fit_glm_ridge_callback") + if (fam$family != "gaussian" || fam$link != "identity" || + isTRUE(fam$for_latent) || isTRUE(fam$for_augdat) || + !identical(refmodel$div_minimizer, divmin) || + !identical(glm_fitter, "fit_glm_ridge_callback") || + any(refmodel$offset != 0)) { + return(NULL) + } + # fit_glm_ridge_callback() uses `regul` and passes only the formal arguments + # of glm_ridge() on: + args_glm_ridge <- intersect(names(divmin_args), methods::formalArgs(glm_ridge)) + if (!all(args_glm_ridge %in% "thresh")) { + return(NULL) + } + return(divmin_args$regul %||% 1e-4) +} + +# Term labels with the variables of an interaction in sorted order (as terms() +# treats them). +norm_label_gauss <- function(lbl) { + vapply(strsplit(lbl, ":", fixed = TRUE), function(parts) { + paste(sort(trimws(parts)), collapse = ":") + }, character(1)) +} + +# The model matrix of the terms in `term_strings` (a string can consist of +# several terms, for example "x1 + x2 + x1:x2") for the training data of +# `refmodel`, and its cross-products. Returns `NULL` for group-level or smooth +# terms, a factor in an interaction term (its dummy coding depends on the other +# terms of a submodel), missing values, or a zero-variance column. +gauss_crossprod_design <- function(refmodel, term_strings, regul) { + terms_all <- setdiff(unique(trimws(unlist( + strsplit(term_strings, "+", fixed = TRUE) + ))), "1") + if (length(terms_all) == 0) { + return(NULL) + } + trms <- extract_terms_response(make_formula(terms_all)) + if (length(trms$group_terms) > 0 || length(trms$additive_terms) > 0) { + return(NULL) + } + tt <- terms(as.formula(paste("~", paste(terms_all, collapse = " + ")))) + data <- refmodel$fetch_data() + mf <- model.frame(tt, data = data, drop.unused.levels = TRUE) + if (nrow(mf) != nrow(data)) { + return(NULL) + } + fac <- attr(tt, "factors") + is_categ <- vapply(rownames(fac), function(v) { + is.factor(mf[[v]]) || is.character(mf[[v]]) || is.logical(mf[[v]]) + }, logical(1)) + if (any(attr(tt, "order") > 1 & + colSums(fac[is_categ, , drop = FALSE]) > 0)) { + return(NULL) + } + x <- model.matrix(tt, data = mf) + assign_x <- attr(x, "assign") + keep <- assign_x > 0 + x <- x[, keep, drop = FALSE] + assign_x <- assign_x[keep] + cols_by_term <- split(seq_along(assign_x), + norm_label_gauss(attr(tt, "term.labels"))[assign_x]) + wobs <- refmodel$wobs + sum_w <- sum(wobs) + n <- nrow(x) + x_means <- colSums(x * wobs) / sum_w + x_c <- sweep(x, 2, x_means) + g_mat <- crossprod(x_c, x_c * wobs) + sd_x <- sqrt(n / (n - 1) * diag(g_mat) / sum_w) + if (n < 2 || any(!is.finite(sd_x)) || any(sd_x == 0)) { + return(NULL) + } + lambda_x <- regul * sd_x^2 + a_full <- g_mat + diag(a_full) <- diag(a_full) + lambda_x + # For new data, the terms of the model frame keep the data-dependent bases + # of the training data (attribute "predvars", for example of poly()): + return(nlist(tt_pred = attr(mf, "terms"), xlev = .getXlevels(tt, mf), keep, + x_c, x_means, cols_by_term, g_mat, a_full, lambda_x, wobs, + sum_w)) +} + +# The term labels of a term string (see norm_label_gauss()). +labels_gauss <- function(term_string) { + if (identical(term_string, "1")) { + return(character(0)) + } + return(unique(norm_label_gauss(trimws( + strsplit(term_string, "+", fixed = TRUE)[[1]] + )))) +} + +# The performance evaluation of perf_eval() along the predictor ranking for +# Gaussian submodels, from cross-products (see the top of this file): for each +# size in `nterms`, the projected dispersion, the cross-entropy and the +# predictions for the test data, without the draw-wise submodel fits. Returns +# the list `out_by_size` of perf_eval() (elements `sub_summary`, `ce` and +# `dis_sub`), or `NULL` if a condition does not hold. +perf_eval_gauss <- function(search_path, nterms, refmodel, p_ref, + refmodel_fulldata, indices_test, newdata_test, + offset_test, wobs_test, y_test, y_oscale_test, + divmin_args) { + regul <- gauss_crossprod_regul(refmodel, divmin_args) + ranking <- search_path$predictor_ranking + if (is.null(regul) || length(nterms) == 0 || any(nterms > length(ranking)) || + any(offset_test != 0)) { + return(NULL) + } + des <- gauss_crossprod_design(refmodel, c("1", ranking), regul) + if (is.null(des)) { + return(NULL) + } + data_test <- newdata_test %||% refmodel_fulldata$fetch_data(obs = indices_test) + x_test <- tryCatch({ + mf_test <- model.frame(des$tt_pred, data = data_test, xlev = des$xlev, + na.action = na.pass) + model.matrix(des$tt_pred, data = mf_test)[, des$keep, drop = FALSE] + }, error = function(e) NULL) + if (is.null(x_test) || nrow(x_test) != nrow(data_test) || anyNA(x_test)) { + return(NULL) + } + x_test_c <- sweep(x_test, 2, des$x_means) + + # Quantities of the projected reference draws: + wobs <- des$wobs + sum_w <- des$sum_w + mu_ref <- p_ref$mu_offs + ybar <- colSums(mu_ref * wobs) / sum_w + sum_wr2 <- colSums(sweep(mu_ref, 2, ybar)^2 * wobs) + sum_wvar <- colSums(p_ref$var * wobs) + c_mat <- crossprod(des$x_c, mu_ref * wobs) + nprj <- ncol(mu_ref) + + # Columns of the submodel of each size (the submodels are nested): + labels_size <- lapply(seq_along(ranking), function(j) { + unique(unlist(lapply(ranking[seq_len(j)], labels_gauss))) + }) + if (!all(unlist(labels_size) %in% names(des$cols_by_term))) { + return(NULL) + } + + out_size <- function(cols, r_chol, y_mat) { + if (length(cols) == 0) { + quad <- rep(0, nprj) + mu_test <- matrix(ybar, nrow = nrow(x_test), ncol = nprj, byrow = TRUE) + } else { + beta <- backsolve(r_chol, y_mat) + quad <- colSums(y_mat^2) + colSums(des$lambda_x[cols] * beta^2) + mu_test <- sweep(x_test_c[, cols, drop = FALSE] %*% beta, 2, ybar, "+") + } + rownames(mu_test) <- rownames(x_test) + dis <- sqrt((sum_wvar + sum_wr2 - quad) / sum_w) + ce <- weighted.mean((-ybar^2 * sum_w - quad) / sum_w, p_ref$wdraws_prj) + sub_summary <- weighted_summary_means( + y_wobs_test = data.frame(y = y_test, y_oscale = y_oscale_test, + wobs = wobs_test), + data_aux_test = data_test, + family = refmodel$family, + wdraws = p_ref$wdraws_prj, + mu = mu_test, + dis = dis, + cl_ref = p_ref$cl, + wdraws_ref = p_ref$wdraws_orig + ) + return(list(sub_summary = sub_summary, ce = ce, dis_sub = dis)) + } + + # Grow the Cholesky factor along the ranking (block Cholesky): + out_by_size <- vector("list", length(nterms)) + cols <- integer(0) + r_chol <- NULL + y_mat <- NULL + for (size_j in seq(0, max(nterms))) { + if (size_j > 0) { + cols_new <- setdiff( + unlist(des$cols_by_term[labels_size[[size_j]]], use.names = FALSE), + cols + ) + if (length(cols_new) > 0) { + a_nn <- des$a_full[cols_new, cols_new, drop = FALSE] + rhs_n <- c_mat[cols_new, , drop = FALSE] + if (length(cols) > 0) { + r_cn <- forwardsolve(t(r_chol), + des$a_full[cols, cols_new, drop = FALSE]) + a_nn <- a_nn - crossprod(r_cn) + rhs_n <- rhs_n - crossprod(r_cn, y_mat) + } + r_n <- chol(a_nn) + y_n <- forwardsolve(t(r_n), rhs_n) + if (length(cols) == 0) { + r_chol <- r_n + y_mat <- y_n + } else { + r_chol <- rbind( + cbind(r_chol, r_cn), + cbind(matrix(0, length(cols_new), length(cols)), r_n) + ) + y_mat <- rbind(y_mat, y_n) + } + cols <- c(cols, cols_new) + } + } + for (k in which(nterms == size_j)) { + out_by_size[[k]] <- out_size(cols, r_chol, y_mat) + } + } + return(out_by_size) +} diff --git a/R/projfun.R b/R/projfun.R index aa79dbf31..61625eb47 100644 --- a/R/projfun.R +++ b/R/projfun.R @@ -106,7 +106,20 @@ perf_eval <- function(search_path, txt_clust_draws(p_ref[["clust_used"]], p_ref[["nprjdraws"]]), " (`refit_prj = ", refit_prj, "`) ...", verbose = verbose && !return_submodls) - out_by_size <- lapply(nterms, function(size_j) { + out_by_size <- NULL + if (refit_prj && !is.null(reweighting_args) && !return_submodls && + !return_preds && getOption("projpred.search_gauss_crossprod", TRUE)) { + # In the LOO CV with validated search, compute the evaluation of Gaussian + # submodels from cross-products (`NULL` if this case does not apply): + out_by_size <- perf_eval_gauss( + search_path = search_path, nterms = nterms, refmodel = refmodel, + p_ref = p_ref, refmodel_fulldata = refmodel_fulldata, + indices_test = indices_test, newdata_test = newdata_test, + offset_test = offset_test, wobs_test = wobs_test, y_test = y_test, + y_oscale_test = y_oscale_test, divmin_args = list(...) + ) + } + if (is.null(out_by_size)) 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, ...) diff --git a/R/projpred-package.R b/R/projpred-package.R index f571a7927..011696408 100644 --- a/R/projpred-package.R +++ b/R/projpred-package.R @@ -267,9 +267,10 @@ #' request [#353](https://github.com/stan-dev/projpred/pull/353)). Currently, #' this is an experimental feature. #' * `projpred.search_gauss_crossprod`: A single logical value indicating -#' whether to compute the forward search for Gaussian submodels with the -#' identity link from cross-products of the model matrix (default `TRUE`). With -#' `FALSE`, \pkg{projpred} projects onto each candidate submodel draw by draw. +#' whether to compute the forward search and the evaluation in the LOO CV with +#' validated search for Gaussian submodels with the identity link from +#' cross-products of the model matrix (default `TRUE`). With `FALSE`, +#' \pkg{projpred} projects onto each candidate submodel draw by draw. #' #' @name internal-global-options #' @keywords internal diff --git a/R/search.R b/R/search.R index 542eabb5d..e7ba539c7 100644 --- a/R/search.R +++ b/R/search.R @@ -173,85 +173,26 @@ search_forward <- function(p_ref, refmodel, nterms_max, verbose = TRUE, # its own search. search_forward_gauss <- function(p_ref, refmodel, nterms_max, verbose, search_terms, ...) { - fam <- refmodel$family - glm_fitter <- getOption("projpred.glm_fitter", "fit_glm_ridge_callback") - # Custom search terms are possible, because the candidates come from - # select_possible_terms_size(), as in search_forward(). - if (fam$family != "gaussian" || - fam$link != "identity" || isTRUE(fam$for_latent) || - isTRUE(fam$for_augdat) || !identical(refmodel$div_minimizer, divmin) || - !identical(glm_fitter, "fit_glm_ridge_callback") || - any(refmodel$offset != 0)) { - return(NULL) - } # The arguments for the divergence minimizer, as in proj_to_submodl(): dot_args <- list(...) divmin_args <- dot_args$search_control %||% dot_args[setdiff(names(dot_args), "search_control")] - # fit_glm_ridge_callback() uses `regul` and passes only the formal arguments - # of glm_ridge() on. Of these, only `thresh` (not used for this family) is - # allowed here: - args_glm_ridge <- intersect(names(divmin_args), methods::formalArgs(glm_ridge)) - if (!all(args_glm_ridge %in% "thresh")) { - return(NULL) - } - regul <- divmin_args$regul %||% 1e-4 - - # Model matrix of all search terms (a search term can consist of several - # terms, for example "x1 + x2 + x1:x2"): - terms_all <- setdiff(unique(trimws(unlist( - strsplit(search_terms, "+", fixed = TRUE) - ))), "1") - if (length(terms_all) == 0) { - return(NULL) - } - trms <- extract_terms_response(make_formula(terms_all)) - if (length(trms$group_terms) > 0 || length(trms$additive_terms) > 0) { + regul <- gauss_crossprod_regul(refmodel, divmin_args) + if (is.null(regul)) { return(NULL) } - tt <- terms(as.formula(paste("~", paste(terms_all, collapse = " + ")))) - data <- refmodel$fetch_data() - mf <- model.frame(tt, data = data, drop.unused.levels = TRUE) - if (nrow(mf) != nrow(data)) { - return(NULL) - } - # The dummy coding of a factor in an interaction depends on the other terms - # of the submodel, so these cases are not covered: - fac <- attr(tt, "factors") - is_categ <- vapply(rownames(fac), function(v) { - is.factor(mf[[v]]) || is.character(mf[[v]]) || is.logical(mf[[v]]) - }, logical(1)) - if (any(attr(tt, "order") > 1 & - colSums(fac[is_categ, , drop = FALSE]) > 0)) { - return(NULL) - } - x <- model.matrix(tt, data = mf) - assign_x <- attr(x, "assign") - x <- x[, assign_x > 0, drop = FALSE] - assign_x <- assign_x[assign_x > 0] - norm_label <- function(lbl) { - vapply(strsplit(lbl, ":", fixed = TRUE), function(parts) { - paste(sort(trimws(parts)), collapse = ":") - }, character(1)) - } - cols_by_term <- split(seq_along(assign_x), - norm_label(attr(tt, "term.labels"))[assign_x]) - - # Cross-products: - wobs <- refmodel$wobs - sum_w <- sum(wobs) - n <- nrow(x) - x_c <- sweep(x, 2, colSums(x * wobs) / sum_w) - g_mat <- crossprod(x_c, x_c * wobs) - c_mat <- crossprod(x_c, p_ref$mu * wobs) - sd_x <- sqrt(n / (n - 1) * diag(g_mat) / sum_w) - if (n < 2 || any(!is.finite(sd_x)) || any(sd_x == 0)) { + # Custom search terms are possible, because the candidates come from + # select_possible_terms_size(), as in search_forward(): + des <- gauss_crossprod_design(refmodel, search_terms, regul) + if (is.null(des)) { return(NULL) } + cols_by_term <- des$cols_by_term + c_mat <- crossprod(des$x_c, p_ref$mu * des$wobs) + lambda_x <- des$lambda_x + a_full <- des$a_full wdraws <- p_ref$wdraws_prj - lambda_x <- regul * sd_x^2 - a_full <- g_mat - diag(a_full) <- diag(a_full) + lambda_x + norm_label <- norm_label_gauss # The term labels of a candidate string. The same strings come back in every # step, so the labels are computed once: diff --git a/man/internal-global-options.Rd b/man/internal-global-options.Rd index 001878a63..088914b4e 100644 --- a/man/internal-global-options.Rd +++ b/man/internal-global-options.Rd @@ -41,9 +41,10 @@ projections (see GitHub issue request \href{https://github.com/stan-dev/projpred/pull/353}{#353}). Currently, this is an experimental feature. \item \code{projpred.search_gauss_crossprod}: A single logical value indicating -whether to compute the forward search for Gaussian submodels with the -identity link from cross-products of the model matrix (default \code{TRUE}). With -\code{FALSE}, \pkg{projpred} projects onto each candidate submodel draw by draw. +whether to compute the forward search and the evaluation in the LOO CV with +validated search for Gaussian submodels with the identity link from +cross-products of the model matrix (default \code{TRUE}). With \code{FALSE}, +\pkg{projpred} projects onto each candidate submodel draw by draw. } } \keyword{internal} diff --git a/tests/testthat/test_search_gauss.R b/tests/testthat/test_search_gauss.R index 78c886338..53fe735d0 100644 --- a/tests/testthat/test_search_gauss.R +++ b/tests/testthat/test_search_gauss.R @@ -51,6 +51,43 @@ test_that(paste( } }) +test_that(paste( + "the LOO CV with validated search gives the same results with and without", + "the cross-product computations" +), { + # Without poly(): the draw-wise path cannot predict a poly() term for one + # left-out observation. + refm <- refm_search_gauss(y ~ x1 + x2 + x3 + x4 + f + x1:x2, n = 40) + args_cv <- list(refm, cv_method = "LOO", validate_search = TRUE, nloo = 10, + nclusters = 3, ndraws_pred = 20, nterms_max = 4, seed = 5821, + verbose = FALSE) + cv_crossprod <- suppressWarnings(do.call(cv_varsel, args_cv)) + opts_old <- options(projpred.search_gauss_crossprod = FALSE) + cv_drawwise <- suppressWarnings(do.call(cv_varsel, args_cv)) + options(opts_old) + expect_identical(ranking(cv_crossprod), ranking(cv_drawwise)) + expect_equal(cv_crossprod$summaries$sub, cv_drawwise$summaries$sub, + tolerance = 1e-10) + expect_equal(performances(cv_crossprod, deltas = FALSE), + performances(cv_drawwise, deltas = FALSE), tolerance = 1e-10) +}) + +test_that("perf_eval_gauss() applies in the LOO CV with validated search", { + # Without poly(): the draw-wise path cannot predict a poly() term for one + # left-out observation. + refm <- refm_search_gauss(y ~ x1 + x2 + x3 + x4 + f + x1:x2, n = 40) + n_used <- 0L + trace("perf_eval_gauss", exit = function() { + if (!is.null(returnValue())) n_used <<- n_used + 1L + }, print = FALSE, where = asNamespace("projpred")) + on.exit(untrace("perf_eval_gauss", where = asNamespace("projpred"))) + cv <- suppressWarnings(cv_varsel( + refm, cv_method = "LOO", validate_search = TRUE, nloo = 5, nclusters = 3, + ndraws_pred = 20, nterms_max = 3, seed = 5821, verbose = FALSE + )) + expect_identical(n_used, 5L) +}) + test_that("search_forward_gauss() applies to Gaussian GLM-type submodels", { refm <- refm_search_gauss(fml_search_gauss) p_ref <- get_refdist(refm, nclusters = 3) From e44a2b3f71eee799a1671c5050b1579d2af4e414 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Fri, 23 Oct 2026 13:22:55 +0300 Subject: [PATCH 3/8] Compute the cross-product criterion for one-column candidates with matrix operations --- NEWS.md | 2 ++ R/search.R | 53 +++++++++++++++++++++++++++++++++++++++++++++++------ 2 files changed, 49 insertions(+), 6 deletions(-) diff --git a/NEWS.md b/NEWS.md index 94ffe21cd..f37180a88 100644 --- a/NEWS.md +++ b/NEWS.md @@ -7,6 +7,8 @@ If you read this from a place other than ) + ## Minor changes ## Bug fixes and maintenance diff --git a/R/search.R b/R/search.R index e7ba539c7..0b3b84b56 100644 --- a/R/search.R +++ b/R/search.R @@ -248,6 +248,37 @@ search_forward_gauss <- function(p_ref, refmodel, nterms_max, verbose, return(list(crit = -sum(wdraws * quad), r_cn = r_cn, r_n = r_n, y_n = y_n)) } + # crit_cand() for many candidates that each add one column (`cols_single`) + # at once, with matrix operations instead of one call per candidate. For + # candidate j, beta_c = b0 - u_j beta_n[j, ]' with b0 = R^-1 y_chosen and + # u_j = R^-1 r_cj, so the penalty term of the chosen columns is + # sum(lambda b0^2) - 2 beta_n[j, ] p[j, ] + beta_n[j, ]^2 l[j], with + # p = (u * lambda)' b0 and l = colSums(lambda u^2). + crit_single <- function(cols_single, cols_chosen, r_chosen, y_chosen) { + a_jj <- a_full[cbind(cols_single, cols_single)] + rhs <- c_mat[cols_single, , drop = FALSE] + if (length(cols_chosen) > 0) { + r_cj <- forwardsolve(t(r_chosen), + a_full[cols_chosen, cols_single, drop = FALSE]) + a_jj <- a_jj - colSums(r_cj^2) + rhs <- rhs - crossprod(r_cj, y_chosen) + } + y_n <- rhs / sqrt(a_jj) + beta_n <- y_n / sqrt(a_jj) + quad <- y_n^2 + lambda_x[cols_single] * beta_n^2 + if (length(cols_chosen) > 0) { + lambda_c <- lambda_x[cols_chosen] + b0 <- backsolve(r_chosen, y_chosen) + u <- backsolve(r_chosen, r_cj) + q0 <- colSums(y_chosen^2) + colSums(lambda_c * b0^2) + quad <- quad + + matrix(q0, nrow = length(cols_single), ncol = length(q0), byrow = TRUE) - + 2 * beta_n * crossprod(u * lambda_c, b0) + + colSums(lambda_c * u^2) * beta_n^2 + } + return(-drop(quad %*% wdraws)) + } + chosen <- character() labels_chosen <- character() cols_chosen <- integer(0) @@ -264,15 +295,25 @@ search_forward_gauss <- function(p_ref, refmodel, nterms_max, verbose, if (!all(unlist(labels_new) %in% names(cols_by_term))) { return(NULL) } - res_cands <- lapply(labels_new, function(lbls) { - crit_cand(unlist(cols_by_term[lbls], use.names = FALSE), cols_chosen, - r_chosen, y_chosen) + cols_cands <- lapply(labels_new, function(lbls) { + unlist(cols_by_term[lbls], use.names = FALSE) }) - imin <- which.min(vapply(res_cands, `[[`, numeric(1), "crit")) + crits <- rep(NA_real_, length(cands)) + is_single <- lengths(cols_cands) == 1L + if (any(is_single)) { + crits[is_single] <- crit_single(unlist(cols_cands[is_single]), + cols_chosen, r_chosen, y_chosen) + } + for (k in which(!is_single)) { + crits[k] <- crit_cand(cols_cands[[k]], cols_chosen, r_chosen, + y_chosen)$crit + } + imin <- which.min(crits) chosen <- c(chosen, cands[imin]) - cols_new <- unlist(cols_by_term[labels_new[[imin]]], use.names = FALSE) + cols_new <- cols_cands[[imin]] if (length(cols_new) > 0) { - best <- res_cands[[imin]] + # The factors of the chosen candidate, for the next step: + best <- crit_cand(cols_new, cols_chosen, r_chosen, y_chosen) if (length(cols_chosen) == 0) { r_chosen <- best$r_n y_chosen <- best$y_n From f1b025f624ce3d3f4bd40d3640d96dfa68b31ea9 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Sun, 25 Oct 2026 17:14:20 +0200 Subject: [PATCH 4/8] Add the PR number to the NEWS entry --- NEWS.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index f37180a88..3ec30ca86 100644 --- a/NEWS.md +++ b/NEWS.md @@ -7,7 +7,7 @@ If you read this from a place other than ) +* Faster forward search and LOO CV with validated search for Gaussian submodels with the identity link, computed from cross-products of the model matrix. For example, `varsel()` with forward search for a Gaussian model with 50 observations and 200 predictor terms is about 6 times faster. (GitHub: #579) ## Minor changes From 4f19ff7ee166c4bf58138908bbf69ae339dfb59d Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Fri, 23 Oct 2026 14:04:40 +0300 Subject: [PATCH 5/8] Skip the submodel fits along the ranking in the validated search when the evaluation refits --- NEWS.md | 2 +- R/cv_varsel.R | 8 ++++-- R/search.R | 26 ++++++++++++------ R/varsel.R | 5 ++-- tests/testthat/test_search_gauss.R | 44 ++++++++++++++++++++++++++++++ 5 files changed, 72 insertions(+), 13 deletions(-) diff --git a/NEWS.md b/NEWS.md index 3ec30ca86..f7a7cf8ea 100644 --- a/NEWS.md +++ b/NEWS.md @@ -7,7 +7,7 @@ If you read this from a place other than Date: Sun, 25 Oct 2026 22:23:09 +0200 Subject: [PATCH 6/8] Use the Gaussian cross-product computations also for the latent projection --- NEWS.md | 2 +- R/gauss_crossprod.R | 10 +++--- R/projpred-package.R | 7 ++-- man/internal-global-options.Rd | 7 ++-- tests/testthat/test_search_gauss.R | 57 ++++++++++++++++++++++++++++++ 5 files changed, 71 insertions(+), 12 deletions(-) diff --git a/NEWS.md b/NEWS.md index f7a7cf8ea..375f9b0d0 100644 --- a/NEWS.md +++ b/NEWS.md @@ -7,7 +7,7 @@ If you read this from a place other than Date: Wed, 28 Oct 2026 20:04:42 +0200 Subject: [PATCH 7/8] Include a poly() term in the validated LOO tests of the cross-product computations --- tests/testthat/test_search_gauss.R | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/tests/testthat/test_search_gauss.R b/tests/testthat/test_search_gauss.R index 3f0ef7f3e..d2f43a687 100644 --- a/tests/testthat/test_search_gauss.R +++ b/tests/testthat/test_search_gauss.R @@ -55,9 +55,7 @@ test_that(paste( "the LOO CV with validated search gives the same results with and without", "the cross-product computations" ), { - # Without poly(): the draw-wise path cannot predict a poly() term for one - # left-out observation. - refm <- refm_search_gauss(y ~ x1 + x2 + x3 + x4 + f + x1:x2, n = 40) + refm <- refm_search_gauss(fml_search_gauss, n = 40) args_cv <- list(refm, cv_method = "LOO", validate_search = TRUE, nloo = 10, nclusters = 3, ndraws_pred = 20, nterms_max = 4, seed = 5821, verbose = FALSE) @@ -73,9 +71,7 @@ test_that(paste( }) test_that("perf_eval_gauss() applies in the LOO CV with validated search", { - # Without poly(): the draw-wise path cannot predict a poly() term for one - # left-out observation. - refm <- refm_search_gauss(y ~ x1 + x2 + x3 + x4 + f + x1:x2, n = 40) + refm <- refm_search_gauss(fml_search_gauss, n = 40) n_used <- 0L trace("perf_eval_gauss", exit = function() { if (!is.null(returnValue())) n_used <<- n_used + 1L From cde0579a12c71888edfff99ffc2e95f85e9fefbe Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Wed, 28 Oct 2026 20:39:46 +0200 Subject: [PATCH 8/8] Restore the option of the cross-product computations also if a test fails --- tests/testthat/test_search_gauss.R | 3 +++ 1 file changed, 3 insertions(+) diff --git a/tests/testthat/test_search_gauss.R b/tests/testthat/test_search_gauss.R index d2f43a687..98aad939c 100644 --- a/tests/testthat/test_search_gauss.R +++ b/tests/testthat/test_search_gauss.R @@ -36,6 +36,7 @@ test_that(paste( seed = 3317, verbose = FALSE )) opts_old <- options(projpred.search_gauss_crossprod = FALSE) + on.exit(options(opts_old), add = TRUE) vs_drawwise <- suppressWarnings(varsel( refm, nclusters = nclusters_crr, nclusters_pred = 7, nterms_max = 6, seed = 3317, verbose = FALSE @@ -61,6 +62,7 @@ test_that(paste( verbose = FALSE) cv_crossprod <- suppressWarnings(do.call(cv_varsel, args_cv)) opts_old <- options(projpred.search_gauss_crossprod = FALSE) + on.exit(options(opts_old), add = TRUE) cv_drawwise <- suppressWarnings(do.call(cv_varsel, args_cv)) options(opts_old) expect_identical(ranking(cv_crossprod), ranking(cv_drawwise)) @@ -182,6 +184,7 @@ test_that(paste( vs_crossprod <- suppressWarnings(do.call(varsel, args_vs)) cv_crossprod <- suppressWarnings(do.call(cv_varsel, args_cv)) opts_old <- options(projpred.search_gauss_crossprod = FALSE) + on.exit(options(opts_old), add = TRUE) vs_drawwise <- suppressWarnings(do.call(varsel, args_vs)) cv_drawwise <- suppressWarnings(do.call(cv_varsel, args_cv)) options(opts_old)