Predictive Stacking: Optimal Model Ensemble

• Predictive stacking is a rigorous statistical framework that convexly combines multiple probabilistic models to maximize out-of-sample predictive performance.
• It learns optimal stacking weights through cross-validation techniques such as LOO or K-fold, and supports implementations from classical linear methods to neural networks.
• The method enhances robustness and computational efficiency in high-dimensional, spatial, temporal, and transfer learning applications while quantifying model uncertainty.

Predictive stacking is a rigorous statistical framework for optimally combining multiple predictive models—often Bayesian or otherwise probabilistic—so as to maximize out-of-sample predictive performance according to a chosen proper scoring rule, typically the log-score. The central idea is to form a convex mixture of candidate predictive distributions, learning the stacking weights by minimizing a cross-validated loss. Predictive stacking has established itself as a principled alternative to Bayesian model averaging, particularly in settings where the true data-generating process is not among the candidates (the "M-open" scenario). It has achieved widespread adoption in spatial modeling, regression, machine learning, time-series, and scientific analyses where model uncertainty, computational efficiency, and robust uncertainty quantification are critical.

1. Formal Definition and Mathematical Formulation

Predictive stacking operates by considering a set of $M$ candidate models $\{\mathcal{M}_1, \ldots, \mathcal{M}_M\}$, each yielding a posterior predictive density $p_m(\tilde{y}|y)$. The stacked predictive density is defined as

$$p_{\mathrm{stack}}(\tilde{y}|y) = \sum_{m=1}^M w_m\, p_m(\tilde{y}|y)\,,$$

where $w = (w_1, \ldots, w_M)$ are non-negative weights summing to unity ($w_m \geq 0$, $\sum_m w_m = 1$).

The key aspect of stacking is the estimation of weights $w$. These are learned by maximizing a proper scoring rule—typically the average leave-one-out (LOO) log-likelihood of the stacked predictive density: $$w^* = \arg\max_{w \in \Delta_{M-1}} \frac{1}{n} \sum_{i=1}^n \log\left( \sum_{m=1}^M w_m\, p_m(y_i|y_{-i}) \right)\,,$$ where $y_{-i}$ denotes the data with the $i$-th observation removed, and $\Delta_{M-1}$ is the probability simplex (Yao et al., 2017).

Cross-validation (often LOO or $K$-fold) is essential to prevent overfitting and to provide an unbiased estimate of predictive utility for candidate model combinations. This optimization problem is strictly convex when the $p_m$ are linearly independent, ensuring a unique solution. Variants include stacking of means or medians (minimizing squared or absolute error) versus full distributional stacking (maximizing log-score/KL), as well as feature-adaptive weights via Bayesian hierarchical or neural-network structures (Yao et al., 2021, Coscrato et al., 2019).

2. Algorithmic Realizations and Extensions

Stacking can be implemented via several architectures, which differ both in the form of the base learners and in the meta-level model:

• Classical weighted linear stacking learns $w$ by (penalized) least-squares on cross-validated out-of-fold predictions, optionally with constraints such as non-negativity or sum-to-one (Chen et al., 2023).
• Constrained and unconstrained stacking: Constraints on $w$ (e.g., $w_j \geq 0$, $\sum_j w_j = 1$) can be enforced or relaxed. Empirical and theoretical results indicate that unconstrained stacking can further reduce CV error, especially in M-open settings (Le et al., 2016).
• Distributional stacking (Bayesian predictive stacking) combines full predictive densities:

$p_{\mathrm{stack}}(y^*) = \sum_k w_k p_k(y^*|y)\,,$

with $w$ selected by maximizing LOO or K-fold cross-validated log-score (Yao et al., 2017).

• Neural network stacking (NNS/CNNS) uses a neural network to parameterize $w(x)$ as a smooth function of input features, optimizing a regularized MSE on meta-features from base learners (Coscrato et al., 2019).
• Bayesian hierarchical stacking treats $w(x)$ as a function of discrete or continuous covariates, with a hierarchical prior to control the degree of pooling and adaptation, fit by MCMC or variational Bayes (Yao et al., 2021).
• Double stacking for transfer learning partitions data into shards, stacks within each, then stacks across shards for scalable compositional inference (Presicce et al., 2024).
• Hybrid/Alternative pooling: Log-linear pooling ("locking") and quantum-inspired superposition ("quacking") replace or generalize convex mixtures, weights are fit via alternative proper scoring rules such as the Hyvärinen score (Yao et al., 2023).

