Module inference

Expand description

Statistical inference for Boulevard boosting: confidence intervals for the regression function f(x), prediction intervals for new labels, and reproduction intervals, with asymptotic (central-limit) guarantees under the assumptions below and validated only in the regimes listed in Validation. Opt-in: train with BoosterKind::Boulevard, then fit a BoulevardInference on the training rows. A Boulevard EBM (crate::ebm, ebm_boulevard) gets bands on its shape functions from EbmInference, built on the same solvers.

Unlike crate::conformal, whose intervals have finite-sample marginal coverage of the label, these intervals are about f itself and hold pointwise (conditionally on x), but only asymptotically and under the assumptions below.

§The algorithms

Boulevard (Zhou & Hooker, Boulevard: Regularized Stochastic Gradient Boosted Trees and Their Limiting Distribution, JMLR 23, 2022) averages its trees instead of summing them: after b rounds the ensemble is f_b = (λ / b) Σ_{i ≤ b} t_i, each tree fitted to the residuals of the current average on a row subsample. Fang, Tan & Hooker (Statistical Inference for Gradient Boosting Regression, NeurIPS 2025) add two variants that recover more of the signal, both implemented here:

BRAT-D (their Algorithm 1; num_parallel_tree = 1): round b drops each earlier tree independently with probability p (Boulevard::dropout; p = 0 is Zhou & Hooker’s Boulevard) and fits the new tree to y − μ − (λ / (b−1)) Σ_{kept} t_s(x), dividing by every earlier tree, not only the kept ones. The model predicts μ + ((1 + λq) / B) Σ t_b with q = 1 − p and learning rate λ = eta ∈ (0, 1].
BRAT-P (Algorithm 2; num_parallel_tree = K ≥ 2): the first iteration boosts K trees in sequence; afterwards tree k of round b fits y − μ − Σ_{l ≠ k} ā_l(x), where ā_l averages slot l’s earlier trees, so a round’s trees grow in parallel. The model predicts μ + (1/B) Σ_{b,k} t_{b,k}.

μ is the intercept (the label mean unless base_score is set). With Boulevard::truncation M > 0 the subtracted ensemble part is clipped to [−M, M], the Γ_M of the convergence proofs. The trained trees are ordinary RegTrees whose leaves already carry the final scale, so prediction, SHAP, slicing, and every export work as for a gbtree model; the native binary and JSON formats also keep the BoulevardInfo inference reads.

§The variance

As the number of rounds grows, both algorithms converge to a kernel ridge regression f̂(x) = μ + s k(x)ᵀ (c I + K)⁻¹ (y − μ 1) in the leaf kernel of the ensemble: K_ij averages, over the trees, 1 / (n_ℓ + κ) when rows i and j share a leaf ℓ holding n_ℓ training rows (κ = lambda / subsample, so a leaf’s value Σ z / (m + lambda) over its m ≈ ξ n_ℓ sampled rows is matched), and k(x) is the same average between x and the training rows. BRAT-D has c = 1 / (λq) and s = (1 + λq) / (λq); BRAT-P c = 1 / (K−1) and s = K / (K−1). The estimate is linear in y with weights w(x) = s u(x) + γ(x) 1, where u = (c I + K)⁻¹ k(x) and, for a label-mean intercept, γ = (1 − s 1ᵀu) / n, so under the regression model y = f(x) + ε with independent noise of variance σ²,

f̂(x) ≈ N(f(x), σ² ‖w(x)‖²)     (Fang, Tan & Hooker, Theorem 2)

The intervals (BoulevardInference::confidence_intervals and siblings) plug in an estimate σ̂² (NoiseVariance) and the normal quantile. The kernel is estimated by the trained trees themselves (the paper’s equations (1)–(2)), with each tree’s row sample replaced by its expectation, as the authors’ reference implementation does. KernelSolver::Exact factors the n × n system (O(n³) time, O(n²) memory); KernelSolver::Nystrom uses the paper’s Appendix A Nyström approximation from s uniformly sampled landmark rows (O(n s²) time, O(n s) memory).

§Assumptions

The asymptotic guarantees rest on the papers’ conditions; the estimates are computed regardless, so read them as approximations when these fail:

Regression with squared error, y = f(x) + ε with independent, homoscedastic, sub-Gaussian noise. booster = boulevard refuses every other objective, row weights, base margins, row sampling that depends on the labels or gradients (class-balanced bagging, gradient-based sampling), and every option that makes leaf values nonlinear in the labels (see TrainingParams::validate).
Structure–value isolation (the tree structures independent of the labels the leaves average): not true of trees grown greedily on the same labels. honest_refit provides it, refitting every leaf on an independent sample through the same Boulevard recursion (the NeurIPS paper’s “integrity”); fit the inference on that sample.
Non-adaptivity (tree structures eventually drawn from a fixed distribution), bounded leaf diameters and a minimal leaf size growing with n (set min_child_weight, which counts rows here), a row subsample, and, for BRAT-P, balanced splits. Both papers report that the intervals behave well in practice without enforcing all of them.
Enough rounds that the ensemble is near its limit: the variance is the limit’s.
Every variance here is conditional on the tree structures (the kernel is treated as fixed, as in both papers). How the structures, and with them the fit’s bias, vary between training samples is not included. In one dimension with small bias this is negligible; in several dimensions the true variance can be a multiple of the estimate and the intervals under-cover (see Validation).
Prediction intervals additionally need Gaussian noise: the estimate’s error is asymptotically normal, but the new label’s own noise is not averaged, so ± z σ̂ covers 1 − α of it only when that noise is normal (uniform noise, for one, is covered with probability 1 at 95%). Otherwise use calibrated_prediction_intervals, which rescales the widths by an empirical quantile, or crate::conformal.

