From 77e09732f8adea07333bdcdd25db8d00f2895e70 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Wed, 21 Oct 2026 11:45:33 +0300 Subject: [PATCH 1/2] Compute stacking weights with an active-set Newton method --- NAMESPACE | 1 - R/helpers.R | 21 -- R/loo_model_weights.R | 233 +++++++++++++++++------ man/loo_model_weights.Rd | 44 +++-- tests/testthat/_snaps/model_weighting.md | 2 +- tests/testthat/test_model_weighting.R | 105 +++++++--- 6 files changed, 282 insertions(+), 124 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index c9c5383d..00ced50a 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -171,7 +171,6 @@ importFrom(parallel,makePSOCKcluster) importFrom(parallel,mclapply) importFrom(parallel,parLapply) importFrom(parallel,stopCluster) -importFrom(stats,constrOptim) importFrom(stats,nobs) importFrom(stats,qnorm) importFrom(stats,quantile) diff --git a/R/helpers.R b/R/helpers.R index 37d407a8..9f4e95ab 100644 --- a/R/helpers.R +++ b/R/helpers.R @@ -51,27 +51,6 @@ difference_of_squares <- function(x, y) { (x - y) * (x + y) } -#' More stable version of `(exp(a) - exp(b)) / exp(c)` -#' -#' @noRd -#' @param a,b,c Numeric vectors of the same length. -#' @return A numeric vector equal to `(exp(a) - exp(b)) / exp(c)`. Elements -#' with `a == b` are returned as an exact zero regardless of `c`; elsewhere -#' `NA` and `NaN` inputs propagate. -#' -exp_diff_over_exp <- function(a, b, c) { - # `a >= b` is NA if `a` or `b` is NA or NaN, and R silently ignores NA - # indices in `[<-`. Seed the result from the inputs and index with which() - # so that missing values propagate instead of leaving a zero behind. - out <- a + b + c - larger <- which(a >= b) - smaller <- which(a < b) - out[larger] <- exp(a[larger] - c[larger]) * -expm1(b[larger] - a[larger]) - out[smaller] <- exp(b[smaller] - c[smaller]) * expm1(a[smaller] - b[smaller]) - out[which(a == b)] <- 0 - out -} - #' Compute point estimates and standard errors from pointwise vectors #' #' @noRd diff --git a/R/loo_model_weights.R b/R/loo_model_weights.R index 9441fd09..b9781345 100644 --- a/R/loo_model_weights.R +++ b/R/loo_model_weights.R @@ -30,12 +30,9 @@ #' @param alpha Positive scalar shape parameter in the Dirichlet distribution #' used for the Bayesian bootstrap. The default is `alpha=1`, which #' corresponds to a uniform distribution on the simplex space. -#' @param optim_method If `method="stacking"`, a string passed to the `method` -#' argument of [stats::constrOptim()] to specify the optimization algorithm. -#' The default is `optim_method="BFGS"`, but other options are available (see -#' [stats::optim()]). -#' @param optim_control If `method="stacking"`, a list of control parameters for -#' optimization passed to the `control` argument of [stats::constrOptim()]. +#' @param optim_method,optim_control Deprecated and ignored. Earlier versions +#' passed them to [stats::constrOptim()]. The stacking weights are now +#' computed with an active-set Newton method, see Details. #' @param r_eff_list Optionally, a list of relative effective sample size #' estimates for the likelihood `(exp(log_lik))` of each observation in #' each model. See [psis()] and [relative_eff()] helper @@ -69,11 +66,28 @@ #' In general, we recommend stacking for averaging predictive distributions, #' while pseudo-BMA+ can serve as a computationally easier alternative. #' +#' `stacking_weights()` maximizes the log score with an active-set method. It +#' starts with the model with the largest total log predictive density, finds +#' the optimal weights of the models in the active set with Newton steps, and +#' adds the model that improves the log score most, until no model can improve +#' it. A model whose weight reaches zero leaves the active set. At the returned +#' weights, the log score is at most \eqn{N \times 10^{-8}} below its maximum. +#' The optimal weights are usually sparse, so the method is fast also for many +#' models. The method is a variant of the constrained Newton method for mixture +#' proportions (Wang, 2007), and the stopping rule uses the optimality condition +#' of Lindsay (1983). +#' +#' Wang, Y. (2007). On fast computation of the non-parametric maximum +#' likelihood estimate of a mixing distribution. *Journal of the Royal +#' Statistical Society: Series B*, 69(2), 185--198. +#' +#' Lindsay, B. G. (1983). The geometry of mixture likelihoods: a general +#' theory. *The Annals of Statistics*, 11(1), 86--94. +#' #' @seealso #' * The __loo__ package [vignettes](https://mc-stan.org/loo/articles/), particularly #' [Bayesian Stacking and Pseudo-BMA weights using the __loo__ package](https://mc-stan.org/loo/articles/loo2-weights.html). #' * [loo()] for details on leave-one-out ELPD estimation. -#' * [constrOptim()] for the choice of optimization methods and control-parameters. #' * [relative_eff()] for computing `r_eff`. #' #' @template loo-and-psis-references @@ -124,8 +138,7 @@ #' wts1 <- loo_model_weights( #' log_lik_list, #' method = "stacking", -#' r_eff_list = r_eff_list, -#' optim_control = list(reltol=1e-10) +#' r_eff_list = r_eff_list #' ) #' print(wts1) #' @@ -134,11 +147,7 @@ #' loo(log_lik_list[[j]], r_eff = r_eff_list[[j]]) #' }) #' -#' wts2 <- loo_model_weights( -#' loo_list, -#' method = "stacking", -#' optim_control = list(reltol=1e-10) -#' ) +#' wts2 <- loo_model_weights(loo_list, method = "stacking") #' all.equal(wts1, wts2) #' #' # can provide names to be used in the results @@ -172,8 +181,8 @@ loo_model_weights.default <- function(x, ..., method = c("stacking", "pseudobma"), - optim_method = "BFGS", - optim_control = list(), + optim_method = NULL, + optim_control = NULL, BB = TRUE, BB_n = 1000, alpha = 1, @@ -244,12 +253,10 @@ loo_model_weights.default <- #' approximately using [loo()] or by running exact leave-one-out or K-fold #' cross-validation. #' -#' @importFrom stats constrOptim -#' stacking_weights <- function(lpd_point, - optim_method = "BFGS", - optim_control = list()) { + optim_method = NULL, + optim_control = NULL) { stopifnot(is.matrix(lpd_point)) N <- nrow(lpd_point) @@ -262,43 +269,17 @@ stacking_weights <- stop("Each observation must have a finite predictive density for at least one model.") } - negative_log_score_loo <- function(w) { - # objective function: log score - stopifnot(length(w) == K - 1) - w_full <- c(w, 1 - sum(w)) - # avoid over- and underflows using log weights and rowLogSumExps - sum <- sum(matrixStats::rowLogSumExps(sweep(lpd_point[1:N,], 2, log(w_full), '+'))) - return(-as.numeric(sum)) - } - - gradient <- function(w) { - # gradient of the objective function - stopifnot(length(w) == K - 1) - w_full <- c(w, 1 - sum(w)) - mixture_lpd <- matrixStats::rowLogSumExps( - sweep(lpd_point, 2, log(w_full), "+") + if (!is.null(optim_method) || !is.null(optim_control)) { + warning( + "The arguments 'optim_method' and 'optim_control' are deprecated and ", + "ignored. The stacking weights are now computed with an active-set ", + "Newton method.", + call. = FALSE ) - -colSums(exp_diff_over_exp( - lpd_point[, seq_len(K - 1), drop = FALSE], - matrix(lpd_point[, K], N, K - 1), - matrix(mixture_lpd, N, K - 1) - )) } - ui <- rbind(rep(-1, K - 1), diag(K - 1)) # K-1 simplex constraint matrix - ci <- c(-1, rep(0, K - 1)) - w <- constrOptim( - theta = rep(1 / K, K - 1), - f = negative_log_score_loo, - grad = gradient, - ui = ui, - ci = ci, - method = optim_method, - control = optim_control - )$par - structure( - c(w, 1 - sum(w)), + stacking_weights_active_set(lpd_point), names = paste0("model", 1:K), class = c("stacking_weights") ) @@ -346,6 +327,152 @@ pseudobma_weights <- } +#' Stacking weights with an active-set Newton method +#' +#' Maximizes `sum_i log(sum_k w_k exp(lpd_point[i, k]))` over the simplex. The +#' problem is concave in `w`, and the optimal weights are usually sparse. The +#' method optimizes the weights of the active models with Newton steps and adds +#' the model with the largest `g_k = mean_i exp(lpd_point[i, k]) / p_i`, where +#' `p_i` is the predictive density of the combination, until `max_k g_k` is at +#' most `exp(tol)`. Then the log score is at most `N * tol` below its maximum, +#' by Jensen's inequality. +#' +#' @noRd +#' @param lpd_point N by K matrix of pointwise log predictive densities. Each +#' row must have at least one finite value. +#' @param tol Tolerance for `max_k log(g_k)`. +#' @return Vector of K weights that sums to 1. +#' +stacking_weights_active_set <- function(lpd_point, tol = 1e-8) { + N <- nrow(lpd_point) + K <- ncol(lpd_point) + + # start with the best model, and add the best model of each observation + # that no model in the set can predict, so that every row has a finite value + active <- which.max(colSums(pmax(lpd_point, -1e300))) + uncovered <- which(rowSums(is.finite(lpd_point[, active, drop = FALSE])) == 0) + while (length(uncovered) > 0) { + active <- c(active, which.max(lpd_point[uncovered[1], ])) + uncovered <- which( + rowSums(is.finite(lpd_point[, active, drop = FALSE])) == 0 + ) + } + w <- rep(1 / length(active), length(active)) + + neg_log_score <- function(w, lpd) { + -sum(matrixStats::rowLogSumExps(sweep(lpd, 2, log(w), "+"))) + } + + for (outer in seq_len(K + 100)) { + w_res <- stacking_newton(w, lpd_point[, active, drop = FALSE], + neg_log_score, tol = tol) + keep <- w_res > 0 + active <- active[keep] + w <- w_res[keep] / sum(w_res[keep]) + + # log g_k for all models + log_p <- matrixStats::rowLogSumExps( + sweep(lpd_point[, active, drop = FALSE], 2, log(w), "+") + ) + log_g <- matrixStats::colLogSumExps(lpd_point - log_p) - log(N) + k_new <- which.max(log_g) + if (log_g[k_new] <= tol || k_new %in% active) { + # optimal, or no model can be added within the numerical precision + break + } + active <- c(active, k_new) + w <- c(w * (1 - 1e-3), 1e-3) + } + + out <- numeric(K) + out[active] <- w + out +} + +#' Newton steps for the stacking weights of a fixed set of models +#' +#' @noRd +#' @param w Starting weights, positive and summing to 1. +#' @param lpd N by |active| matrix of log predictive densities. Each row has at +#' least one finite value. +#' @param neg_log_score Function of `w` and `lpd` that returns the negative log +#' score. +#' @param tol Tolerance of `stacking_weights_active_set()`. Newton stops when +#' every `g_k` of the set is within `tol / 10` of 1. +#' @return Weights, non-negative and summing to 1. Models whose weight reaches +#' zero leave the set, and their weights in the result are zero. +#' +stacking_newton <- function(w, lpd, neg_log_score, tol, max_iter = 100) { + out <- numeric(length(w)) + in_set <- seq_along(w) + for (iter in seq_len(max_iter)) { + n_active <- length(w) + if (n_active == 1) { + break + } + lpd_set <- lpd[, in_set, drop = FALSE] + # scale each row by its maximum over the current set, so that the + # densities are at most 1 and the largest one in each row is 1 + dens <- exp(lpd_set - matrixStats::rowMaxs(lpd_set)) + p <- drop(dens %*% w) + grad <- -colSums(dens / p) + hess <- crossprod(dens / p) + # a small ridge keeps the system solvable for identical models + hess <- hess + diag(1e-12 * mean(diag(hess)), n_active) + kkt <- rbind(cbind(hess, 1), c(rep(1, n_active), 0)) + step <- tryCatch( + solve(kkt, c(-grad, 0))[seq_len(n_active)], + error = function(e) NULL + ) + if (is.null(step)) { + break + } + # at the optimum of this set, g_k = -grad_k / N is 1 for every model + if (max(abs(-grad / nrow(lpd) - 1)) < tol / 10) { + break + } + f0 <- neg_log_score(w, lpd_set) + decrement <- -sum(grad * step) + # close to the optimum the change of the log score is below its rounding + # error, and the line search cannot see it; take the Newton step there + near_optimum <- decrement < 1e-10 * max(1, abs(f0)) + # largest step that keeps the weights non-negative + ratio <- ifelse(step < 0, -w / step, Inf) + t_max <- min(1, ratio) + if (t_max < 1e-10) { + # a weight is almost zero and the step would make it negative: the + # optimum of this set has that weight at zero + w[ratio < 1e-10] <- 0 + } else { + t <- t_max + repeat { + w_new <- pmax(w + t * step, 0) + f_new <- neg_log_score(w_new, lpd_set) + # a step to the boundary sets a weight to zero; accept it if the log + # score does not get worse beyond rounding + to_boundary <- t == t_max && t_max < 1 && + f_new <= f0 + 1e-12 * max(1, abs(f0)) + if (near_optimum || f_new <= f0 - 1e-4 * t * decrement || + to_boundary || t < 1e-10) { + break + } + t <- t / 2 + } + if (t < 1e-10) { + break + } + w <- w_new + } + w <- w / sum(w) + # models whose weight reached zero leave the set + keep <- w > 0 + in_set <- in_set[keep] + w <- w[keep] / sum(w[keep]) + } + out[in_set] <- w + out +} + #' @export print.stacking_weights <- function(x, digits = 3, ...) { cat("Method: stacking\n------\n") diff --git a/man/loo_model_weights.Rd b/man/loo_model_weights.Rd index e7323912..0bca1157 100644 --- a/man/loo_model_weights.Rd +++ b/man/loo_model_weights.Rd @@ -13,8 +13,8 @@ loo_model_weights(x, ...) x, ..., method = c("stacking", "pseudobma"), - optim_method = "BFGS", - optim_control = list(), + optim_method = NULL, + optim_control = NULL, BB = TRUE, BB_n = 1000, alpha = 1, @@ -22,7 +22,7 @@ loo_model_weights(x, ...) cores = getOption("mc.cores", 1) ) -stacking_weights(lpd_point, optim_method = "BFGS", optim_control = list()) +stacking_weights(lpd_point, optim_method = NULL, optim_control = NULL) pseudobma_weights(lpd_point, BB = TRUE, BB_n = 1000, alpha = 1) } @@ -47,13 +47,9 @@ to use for obtaining the weights. \code{"stacking"} refers to stacking of predictive distributions and \code{"pseudobma"} refers to pseudo-BMA+ weighting (or plain pseudo-BMA weighting if argument \code{BB} is \code{FALSE}).} -\item{optim_method}{If \code{method="stacking"}, a string passed to the \code{method} -argument of \code{\link[stats:constrOptim]{stats::constrOptim()}} to specify the optimization algorithm. -The default is \code{optim_method="BFGS"}, but other options are available (see -\code{\link[stats:optim]{stats::optim()}}).} - -\item{optim_control}{If \code{method="stacking"}, a list of control parameters for -optimization passed to the \code{control} argument of \code{\link[stats:constrOptim]{stats::constrOptim()}}.} +\item{optim_method, optim_control}{Deprecated and ignored. Earlier versions +passed them to \code{\link[stats:constrOptim]{stats::constrOptim()}}. The stacking weights are now +computed with an active-set Newton method, see Details.} \item{BB}{Logical used when \code{"method"}=\code{"pseudobma"}. If \code{TRUE} (the default), the Bayesian bootstrap will be used to adjust @@ -125,6 +121,24 @@ regularizes the weights away from the extremes of 0 and 1. In general, we recommend stacking for averaging predictive distributions, while pseudo-BMA+ can serve as a computationally easier alternative. + +\code{stacking_weights()} maximizes the log score with an active-set method. It +starts with the model with the largest total log predictive density, finds +the optimal weights of the models in the active set with Newton steps, and +adds the model that improves the log score most, until no model can improve +it. A model whose weight reaches zero leaves the active set. At the returned +weights, the log score is at most \eqn{N \times 10^{-8}} below its maximum. +The optimal weights are usually sparse, so the method is fast also for many +models. The method is a variant of the constrained Newton method for mixture +proportions (Wang, 2007), and the stopping rule uses the optimality condition +of Lindsay (1983). + +Wang, Y. (2007). On fast computation of the non-parametric maximum +likelihood estimate of a mixing distribution. \emph{Journal of the Royal +Statistical Society: Series B}, 69(2), 185--198. + +Lindsay, B. G. (1983). The geometry of mixture likelihoods: a general +theory. \emph{The Annals of Statistics}, 11(1), 86--94. } \examples{ \dontrun{ @@ -171,8 +185,7 @@ r_eff_list <- lapply(model_list, function(x) { wts1 <- loo_model_weights( log_lik_list, method = "stacking", - r_eff_list = r_eff_list, - optim_control = list(reltol=1e-10) + r_eff_list = r_eff_list ) print(wts1) @@ -181,11 +194,7 @@ loo_list <- lapply(1:length(log_lik_list), function(j) { loo(log_lik_list[[j]], r_eff = r_eff_list[[j]]) }) -wts2 <- loo_model_weights( - loo_list, - method = "stacking", - optim_control = list(reltol=1e-10) -) +wts2 <- loo_model_weights(loo_list, method = "stacking") all.equal(wts1, wts2) # can provide names to be used in the results @@ -231,7 +240,6 @@ stacking to average Bayesian predictive distributions. \item The \strong{loo} package \href{https://mc-stan.org/loo/articles/}{vignettes}, particularly \href{https://mc-stan.org/loo/articles/loo2-weights.html}{Bayesian Stacking and Pseudo-BMA weights using the \strong{loo} package}. \item \code{\link[=loo]{loo()}} for details on leave-one-out ELPD estimation. -\item \code{\link[=constrOptim]{constrOptim()}} for the choice of optimization methods and control-parameters. \item \code{\link[=relative_eff]{relative_eff()}} for computing \code{r_eff}. } } diff --git a/tests/testthat/_snaps/model_weighting.md b/tests/testthat/_snaps/model_weighting.md index ba0fad93..dc174574 100644 --- a/tests/testthat/_snaps/model_weighting.md +++ b/tests/testthat/_snaps/model_weighting.md @@ -1,6 +1,6 @@ # loo_model_weights (stacking and pseudo-BMA) gives expected result - WAoAAAACAAQFAgACAwAAAAAOAAAAAz/KEXnf1DM9P+l7oXTYyi4+YzIpAwAAAA== + WAoAAAACAAQGAQACAwAAAAAOAAAAAz/KD9BbGOhwP+l8C+k5xeQAAAAAAAAAAA== --- diff --git a/tests/testthat/test_model_weighting.R b/tests/testthat/test_model_weighting.R index eb64af92..252dc2df 100644 --- a/tests/testthat/test_model_weighting.R +++ b/tests/testthat/test_model_weighting.R @@ -20,31 +20,6 @@ loo_list <- lapply(1:length(ll_list), function(j) { tol <- 0.01 # absolute tolerance of weights -test_that("stacking gradient is stable for similar model predictions", { - a <- -30 + 1e-14 - b <- -30 - - expect_equal(exp_diff_over_exp(a, b, b), 1.065814103640156e-14) - expect_equal(exp_diff_over_exp(b, a, b), -1.065814103640156e-14) - expect_equal( - exp_diff_over_exp(c(-Inf, 0, Inf), c(-Inf, 0, Inf), c(0, 0, 0)), - c(0, 0, 0) - ) -}) - -test_that("exp_diff_over_exp propagates missing values", { - expect_equal(exp_diff_over_exp(NaN, 0, 0), NaN) - expect_equal(exp_diff_over_exp(0, NaN, 0), NaN) - expect_equal(exp_diff_over_exp(0, -1, NaN), NaN) - expect_equal(exp_diff_over_exp(NA_real_, 0, 0), NA_real_) - # equal numerators short-circuit to an exact zero, whatever the denominator - expect_equal(exp_diff_over_exp(0, 0, NaN), 0) - expect_equal( - exp_diff_over_exp(c(NaN, 0), c(0, -1), c(0, 0)), - c(NaN, exp_diff_over_exp(0, -1, 0)) - ) -}) - test_that("stacking handles equal infinite log predictive densities", { lpd <- matrix( c(-Inf, 0, -Inf, -1, -1, -1), @@ -146,6 +121,23 @@ test_that("loo_model_weights (stacking and pseudo-BMA) gives expected result", { expect_identical(w3, w3_b) }) +# Upper bound of how far the log score of the weights `w` is below its +# maximum: N * log(max_k g_k), with g_k = mean_i p_ik / sum_j w_j p_ij +stacking_gap_bound <- function(w, lpd) { + log_p <- matrixStats::rowLogSumExps(sweep(lpd, 2, log(w), "+")) + nrow(lpd) * max(matrixStats::colLogSumExps(lpd - log_p) - log(nrow(lpd))) +} + +# EM for mixture weights with fixed components: a slow independent solver +stacking_em <- function(lpd, iter = 20000) { + w <- rep(1 / ncol(lpd), ncol(lpd)) + for (i in seq_len(iter)) { + a <- sweep(lpd, 2, log(w), "+") + w <- colMeans(exp(a - matrixStats::rowLogSumExps(a))) + } + w +} + test_that("stacking_weights gives expected result", { set.seed(10) lpd_point <- matrix(rnorm(2000, -1), 400, 5) @@ -155,12 +147,65 @@ test_that("stacking_weights gives expected result", { expect_s3_class(actual, "stacking_weights") expect_named(actual, paste0("model", 1:5)) expect_equal(sum(actual), 1) - expect_equal( - as.numeric(actual), - c(0.158382400055713, 0.176065699214342, 0.19918082992972, - 0.264490615659226, 0.201880455141), - tolerance = 1e-6 + expect_lt(stacking_gap_bound(as.numeric(actual), lpd_point), 400 * 1e-8) + expect_equal(as.numeric(actual), stacking_em(lpd_point), tolerance = 1e-5) +}) + +test_that("stacking_weights finds the optimum for many models", { + set.seed(119) + N <- 100 + K <- 500 + y <- rnorm(N) + mu <- rnorm(K, 0, 0.5) + s <- exp(rnorm(K, 0, 0.3)) + lpd <- sapply(seq_len(K), function(k) dnorm(y, mu[k], s[k], log = TRUE)) + + w <- as.numeric(stacking_weights(lpd)) + expect_equal(sum(w), 1) + expect_true(all(w >= 0)) + expect_lt(stacking_gap_bound(w, lpd), N * 1e-8) + # the optimal weights are sparse + expect_lt(sum(w > 0), 20) +}) + +test_that("stacking_weights handles -Inf, duplicate and identical models", { + set.seed(2) + N <- 60 + lpd <- matrix(rnorm(3 * N, -1), N, 3) + # each model is impossible for a third of the observations, so the + # combination needs all three + lpd[1:20, 1] <- -Inf + lpd[21:40, 2] <- -Inf + lpd[41:60, 3] <- -Inf + w <- as.numeric(stacking_weights(lpd)) + expect_true(all(w > 0)) + expect_lt(stacking_gap_bound(w, lpd), N * 1e-8) + + # duplicate models share the weight that one of them would get + lpd2 <- matrix(rnorm(4 * N, -1), N, 4) + w2 <- as.numeric(stacking_weights(lpd2)) + w_dup <- as.numeric(stacking_weights(cbind(lpd2, lpd2[, 1]))) + expect_equal(c(w_dup[1] + w_dup[5], w_dup[2:4]), w2, tolerance = 1e-6) + + # identical models: all weight to one of them + lpd3 <- matrix(rep(rnorm(N, -1), 3), N, 3) + expect_equal(sort(as.numeric(stacking_weights(lpd3))), c(0, 0, 1)) +}) + +test_that("optim_method and optim_control are deprecated", { + set.seed(3) + lpd <- matrix(rnorm(40, -1), 20, 2) + w <- stacking_weights(lpd) + expect_warning( + w2 <- stacking_weights(lpd, optim_control = list(reltol = 1e-10)), + "deprecated and ignored" ) + expect_identical(w2, w) + expect_warning( + loo_model_weights(ll_list, optim_method = "Nelder-Mead"), + "deprecated and ignored" + ) + expect_no_warning(loo_model_weights(ll_list)) }) test_that("Bayesian bootstrap gives expected result", { From 07446dbffc952bd613123d21555508b5bcbb3212 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Wed, 21 Oct 2026 12:06:23 +0300 Subject: [PATCH 2/2] chore: add news --- NEWS.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/NEWS.md b/NEWS.md index 877c1afe..4785be21 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,7 @@ # loo (development version) +* Faster and more accurate stacking weights with an active-set Newton method; `optim_method` +and `optim_control` are deprecated. by @avehtari in #412 * Improve numerical stability in `loo()`, `psis()`, model weighting, subsampling, and moment matching in #395 * Fix `loo_compare()` when used with subsampling: compute model comparison by