DEV Community: Kasiuk Vadim

Uncertainty Estimates of Predictions via a General Bias-Variance Decomposition

Kasiuk Vadim — Thu, 14 May 2026 22:06:12 +0000

A General Bias‑Variance Decomposition for Proper Scoring Rules – Finally!

Or: Why your ensemble works, how to build confidence regions in logit space, and what Bregman information really does for uncertainty estimation.

If you’ve ever trained a classifier, you’ve heard the mantra:

Bias‑variance trade‑off.

But look closely – the classical decomposition works for squared error only.

What about log‑loss? Brier score? CRPS?

For years, we had no general, closed‑form bias‑variance decomposition for strictly proper scoring rules.

Until now.

In their AISTATS 2023 paper, Gruber & Buettner (PDF) finally fill this gap.

And they give us practical tools:

Explain ensembles via a law of total Bregman variance.
Build confidence regions directly in logit space.
Detect out‑of‑distribution inputs better than raw softmax confidence.

Let’s dive in.

The problem: Uncertainty under domain drift

Your model says “cat” with 0.99 probability – but the image is heavily corrupted.

You know from Ovadia et al. (2019) that softmax confidence is not reliable under dataset shift.

What we need is a variance‑based uncertainty measure that works for any proper loss.

And we need a theory that explains why – for example – ensembling always helps.

Missing piece: A general bias‑variance decomposition for strictly proper scoring rules.

Background: Bregman divergences & proper scoring rules

Bregman divergence

Given a differentiable convex function $ϕ\phi$ , the Bregman divergence is

d_\phi(x, y) = \phi(y) - \phi(x) - \langle \nabla \phi(x), y-x \rangle.

Example: $ϕ(x)=x2\phi(x)=x^2$ gives $dϕ(x,y)=(x−y)2d_\phi(x,y)=(x-y)^2$ (squared error).

Example: $ϕ(x)=xln⁡x\phi(x)=x\ln x$ gives the KL divergence.

Strictly proper scoring rule

A scoring rule $S (P, y)$ is strictly proper if the expected score is maximised only when $P$ equals the true data distribution $Q$ .

Common examples:

Log score: $S(P,y)=log⁡p(y)S(P,y)=\log p(y)$
Brier score: $S(P,y)=−∣δy−P∣2S(P,y)=-|\delta_y - P|^2$
CRPS (continuous ranked probability score)

Every strictly proper scoring rule corresponds to a Bregman divergence generated by the negative entropy $G$ (Ovcharov, 2018).

The main result: A general bias‑variance decomposition

Let $f^\hat{f}$ be a random prediction (e.g., from different training sets), and $\sim Q$ the true outcome.

Let $S$ be a strictly proper scoring rule with negative entropy $G$ , and $G^*$ its convex conjugate.

Theorem (Gruber & Buettner, 2023)