§Validation

Simulations (50 training samples each, 100 fixed test points, coverage of the true f, σ̂² from a holdout of n/2 rows):

setting	`n`	90% CI	95% CI	95% PI
`f = sin 2πx + x²/2` (1-d), BRAT-D `λ = 0.6`, `p = 0.6`, `ξ = 0.6`, depth 8, 200 trees, `honest_refit`	500	0.899	0.949	0.959
same	1000	0.904	0.953	0.955
same	2000	0.901	0.952	0.955
same, Nyström `s = 1000`	4000	0.909	0.957	0.953
same, `p = 0` (Zhou & Hooker’s Boulevard)	1000	0.903	0.952	0.955
same, BRAT-P `K = 4`, 100 rounds	1000	0.883	0.936	0.956
same, without `honest_refit`	1000	0.756	0.845	0.957
`f = 4x₁ − x₂²` on `[0, 1]³` (the NeurIPS paper’s §6 test function), `λ = 1`, `p = 0.95`, `ξ = 1`, depth 6, 100 trees, `honest_refit`	1000	0.651	0.726	—
same	2000	0.626	0.711	—
same, without `honest_refit`	1000	0.543	0.632	—
Friedman #1 (5-d), BRAT-D `p = 0.6`, depth 6, `honest_refit` (20 samples)	2000	0.122	0.148	0.959

In the 1-d setting the estimated variance matches the across-sample variance of the estimate to within 5% (BRAT-P: 8% low). In the 3-d setting it is about half of it: 1.06 × when the tree structures are held fixed and only the refit sample is redrawn, 1.96 × when they are retrained, so the missing term is the structures’ sample-to-sample variation. In 5-d the fit’s bias dominates. The authors’ reference package (boulevard-boosting 0.1.0a1) gives the same coverage, interval widths, and MSE as this module in the settings compared (within seed-to-seed noise). Treat the confidence intervals as validated for low-dimensional smooth signals with honest refits only.

The NeurIPS paper’s variable-importance test (§4) is not provided: with the variances above its statistic is anti-conservative (82–94% rejection of a true null at a nominal 5% in the §6 setup with honest refits; 95% with the reference package’s own weights), and the paper’s reported size comes from a different regime (both fits on the same training sample, noise standard deviation 0.01, n ≤ 200, depth 8).

§Example

use hessboost::config::{BoosterKind, Boulevard};
use hessboost::inference::{BoulevardInference, KernelSolver, NoiseVariance};
use hessboost::prelude::*;

let n = 300;
let x: Vec<f32> = (0..n).map(|i| ((i * 37) % n) as f32 / n as f32).collect();
let y: Vec<f32> = x
    .iter()
    .enumerate()
    .map(|(i, v)| (6.0 * v).sin() + 0.1 * ((i * 7919 % 101) as f32 / 50.0 - 1.0))
    .collect();
let all = DMatrix::from_dense(&x, n, 1)?.with_labels(&y)?;
let (fit_rows, cal_rows): (Vec<usize>, Vec<usize>) = (0..n).partition(|i| i % 3 != 0);
let (dtrain, dcal) = (all.select_rows(&fit_rows)?, all.select_rows(&cal_rows)?);

let params = TrainingParams::builder()
    .booster(BoosterKind::Boulevard(Boulevard::builder().dropout(0.5).build()?))
    .eta(0.8)
    .subsample(0.8)
    .max_depth(3)
    .min_child_weight(5.0)
    .build()?;
let model = train(&params, &dtrain, 100)?;

let inference = BoulevardInference::fit(
    &model,
    &dtrain,
    NoiseVariance::Holdout(&dcal),
    KernelSolver::Exact,
)?;
let ci = inference.confidence_intervals(&dcal, 0.05)?;
let pi = inference.prediction_intervals(&dcal, 0.05)?;
assert!(ci.iter().zip(&pi).all(|(c, p)| p.lower < c.lower && c.upper < p.upper));

Structs§

BoulevardInference: The variance machinery of one Boulevard model: its leaf kernel over the training rows, the factored ridge system, and a noise estimate. Fit once with fit, then query any rows.
BoulevardInfo: How a booster = boulevard model was trained: the settings its inference and honest_refit read, recorded by training (BoostedModel::boulevard). Whether it is BRAT-D or BRAT-P follows from the model’s num_parallel_tree (1 or more).
EbmInference: The variance machinery of a Boulevard EBM (Ebm::boulevard): the additive term kernels over the training rows, their factored ridge systems, and a noise estimate. Fit once with fit, then query shape-function bands (term_bands) or intervals for f(x).
TermBands: Pointwise confidence bands of one term’s shape function on its grid (EbmInference::term_bands), aligned with TermShape::values.

Enums§

KernelSolver: How BoulevardInference::fit solves the kernel ridge systems.
NoiseVariance: Where BoulevardInference::fit takes the noise variance σ² from.

Constants§

MAX_EXACT_ROWS: Largest training set KernelSolver::Exact factors (its n × n system takes 8 n² bytes: 512 MiB here).

Functions§

honest_refit: Refit every leaf of the Boulevard model model on values, labelled rows independent of its training data, keeping every tree’s structure: the Boulevard recursion the model was trained with (BRAT-D with its dropout, or BRAT-P, with the same learning rate, row subsample ratio, L2 penalty, and truncation) is rerun on values, each tree’s leaves set to Σ z / (m + lambda) over the m rows of a fresh row sample reaching them (0 for a leaf none reaches), and the result scaled as training scales it. A label-mean intercept is re-estimated on values.