← All posts

Learning models that hold up on a population you never labelled

The Annals of Statistics, 2026 · with Peter Bühlmann, Zijian Guo · arXiv:2309.02211

The problem

You have labelled data from several hospitals, or several regions, or several years, and you want a model that works on a different population. Crucially, you are not flying blind about that population: you have plenty of its covariates — patient records, weather readings, sensor logs — you simply have none of its outcome labels. This is multisource unsupervised domain adaptation, and that unlabelled target sample is what makes the problem tractable.

The default move, pooling everything and minimizing average error, optimizes for whichever source you happened to sample most, and can fail exactly where the target differs most. What we want instead is a model that degrades gracefully across any plausible target consistent with what we have seen.

The obstacle is that the target's label mechanism \(\mathbf Q_{Y\mid X}\) is fundamentally unidentifiable — no amount of unlabelled data reveals it. So rather than guess it, we hold the observable covariate law \(\mathbf Q_X\) fixed and hedge over a set of plausible conditionals.

Scoring a model by what it explains

We score a model \(f\) on a distribution \(\mathbf T\) by its reward — how much better it does than predicting nothing at all:

\[ \mathbf R_{\mathbf T}(f) \;=\; \mathbb E_{(X,Y)\sim \mathbf T}\big[\, Y^2 - (Y - f(X))^2 \,\big]. \]

Then we take the uncertainty set to be every target that shares the observed covariates \(\mathbf Q_X\) but whose conditional is some mixture of the sources',

\[ \mathcal C(\mathbf Q_X, \mathcal H) \;=\; \Big\{\, \mathbf T = (\mathbf Q_X,\, \mathbf T_{Y\mid X}) \ :\ \mathbf T_{Y\mid X} = \textstyle\sum_{l=1}^{L} q_l\, \mathbf P^{(l)}_{Y\mid X},\ \ q \in \mathcal H \,\Big\}, \]

for a convex \(\mathcal H \subseteq \Delta^L\) — take \(\mathcal H = \Delta^L\), the whole simplex, if you know nothing. The robust model maximizes the worst-case reward over that set:

\[ f^\star_{\mathcal H} \;=\; \arg\max_{f\in\mathcal F}\ \min_{\mathbf T\in\mathcal C(\mathbf Q_X, \mathcal H)}\ \mathbf R_{\mathbf T}(f). \]

The characterization that makes it practical

The headline result (Theorem 2.1) is that this minimax model is nothing more exotic than a weighted average of the individual source models \(f^{(l)}\), with weights solving a small convex quadratic program:

\[ f^\star_{\mathcal H} \;=\; \sum_{l=1}^{L} q^\star_l\, f^{(l)}, \qquad q^\star \;=\; \arg\min_{q\in\mathcal H}\ q^\top \Gamma\, q, \qquad \Gamma_{k,l} = \mathbb E_{\mathbf Q_X}\big[f^{(k)}(X)\, f^{(l)}(X)\big]. \]

The reason is a one-line geometric fact: \(q^\top\Gamma q = \mathbb E_{\mathbf Q_X}\big[\big(\sum_l q_l f^{(l)}(X)\big)^2\big]\) is the squared \(L^2(\mathbf Q_X)\) distance from the aggregated model to the null model \(f \equiv 0\). So \(f^\star\) is simply the point of the convex hull of the source models lying closest to zero — the least aggressive model that still guards against every mixture in the set.

Because \(f^\star\) is an aggregate of per-source models, you can fit each \(f^{(l)}\) with any machine-learning algorithm you like — random forests, gradient boosting, neural networks — and then combine them by solving the tiny QP. That also gives the method a federated reading: sources never share raw data, only their fitted models, and the optimal weights \(q^\star\) tell you how much each source matters for the target at hand.

Why the reward, and not squared error or regret?

This is the part I find most instructive. Three natural losses give three different robust models, and only one of them is both indifferent to noise levels and computationally tractable.

The squared-error failure is worth seeing concretely. With two sources and \(\mathcal H = \Delta^2\), its robust model is \(q_1 f^{(1)} + (1-q_1) f^{(2)}\) with

\[ q_1 \;=\; 0 \vee \left( \frac12 + \frac{\sigma_1^2 - \sigma_2^2}{2\,\mathbb E_{\mathbf Q}\big[(f^{(1)}(X)-f^{(2)}(X))^2\big]} \right) \wedge 1 . \]

Read what this says. When the two sources have equal noise, \(q_1 = \tfrac12\) and you get the plain average. But as source 1 grows noisier, \(q_1 \to 1\) — the robust model collapses onto the model of the source with the most noise, purely because that source is the hardest to predict. Worst-case squared error confuses irreducible noise with poor fit. The reward, which subtracts off \(\mathbb E[Y^2]\), never makes that mistake.

All three coincide when there is only one source; they diverge only in the multisource regime. Regret's robust model is a Chebyshev-centre problem, NP-hard for general convex \(\mathcal H\) — and for a ball-shaped prior it ignores the ball's radius entirely.

Estimating the weights is the delicate part