Efficient cross-validation is crucial; Pareto-smoothed importance sampling (PSIS-LOO) is often employed for Bayesian models due to its computational scalability and accuracy diagnostics (Yao et al., 2017). In geospatial and spatiotemporal settings, closed-form conjugate posteriors and block updates support large-scale, parallelizable stacking (Zhang et al., 2023, Presicce et al., 2024).

3. Theoretical Properties and Consistency

Stacking possesses robust theoretical guarantees:

• Asymptotic optimality: Under regularity conditions, cross-validated stacking weights converge to minimize posterior expected loss. For the log-score, stacking converges to the mixture $w^*$ such that $p_{\mathrm{stack}}$ is closest (in Kullback-Leibler divergence) to the true data-generating process within the convex hull of candidate predictive densities (Le et al., 2016, Yao et al., 2017).
• Risk reduction: In regression with nested linear models, stacking with appropriate $l_1$ or $l_0$ complexity penalization yields strictly smaller expected risk than the best single model selected by standard criteria (e.g. AIC, BIC). The optimal stacked estimator is computed via a weighted isotonic regression (Chen et al., 2023).
• Stability: Stacking improves hypothesis stability over both bagging and subbagging, especially when using regularized or stable meta-learners. The stability of the ensemble is the product of stabilities of the base learners and the combiner, and is further improved via bootstrapping or subsampling strategies (Arsov et al., 2019).
• Exchangeability and conformal prediction: Stacked ensembles can be conformalized to achieve finite-sample valid predictive intervals, exploiting stability and cross-fitting to approximate marginal coverage without separate calibration (F, 18 May 2025).

These properties generalize across M-complete and M-open settings; the stacking objective adapts to either scenario and achieves consistency under broad conditions (Le et al., 2016).

4. Application in High-Dimensional and Complex Models

Stacking is especially impactful in modern settings where traditional marginal likelihood-based Bayesian model averaging (BMA) is either computationally prohibitive or statistically suboptimal:

• High-dimensional Gaussian Process (GP) regression: Sketching–stacking pipelines first project high-dimensional features through random or learned matrices, fitting multiple GPs on sketched features, then using predictive stacking to combine the resulting analytically tractable posterior predictions. This strategy bypasses the issues of poorly-identified variables, MCMC convergence failures, and computational non-scalability, dominating both MCMC-based and other ensemble methods in accuracy and speed for spatial tasks such as air pollution mapping (Gailliot et al., 2024).
• Spatial-temporal and non-Gaussian models: Predictive stacking with closed-form conjugate (Diaconis-Ylvisaker or GCM) priors is used to avoid iterative sampling over weakly-identified parameters, while still capturing full predictive uncertainty (Pan et al., 2024, Pan et al., 30 May 2025).
• Meta-learning and transfer learning for massive geospatial data: Double stacking allows split-and-combine analysis, supporting streaming architectures and rapid uncertainty-quantified inference across terascale datasets (Presicce et al., 2024).
• Machine learning pipelines: Stacking meta-learners (e.g., KN, SVM, MLP, tree ensembles) or implementing stacking with deep neural networks, supports both classification and regression with features ranging from images and text to high-dimensional sensor or simulation data (Chatzimparmpas et al., 2020, E et al., 2024).
• Probabilistic time series forecasting: Bayesian regression stacking enables robust ensemble forecasts with principled uncertainty quantification (e.g., Value-at-Risk), supporting group-specific partial pooling via hierarchical stacking (Pavlyshenko, 2022).

5. Practical Implementation: Algorithms and Pseudocode

The stacking workflow involves the following core steps (modulo context and specific variant):

1. Model construction: Fit each candidate predictive model $\mathcal{M}_m$ to the training data. Save either out-of-fold predictions or posterior predictive draws, as appropriate.
2. Cross-validated predictive evaluation: For each data point $i$ and model $m$, compute $p_m(y_i|y_{-i})$ (or leave-one-group-out, or out-of-fold predictions).
3. Weight optimization:
   • Stacking of means: Solve the quadratic program

   $$\min_{w\in\Delta_{M-1}} \sum_i (y_i - \sum_m w_m \hat{y}_{m,-i}(x_i))^2\,,$$

• Stacking of densities: Solve

  $$\max_{w\in\Delta_{M-1}} \frac{1}{n} \sum_{i=1}^n \log \left( \sum_m w_m\,p_m(y_i|y_{-i}) \right)$$

  via constrained convex optimization (e.g., projected gradient, MOSEK, CVXR, SciPy optimizers).