ParseError: KaTeX parse error: Expected group after '^' at position 98: …\underbrace{B{G^̲}[S(\hat{f})]}{…

What does each term mean?

$B_{G^*}[X]$ – Bregman information (generalised variance). For $ϕ(x)=x2\phi(x)=x^2$ , $Bϕ[X]=Var(X)B_\phi[X] = \mathrm{Var}(X)$ .
$d_{G^*, S^{-1}}$ – Bregman divergence in the dual space – that’s the squared bias.

So the classical MSE decomposition ( $error=noise+var+bias2\text{error} = \text{noise} + \text{var} + \text{bias}^2$ ) is a special case of this theorem.

Bregman information – the “variance” term

Definition (Banerjee et al., 2005):

B_\phi[X] = \mathbb{E}[d_\phi(\mathbb{E}[X], X)] = \mathbb{E}[\phi(X)] - \phi(\mathbb{E}[X]).

It measures spread around the mean in the sense of a Bregman divergence.

Figure 2 in the paper shows $Bσ+B_{\sigma_+}$ for the softplus function $σ+(x)=ln⁡(1+ex)\sigma_+(x)=\ln(1+e^x)$ – this controls variance for binary classification in logit space.

[!NOTE]
When $ϕ\phi$ is the squared function, $BϕB_\phi$ is the classical variance.

When $ϕ\phi$ is the log‑sum‑exp function (LSE), $BϕB_\phi$ is the variance in logit space.

Special case: Exponential families

For an exponential family $pθ(y)=exp⁡(⟨θ,T(y)⟩−A(θ))h(y)p_\theta(y) = \exp(\langle \theta, T(y)\rangle - A(\theta))h(y)$ , the decomposition becomes:

E[−ln⁡pθ^(Y)]=A(θ)+BA[θ^]+dA(θ,E[θ^]). \mathbb{E}[-\ln p_{\hat{\theta}}(Y)] = A(\theta) + B_A[\hat{\theta}] + d_A(\theta, \mathbb{E}[\hat{\theta}]).

$BA[θ^]B_A[\hat{\theta}]$ – variance in the natural parameter space (classical variance weighted by the log‑partition function $A$ ).
Perfectly recovers the classical MSE case when $A(θ)=θ2/2A(\theta)=\theta^2/2$ .

Special case: Classification (logit space) – this is huge

Let $z^∈Rk\hat{z} \in \mathbb{R}^k$ be the logits (before softmax).

Let $sm(z)\text{sm}(z)$ be the softmax probabilities.

Use the negative log‑likelihood (log loss) as scoring rule.

Corollary:

E[−ln⁡smY(z^)]=H(Q) + BLSE[z^] + dLSE(sm−1(Q), E[z^]), \mathbb{E}[-\ln \text{sm}Y(\hat{z})] = H(Q) \;+\; B{\mathrm{LSE}}[\hat{z}] \;+\; d_{\mathrm{LSE}}\big(\mathrm{sm}^{-1}(Q),\,\mathbb{E}[\hat{z}]\big),

where $LSE(x)=ln⁡∑exi\mathrm{LSE}(x) = \ln\sum e^{x_i}$ (LogSumExp).

Why is this surprising?

The variance term $BLSE[z^]B_{\mathrm{LSE}}[\hat{z}]$ is computed directly on the logits, without applying softmax.
No normalisation to probabilities needed – numerically stable and conceptually clean.

This is perfect for deep neural networks:

To estimate predictive uncertainty, just compute the Bregman information of the logits over an ensemble or multiple forward passes.

Applications

1. Why ensembles reduce uncertainty

The law of total Bregman information:

B_G[X] = \mathbb{E}[B_G[X \mid Y]] + B_G[\mathbb{E}[X \mid Y]].

For an ensemble that averages over random initialisations $W$ :

As number of ensemble members $\to \infty$ ,

BA[θ^D(n)]→BA[EW[θ^W,D]], B_A[\hat{\theta}D^{(n)}] \to B_A[\mathbb{E}_W[\hat{\theta}{W,D}]],

i.e., the variance due to $W$ disappears.

The expected score strictly improves.

This is the first general theoretical justification for why ensembles are almost always beneficial.

2. Confidence regions via Markov’s inequality

Using Markov’s inequality on the Bregman divergence:

P\Big(d_G(\mathbb{E}[X], X) \ge \frac{1}{\alpha} B_G[X]\Big) \le \alpha.

Thus a $(1−α)(1-\alpha)$ -confidence region is:

ParseError: KaTeX parse error: Expected group as argument to '\big' at position 92: …]}{\alpha} \big}̲.

Figure 3 & 4 in the paper:

Binary classification – confidence intervals on the probability simplex.
Iris dataset – convex confidence regions for three classes.

No need for normality assumptions – works with any proper score.

3. Out‑of‑distribution detection (CIFAR‑10C / ImageNet‑C)

Setup: Train on clean images, test on corrupted versions (CIFAR‑10C).

We want to discard uncertain predictions so that the remaining predictions have high accuracy.

Result (Figure 1 in the paper):

To reach 90% validation accuracy, using max softmax confidence you must discard ≈14% of data.
Using Bregman information $BLSEB_{\mathrm{LSE}}$ you only discard ≈7% of data.

→ Bregman information is a superior uncertainty measure under domain drift.

