From 6b44016c36d103d62dd8d056714596f787a88be5 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Thu, 22 Oct 2026 23:16:10 +0300 Subject: [PATCH 1/3] Share one model matrix across the draw-wise GLM submodel fits --- R/divergence_minimizers.R | 190 ++++++++++++++++++++++++++++++-------- 1 file changed, 151 insertions(+), 39 deletions(-) diff --git a/R/divergence_minimizers.R b/R/divergence_minimizers.R index 36345eda8..7e6a3855f 100644 --- a/R/divergence_minimizers.R +++ b/R/divergence_minimizers.R @@ -77,6 +77,15 @@ divmin <- function( rep(projpred_formulas_no_random, length(formulas)) ) } + # The draw-wise formulas of a non-multilevel non-additive submodel differ only + # in the response, so fit_glm_ridge_callback() can use one model matrix for + # all draws: + if (identical(sdivmin, fit_glm_ridge_callback) && length(formulas) > 1) { + shared_x <- shared_model_matrix(formulas, data = list(...)$data) + if (!is.null(shared_x)) { + sdivmin <- glm_ridge_callback_shared(shared_x) + } + } if (is.null(getOption("projpred.parallel_proj_trigger")) && !is.null(getOption("projpred.prll_prj_trigger"))) { @@ -173,28 +182,22 @@ divmin <- function( # submodels: fit_glm_ridge_callback <- function(formula, data, projpred_var = matrix(nrow = nrow(data)), - regul = 1e-4, thresh_conv = 1e-7, ...) { + regul = 1e-4, thresh_conv = 1e-7, + projpred_shared_x = NULL, ...) { # Preparations: - fr <- model.frame(formula, data = data, drop.unused.levels = TRUE) - da_classes <- attr(attr(fr, "terms"), "dataClasses") - nms_chr_fac <- names(da_classes)[da_classes %in% c("character", "factor")] - resp_nm <- all.vars(attr(fr, "terms"))[attr(attr(fr, "terms"), "response")] - nms_chr_fac <- setdiff(nms_chr_fac, resp_nm) - if (length(nms_chr_fac) > 0) { - xlvls <- lapply(setNames(nm = nms_chr_fac), function(nm_chr_fac) { - levels(as.factor(fr[[nm_chr_fac]])) - }) + if (!is.null(projpred_shared_x)) { + # Model matrix built once by shared_model_matrix(); only the response + # differs between the draws: + x <- projpred_shared_x$x + xlvls <- projpred_shared_x$xlvls + y <- data[[all.vars(formula[[2]])]] + names(y) <- projpred_shared_x$y_names } else { - xlvls <- NULL + xy <- glm_model_matrix(formula, data = data) + x <- xy$x + xlvls <- xy$xlvls + y <- xy$y } - # TODO: In the following model.matrix() call, allow user-specified contrasts - # to be passed to argument `contrasts.arg`. The `contrasts.arg` default - # (`NULL`) uses `options("contrasts")` internally, but it might be more - # convenient to let users specify contrasts directly. At that occasion, - # contrasts should also be tested thoroughly (not done until now). - x <- model.matrix(formula, data = fr) - x <- x[, colnames(x) != "(Intercept)", drop = FALSE] - y <- model.response(fr) # Exclude arguments from `...` which cannot be passed to glm_ridge(): dot_args <- list(...) dot_args <- dot_args[intersect( @@ -221,6 +224,72 @@ fit_glm_ridge_callback <- function(formula, data, return(sub) } +# Model frame preparations for fit_glm_ridge_callback(): predictor model matrix +# (without intercept column), levels of character and factor predictors, and +# response. +glm_model_matrix <- function(formula, data) { + fr <- model.frame(formula, data = data, drop.unused.levels = TRUE) + da_classes <- attr(attr(fr, "terms"), "dataClasses") + nms_chr_fac <- names(da_classes)[da_classes %in% c("character", "factor")] + resp_nm <- all.vars(attr(fr, "terms"))[attr(attr(fr, "terms"), "response")] + nms_chr_fac <- setdiff(nms_chr_fac, resp_nm) + if (length(nms_chr_fac) > 0) { + xlvls <- lapply(setNames(nm = nms_chr_fac), function(nm_chr_fac) { + levels(as.factor(fr[[nm_chr_fac]])) + }) + } else { + xlvls <- NULL + } + # TODO: In the following model.matrix() call, allow user-specified contrasts + # to be passed to argument `contrasts.arg`. The `contrasts.arg` default + # (`NULL`) uses `options("contrasts")` internally, but it might be more + # convenient to let users specify contrasts directly. At that occasion, + # contrasts should also be tested thoroughly (not done until now). + x <- model.matrix(formula, data = fr) + x <- x[, colnames(x) != "(Intercept)", drop = FALSE] + return(list(x = x, xlvls = xlvls, y = model.response(fr), n_fr = nrow(fr))) +} + +# Model matrix shared by the draw-wise formulas `formulas` (which have to differ +# only in the response). Returns `NULL` if sharing is not possible: different +# right-hand sides, a response that is not a plain column of `data`, or missing +# values (which would make the model frames differ between the draws). +shared_model_matrix <- function(formulas, data) { + if (!is.data.frame(data)) { + return(NULL) + } + rhs_1 <- formulas[[1]][[3]] + if (!all(vapply(formulas, function(fml) { + length(fml) == 3 && identical(fml[[3]], rhs_1) + }, logical(1)))) { + return(NULL) + } + resp_nms <- vapply(formulas, function(fml) { + if (is.name(fml[[2]])) as.character(fml[[2]]) else NA_character_ + }, character(1)) + if (anyNA(resp_nms) || !all(resp_nms %in% names(data)) || + anyNA(data[resp_nms])) { + return(NULL) + } + xy <- glm_model_matrix(formulas[[1]], data = data) + if (xy$n_fr != nrow(data)) { + return(NULL) + } + # The response names (as from model.response()) are also shared, so the + # draw-wise responses do not each get their own copy: + return(list(x = xy$x, xlvls = xy$xlvls, y_names = names(xy$y))) +} + +# fit_glm_ridge_callback() with a fixed shared model matrix. The enclosing +# environment holds only `shared_x` (relevant for the parallel case in divmin(), +# where the function is exported). +glm_ridge_callback_shared <- function(shared_x) { + force(shared_x) + function(...) { + fit_glm_ridge_callback(..., projpred_shared_x = shared_x) + } +} + # Alternative to fit_glm_ridge_callback() (may be used via global option # `projpred.glm_fitter`): fit_glm_callback <- function(formula, family, ...) { @@ -1159,6 +1228,12 @@ check_conv_s <- function(fit_s) { # Prediction functions for submodels -------------------------------------- subprd <- function(fits, newdata) { + if (!is.null(newdata) && length(fits) > 1 && inherits(fits[[1]], "subfit")) { + prd <- subprd_subfit_shared(fits, newdata) + if (!is.null(prd)) { + return(prd) + } + } prd_list <- lapply(fits, function(fit) { is_glmmPQL <- inherits(fit, "glmmPQL") is_glmm <- inherits(fit, c("lmerMod", "glmerMod")) @@ -1290,28 +1365,65 @@ predict.subfit <- function(object, newdata = NULL, ...) { if (is.null(beta)) { return(as.matrix(rep(alpha, nrow(newdata)))) } else { - # TODO: In the following model.matrix() call, allow user-specified - # contrasts to be passed to argument `contrasts.arg`. The `contrasts.arg` - # default (`NULL`) uses `options("contrasts")` internally, but it might be - # more convenient to let users specify contrasts directly. At that - # occasion, contrasts should also be tested thoroughly (not done until - # now). - x <- model.matrix(delete.response(terms(object$formula)), data = newdata, - xlev = object$xlvls) - # Need to ensure that the columns of `x` match the coefficients (note that - # `rownames(beta)` is the same as `colnames(object$x)`, see the tests): - nms_to_use <- c("(Intercept)", colnames(object$x)) - nms_to_use <- reorder_ias(nms_to_use, colnames(x)) - if (all(!is.na(nms_to_use)) && all(nms_to_use %in% colnames(x))) { - x <- x[, nms_to_use, drop = FALSE] - } else { - stop("The column names of the new model matrix don't match the ", - "column names of the original model matrix. Please notify the ", - "package maintainer.") - } - return(x %*% rbind(alpha, beta)) + return(subfit_newdata_x(object, newdata) %*% rbind(alpha, beta)) + } + } +} + +# Model matrix (with intercept column) of a `subfit` object for new data, with +# the columns in the order of the coefficients. +subfit_newdata_x <- function(object, newdata) { + # TODO: In the following model.matrix() call, allow user-specified + # contrasts to be passed to argument `contrasts.arg`. The `contrasts.arg` + # default (`NULL`) uses `options("contrasts")` internally, but it might be + # more convenient to let users specify contrasts directly. At that + # occasion, contrasts should also be tested thoroughly (not done until + # now). + x <- model.matrix(delete.response(terms(object$formula)), data = newdata, + xlev = object$xlvls) + # Need to ensure that the columns of `x` match the coefficients (note that + # `rownames(beta)` is the same as `colnames(object$x)`, see the tests): + nms_to_use <- c("(Intercept)", colnames(object$x)) + nms_to_use <- reorder_ias(nms_to_use, colnames(x)) + if (all(!is.na(nms_to_use)) && all(nms_to_use %in% colnames(x))) { + x <- x[, nms_to_use, drop = FALSE] + } else { + stop("The column names of the new model matrix don't match the ", + "column names of the original model matrix. Please notify the ", + "package maintainer.") + } + return(x) +} + +# Predictions of several `subfit` objects with the same predictor terms (the +# draw-wise fits of one submodel) for new data: one model matrix, one matrix +# product. Returns `NULL` if the fits do not share their predictor terms. +subprd_subfit_shared <- function(fits, newdata) { + fit_1 <- fits[[1]] + rhs_1 <- fit_1$formula[[3]] + nms_1 <- rownames(fit_1$beta) + is_shared <- all(vapply(fits, function(fit) { + inherits(fit, "subfit") && identical(fit$formula[[3]], rhs_1) && + identical(rownames(fit$beta), nms_1) && identical(fit$xlvls, fit_1$xlvls) + }, logical(1))) + if (!is_shared) { + return(NULL) + } + alphas <- vapply(fits, function(fit) as.numeric(fit$alpha), numeric(1)) + if (is.null(fit_1$beta) || NROW(fit_1$beta) == 0) { + if (is.null(fit_1$beta)) { + # As predict.subfit() for `beta = NULL`: + return(matrix(alphas, nrow = nrow(newdata), ncol = length(fits), + byrow = TRUE)) } + return(NULL) + } + betas <- vapply(fits, function(fit) as.numeric(fit$beta), + numeric(length(nms_1))) + if (!is.matrix(betas)) { + betas <- matrix(betas, nrow = 1) } + return(subfit_newdata_x(fit_1, newdata) %*% rbind(alphas, betas)) } #' @noRd From 3325739fecc23fcd8a994bfe4346dc36ba494ff5 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Sat, 24 Oct 2026 11:23:15 +0300 Subject: [PATCH 2/3] Solve Gaussian identity-link submodel projections for all draws at once --- NEWS.md | 2 + R/divergence_minimizers.R | 71 +++++++++++- tests/testthat/test_submodel_fit_shared.R | 125 ++++++++++++++++++++++ 3 files changed, 196 insertions(+), 2 deletions(-) create mode 100644 tests/testthat/test_submodel_fit_shared.R diff --git a/NEWS.md b/NEWS.md index 94ffe21cd..9ae4aed91 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/divergence_minimizers.R b/R/divergence_minimizers.R index 7e6a3855f..eab961494 100644 --- a/R/divergence_minimizers.R +++ b/R/divergence_minimizers.R @@ -79,10 +79,18 @@ divmin <- function( } # The draw-wise formulas of a non-multilevel non-additive submodel differ only # in the response, so fit_glm_ridge_callback() can use one model matrix for - # all draws: - if (identical(sdivmin, fit_glm_ridge_callback) && length(formulas) > 1) { + # all draws (and the joint Gaussian solve can fit them, also for a single + # draw): + if (identical(sdivmin, fit_glm_ridge_callback) && length(formulas) >= 1) { shared_x <- shared_model_matrix(formulas, data = list(...)$data) if (!is.null(shared_x)) { + # Gaussian submodels with the identity link: one solve for all draws + # (this fitter cannot throw messages or warnings in this case): + outdmin <- glm_ridge_gauss_multi(formulas, shared_x = shared_x, + dot_args = list(...)) + if (!is.null(outdmin)) { + return(outdmin) + } sdivmin <- glm_ridge_callback_shared(shared_x) } } @@ -290,6 +298,65 @@ glm_ridge_callback_shared <- function(shared_x) { } } +# Fits all draw-wise submodels of fit_glm_ridge_callback() at once, for a +# Gaussian submodel with the identity link. For this family, glm_ridge() solves +# one weighted ridge least-squares problem in one exact step. The draws share +# the model matrix, the observation weights, the standardization and the +# penalty. Thus, one factorization with one right-hand side per draw gives the +# same solutions. Returns `NULL` if one of the checks below fails; then divmin() +# uses the draw-wise path. +glm_ridge_gauss_multi <- function(formulas, shared_x, dot_args) { + family <- dot_args$family + if (is.null(family) || family$family != "gaussian" || + family$link != "identity") { + return(NULL) + } + # fit_glm_ridge_callback() passes only the formal arguments of glm_ridge() on + # and ignores the others. Of these, this function uses only `family`, + # `weights` and `thresh` (not used for this family): + args_glm_ridge <- intersect(names(dot_args), methods::formalArgs(glm_ridge)) + if (!all(args_glm_ridge %in% c("family", "weights", "thresh"))) { + return(NULL) + } + x <- shared_x$x + if (ncol(x) == 0) { + return(NULL) + } + data <- dot_args$data + weights <- dot_args$weights %||% rep(1.0, nrow(x)) + regul <- dot_args$regul %||% 1e-4 + if (length(weights) != nrow(x) || anyNA(weights)) { + return(NULL) + } + resp_nms <- vapply(formulas, function(fml) as.character(fml[[2]]), + character(1)) + y <- as.matrix(data[resp_nms]) + # As in glm_ridge(): + transf <- standardization(x, center = TRUE, scale = TRUE, weights = weights) + if (any(transf$scale == 0)) { + return(NULL) + } + x_std <- cbind(1, t((t(x) - transf$shift) / transf$scale)) + sqrt_w <- sqrt(weights) + x_w <- x_std * sqrt_w + regmat <- diag(c(0, rep(regul, ncol(x))), nrow = ncol(x_std)) + coefs <- solve(crossprod(x_w) + regmat, crossprod(x_w, y * sqrt_w)) + beta <- coefs[-1, , drop = FALSE] / transf$scale + beta0 <- coefs[1, ] - colSums(transf$shift * beta) + # As from glm_ridge_c() (double, also for integer observation weights): + w <- matrix(as.numeric(weights), ncol = 1) + return(lapply(seq_along(formulas), function(s) { + beta_s <- beta[, s, drop = FALSE] + dimnames(beta_s) <- list(colnames(x), NULL) + y_s <- y[, s] + names(y_s) <- shared_x$y_names + structure(list(alpha = unname(beta0[s]), beta = beta_s, w = w, + formula = formulas[[s]], x = x, y = y_s, + xlvls = shared_x$xlvls), + class = "subfit") + })) +} + # Alternative to fit_glm_ridge_callback() (may be used via global option # `projpred.glm_fitter`): fit_glm_callback <- function(formula, family, ...) { diff --git a/tests/testthat/test_submodel_fit_shared.R b/tests/testthat/test_submodel_fit_shared.R new file mode 100644 index 000000000..d5560f0d8 --- /dev/null +++ b/tests/testthat/test_submodel_fit_shared.R @@ -0,0 +1,125 @@ +context("shared model matrix and joint Gaussian solve for submodel fits") + +# Data with one response column per projected draw, as divmin() gets them, and +# the draw-wise formulas. The formulas include a factor and an interaction. +shared_fit_data <- function(n = 40, S = 5, seed = 4721) { + set.seed(seed) + dat <- data.frame(x1 = rnorm(n), x2 = rnorm(n), + f = factor(sample(c("a", "b", "c"), n, replace = TRUE))) + for (s in seq_len(S)) { + dat[[paste0("y", s)]] <- 0.5 * dat$x1 - 0.3 * dat$x2 + (dat$f == "b") + + rnorm(n, sd = 0.3) + } + fmls <- lapply(seq_len(S), function(s) { + as.formula(paste0("y", s, " ~ x1 + x2 + f + x1:x2")) + }) + list(data = dat, formulas = fmls, w = runif(n, 0.5, 2)) +} + +test_that(paste( + "the joint Gaussian solve gives the same fits as the draw-wise", + "fit_glm_ridge_callback()" +), { + d <- shared_fit_data() + shared_x <- shared_model_matrix(d$formulas, data = d$data) + expect_false(is.null(shared_x)) + fits_multi <- glm_ridge_gauss_multi( + d$formulas, shared_x = shared_x, + dot_args = list(data = d$data, family = gaussian(), weights = d$w) + ) + expect_length(fits_multi, length(d$formulas)) + fits_loop <- lapply(d$formulas, fit_glm_ridge_callback, data = d$data, + family = gaussian(), weights = d$w) + for (s in seq_along(d$formulas)) { + expect_s3_class(fits_multi[[s]], "subfit") + expect_identical(names(fits_multi[[s]]), names(fits_loop[[s]])) + expect_equal(unname(fits_multi[[s]]$alpha), unname(fits_loop[[s]]$alpha), + tolerance = 1e-10) + expect_equal(fits_multi[[s]]$beta, fits_loop[[s]]$beta, tolerance = 1e-10) + expect_identical(fits_multi[[s]]$x, fits_loop[[s]]$x) + expect_identical(fits_multi[[s]]$y, fits_loop[[s]]$y) + expect_identical(fits_multi[[s]]$xlvls, fits_loop[[s]]$xlvls) + expect_equal(fits_multi[[s]]$w, fits_loop[[s]]$w) + } +}) + +test_that("the joint Gaussian solve falls back for other families", { + d <- shared_fit_data() + shared_x <- shared_model_matrix(d$formulas, data = d$data) + expect_null(glm_ridge_gauss_multi( + d$formulas, shared_x = shared_x, + dot_args = list(data = d$data, family = gaussian(link = "log")) + )) + expect_null(glm_ridge_gauss_multi( + d$formulas, shared_x = shared_x, + dot_args = list(data = d$data, family = poisson()) + )) +}) + +test_that("shared_model_matrix() returns NULL if the model matrix differs", { + d <- shared_fit_data() + # Different right-hand sides: + fmls_rhs <- d$formulas + fmls_rhs[[2]] <- y2 ~ x1 + x2 + expect_null(shared_model_matrix(fmls_rhs, data = d$data)) + # A response that is not a plain column of the data: + fmls_resp <- d$formulas + fmls_resp[[1]] <- log(abs(y1)) ~ x1 + x2 + f + x1:x2 + expect_null(shared_model_matrix(fmls_resp, data = d$data)) + # A missing predictor value: + dat_na <- d$data + dat_na$x2[3] <- NA + expect_null(shared_model_matrix(d$formulas, data = dat_na)) +}) + +test_that(paste( + "the joint predictions for new data equal the draw-wise predictions of", + "predict.subfit()" +), { + d <- shared_fit_data() + shared_x <- shared_model_matrix(d$formulas, data = d$data) + fits <- glm_ridge_gauss_multi( + d$formulas, shared_x = shared_x, + dot_args = list(data = d$data, family = gaussian(), weights = d$w) + ) + nd <- d$data[c(2, 5, 7, 11), ] + prd_loop <- do.call(cbind, lapply(fits, predict, newdata = nd)) + prd_shared <- subprd_subfit_shared(fits, nd) + expect_equal(unname(prd_shared), unname(prd_loop), tolerance = 1e-12) + expect_equal(unname(subprd(fits, nd)), unname(prd_loop), tolerance = 1e-12) +}) + +test_that("project() uses the joint Gaussian solve for a Gaussian submodel", { + d <- shared_fit_data(n = 60) + dat <- d$data[, c("x1", "x2", "f")] + dat$y <- d$data$y1 + fml <- y ~ x1 + x2 + f + mm <- model.matrix(fml, dat) + ndraws <- 100 + beta <- drop(solve(crossprod(mm), crossprod(mm, dat$y))) + + matrix(rnorm(ncol(mm) * ndraws, sd = 0.05), ncol(mm)) + fit <- list(beta = beta, data = dat) + ref_predfun <- function(fit, newdata = NULL) { + nd <- if (is.null(newdata)) fit$data else newdata + model.matrix(~ x1 + x2 + f, nd) %*% fit$beta + } + refm <- suppressWarnings(init_refmodel( + fit, data = dat, formula = fml, family = gaussian(), + ref_predfun = ref_predfun, dis = rep(0.3, ndraws) + )) + calls <- new.env() + calls$n <- 0L + ns <- environment(glm_ridge_gauss_multi) + suppressMessages(trace( + "glm_ridge_gauss_multi", where = ns, print = FALSE, + exit = bquote(if (!is.null(returnValue(NULL))) { + assign("n", get("n", envir = .(calls)) + 1L, envir = .(calls)) + }) + )) + on.exit(suppressMessages(untrace("glm_ridge_gauss_multi", where = ns)), + add = TRUE) + prj <- project(refm, predictor_terms = c("x1", "f"), nclusters = 5, + seed = 3906, verbose = FALSE) + expect_gt(calls$n, 0L) + expect_length(prj$outdmin, 5) +}) From 3eebda0c2f7ffc19036bbde801fbb3a7675fe7c9 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Sat, 24 Oct 2026 13:11:26 +0300 Subject: [PATCH 3/3] 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 9ae4aed91..96fc28373 100644 --- a/NEWS.md +++ b/NEWS.md @@ -7,7 +7,7 @@ If you read this from a place other than ) +* 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) ## Minor changes