From c29bbb385301d336a864e1b08dfb07a815505560 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Sun, 25 Oct 2026 23:28:33 +0200 Subject: [PATCH 1/3] Faster augmented-data projection for cumulative-link submodels: C++ fits for all draws, predictions from the coefficients and a vectorized cross-entropy --- NEWS.md | 2 + R/RcppExports.R | 4 + R/divergence_minimizers.R | 193 ++++++++++++++++++ R/extend_family.R | 34 ++- src/RcppExports.cpp | 18 ++ src/cumulfit.cpp | 295 +++++++++++++++++++++++++++ tests/testthat/test_augdat_predict.R | 115 +++++++++++ tests/testthat/test_cumul_native.R | 128 ++++++++++++ 8 files changed, 771 insertions(+), 18 deletions(-) create mode 100644 src/cumulfit.cpp create mode 100644 tests/testthat/test_augdat_predict.R create mode 100644 tests/testthat/test_cumul_native.R diff --git a/NEWS.md b/NEWS.md index 94ffe21cd..b5ca7bd4f 100644 --- a/NEWS.md +++ b/NEWS.md @@ -7,6 +7,8 @@ If you read this from a place other than do not edit by hand # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 +cumul_fit_multi_c <- function(x, w_all, ncat, offset, link, thresh, maxit) { + .Call('_projpred_cumul_fit_multi_c', PACKAGE = 'projpred', x, w_all, ncat, offset, link, thresh, maxit) +} + glm_elnet_c <- function(x, pseudo_obs, lambda, alpha, intercept, penalty, thresh, qa_updates_max, pmax, pmax_strict, beta, beta0, w0, as_updates_max = 50L) { .Call('_projpred_glm_elnet_c', PACKAGE = 'projpred', x, pseudo_obs, lambda, alpha, intercept, penalty, thresh, qa_updates_max, pmax, pmax_strict, beta, beta0, w0, as_updates_max) } diff --git a/R/divergence_minimizers.R b/R/divergence_minimizers.R index 36345eda8..9bb59015d 100644 --- a/R/divergence_minimizers.R +++ b/R/divergence_minimizers.R @@ -694,6 +694,18 @@ divmin_augdat <- function( projpred_ws_aug <- unclass(projpred_ws_aug) attr(projpred_ws_aug, "ndiscrete") <- NULL + # Cumulative-link submodels without group-level terms: all draws in one C++ + # call (see cumul_native_multi()): + if (identical(sdivmin, fit_cumul)) { + outdmin <- cumul_native_multi(formula, data = data, family = family, + ws_aug = projpred_ws_aug, + dot_args = list(...)) + if (!is.null(outdmin)) { + check_conv(outdmin, rep(0L, length(outdmin)), do_check = do_check_conv) + return(outdmin) + } + } + if (is.null(getOption("projpred.parallel_proj_trigger")) && !is.null(getOption("projpred.prll_prj_trigger"))) { warning( @@ -868,6 +880,98 @@ fit_cumul <- function(formula, data, family, weights, ...) { return(fitobj) } +# Fits the cumulative-link submodels of all draws in one C++ call +# (cumul_fit_multi_c()) with the objective of MASS::polr(), instead of one +# fit_cumul() call per draw. The augmented data hold C copies of the n +# observations (one per response category), so this builds the model matrix +# once, and the augmented weights of each draw form an n x C matrix. Newton +# iterations replace the BFGS iterations of polr(), so the estimates agree with +# those of polr() up to its convergence tolerance. The returned `polr` objects +# hold the elements that predict.polr() (with `newdata`), coef() and +# check_conv_s() use, without `fitted.values` and `lp`. Returns `NULL` for other +# cases; divmin_augdat() then calls fit_cumul() per draw. Global option +# `projpred.augdat_native = FALSE` switches this off. +cumul_native_multi <- function(formula, data, family, ws_aug, dot_args) { + if (!getOption("projpred.augdat_native", TRUE)) { + return(NULL) + } + # fit_cumul() passes these arguments on to MASS::polr(); they are not + # reproduced here: + if (length(intersect(names(dot_args), + c(methods::formalArgs(MASS::polr), + methods::formalArgs(stats::optim))))) { + return(NULL) + } + method <- family$link + if (method == "logit") { + method <- "logistic" + } else if (method == "probit_approx") { + method <- "probit" + } + link_code <- match(method, c("logistic", "probit", "loglog", "cloglog", + "cauchit")) + if (is.na(link_code)) { + return(NULL) + } + # As in MASS::polr(): + mfr <- model.frame(formula, data = data) + if (nrow(mfr) != nrow(data) || !is.null(model.offset(mfr))) { + return(NULL) + } + trms <- attr(mfr, "terms") + x <- model.matrix(trms, mfr) + cons <- attr(x, "contrasts") + xint <- match("(Intercept)", colnames(x), nomatch = 0L) + if (xint == 0L) { + return(NULL) + } + x <- x[, -xint, drop = FALSE] + y <- model.response(mfr) + if (!is.factor(y)) { + return(NULL) + } + lev <- levels(y) + ncat <- length(lev) + n <- nrow(x) %/% ncat + ws_aug <- as.matrix(ws_aug) + # The augmented data have to be C copies of the n observations, in the order + # of the response categories: + if (ncat < 3L || n * ncat != nrow(x) || nrow(ws_aug) != nrow(x) || + !identical(as.integer(y), rep(seq_len(ncat), each = n)) || + anyNA(ws_aug) || any(ws_aug < 0)) { + return(NULL) + } + x1 <- x[seq_len(n), , drop = FALSE] + if (ncol(x1) > 0 && !all(x == x1[rep(seq_len(n), ncat), , drop = FALSE])) { + return(NULL) + } + nprj <- ncol(ws_aug) + out <- cumul_fit_multi_c( + unname(x1), matrix(ws_aug, nrow = n), ncat, rep(0, n), link_code, + thresh = 1e-12, maxit = 100L + ) + xlvls <- .getXlevels(trms, mfr) + zeta_nms <- paste(lev[-ncat], lev[-1L], sep = "|") + pc <- ncol(x1) + # As in polr(), also for a submodel without predictors: + beta_nms <- colnames(x1) %||% character() + return(lapply(seq_len(nprj), function(s) { + beta <- setNames(out$beta[, s], beta_nms) + zeta <- setNames(out$zeta[, s], zeta_nms) + wsum <- sum(ws_aug[, s]) + structure( + list(coefficients = beta, zeta = zeta, deviance = 2 * out$loss[s], + lev = lev, terms = trms, df.residual = wsum - pc - (ncat - 1L), + edf = pc + ncat - 1L, n = wsum, nobs = wsum, + call = quote(projpred_cumul_fit()), method = method, + convergence = out$convergence[s], + niter = c(f.evals = out$niter[s], g.evals = out$niter[s]), + contrasts = cons, xlevels = xlvls), + class = "polr" + ) + })) +} + # Use ordinal::clmm() to fit multilevel submodels for the brms::cumulative() # family: fit_cumul_mlvl <- function(formula, data, family, weights, ...) { @@ -1214,6 +1318,12 @@ subprd <- function(fits, newdata) { } subprd_augdat <- function(fits, newdata) { + # `polr` and `multinom` fits of the same submodel: linear predictors from the + # coefficients, for all draws at once (see subprd_augdat_direct()): + prd_direct <- subprd_augdat_direct(fits, newdata) + if (!is.null(prd_direct)) { + return(prd_direct) + } if (inherits(fits[[1]], "clmm")) { y_nm <- extract_terms_response(formula(fits[[1]]))$response } @@ -1269,6 +1379,89 @@ subprd_augdat <- function(fits, newdata) { return(do.call(abind::abind, list(prd_list, rev.along = 0))) } +# Linear predictors (on the link scale used by subprd_augdat()) of `polr` or +# `multinom` fits that share the terms, the factor levels and the contrasts (as +# the draw-wise fits of one submodel do). The model matrix for `newdata` is +# built once, as in predict.polr() and predict.multinom(). For `polr`, the +# linear predictors of the cumulative probabilities are zeta_j - eta, which +# subprd_augdat() computes from the predicted category probabilities; for +# `multinom`, they are the log odds against the first category. This function +# computes both directly, so they differ from those of subprd_augdat() only by +# rounding. Returns `NULL` for other cases. +subprd_augdat_direct <- function(fits, newdata) { + fit_1 <- fits[[1]] + if (inherits(fit_1, "polr")) { + is_polr <- TRUE + } else if (inherits(fit_1, "multinom") && is.matrix(coef(fit_1))) { + is_polr <- FALSE + } else { + return(NULL) + } + same_cls <- function(fit) { + if (is_polr) { + inherits(fit, "polr") && !inherits(fit, "clmm") + } else { + inherits(fit, "multinom") && is.matrix(coef(fit)) + } + } + # The terms are compared without the environment of their formula: + terms_noenv <- function(trms) { + attr(trms, ".Environment") <- NULL + trms + } + trms_1 <- terms_noenv(fit_1$terms) + if (!all(vapply(fits, function(fit) { + same_cls(fit) && identical(terms_noenv(fit$terms), trms_1) && + identical(fit$xlevels, fit_1$xlevels) && + identical(fit$contrasts, fit_1$contrasts) + }, logical(1)))) { + return(NULL) + } + trms <- delete.response(fit_1$terms) + mfr <- model.frame(trms, newdata, na.action = function(x) x, + xlev = fit_1$xlevels) + if (!is.null(model.offset(mfr))) { + return(NULL) + } + x <- model.matrix(trms, mfr, contrasts.arg = fit_1$contrasts) + n <- nrow(x) + nprj <- length(fits) + if (is_polr) { + xint <- match("(Intercept)", colnames(x), nomatch = 0L) + if (xint > 0L) { + x <- x[, -xint, drop = FALSE] + } + if (!identical(colnames(x) %||% character(), + names(coef(fit_1)) %||% character())) { + return(NULL) + } + zetas <- do.call(cbind, lapply(fits, "[[", "zeta")) + if (ncol(x) > 0) { + etas <- x %*% do.call(cbind, lapply(fits, coef)) + } else { + etas <- matrix(0, n, nprj) + } + ncat_m1 <- nrow(zetas) + out <- array(NA_real_, dim = c(n, ncat_m1, nprj)) + for (s in seq_len(nprj)) { + out[, , s] <- outer(-etas[, s], zetas[, s], "+") + } + cat_nms <- utils::head(fit_1$lev, -1) + } else { + coef_1 <- coef(fit_1) + if (!identical(colnames(x), colnames(coef_1))) { + return(NULL) + } + out <- array(NA_real_, dim = c(n, nrow(coef_1), nprj)) + for (s in seq_len(nprj)) { + out[, , s] <- x %*% t(coef(fits[[s]])) + } + cat_nms <- rownames(coef_1) + } + dimnames(out) <- list(rownames(newdata), cat_nms, NULL) + return(out) +} + subprd_augdat_binom <- function(fits, newdata) { augprd <- subprd(fits, newdata = newdata) return(array(augprd, dim = c(nrow(augprd), 1L, ncol(augprd)))) diff --git a/R/extend_family.R b/R/extend_family.R index e3f9c3999..cc026d993 100644 --- a/R/extend_family.R +++ b/R/extend_family.R @@ -420,24 +420,22 @@ extend_family <- function(family, } else if (length(w_obs) != n_obs) { stop("Argument `w_obs` needs to be of length 0, 1, or `n_obs`.") } - return(do.call(rbind, lapply(seq_len(n_obs), function(i_obs) { - do.call(c, lapply(seq_len(n_prjdraws), function(i_prjdraws) { - prbs_ref <- mu_ref_arr[i_obs, , i_prjdraws] - prbs_sub <- mu_sub_arr[i_obs, , i_prjdraws] - # Assign some nonzero value to have a finite log() value (the specific - # value doesn't matter for `prbs_ref` because of the multiplication - # with zero): - prbs_ref[prbs_ref == 0] <- .Machine$double.eps - # Assign some nonzero value to have a finite log() value (for - # `prbs_sub`, the specific value matters, in contrast to `prbs_ref`): - prbs_sub[prbs_sub == 0] <- .Machine$double.eps - # In analogy to `extend_family(binomial())$ce()`, - # `extend_family(poisson())$ce()`, and - # `extend_family(gaussian())$ce()`, multiply by the observation - # weight: - return(w_obs[i_obs] * (-sum(prbs_ref * log(prbs_sub)))) - })) - }))) + # Assign some nonzero value to have a finite log() value (the specific + # value doesn't matter for `mu_ref_arr` because of the multiplication with + # zero): + mu_ref_arr[mu_ref_arr == 0] <- .Machine$double.eps + # Assign some nonzero value to have a finite log() value (for + # `mu_sub_arr`, the specific value matters, in contrast to `mu_ref_arr`): + mu_sub_arr[mu_sub_arr == 0] <- .Machine$double.eps + # Sum over the response categories (the second dimension) for each + # observation and draw: rowSums() of the observations x draws x + # categories array: + ce_unw <- -rowSums(aperm(mu_ref_arr * log(mu_sub_arr), c(1, 3, 2)), + dims = 2) + # In analogy to `extend_family(binomial())$ce()`, + # `extend_family(poisson())$ce()`, and `extend_family(gaussian())$ce()`, + # multiply by the observation weight: + return(matrix(w_obs * ce_unw, nrow = n_obs, ncol = n_prjdraws)) } family$ce <- function(pref, data, psub) { ce_sums <- colSums( diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index c85e42fe8..9bf3972dc 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -11,6 +11,23 @@ Rcpp::Rostream& Rcpp::Rcout = Rcpp::Rcpp_cout_get(); Rcpp::Rostream& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get(); #endif +// cumul_fit_multi_c +List cumul_fit_multi_c(arma::mat x, arma::mat w_all, int ncat, arma::vec offset, int link, double thresh, int maxit); +RcppExport SEXP _projpred_cumul_fit_multi_c(SEXP xSEXP, SEXP w_allSEXP, SEXP ncatSEXP, SEXP offsetSEXP, SEXP linkSEXP, SEXP threshSEXP, SEXP maxitSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< arma::mat >::type x(xSEXP); + Rcpp::traits::input_parameter< arma::mat >::type w_all(w_allSEXP); + Rcpp::traits::input_parameter< int >::type ncat(ncatSEXP); + Rcpp::traits::input_parameter< arma::vec >::type offset(offsetSEXP); + Rcpp::traits::input_parameter< int >::type link(linkSEXP); + Rcpp::traits::input_parameter< double >::type thresh(threshSEXP); + Rcpp::traits::input_parameter< int >::type maxit(maxitSEXP); + rcpp_result_gen = Rcpp::wrap(cumul_fit_multi_c(x, w_all, ncat, offset, link, thresh, maxit)); + return rcpp_result_gen; +END_RCPP +} // glm_elnet_c List glm_elnet_c(arma::mat x, Function pseudo_obs, arma::vec lambda, double alpha, bool intercept, arma::vec penalty, double thresh, int qa_updates_max, int pmax, bool pmax_strict, arma::vec beta, double beta0, arma::vec w0, int as_updates_max); RcppExport SEXP _projpred_glm_elnet_c(SEXP xSEXP, SEXP pseudo_obsSEXP, SEXP lambdaSEXP, SEXP alphaSEXP, SEXP interceptSEXP, SEXP penaltySEXP, SEXP threshSEXP, SEXP qa_updates_maxSEXP, SEXP pmaxSEXP, SEXP pmax_strictSEXP, SEXP betaSEXP, SEXP beta0SEXP, SEXP w0SEXP, SEXP as_updates_maxSEXP) { @@ -58,6 +75,7 @@ END_RCPP } static const R_CallMethodDef CallEntries[] = { + {"_projpred_cumul_fit_multi_c", (DL_FUNC) &_projpred_cumul_fit_multi_c, 7}, {"_projpred_glm_elnet_c", (DL_FUNC) &_projpred_glm_elnet_c, 14}, {"_projpred_glm_ridge_c", (DL_FUNC) &_projpred_glm_ridge_c, 11}, {NULL, NULL, 0} diff --git a/src/cumulfit.cpp b/src/cumulfit.cpp new file mode 100644 index 000000000..3386e42ce --- /dev/null +++ b/src/cumulfit.cpp @@ -0,0 +1,295 @@ +#include +#include +#include + +//[[Rcpp::depends(RcppArmadillo)]] + +using namespace Rcpp; +using namespace arma; + +// Cumulative-link model fits for the augmented-data projection, with the +// objective of MASS::polr(): for observation i with linear predictor eta_i and +// category weights w_ic, +// L(beta, zeta) = -sum_i sum_c w_ic log(F(z1_ic) - F(z2_ic)), +// z1_ic = min(100, zeta_c - eta_i), z2_ic = max(-100, zeta_{c-1} - eta_i), +// with zeta_0 = -Inf and zeta_C = Inf. The objective, with its bounds of +-100, +// and the latent distributions transcribe the negative log-likelihood and the +// `pfun` / `dfun` of MASS::polr() (R package MASS, GPL-2 | GPL-3). The fit is +// new: polr() minimizes with optim() (BFGS), and this code uses Newton steps +// with the analytic gradient and Hessian, and a backtracking line search that +// keeps the thresholds increasing. + +// Link codes (as polr()'s `method`): +enum { CL_LOGISTIC = 1, CL_PROBIT = 2, CL_LOGLOG = 3, CL_CLOGLOG = 4, + CL_CAUCHIT = 5 }; + +// Distribution function F, survival function S = 1 - F (computed directly, so +// that upper-tail probabilities do not cancel), density f and derivative of the +// density f' of the latent variable, as in polr()'s `pfun` and `dfun`: +static inline void cl_funs(int link, double z, double& F, double& S, + double& f, double& fd) { + switch (link) { + case CL_LOGISTIC: { + F = R::plogis(z, 0.0, 1.0, 1, 0); + S = R::plogis(z, 0.0, 1.0, 0, 0); + f = R::dlogis(z, 0.0, 1.0, 0); + fd = f * (S - F); + break; + } + case CL_PROBIT: { + F = R::pnorm(z, 0.0, 1.0, 1, 0); + S = R::pnorm(z, 0.0, 1.0, 0, 0); + f = R::dnorm(z, 0.0, 1.0, 0); + fd = -z * f; + break; + } + case CL_LOGLOG: { + // MASS::pgumbel(): exp(-exp(-z)) + double e = std::exp(-z); + F = std::exp(-e); + S = -std::expm1(-e); + f = std::exp(-z - e); + fd = f * (e - 1.0); + break; + } + case CL_CLOGLOG: { + // MASS:::pGumbel(): 1 - exp(-exp(z)) + double e = std::exp(z); + F = -std::expm1(-e); + S = std::exp(-e); + f = std::exp(z - e); + fd = f * (1.0 - e); + break; + } + case CL_CAUCHIT: { + F = R::pcauchy(z, 0.0, 1.0, 1, 0); + S = R::pcauchy(z, 0.0, 1.0, 0, 0); + f = R::dcauchy(z, 0.0, 1.0, 0); + fd = -2.0 * z * M_PI * f * f; + break; + } + default: + Rcpp::stop("Unknown link code."); + } +} + +static inline double cl_quantile(int link, double p) { + switch (link) { + case CL_LOGISTIC: return R::qlogis(p, 0.0, 1.0, 1, 0); + case CL_PROBIT: return R::qnorm(p, 0.0, 1.0, 1, 0); + case CL_LOGLOG: return -std::log(-std::log(p)); + case CL_CLOGLOG: return std::log(-std::log1p(-p)); + case CL_CAUCHIT: return R::qcauchy(p, 0.0, 1.0, 1, 0); + default: Rcpp::stop("Unknown link code."); + } + return 0.0; +} + +// Objective (and optionally gradient and Hessian) at (beta, zeta). Returns +// Inf if a category probability with positive weight is not positive. +static double cl_eval(const mat& x, const mat& w, const vec& offset, + const vec& beta, const vec& zeta, int link, + bool deriv, vec& grad, mat& hess) { + const uword n = w.n_rows; + const uword ncat = w.n_cols; + const uword q = ncat - 1; + const uword p = x.n_cols; + vec eta = offset; + if (p > 0) eta += x * beta; + vec g_eta, h_etaeta; + mat h_etazeta; + vec g_zeta; + mat h_zetazeta; + if (deriv) { + g_eta.zeros(n); + h_etaeta.zeros(n); + h_etazeta.zeros(n, q); + g_zeta.zeros(q); + h_zetazeta.zeros(q, q); + } + double loss = 0.0; + for (uword i = 0; i < n; ++i) { + for (uword c = 0; c < ncat; ++c) { + double wic = w(i, c); + if (wic == 0.0) continue; + // Upper and lower threshold of category c (0-based): zeta_c and + // zeta_{c-1}: + double z1 = (c < q) ? std::min(100.0, zeta(c) - eta(i)) : 100.0; + double z2 = (c > 0) ? std::max(-100.0, zeta(c - 1) - eta(i)) : -100.0; + double F1, S1, f1, fd1, F2, S2, f2, fd2; + cl_funs(link, z1, F1, S1, f1, fd1); + cl_funs(link, z2, F2, S2, f2, fd2); + if (c == q) { + f1 = 0.0; + fd1 = 0.0; + } + if (c == 0) { + f2 = 0.0; + fd2 = 0.0; + } + // F(z1) - F(z2), from the upper tails if both arguments are positive: + double pr = (z2 > 0.0) ? (S2 - S1) : (F1 - F2); + if (!(pr > 0.0)) return std::numeric_limits::infinity(); + loss -= wic * std::log(pr); + if (!deriv) continue; + double a = f1 / pr, b = f2 / pr; + double ad = fd1 / pr, bd = fd2 / pr; + // Derivatives of -w * log(pr): + g_eta(i) -= wic * (-(a - b)); + h_etaeta(i) -= wic * ((ad - bd) - (a - b) * (a - b)); + if (c < q) { + g_zeta(c) -= wic * a; + h_zetazeta(c, c) -= wic * (ad - a * a); + h_etazeta(i, c) -= wic * (-ad + a * (a - b)); + } + if (c > 0) { + g_zeta(c - 1) -= wic * (-b); + h_zetazeta(c - 1, c - 1) -= wic * (-bd - b * b); + h_etazeta(i, c - 1) -= wic * (bd - b * (a - b)); + } + if (c > 0 && c < q) { + h_zetazeta(c, c - 1) -= wic * (a * b); + h_zetazeta(c - 1, c) -= wic * (a * b); + } + } + } + if (deriv) { + grad.set_size(p + q); + hess.set_size(p + q, p + q); + if (p > 0) { + grad.head(p) = x.t() * g_eta; + hess.submat(0, 0, p - 1, p - 1) = x.t() * (x.each_col() % h_etaeta); + mat hbz = x.t() * h_etazeta; + hess.submat(0, p, p - 1, p + q - 1) = hbz; + hess.submat(p, 0, p + q - 1, p - 1) = hbz.t(); + } + grad.tail(q) = g_zeta; + hess.submat(p, p, p + q - 1, p + q - 1) = h_zetazeta; + } + return loss; +} + +// Fits one weight matrix `w` (n x C). Writes the estimates into `beta` and +// `zeta`; returns the objective. `conv` is 0 on convergence. +static double cl_fit_one(const mat& x, const mat& w, const vec& offset, + int link, double thresh, int maxit, vec& beta, + vec& zeta, int& conv, int& niter) { + const uword p = x.n_cols; + const uword ncat = w.n_cols; + const uword q = ncat - 1; + // Start: beta = 0 and the thresholds of the model without predictors: + beta.zeros(p); + zeta.set_size(q); + rowvec wc = sum(w, 0); + double wtot = accu(wc); + double cum = 0.0; + for (uword c = 0; c < q; ++c) { + cum += wc(c); + double pc = std::min(std::max(cum / wtot, 1e-8), 1.0 - 1e-8); + zeta(c) = cl_quantile(link, pc); + if (c > 0 && zeta(c) <= zeta(c - 1)) zeta(c) = zeta(c - 1) + 1e-6; + } + vec grad; + mat hess; + double loss = cl_eval(x, w, offset, beta, zeta, link, true, grad, hess); + conv = 1; + niter = 0; + vec theta(p + q), dir, theta_new; + for (int it = 0; it < maxit; ++it) { + niter = it + 1; + // Newton direction (damped if the Hessian is not positive definite): + mat hreg = hess; + double damp = 0.0; + bool ok = false; + double dscale = std::max(1e-10, max(abs(hess.diag()))); + for (int k = 0; k < 30 && !ok; ++k) { + mat r; + if (damp > 0.0) hreg = hess + damp * eye(p + q, p + q); + ok = chol(r, hreg); + if (ok) { + dir = -solve(trimatu(r), solve(trimatl(r.t()), grad)); + } else { + damp = (damp == 0.0) ? 1e-8 * dscale : 10.0 * damp; + } + } + if (!ok) break; + double decrement = -dot(grad, dir); + if (decrement < thresh * std::max(1.0, std::abs(loss))) { + conv = 0; + break; + } + if (p > 0) theta.head(p) = beta; + theta.tail(q) = zeta; + // Backtracking line search (Armijo), keeping the thresholds increasing: + double t = 1.0; + bool accepted = false; + vec grad_new; + mat hess_new; + for (int ls = 0; ls < 60; ++ls) { + theta_new = theta + t * dir; + vec zeta_new = theta_new.tail(q); + bool ordered = true; + for (uword c = 1; c < q; ++c) { + if (!(zeta_new(c) > zeta_new(c - 1))) { + ordered = false; + break; + } + } + if (ordered) { + vec beta_new = (p > 0) ? vec(theta_new.head(p)) : vec(); + double loss_new = cl_eval(x, w, offset, beta_new, zeta_new, link, + false, grad_new, hess_new); + if (std::isfinite(loss_new) && + loss_new <= loss - 1e-4 * t * decrement) { + accepted = true; + break; + } + } + t *= 0.5; + } + if (!accepted) { + // No decrease possible along the direction: at the optimum up to + // rounding if the decrement is tiny, otherwise not converged. + if (decrement < 1e3 * thresh * std::max(1.0, std::abs(loss))) conv = 0; + break; + } + if (p > 0) beta = theta_new.head(p); + zeta = theta_new.tail(q); + loss = cl_eval(x, w, offset, beta, zeta, link, true, grad, hess); + } + return loss; +} + +// Fits the cumulative-link model for S weight matrices, given as the n x C +// column blocks of `w_all` (n x (C * S)). The model matrix `x` (n x p, without +// intercept) and the offsets are shared. +// [[Rcpp::export]] +List cumul_fit_multi_c(arma::mat x, arma::mat w_all, int ncat, + arma::vec offset, int link, double thresh, + int maxit) { + const uword n = x.n_rows; + const uword p = x.n_cols; + const uword q = ncat - 1; + const uword S = w_all.n_cols / ncat; + if (w_all.n_rows != n || S * ncat != w_all.n_cols) { + Rcpp::stop("Dimensions of `w_all` do not match."); + } + mat beta_out(p, S); + mat zeta_out(q, S); + vec loss_out(S); + IntegerVector conv_out(S), niter_out(S); + for (uword s = 0; s < S; ++s) { + mat w = w_all.cols(s * ncat, (s + 1) * ncat - 1); + vec beta, zeta; + int conv, niter; + loss_out(s) = cl_fit_one(x, w, offset, link, thresh, maxit, beta, zeta, + conv, niter); + if (p > 0) beta_out.col(s) = beta; + zeta_out.col(s) = zeta; + conv_out[s] = conv; + niter_out[s] = niter; + } + return List::create(Named("beta") = beta_out, Named("zeta") = zeta_out, + Named("loss") = loss_out, Named("convergence") = conv_out, + Named("niter") = niter_out); +} diff --git a/tests/testthat/test_augdat_predict.R b/tests/testthat/test_augdat_predict.R new file mode 100644 index 000000000..4b8b2e749 --- /dev/null +++ b/tests/testthat/test_augdat_predict.R @@ -0,0 +1,115 @@ +context("augdat_predict") + +# Needed to clean up the workspace afterwards (i.e, after this test file): +ls_bu <- ls() + +if (exists(".Random.seed", envir = .GlobalEnv)) { + rng_old <- get(".Random.seed", envir = .GlobalEnv) +} +set.seed(97531) +n_ap <- 120 +dat_ap <- data.frame(a = rnorm(n_ap), b = rnorm(n_ap), + f = factor(sample(c("u", "v", "w"), n_ap, replace = TRUE))) +eta_ap <- with(dat_ap, 0.8 * a - 0.5 * b + 0.6 * (f == "v")) +dat_ap$yo <- cut(eta_ap + rlogis(n_ap), c(-Inf, -1, 0, 1, Inf), + labels = paste0("c", 1:4), ordered_result = TRUE) +dat_ap$yc <- factor(apply(cbind(0, eta_ap, -eta_ap, 0.3), 1, function(e) { + sample(4, 1, prob = exp(e)) +})) +newdat_ap <- data.frame(a = c(-1, 0.5, 2), b = c(0, 1, -1), + f = factor(c("w", "u", "v"), levels = c("u", "v", "w"))) +rownames(newdat_ap) <- c("r1", "r2", "r3") + +# subprd_augdat() without subprd_augdat_direct(), as on master: +lpreds_via_probs <- function(fits, newdata) { + arrs <- lapply(fits, function(fit) { + prbs <- predict(fit, newdata = newdata, type = "probs") + if (inherits(fit, "polr")) { + cprbs <- t(apply(prbs, 1, cumsum))[, -ncol(prbs), drop = FALSE] + linkfun_raw(cprbs, link_nm = fit$method) + } else { + log(sweep(prbs[, -1, drop = FALSE], 1, prbs[, 1], "/")) + } + }) + do.call(abind::abind, list(arrs, rev.along = 0)) +} + +test_that("subprd_augdat_direct() gives the linear predictors of polr fits", { + for (method in c("logistic", "probit", "cloglog")) { + fits <- lapply(1:3, function(s) { + wt <- runif(n_ap, 0.5, 1.5) + suppressWarnings(MASS::polr(yo ~ a + b + f, data = cbind(dat_ap, wt = wt), + weights = wt, method = method)) + }) + prd <- subprd_augdat_direct(fits, newdat_ap) + expect_equal(prd, lpreds_via_probs(fits, newdat_ap), tolerance = 1e-10, + info = method) + expect_identical(dimnames(prd), list(rownames(newdat_ap), c("c1", "c2", "c3"), + NULL)) + } + # Without predictors: + fits_icpt <- lapply(1:2, function(s) MASS::polr(yo ~ 1, data = dat_ap)) + expect_equal(subprd_augdat_direct(fits_icpt, newdat_ap), + lpreds_via_probs(fits_icpt, newdat_ap), tolerance = 1e-10) +}) + +test_that("subprd_augdat_direct() gives the linear predictors of multinom fits", { + fits <- lapply(1:3, function(s) { + wt <- runif(n_ap, 0.5, 1.5) + nnet::multinom(yc ~ a + b + f, data = cbind(dat_ap, wt = wt), weights = wt, + trace = FALSE) + }) + prd <- subprd_augdat_direct(fits, newdat_ap) + expect_equal(prd, lpreds_via_probs(fits, newdat_ap), tolerance = 1e-10) + expect_identical(dimnames(prd), list(rownames(newdat_ap), c("2", "3", "4"), + NULL)) +}) + +test_that("subprd_augdat_direct() returns `NULL` for other cases", { + fit_p1 <- MASS::polr(yo ~ a, data = dat_ap) + fit_p2 <- MASS::polr(yo ~ b, data = dat_ap) + fit_m <- nnet::multinom(yc ~ a, data = dat_ap, trace = FALSE) + expect_null(subprd_augdat_direct(list(fit_p1, fit_p2), newdat_ap)) + expect_null(subprd_augdat_direct(list(fit_p1, fit_m), newdat_ap)) + expect_null(subprd_augdat_direct(list(lm(a ~ b, data = dat_ap)), newdat_ap)) +}) + +test_that("`ce_ptwise` of the augmented-data families sums over the categories", { + fam <- extend_family( + structure(list(family = "cumulative_rstanarm", link = "logit", + cats = paste0("c", 1:4)), class = "family"), + augdat_y_unqs = paste0("c", 1:4), augdat_link = augdat_link_cumul, + augdat_ilink = augdat_ilink_cumul, augdat_args_link = list(link = "logit"), + augdat_args_ilink = list(link = "logit") + ) + n_obs <- 6 + n_cat <- 4 + n_draws <- 3 + rand_prbs <- function() { + arr <- array(runif(n_obs * n_cat * n_draws), c(n_obs, n_cat, n_draws)) + sweep(arr, c(1, 3), apply(arr, c(1, 3), sum), "/") + } + prbs_ref <- rand_prbs() + prbs_sub <- rand_prbs() + prbs_ref[1, 2, 1] <- 0 + prbs_sub[2, 3, 2] <- 0 + w_obs <- runif(n_obs) + ce_expected <- matrix(NA_real_, n_obs, n_draws) + for (i in seq_len(n_obs)) { + for (s in seq_len(n_draws)) { + pr <- prbs_ref[i, , s] + ps <- prbs_sub[i, , s] + pr[pr == 0] <- .Machine$double.eps + ps[ps == 0] <- .Machine$double.eps + ce_expected[i, s] <- w_obs[i] * (-sum(pr * log(ps))) + } + } + expect_identical( + fam$ce_ptwise(arr2augmat(prbs_ref), arr2augmat(prbs_sub), w_obs = w_obs), + ce_expected + ) +}) + +if (exists("rng_old")) assign(".Random.seed", rng_old, envir = .GlobalEnv) +# Clean up the workspace: +rm(list = setdiff(ls(), ls_bu)) diff --git a/tests/testthat/test_cumul_native.R b/tests/testthat/test_cumul_native.R new file mode 100644 index 000000000..a588f77d9 --- /dev/null +++ b/tests/testthat/test_cumul_native.R @@ -0,0 +1,128 @@ +context("cumul_native") + +# Needed to clean up the workspace afterwards (i.e, after this test file): +ls_bu <- ls() + +if (exists(".Random.seed", envir = .GlobalEnv)) { + rng_old <- get(".Random.seed", envir = .GlobalEnv) +} +set.seed(8642) +n_cn <- 200 +ncat_cn <- 4 +x_cn <- matrix(rnorm(n_cn * 3), n_cn, 3, dimnames = list(NULL, c("a", "b", "c"))) +eta_cn <- drop(x_cn %*% c(1, -0.5, 0.3)) +links_cn <- c("logistic", "probit", "loglog", "cloglog", "cauchit") +pfun_cn <- list(logistic = plogis, probit = pnorm, + loglog = function(q) exp(-exp(-q)), + cloglog = function(q) 1 - exp(-exp(q)), cauchit = pcauchy) +# Category weights (n x C) as from the projected draws of a reference model: +wts_cn <- function(method, eta = eta_cn) { + cp <- cbind(0, pfun_cn[[method]](outer(-eta, c(-1, 0, 1), "+")), 1) + pmax(t(apply(cp, 1, diff)), 1e-10) * sample(1:3, length(eta), replace = TRUE) +} +# Augmented data (C copies of the observations) with the weights in column +# `wt`: +augdat_cn <- function(w) { + data.frame(x_cn[rep(seq_len(n_cn), ncat_cn), ], + y = factor(rep(seq_len(ncat_cn), each = n_cn)), + wt = as.vector(w)) +} + +test_that("cumul_fit_multi_c() finds the optimum of the polr() objective", { + for (method in setdiff(links_cn, "cloglog")) { + w <- wts_cn(method) + fit_polr <- suppressWarnings(MASS::polr( + y ~ a + b + c, data = augdat_cn(w), weights = wt, method = method, + start = c(1, -0.5, 0.3, -1, 0, 1), + control = list(reltol = 1e-15, maxit = 5000) + )) + fit_nat <- cumul_fit_multi_c(x_cn, w, ncat_cn, rep(0, n_cn), + match(method, links_cn), 1e-12, 100L) + expect_identical(fit_nat$convergence, 0L, info = method) + expect_equal(drop(fit_nat$beta), unname(coef(fit_polr)), tolerance = 1e-5, + info = method) + expect_equal(drop(fit_nat$zeta), unname(fit_polr$zeta), tolerance = 1e-5, + info = method) + expect_lt(2 * fit_nat$loss, fit_polr$deviance + 1e-6) + } + # cloglog (polr() often fails here because its probabilities of the highest + # category round to zero): the gradient is zero at the optimum. + w <- wts_cn("cloglog") + fit_nat <- cumul_fit_multi_c(x_cn, w, ncat_cn, rep(0, n_cn), 4L, 1e-12, 100L) + expect_identical(fit_nat$convergence, 0L) + # Survival function of the cloglog latent variable (1 - pfun, without + # cancellation): + sfun <- function(q) exp(-exp(q)) + obj <- function(theta) { + eta <- drop(x_cn %*% theta[1:3]) + g <- c(-Inf, theta[4:6], Inf) + -sum(vapply(seq_len(ncat_cn), function(k) { + z1 <- pmin(100, g[k + 1] - eta) + z2 <- pmax(-100, g[k] - eta) + pr <- ifelse(z2 > 0, sfun(z2) - sfun(z1), + pfun_cn$cloglog(z1) - pfun_cn$cloglog(z2)) + sum(w[, k] * log(pr)) + }, numeric(1))) + } + theta_hat <- c(fit_nat$beta, fit_nat$zeta) + grad_num <- vapply(seq_along(theta_hat), function(j) { + h <- 1e-6 + e <- replace(numeric(length(theta_hat)), j, h) + (obj(theta_hat + e) - obj(theta_hat - e)) / (2 * h) + }, numeric(1)) + expect_lt(max(abs(grad_num)), 1e-3) + # Several weight matrices in one call give the same fits as one at a time: + w2 <- wts_cn("logistic") + fit_both <- cumul_fit_multi_c(x_cn, cbind(w, w2), ncat_cn, rep(0, n_cn), 1L, + 1e-12, 100L) + fit_2 <- cumul_fit_multi_c(x_cn, w2, ncat_cn, rep(0, n_cn), 1L, 1e-12, 100L) + expect_identical(fit_both$beta[, 2], drop(fit_2$beta)) + expect_identical(fit_both$zeta[, 2], drop(fit_2$zeta)) +}) + +test_that("cumul_native_multi() returns usable `polr` objects", { + method <- "logistic" + ws <- cbind(as.vector(wts_cn(method)), as.vector(wts_cn(method))) + daug <- augdat_cn(ws[, 1])[c("a", "b", "c", "y")] + fam <- structure(list(family = "cumulative", link = "logit"), + class = "family") + fits <- cumul_native_multi(y ~ a + b, data = daug, family = fam, + ws_aug = ws, dot_args = list()) + expect_length(fits, 2) + newdat <- data.frame(a = c(-1, 0, 2), b = c(0.5, 0, -1)) + for (s in 1:2) { + daug$wt <- ws[, s] + fit_polr <- suppressWarnings(MASS::polr( + y ~ a + b, data = daug, weights = wt, method = method, + control = list(reltol = 1e-15, maxit = 5000) + )) + expect_s3_class(fits[[s]], "polr") + expect_identical(fits[[s]]$convergence, 0L) + expect_identical(names(coef(fits[[s]])), names(coef(fit_polr))) + expect_identical(names(fits[[s]]$zeta), names(fit_polr$zeta)) + expect_equal(coef(fits[[s]]), coef(fit_polr), tolerance = 1e-5) + expect_equal(fits[[s]]$zeta, fit_polr$zeta, tolerance = 1e-5) + expect_equal(predict(fits[[s]], newdata = newdat, type = "probs"), + predict(fit_polr, newdata = newdat, type = "probs"), + tolerance = 1e-5) + expect_true(check_conv_s(fits[[s]])) + } + # Submodel without predictors: named empty coefficient vector, as from polr(): + fits_icpt <- cumul_native_multi(y ~ 1, data = daug, family = fam, + ws_aug = ws, dot_args = list()) + expect_identical(coef(fits_icpt[[1]]), setNames(numeric(), character())) + # Cases that use fit_cumul() per draw: + expect_null(cumul_native_multi(y ~ a + b, data = daug, family = fam, + ws_aug = ws, dot_args = list(Hess = TRUE))) + opts_old <- options(projpred.augdat_native = FALSE) + expect_null(cumul_native_multi(y ~ a + b, data = daug, family = fam, + ws_aug = ws, dot_args = list())) + options(opts_old) + daug_perm <- daug[c(2:nrow(daug), 1), ] + expect_null(cumul_native_multi(y ~ a + b, data = daug_perm, family = fam, + ws_aug = ws, dot_args = list())) +}) + +if (exists("rng_old")) assign(".Random.seed", rng_old, envir = .GlobalEnv) +# Clean up the workspace: +rm(list = setdiff(ls(), ls_bu)) From b76723e4df8ec621dac90803c06075a09985449c Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Sun, 25 Oct 2026 22:16:16 +0200 Subject: [PATCH 2/3] Fit the multinomial-logit submodels of the augmented-data projection for all draws in one C++ call --- NEWS.md | 2 + R/RcppExports.R | 4 + R/divergence_minimizers.R | 108 +++++++++++++++++++ src/RcppExports.cpp | 17 +++ src/categfit.cpp | 165 +++++++++++++++++++++++++++++ tests/testthat/test_categ_native.R | 94 ++++++++++++++++ 6 files changed, 390 insertions(+) create mode 100644 src/categfit.cpp create mode 100644 tests/testthat/test_categ_native.R diff --git a/NEWS.md b/NEWS.md index b5ca7bd4f..ad3631c41 100644 --- a/NEWS.md +++ b/NEWS.md @@ -9,6 +9,8 @@ If you read this from a place other than do not edit by hand # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 +categ_fit_multi_c <- function(x, w_all, ncat, icpt, thresh, maxit) { + .Call('_projpred_categ_fit_multi_c', PACKAGE = 'projpred', x, w_all, ncat, icpt, thresh, maxit) +} + cumul_fit_multi_c <- function(x, w_all, ncat, offset, link, thresh, maxit) { .Call('_projpred_cumul_fit_multi_c', PACKAGE = 'projpred', x, w_all, ncat, offset, link, thresh, maxit) } diff --git a/R/divergence_minimizers.R b/R/divergence_minimizers.R index 9bb59015d..b9421c307 100644 --- a/R/divergence_minimizers.R +++ b/R/divergence_minimizers.R @@ -706,6 +706,18 @@ divmin_augdat <- function( } } + # Multinomial-logit submodels without group-level terms: all draws in one C++ + # call (see categ_native_multi()): + if (identical(sdivmin, fit_categ)) { + outdmin <- categ_native_multi(formula, data = data, family = family, + ws_aug = projpred_ws_aug, + dot_args = list(...)) + if (!is.null(outdmin)) { + check_conv(outdmin, rep(0L, length(outdmin)), do_check = do_check_conv) + return(outdmin) + } + } + if (is.null(getOption("projpred.parallel_proj_trigger")) && !is.null(getOption("projpred.prll_prj_trigger"))) { warning( @@ -1075,6 +1087,102 @@ fit_categ <- function(formula, data, family, weights, ...) { return(fitobj) } +# Fits the multinomial-logit submodels of all draws in one C++ call +# (categ_fit_multi_c()) with the objective of nnet::multinom(), instead of one +# fit_categ() call per draw. Like cumul_native_multi(), it builds the model +# matrix once and uses the augmented weights as an n x C matrix. Newton +# iterations replace the BFGS iterations of nnet::multinom() (by default at +# most 100), so the estimates agree with those of nnet::multinom() up to its +# convergence tolerance. The returned objects are copies of one `multinom` +# object (from nnet::multinom() with `maxit = 0`) with the weights, the +# objective and the fitted values of each draw. Returns `NULL` for other cases; +# divmin_augdat() then calls fit_categ() per draw. Global option +# `projpred.augdat_native = FALSE` switches this off. +categ_native_multi <- function(formula, data, family, ws_aug, dot_args) { + if (!getOption("projpred.augdat_native", TRUE) || family$link != "logit") { + return(NULL) + } + # fit_categ() passes these arguments on to nnet::multinom(); they are not + # reproduced here: + if (length(intersect(names(dot_args), + c(methods::formalArgs(nnet::multinom), + methods::formalArgs(nnet::nnet))))) { + return(NULL) + } + mfr <- model.frame(formula, data = data) + if (nrow(mfr) != nrow(data) || !is.null(model.offset(mfr))) { + return(NULL) + } + trms <- attr(mfr, "terms") + x <- model.matrix(trms, mfr) + y <- model.response(mfr) + if (!is.factor(y)) { + return(NULL) + } + lev <- levels(y) + ncat <- length(lev) + n <- nrow(x) %/% ncat + ws_aug <- as.matrix(ws_aug) + # The augmented data have to be C copies of the n observations, in the order + # of the response categories: + if (ncat < 3L || n * ncat != nrow(x) || nrow(ws_aug) != nrow(x) || + !identical(as.integer(y), rep(seq_len(ncat), each = n)) || + anyNA(ws_aug) || any(ws_aug < 0)) { + return(NULL) + } + x1 <- x[seq_len(n), , drop = FALSE] + if (!all(x == x1[rep(seq_len(n), ncat), , drop = FALSE])) { + return(NULL) + } + # The template object (no iterations, so its weights are the zero start + # values): + tmpl <- do.call(nnet::multinom, list(formula = formula, data = data, + weights = ws_aug[, 1], maxit = 0, + trace = FALSE)) + if (!identical(tmpl$vcoefnames, colnames(x1)) || tmpl$n[3L] != ncat || + length(tmpl$wts) != ncat * (ncol(x1) + 1L)) { + return(NULL) + } + icpt <- match("(Intercept)", colnames(x1), nomatch = 0L) - 1L + nprj <- ncol(ws_aug) + out <- categ_fit_multi_c(unname(x1), matrix(ws_aug, nrow = n), ncat, icpt, + thresh = 1e-12, maxit = 100L) + y_ind <- diag(ncat)[as.integer(y), , drop = FALSE] + dimnames(y_ind) <- dimnames(tmpl$fitted.values) + idx_aug <- rep(seq_len(n), ncat) + idx_rows <- seq_len(n) + # A loop instead of lapply() with an anonymous function, which the byte + # compiler would compile again for each projection: + outdmin <- vector("list", nprj) + for (s in seq_len(nprj)) { + coef_s <- out$coef[, , s] + dim(coef_s) <- dim(out$coef)[1:2] + # Weights as nnet::multinom() stores them (see nnet:::coef.multinom()): + # one row per category, the first column and the first row (reference + # category) zero: + wts_mat <- matrix(0, ncat, ncol(x1) + 1L) + wts_mat[-1L, -1L] <- t(coef_s) + eta <- cbind(0, x1 %*% coef_s) + eta_max <- eta[cbind(idx_rows, max.col(eta, ties.method = "first"))] + prbs <- exp(eta - eta_max) + prbs <- prbs / rowSums(prbs) + fitted_s <- prbs[idx_aug, , drop = FALSE] + dimnames(fitted_s) <- dimnames(tmpl$fitted.values) + fit_s <- tmpl + fit_s$wts <- as.vector(t(wts_mat)) + fit_s$value <- out$loss[s] + fit_s$deviance <- 2 * out$loss[s] + fit_s$AIC <- fit_s$deviance + 2 * fit_s$edf + fit_s$convergence <- out$convergence[s] + fit_s$weights <- ws_aug[, s] + fit_s$fitted.values <- fitted_s + fit_s$residuals <- y_ind - fitted_s + fit_s$call <- quote(projpred_categ_fit()) + outdmin[[s]] <- fit_s + } + return(outdmin) +} + # Use mclogit::mblogit() to fit multilevel submodels for the brms::categorical() # family: fit_categ_mlvl <- function(formula, projpred_formula_no_random, diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 9bf3972dc..a745a35bc 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -11,6 +11,22 @@ Rcpp::Rostream& Rcpp::Rcout = Rcpp::Rcpp_cout_get(); Rcpp::Rostream& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get(); #endif +// categ_fit_multi_c +List categ_fit_multi_c(arma::mat x, arma::mat w_all, int ncat, int icpt, double thresh, int maxit); +RcppExport SEXP _projpred_categ_fit_multi_c(SEXP xSEXP, SEXP w_allSEXP, SEXP ncatSEXP, SEXP icptSEXP, SEXP threshSEXP, SEXP maxitSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< arma::mat >::type x(xSEXP); + Rcpp::traits::input_parameter< arma::mat >::type w_all(w_allSEXP); + Rcpp::traits::input_parameter< int >::type ncat(ncatSEXP); + Rcpp::traits::input_parameter< int >::type icpt(icptSEXP); + Rcpp::traits::input_parameter< double >::type thresh(threshSEXP); + Rcpp::traits::input_parameter< int >::type maxit(maxitSEXP); + rcpp_result_gen = Rcpp::wrap(categ_fit_multi_c(x, w_all, ncat, icpt, thresh, maxit)); + return rcpp_result_gen; +END_RCPP +} // cumul_fit_multi_c List cumul_fit_multi_c(arma::mat x, arma::mat w_all, int ncat, arma::vec offset, int link, double thresh, int maxit); RcppExport SEXP _projpred_cumul_fit_multi_c(SEXP xSEXP, SEXP w_allSEXP, SEXP ncatSEXP, SEXP offsetSEXP, SEXP linkSEXP, SEXP threshSEXP, SEXP maxitSEXP) { @@ -75,6 +91,7 @@ END_RCPP } static const R_CallMethodDef CallEntries[] = { + {"_projpred_categ_fit_multi_c", (DL_FUNC) &_projpred_categ_fit_multi_c, 6}, {"_projpred_cumul_fit_multi_c", (DL_FUNC) &_projpred_cumul_fit_multi_c, 7}, {"_projpred_glm_elnet_c", (DL_FUNC) &_projpred_glm_elnet_c, 14}, {"_projpred_glm_ridge_c", (DL_FUNC) &_projpred_glm_ridge_c, 11}, diff --git a/src/categfit.cpp b/src/categfit.cpp new file mode 100644 index 000000000..f2afbc36b --- /dev/null +++ b/src/categfit.cpp @@ -0,0 +1,165 @@ +#include +#include +#include + +//[[Rcpp::depends(RcppArmadillo)]] + +using namespace Rcpp; +using namespace arma; + +// Multinomial-logit fits for the augmented-data projection, with the objective +// of nnet::multinom() (without weight decay): for observation i with category +// weights w_ic, c = 0, ..., C - 1, and linear predictors eta_i0 = 0 and +// eta_ik = x_i' b_k (k = 1, ..., C - 1), +// L(b) = -sum_i sum_c w_ic log(pi_ic), pi_ic = exp(eta_ic) / sum_c exp(eta_ic). +// The objective is convex. nnet::multinom() minimizes it with the BFGS method +// of nnet::nnet(); this code uses Newton steps with the analytic gradient and +// Hessian and a backtracking line search. + +// Objective (and optionally gradient and Hessian) at the coefficients `b` +// (p x K, K = C - 1; stored column-wise in the parameter vector). +static double ml_eval(const mat& x, const mat& w, const vec& nw, const mat& b, + bool deriv, vec& grad, mat& hess) { + const uword n = x.n_rows; + const uword p = x.n_cols; + const uword K = b.n_cols; + mat eta = x * b; // n x K + mat pr(n, K); + double loss = 0.0; + for (uword i = 0; i < n; ++i) { + double m = 0.0; + for (uword k = 0; k < K; ++k) m = std::max(m, eta(i, k)); + double den = std::exp(-m); + for (uword k = 0; k < K; ++k) den += std::exp(eta(i, k) - m); + double lden = m + std::log(den); + // log(pi_i0) = -lden, log(pi_ik) = eta_ik - lden: + if (w(i, 0) != 0.0) loss -= w(i, 0) * (-lden); + for (uword k = 0; k < K; ++k) { + double lp = eta(i, k) - lden; + pr(i, k) = std::exp(lp); + if (w(i, k + 1) != 0.0) loss -= w(i, k + 1) * lp; + } + } + if (deriv) { + grad.set_size(p * K); + hess.set_size(p * K, p * K); + for (uword k = 0; k < K; ++k) { + vec r = nw % pr.col(k) - w.col(k + 1); + grad.subvec(k * p, (k + 1) * p - 1) = x.t() * r; + for (uword l = k; l < K; ++l) { + vec d = -nw % pr.col(k) % pr.col(l); + if (l == k) d += nw % pr.col(k); + mat h = x.t() * (x.each_col() % d); + hess.submat(k * p, l * p, (k + 1) * p - 1, (l + 1) * p - 1) = h; + if (l != k) { + hess.submat(l * p, k * p, (l + 1) * p - 1, (k + 1) * p - 1) = h.t(); + } + } + } + } + return loss; +} + +// Fits one weight matrix `w` (n x C). The model matrix `x` contains the +// intercept column (if any) at position `icpt` (-1 if none). Writes the +// estimates into `b` (p x K); returns the objective. `conv` is 0 on +// convergence. +static double ml_fit_one(const mat& x, const mat& w, int icpt, double thresh, + int maxit, mat& b, int& conv, int& niter) { + const uword p = x.n_cols; + const uword K = w.n_cols - 1; + vec nw = sum(w, 1); + // Start: the intercepts of the model without predictors, otherwise zero: + b.zeros(p, K); + if (icpt >= 0) { + rowvec wc = sum(w, 0); + for (uword k = 0; k < K; ++k) { + if (wc(0) > 0.0 && wc(k + 1) > 0.0) { + b(icpt, k) = std::log(wc(k + 1) / wc(0)); + } + } + } + vec grad; + mat hess; + double loss = ml_eval(x, w, nw, b, true, grad, hess); + conv = 1; + niter = 0; + const uword np = p * K; + for (int it = 0; it < maxit; ++it) { + niter = it + 1; + vec dir; + mat hreg = hess; + double damp = 0.0; + bool ok = false; + double dscale = std::max(1e-10, max(abs(hess.diag()))); + for (int k = 0; k < 30 && !ok; ++k) { + mat r; + if (damp > 0.0) hreg = hess + damp * eye(np, np); + ok = chol(r, hreg); + if (ok) { + dir = -solve(trimatu(r), solve(trimatl(r.t()), grad)); + } else { + damp = (damp == 0.0) ? 1e-8 * dscale : 10.0 * damp; + } + } + if (!ok) break; + double decrement = -dot(grad, dir); + if (decrement < thresh * std::max(1.0, std::abs(loss))) { + conv = 0; + break; + } + mat bdir = reshape(dir, p, K); + double t = 1.0; + bool accepted = false; + mat b_new; + vec grad_new; + mat hess_new; + for (int ls = 0; ls < 60; ++ls) { + b_new = b + t * bdir; + double loss_new = ml_eval(x, w, nw, b_new, false, grad_new, hess_new); + if (std::isfinite(loss_new) && loss_new <= loss - 1e-4 * t * decrement) { + accepted = true; + break; + } + t *= 0.5; + } + if (!accepted) { + if (decrement < 1e3 * thresh * std::max(1.0, std::abs(loss))) conv = 0; + break; + } + b = b_new; + loss = ml_eval(x, w, nw, b, true, grad, hess); + } + return loss; +} + +// Fits the multinomial-logit model for S weight matrices, given as the n x C +// column blocks of `w_all` (n x (C * S)). The model matrix `x` (n x p) is +// shared; `icpt` is the 0-based index of its intercept column (-1 if none). +// Returns the coefficients as a p x K x S array. +// [[Rcpp::export]] +List categ_fit_multi_c(arma::mat x, arma::mat w_all, int ncat, int icpt, + double thresh, int maxit) { + const uword n = x.n_rows; + const uword p = x.n_cols; + const uword K = ncat - 1; + const uword S = w_all.n_cols / ncat; + if (w_all.n_rows != n || S * ncat != w_all.n_cols) { + Rcpp::stop("Dimensions of `w_all` do not match."); + } + cube b_out(p, K, S); + vec loss_out(S); + IntegerVector conv_out(S), niter_out(S); + for (uword s = 0; s < S; ++s) { + mat w = w_all.cols(s * ncat, (s + 1) * ncat - 1); + mat b; + int conv, niter; + loss_out(s) = ml_fit_one(x, w, icpt, thresh, maxit, b, conv, niter); + b_out.slice(s) = b; + conv_out[s] = conv; + niter_out[s] = niter; + } + return List::create(Named("coef") = b_out, Named("loss") = loss_out, + Named("convergence") = conv_out, + Named("niter") = niter_out); +} diff --git a/tests/testthat/test_categ_native.R b/tests/testthat/test_categ_native.R new file mode 100644 index 000000000..570614ef4 --- /dev/null +++ b/tests/testthat/test_categ_native.R @@ -0,0 +1,94 @@ +context("categ_native") + +# Needed to clean up the workspace afterwards (i.e, after this test file): +ls_bu <- ls() + +if (exists(".Random.seed", envir = .GlobalEnv)) { + rng_old <- get(".Random.seed", envir = .GlobalEnv) +} +set.seed(24680) +n_ca <- 150 +ncat_ca <- 4 +dat_ca <- data.frame(a = rnorm(n_ca), b = rnorm(n_ca), + f = factor(sample(c("u", "v"), n_ca, replace = TRUE))) +eta_ca <- with(dat_ca, cbind(0, 0.8 * a, -0.6 * b + 0.5 * (f == "v"), + 0.4 * a + 0.3 * b)) +# Category weights (n x C) as from the projected draws of a reference model: +wts_ca <- function() { + prb <- exp(eta_ca) / rowSums(exp(eta_ca)) + prb * sample(1:3, n_ca, replace = TRUE) +} +# Augmented data (C copies of the observations) with the weights in column +# `wt`: +augdat_ca <- function(w) { + cbind(dat_ca[rep(seq_len(n_ca), ncat_ca), ], + y = factor(rep(seq_len(ncat_ca), each = n_ca)), + wt = as.vector(w)) +} +fit_tight <- function(fml, daug) { + nnet::multinom(fml, data = daug, weights = wt, trace = FALSE, + maxit = 10000, reltol = 1e-15, abstol = 0) +} + +test_that("categ_fit_multi_c() finds the optimum of the multinom() objective", { + w <- wts_ca() + daug <- augdat_ca(w) + fit_mn <- fit_tight(y ~ a + b + f, daug) + x <- model.matrix(~ a + b + f, data = dat_ca) + fit_nat <- categ_fit_multi_c(x, w, ncat_ca, 0L, 1e-12, 100L) + expect_identical(fit_nat$convergence, 0L) + expect_equal(t(fit_nat$coef[, , 1]), unname(coef(fit_mn)), tolerance = 1e-5) + expect_lt(fit_nat$loss, fit_mn$value + 1e-8) + # Several weight matrices in one call give the same fits as one at a time: + w2 <- wts_ca() + fit_both <- categ_fit_multi_c(x, cbind(w, w2), ncat_ca, 0L, 1e-12, 100L) + fit_2 <- categ_fit_multi_c(x, w2, ncat_ca, 0L, 1e-12, 100L) + expect_identical(fit_both$coef[, , 2], fit_2$coef[, , 1]) +}) + +test_that("categ_native_multi() returns usable `multinom` objects", { + ws <- cbind(as.vector(wts_ca()), as.vector(wts_ca())) + daug <- augdat_ca(ws[, 1])[c("a", "b", "f", "y")] + fam <- structure(list(family = "categorical", link = "logit"), + class = "family") + fits <- categ_native_multi(y ~ a + f, data = daug, family = fam, + ws_aug = ws, dot_args = list()) + expect_length(fits, 2) + newdat <- data.frame(a = c(-1, 0.5, 2), + f = factor(c("v", "u", "v"), levels = c("u", "v"))) + for (s in 1:2) { + daug$wt <- ws[, s] + fit_mn <- fit_tight(y ~ a + f, daug) + expect_s3_class(fits[[s]], c("multinom", "nnet")) + expect_identical(fits[[s]]$convergence, 0L) + expect_identical(dimnames(coef(fits[[s]])), dimnames(coef(fit_mn))) + expect_identical(fits[[s]]$lev, fit_mn$lev) + expect_equal(coef(fits[[s]]), coef(fit_mn), tolerance = 1e-5) + expect_equal(fits[[s]]$value, fit_mn$value, tolerance = 1e-8) + expect_equal(predict(fits[[s]], newdata = newdat, type = "probs"), + predict(fit_mn, newdata = newdat, type = "probs"), + tolerance = 1e-5) + expect_equal(unname(fits[[s]]$fitted.values), unname(fit_mn$fitted.values), + tolerance = 1e-5) + expect_true(check_conv_s(fits[[s]])) + } + # Submodel without predictors: + fits_icpt <- categ_native_multi(y ~ 1, data = daug, family = fam, + ws_aug = ws, dot_args = list()) + daug$wt <- ws[, 1] + expect_equal(coef(fits_icpt[[1]]), coef(fit_tight(y ~ 1, daug)), + tolerance = 1e-5) + # Cases that use fit_categ() per draw: + expect_null(categ_native_multi(y ~ a, data = daug, family = fam, + ws_aug = ws, dot_args = list(Hess = TRUE))) + opts_old <- options(projpred.augdat_native = FALSE) + expect_null(categ_native_multi(y ~ a, data = daug, family = fam, + ws_aug = ws, dot_args = list())) + options(opts_old) + expect_null(categ_native_multi(y ~ a, data = daug[c(2:nrow(daug), 1), ], + family = fam, ws_aug = ws, dot_args = list())) +}) + +if (exists("rng_old")) assign(".Random.seed", rng_old, envir = .GlobalEnv) +# Clean up the workspace: +rm(list = setdiff(ls(), ls_bu)) From ea2cfb29b7502e6230d53a7ed16e33ab3a3d13a8 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Mon, 26 Oct 2026 21:13:19 +0200 Subject: [PATCH 3/3] Add the PR number to the NEWS entries --- NEWS.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/NEWS.md b/NEWS.md index ad3631c41..c46761f2d 100644 --- a/NEWS.md +++ b/NEWS.md @@ -7,9 +7,9 @@ If you read this from a place other than