Limitations (real talk)

Computational cost: Estimating Bregman information requires multiple predictions per input (ensemble, MC dropout, or multi‑epoch sampling).
Proper scoring rules only: Doesn’t directly apply to 0‑1 loss (accuracy). But for probabilistic forecasting that’s fine – use log‑loss.
Not Bayesian: It gives a frequentist variance measure, not a full posterior.

Future work: extend to Bayesian neural networks and large language models (uncertainty for hallucinations).

Take‑away

First general closed‑form bias‑variance decomposition for strictly proper scoring rules.
Bregman information emerges as the universal variance term – generalising classical variance.
Logit‑space formulation makes it practical for deep learning.
Demonstrated benefits: ensembling theory, confidence regions, OOD detection.

Code available: GitHub – MLO‑lab/Uncertainty_Estimates_via_BVD

If you liked this, check out my previous post on Bayesian Neural Networks under covariate shift.

And let me know: how do you estimate uncertainty in your models today?

Bayesian Neural Networks Under Covariate Shift: When Theory Fails Practice

Kasiuk Vadim — Sun, 14 Dec 2025 19:25:54 +0000

October 22, 2025 | Machine Learning | Bayesian Methods

The Surprising Failure of Bayesian Robustness

If you've been following Bayesian deep learning literature, you've likely encountered the standard narrative: Bayesian methods provide principled uncertainty quantification, which should make them more robust to distribution shifts. The theory sounds compelling—when faced with out-of-distribution data, Bayesian Model Averaging (BMA) should account for multiple plausible explanations, leading to calibrated uncertainty and better generalization.

But what if this narrative is fundamentally flawed? What if, in practice, Bayesian Neural Networks (BNNs) with exact inference are actually less robust to distribution shift than their classical counterparts?

This is exactly what Izmailov et al. discovered in their NeurIPS 2021 paper, "Dangers of Bayesian Model Averaging under Covariate Shift." Their findings are both surprising and important—they challenge core assumptions about Bayesian methods and have significant implications for real-world applications.

The Counterintuitive Result

Let's start with the most striking finding:

Yes, you read that correctly. On severely corrupted CIFAR-10-C data, a Bayesian Neural Network using Hamiltonian Monte Carlo (HMC) achieves only 44% accuracy, while a simple Maximum a-Posteriori (MAP) estimate achieves 69% accuracy. That's a 25 percentage point gap in favor of the simpler method!

This is particularly surprising because on clean, in-distribution data, the BNN actually outperforms MAP by 5%. So we have a method that's better on standard benchmarks but catastrophically fails under distribution shift.

Why Does This Happen? The "Dead Pixels" Analogy

The authors provide an elegant explanation through what they call the "dead pixels" phenomenon. Consider MNIST digits—they always have black pixels in the corners (intensity = 0). These are "dead pixels" that never activate during training.

The Bayesian Problem

For a BNN with independent Gaussian priors on weights:

