AURORA

A Gaussian process is a model that tells you not just a prediction but how sure it is, which matters anywhere a wrong answer is expensive. The catch is cost: an exact GP scales with the cube of the training set. The usual fix is an approximation, like Random Fourier Features or Nyström, applied everywhere.

AURORA asked whether it has to be applied everywhere. It scores each part of the input space by uncertainty and how sparse the data is there, splits it into low, medium and high importance, and routes each to an approximation of matching cost: Random Fourier Features, then Nyström, then the exact GP. The December 2025 evaluation reported a 34–47% calibration improvement over uniform approximations, 41.1% on average, across five datasets.

When I went back to it in September 2026, the baseline gave it away. On the protein dataset, its error (RMSE 6.1346) was the same as the spread of the target itself (6.1182). In other words it explained nothing: it was predicting the mean.

The cause was a units bug. The kernel's lengthscale was estimated as the median distance between raw inputs, then applied to standardised inputs. On protein that's wrong by a factor of 155,000, which flattened the kernel into a constant. The reported gain was the distance between AURORA and a model that wasn't really modelling anything. With the lengthscale fitted in the space the kernel works in, and the noise fitted rather than hard-coded:

Checking the fix turned up a second problem, and it's the most useful thing in the project for anyone else. The code's “exact” GP used GPyTorch's fast predictive-variance setting. Below 800 training points it does nothing at all. Above 800, GPyTorch switches to iterative solvers and the variance is approximated, which is exactly the quantity calibration is measured from. Predictive standard deviations moved by a median of about 100%, and calibration error at 4,000 points came out three to ten times too high.

The first time I measured this, I blamed the setting itself. A later check of my own experiment found it had reused a cached covariance between runs, which muddled the picture. Measured cleanly, the thing that controls the error is the 800-point threshold, not the flag: raising GPyTorch's Cholesky size limit above the training set removes it, while turning the fast-variance setting off alone leaves about 5% error from the iterative solver. The threshold is a library default that never appears in user code, so anyone measuring calibration on an “exact” GP above 800 points can be badly off without knowing it.

A negative result on five benchmarks doesn't say much if those benchmarks never had regional structure to begin with. So I built datasets where the structure is true by construction, with noise varying fiftyfold between regions, and gave the method every advantage: the true partition and an oracle choosing the best model for each region. It still gained nothing, +0.0003 in calibration error on average.

The reason is the part that generalises. Sparse approximations like these change only one term of a GP's predictive variance and leave the noise term, which each model fits for itself, nearly the same. So they give almost the same uncertainty everywhere. On the constructed data, four models' per-region predictive standard deviations agree within 3% while the true noise varies 50 times. Varying Nyström's rank on five benchmarks changes accuracy by 21–155% but the average predictive uncertainty by at most 12%. Better approximations make better predictions; they don't make more honest uncertainty.

If the approximation is the wrong knob, the noise should be the right one: give each region its own noise level instead of one for the whole space. I ran that fairly on eleven datasets, with the same tuning budget for every method, and a single global noise parameter won calibration on nine. I wrote it up as a negative result, attacked it, and even found what looked like a crossover near 2,500 training points.

A pre-submission audit of the code against the paper found the reason, and it was mine. The out-of-fold squared residual I used as the noise estimate already contains the model's own uncertainty, and the regional arm then added that uncertainty again, making it under-confident. The heteroscedastic baseline took a geometric mean of squared residuals, which for Gaussian noise comes out at 0.28 of the true variance. With the estimator corrected, the result is the opposite: regional noise improves calibration on airfoil, sarcos and calhousing, holding up under multiple-comparison correction in two separate sets of ten splits, and worsens calibration on none. It does worsen likelihood on the smallest dataset, yacht, with 184 training points. The crossover was an artifact too: the corrected gap favours regional noise at every size from 300 to 4,000.

A second claim-by-claim check, three days later, tightened the rest. Results now carry multiple-comparison corrections; claims are scoped to the approximations actually tested; and run-to-run variation in the GPyTorch measurements is reported, since they aren't bit-for-bit reproducible (the NumPy results are). Along the way I found defects in my own audit that are worth knowing about:

Five of the six explanations I proposed along the way didn't survive their own predictions; the one that did rests on a direct measurement. Five headline numbers turned on how a baseline or estimator was configured, three in the original work and two in mine, and both of mine would have reversed a conclusion. The paper reports all of it, including a section on which claims failed and what caught them.

It's titled “Predictive Variance Decides: Region-Aware Gaussian Processes Should Adapt the Noise, Not the Approximation”, with Alekhya Bulusu. It's written, compiled and signed off, and the arXiv submission is the last step.