In practice \(\Gamma\) must be estimated, and plugging in cross-fitted ML estimates \(\widehat f^{(l)}\) leaves a first-order bias. We subtract an estimate of it, reweighting each source's data by the density ratio \(\omega^{(l)} = d\mathbf Q_X / d\mathbf P^{(l)}_X\) to account for covariate shift:

\[ \widehat\Gamma_{k,l} \;=\; \frac1N \sum_{j=1}^{N} \widehat f^{(k)}\!\big(X_j^{\mathbf Q}\big)\,\widehat f^{(l)}\!\big(X_j^{\mathbf Q}\big) \;-\; \widehat{\mathcal D}_{k,l} \;-\; \widehat{\mathcal D}_{l,k}, \] \[ \widehat{\mathcal D}_{k,l} \;=\; \frac{1}{n_l}\sum_{i=1}^{n_l} \widehat\omega^{(l)}\!\big(X_i^{(l)}\big)\, \widehat f^{(k)}\!\big(X_i^{(l)}\big)\, \Big(\widehat f^{(l)}\!\big(X_i^{(l)}\big) - Y_i^{(l)}\Big). \]

Notice that the first term is where the unlabelled target sample \(\{X_j^{\mathbf Q}\}_{j\le N}\) earns its keep: \(\Gamma\) is an expectation under \(\mathbf Q_X\), so we can estimate it without a single target label. The correction provably shrinks the dominant error term whenever the base learners converge more slowly than \(n^{-1/2}\) — which is to say, for essentially every flexible ML method. And it is never worse than the plug-in, even when the density-ratio model is misspecified: an inconsistent \(\widehat\omega\) costs you the improvement, not the validity. In simulation the gap widens with sample size, exactly as the theory predicts.

Approximate median errors over 200 rounds, \(L=3\), \(N=20{,}000\), read from the box plots in the paper's Figure 6 to within about ±0.03. What matters is the pattern, not the third digit: the bias-corrected estimator with estimated density ratios sits essentially on top of the oracle that knows the true ratios, and its margin over the plug-in grows with \(n_l\) rather than shrinking.

Does it work? Air quality across Beijing

Measuring PM2.5 requires specialized equipment; measuring the weather does not. So a natural task is predicting PM2.5 from temperature, pressure, dew point, precipitation, wind speed and wind direction — and the relationship genuinely differs by district, because urbanization and industrial activity do. We take 5 monitoring sites as labelled sources and treat the remaining 7 sites as targets, using only their weather covariates. The score is the worst-case reward across those 7 targets, computed separately for each of the 16 year × season combinations from 2013 to 2016.

The single exception is winter 2015, where Group DRO edges ahead. ERM and importance weighting are weakest almost everywhere: they ignore the shift in the conditional outcome model. Group DRO cannot use the target's covariates at all — which is precisely the information DRoL exploits.

The weights turn out to be as useful as the predictions. They answer, per target, the question "which of my sources actually resembles this deployment population?" — an unusual thing to get out of a robustness method, and a directly interpretable one.

Aggregation weights \(q^\star\) read from the paper's Figure 8; rows need not sum to exactly 1 because of rounding, and dashes are weights at or near zero. Nongzhanguan concentrates almost four-fifths of its weight on Wanliu alone, while the other three targets hedge across three or four sources — a fact about the data that a pooled model would simply average away.

Finally, a little labelled target data goes a long way. Labelling just 5% of a target site (about 100 observations) and using it to centre the prior set \(\mathcal H = \{q : \|q - \widehat q^{\,\mathrm{label}}\|_2 \le \rho\}\) lifts the reward well above the default \(\mathcal H = \Delta^L\): at Nongzhanguan, across the four summers, the reward reaches ≈ 0.62–0.77 against ≈ 0.38–0.63 for the default, while fitting on those 100 points alone reaches only ≈ 0.25–0.45. A wrong prior with a tight radius, on the other hand, hurts — in simulation a uniform prior at \(\rho = 0\) scores ≈ 1.1 against the default's ≈ 1.5, and only recovers as \(\rho\) grows and the bad prior's influence washes out. The radius is the dial trading robustness against predictiveness, and it should reflect how much you actually know.

Why it matters

A model fit to pooled multisource data optimizes for the population you happened to sample most. This framework instead produces a model that performs well for every mixture of the sources, computes through standard quadratic programming, respects privacy constraints, and reports interpretable source-importance weights. It turns "will this model transfer?" from a hope into something you can optimize for directly — and it carries the classical maximin-effects idea beyond linear models to arbitrary machine learning.

Two honest caveats. First, the guarantee is worst-case over the mixture class, so with the full simplex and a genuinely irrelevant source the predictions can be conservative; the prior set \(\mathcal H\) exists precisely to buy that back. Second, what the theory delivers is high-probability convergence rates for the fitted model and its weights — not confidence intervals. Uncertainty quantification for this kind of worst-case estimator is genuinely hard, because the optimal weights can sit on the boundary of the simplex and destroy asymptotic normality. That problem is the subject of a later paper.

© 2026 Zhenyu (Zach) Wang 王振宇