Weights connected to dead pixels don't affect the training loss (always multiplied by zero)
Therefore, the posterior equals the prior for these weights (they're not updated)
At test time with noise, dead pixels might activate
Random weights from the prior get multiplied by non-zero values
Noise propagates through the network → poor predictions

The MAP Solution

For MAP estimation with regularization:

Weights connected to dead pixels get pushed to zero by the regularizer
At test time, even if dead pixels activate, zero weights ignore them
Noise doesn't propagate → robust predictions

Formally, this is captured by Lemma:

If feature $x^i_k = 0$ for all training examples and the prior factorizes, then:

p(w^1_{ij}|\mathcal{D}) = p(w^1_{ij})

The posterior equals the prior, and these weights remain random.

The General Problem: Linear Dependencies

The dead pixels example is just a special case. The real issue is any linear dependency in the training data.

Proposition states that if training data lies in an affine subspace:

\sum_{j=1}^m x_i^j c_j = c_0 \quad \forall i

then:

The posterior of the weight projection $wjc=∑i=1mciwij1−c0bj1w_j^c = \sum_{i=1}^m c_i w^1_{ij} - c_0 b^1_j$ equals the prior
MAP sets $w_j^c = 0$
BMA predictions are sensitive to test data outside the subspace

This explains why certain corruptions hurt BNNs more than others:

The Brilliant Solution: EmpCov Prior

The authors' solution is both simple and elegant: align the prior with the data covariance structure.

The Empirical Covariance (EmpCov) prior for first-layer weights:

p(w^1) = \mathcal{N}\left(0, \alpha\Sigma + \epsilon I\right)

where

Σ=1n−1∑i=1nxixi⊤\Sigma = \frac{1}{n-1} \sum_{i=1}^n x_i x_i^\top

is the empirical data covariance.

How It Works

Eigenvectors of prior = Principal components of data
**Prior variance along PC $pi:ασi2+ϵp_i: \alpha\sigma_i^2 + \epsilon$
For zero-variance direction ( $σi2=0\sigma_i^2 = 0$ ): variance = $ϵ\epsilon$ (tiny)
Result: BNN can't sample large random weights along unimportant directions

The improvements are substantial:

Corruption/Shift	BNN (Gaussian)	BNN (EmpCov)	Improvement
Gaussian noise	21.3%	52.8%	+31.5 pp
Shot noise	24.1%	54.2%	+30.1 pp
MNIST→SVHN	31.2%	45.8%	+14.6 pp

Why Do Other Methods Work Better?

Here's the interesting part: many approximate Bayesian methods don't suffer from this problem. Why?

Method	Why Robust?	Connection to MAP
Deep Ensembles	Average of MAP solutions	Direct
SWAG	Gaussian around SGD trajectory	Indirect
MC Dropout	Implicit regularization	Indirect
Variational Inference	Often collapses to MAP-like solutions	Indirect
BNN (HMC)	Samples exact posterior	None

The common theme: most approximate methods are biased toward MAP solutions, which are robust due to regularization. HMC is unique in sampling the exact posterior, including problematic directions where posterior = prior.

Practical Implications

For Practitioners

Don't assume BNNs are robust: Test on corrupted/out-of-distribution data
Consider deep ensembles: They're often more reliable under shift
If using BNNs: Implement data-aware priors like EmpCov
Benchmark properly: Always include distribution shift evaluations

For Researchers

Re-evaluate Bayesian assumptions: The theory-practice gap needs addressing
Design better priors: Data-dependent priors are crucial
Study intermediate layers: The problem might not be limited to the first layer
Explore hybrid approaches: Combine BNNs with domain adaptation techniques

The Bigger Picture

This paper represents a paradigm shift in how we think about Bayesian methods:

BMA ≠ Automatic Robustness: Averaging over the posterior can actually hurt generalization under shift
Regularization Matters More: MAP's explicit regularization provides unexpected benefits
Context Matters: BNNs are great for calibrated in-distribution uncertainty but not for shift robustness

As the authors note, this problem affects "virtually every real-world application of Bayesian neural networks, since train and test rarely come from exactly the same distribution."

Conclusion

The "Dangers of Bayesian Model Averaging under Covariate Shift" paper is a must-read for anyone working with Bayesian methods or robustness. It:

Identifies a critical failure mode of BNNs under distribution shift
Provides theoretical understanding through linear dependencies
Offers practical solutions with data-aware priors
Challenges conventional wisdom about Bayesian robustness

The key takeaway: Bayesian methods are powerful tools, but they're not magic. Understanding their limitations—especially under distribution shift—is crucial for safe deployment in real-world applications.

As machine learning systems get deployed in increasingly diverse and unpredictable environments, papers like this remind us that robustness needs to be explicitly designed and tested, not just assumed from theoretical principles.

Reference: Izmailov, P., Nicholson, P., Lotfi, S., & Wilson, A. G. (2021). Dangers of Bayesian Model Averaging under Covariate Shift. Advances in Neural Information Processing Systems, 34.

Code: Available at GitHub

This post is based on the NeurIPS 2021 paper "Dangers of Bayesian Model Averaging under Covariate Shift." All credit goes to the original authors for their insightful work. Any errors in interpretation are mine.