1. Prediction: For new input $x^*$, combine the pointwise predictions or predictive densities using the learned $w^*$ to produce the final ensemble output.

More advanced formulations (e.g., hierarchical, feature-adaptive, or double stacking) modify steps 2–3 to incorporate group/covariate hierarchies, nonstationarity, or recursively combine stacking solutions at different data partitions (Yao et al., 2021, Presicce et al., 2024).

6. Empirical Performance, Limitations, and Diagnostics

Stacking consistently demonstrates superior or equal predictive performance relative to BMA, model selection, or single best models—especially in misspecified settings and with large/complex data:

• In spatiotemporal applications (e.g., air pollution, wireless sensor grid, vegetation index), stacking provides prediction error and uncertainty coverage competitive with full MCMC inference at a fraction (less than 1/1000) of the computational cost (Zhang et al., 2023, Presicce et al., 2024).
• In large UCI regression datasets, deep neural stacking achieves lower MSE than both standard stacking and pure meta-learners, provided that sufficient sample size and model diversity are present (Coscrato et al., 2019).
• In transfer learning across massive distributed data, stacking allows transparent assimilation of local posteriors into global inference, with theoretical guarantees on consistency and uncertainty quantification (Presicce et al., 2024).
• Empirical studies confirm that stacking meta-learners enhances both stability and interpretability, provided proper cross-validation and regularization are used (Arsov et al., 2019, Chatzimparmpas et al., 2020).

Limitations and Considerations

• Stacking cannot rescue fundamentally misspecified candidate sets; predictive performance is limited by the capacity of underlying models (Wakayama et al., 2024).
• Weight regularization (e.g., Dirichlet priors, log-barriers, $l_2$ shrinkage) becomes crucial as the number of candidates increases or as the validation data shrinks (Yao et al., 2017).
• Diagnostics via PSIS $k$ values, LOO log-score variance, or the stability of stacking weights across partitions are recommended.

7. Software Ecosystem and Representative Methods

A range of toolkits and packages support predictive stacking workflows:

• loo (R), BayesBlend (Python): Comprehensive implementations of stacking of predictive distributions, pseudo-BMA(+), and hierarchical stacking with PSIS-LOO support (Yao et al., 2017, Haines et al., 2024).
• Stan: MCMC-based routines for Bayesian stacking and hierarchical stacking, including user-adjustable priors and support for covariate-dependent weights (Yao et al., 2021).
• StackGenVis: Visual analytics for interactive metric selection, feature and model pruning, and interpretability (Chatzimparmpas et al., 2020).
• Domain-specific packages: Spatio-temporal and geo-statistical stacking are available in recent spatial Bayesian toolkits (Zhang et al., 2023, Presicce et al., 2024).

Direct code-level APIs and workflows are published in (Yao et al., 2017, Haines et al., 2024), while algorithmic details (including variational and convex optimization) are available in the accompanying technical appendices.

---

References:

• (Yao et al., 2017): Yao et al., "Using stacking to average Bayesian predictive distributions"
• (Yao et al., 2021): Yao et al., "Bayesian hierarchical stacking: Some models are (somewhere) useful"
• (Zhang et al., 2023): Zhang et al., "Bayesian Geostatistics Using Predictive Stacking"
• (Presicce et al., 2024): Banerjee et al., "Bayesian Transfer Learning for Artificially Intelligent Geospatial Systems: A Predictive Stacking Approach"
• (Pan et al., 2024): Banerjee et al., "Bayesian Inference for Spatial-Temporal Non-Gaussian Data Using Predictive Stacking"
• (Chen et al., 2023): Chen–Klusowski–Tan, "Error Reduction from Stacked Regressions"
• (Arsov et al., 2019): Bousquet & Kock, "Stacking and stability"
• (Coscrato et al., 2019): Loureiro et al., "The NN-Stacking: Feature weighted linear stacking through neural networks"
• (Chatzimparmpas et al., 2020): Wanner et al., "StackGenVis: Alignment of Data, Algorithms, and Models for Stacking Ensemble Learning Using Performance Metrics"
• (Haines et al., 2024): Haines & Goold, "BayesBlend: Easy Model Blending using Pseudo-Bayesian Model Averaging, Stacking and Hierarchical Stacking in Python"
• (F, 18 May 2025): Marques, "Stacked conformal prediction"

The cited literature provides comprehensive mathematical, algorithmic, and empirical foundations for predictive stacking and its contemporary variants.