GPGam: Additive Gaussian Processes over Binned Features
1. Overview
A generalized additive model predicts by adding up one function per feature,
y = f₁(x₁) + f₂(x₂) + ... + fₔ(xₔ). Each
f can take any shape, so the model can say that income matters
greatly up to a point and hardly at all beyond it, or that risk rises at both
ends of a range. A linear model gets one number per feature and has to spend it
on a slope. What makes the form worth having is that the sum can be read: one
curve per feature, that curve is the whole story for that feature, and reading
all of them leaves nothing out. With a boosted ensemble you can compute an
explanation afterwards, but the explanation is a summary of the model rather
than the model itself, and two reasonable summary methods can disagree about the
same prediction. Here the thing you read and the thing that predicts are one
object.
GPGam makes every one of those curves a Gaussian process over the
quantile bins of its feature. Binning is what makes this practical: once the
features are binned, the exact GP marginal likelihood depends on the data only
through the bin co-occurrence counts, the bin sums, and y'y, so
every optimizer step costs the same whether the data has a thousand rows or a
hundred thousand. That one likelihood then settles every choice a GAM usually
leaves to its user: how smooth each curve should be, which features earn a curve
at all, which pairs of features get an interaction surface, and how finely to
grid it. There is nothing to tune, nothing is chosen by cross-validation, and no
seed enters anywhere, so two fits on the same data give the same model.
Across four regression benchmarks covering 113 datasets, GPGam is the strongest interpretable model on three of them and second by a small margin on the fourth, ahead of explainable boosting machines, rule sets, shrunk trees and linear models. Only TabPFN, a black box that cannot be fit at all on several of the wider datasets, ranks higher overall. The model was not designed by hand: it was found by an autoresearch loop, and the section on how it was found shows what that loop added and later took away.
The cost of the form is additivity. A pure additive model cannot represent a feature whose effect depends on another feature. GPGam is a GA2M, which buys back part of that with a few two-feature surfaces and picks which pairs to spend them on. What it still cannot represent is an effect that depends on three or more features at once, and the results section is specific about where that limit shows up.
2. How GPGam works
Every curve in GPGam is a Gaussian process over the quantile bins of its feature, and every choice about those curves is made by maximizing one quantity: the exact GP marginal likelihood. Nothing is chosen by cross-validation, there is no validation split, no bagging and no seed, so two fits on the same data give the same model.
Binning makes the exact likelihood cheap
Fitting a GAM as a Gaussian process normally means an n×n kernel matrix, which is hopeless past a few thousand rows. Every component of a GA2M is a function of one feature, or of two on a grid. So once the features are binned, each component is just a function on a small grid, holding one value per bin.
Write Z for the indicator matrix recording which bin each row falls
in. Then the exact marginal likelihood touches the data only through
three quantities:
C = Z.T @ Z # bin co-occurrence counts
b = Z.T @ y # bin sums
y.T @ yThat is P×P numbers plus P more, where
P is the total number of bins. A budget holds P to a
few thousand. One pass over the data computes all of it, and every optimizer
step after that uses those numbers alone. So each step costs O(P³)
no matter how many rows there are. A hundred thousand rows cost what a
thousand rows cost.
What the likelihood decides
Each feature gets a prior built from two covariance kernels on its bin grid. A Matérn ½ kernel at a short lengthscale produces rough shapes that can turn sharply. A squared exponential kernel at a longer lengthscale produces smooth ones. Their amplitudes are free parameters, so maximizing the likelihood picks the blend for each feature separately.
The figure below shows what that means in practice. It is one model fit to one synthetic dataset of 4,000 rows whose five features differ only in how smooth their effect is: a slow wave, a rapid wave, a sharp step, a straight line, and a feature with no effect at all. The same model, with no tuning and no per-feature settings, has to handle all five.
- Slow wave and straight line. Almost all the weight goes to the smooth kernel: 722 times as much for the wave, 339 times for the line. There is nothing sharp to fit, so the rough kernel is left unused.
- Rapid wave. The balance flips. The rough kernel takes 83 times the weight of the smooth one, which is what lets the curve turn ten times across the range instead of smoothing the oscillation away.
- Sharp step. A mixture, and the interesting case. Most of the weight is smooth, because the function is flat almost everywhere, but the rough kernel gets 47 times more weight here than it does for the slow wave, and that is what buys the jump. The fitted step spans 1.41 against a true 1.40, and turns at 0.538 against a true 0.55.
- No effect. Both amplitudes collapse toward zero and the curve is flat. Nothing removed this feature; the same likelihood that fit the others simply found no reason to give it any amplitude.
Those same amplitudes are how the model selects features. A feature that explains nothing has its amplitude driven to zero and drops out, which is automatic relevance determination, obtained free from the same objective. Interactions are screened on the residual and then fit in blocks that share shrinkage, with the grid resolution for each block chosen by comparing marginal likelihoods.
3. How well it does, and how it was found
Benchmarks
Four regression benchmark suites were used, covering 113 datasets that range from 50 to 72,000 rows and up to 1,024 features. GPGam is the strongest interpretable model on three of them, and ties explainable boosting machines on the fourth.
| model | imodels-65 | classic-7 | TabArena-13 | CTR23-28 |
|---|---|---|---|---|
| 65 sets | 7 sets | 13 sets | 28 sets | |
| Interpretable | ||||
| GPGam | 3.58 | 2.57 | 2.85 | 3.36 |
| EBM | 3.91 | 3.57 | 2.77 | 3.46 |
| FIGS | 7.83 | 7.29 | 6.69 | 6.64 |
| RuleFit | 6.74 | 7.86 | 7.38 | 7.61 |
| Ridge | 6.83 | 9.43 | 8.77 | 7.75 |
| Hierarchical shrinkage | 7.71 | 8.43 | 8.23 | 7.96 |
| Decision tree | 8.14 | 8.29 | 8.46 | 8.71 |
| Not interpretable | ||||
| TabPFN | 2.52 | 1.86 | 3.54 | 2.46 |
| Gradient boosting | 4.78 | 4.71 | 4.69 | 4.61 |
| Random forest | 5.15 | 3.29 | 3.46 | 4.68 |
| Neural net | 8.80 | 8.71 | 9.15 | 8.75 |
Mean rank across the datasets of each suite, lower is better, over the same eleven models everywhere. Ranks are only comparable down a column, not across rows of different pools. Bold marks the best model within its section. TabPFN could not be fit on two TabArena datasets, one with 1,024 features and one with 81, and is ranked last on those.
What was subsampled. The four columns do not measure the same thing, and it is worth being precise about how. The imodels suite cuts every dataset down to at most 1,000 training rows and 50 features at a fixed seed, so that column measures the small data regime rather than the datasets at their real size. The other three suites use every row and column, from 517 rows up to 72,000, and up to 1,024 features.
TabPFN is handicapped in the three full-size suites. It is capped at 2,500 training rows by GPU memory, and on a dataset wider than 500 features only the first 500 are passed to it. Where it still ranks first, it does so on a fraction of the data the other models saw; where it ranks last on two TabArena datasets, that is because no fit was possible at all. Its numbers should be read as a lower bound on what the method can do rather than a measurement of it.
Every result here is one 80/20 split at a fixed seed with all models at library defaults, which is not how these benchmarks are meant to be scored. TabArena and CTR23 both define protocols with repeated folds and hyperparameter search, so none of these numbers are comparable to published leaderboard entries. CTR23 is also cut from 35 datasets to 28, dropping the seven it shares with TabArena so the two columns stay independent.
The model fails in the same way on every suite, which is worth knowing before you use it. It loses on smooth simulations and on data with strong interactions, such as wave energy converters, building energy simulation, and robot arm dynamics. In all of those the target depends on three or more inputs at once, and no sum of pairwise pieces can represent that. This is a limit of the model class, not of the fitting procedure.
How the model was found
GPGam was not designed by hand. It came out of the autoresearch loop described
in Agentic-imodels: Evolving agentic
interpretability tools via autoresearch, and the code that produced it is in
the agentic-imodels
repository, under evolve_gam/. The loop was given a benchmark, a
rule that the model stay additive with at most pairwise terms, and one
instruction: beat explainable boosting machines. The loop proposed a change,
measured it on 65 regression datasets, and kept it only if the numbers held.
Every idea that failed was recorded with the cost of failing, which is why the
parameter list is short: anything that could not pay for itself was removed.
Nothing in that instruction mentions Gaussian processes, binning, or marginal likelihood. Those came out of the search.
The two halves of that plot are different activities. Up to v38 the loop was adding capability: v35 is the first version that fit a GAM as a Gaussian process, but it formed a kernel matrix over the training rows and so could only handle small data. Versions 36 to 38 added the binned formulation that made the exact likelihood affordable at any size, then bolted a boosted tree stage onto the residual. The file roughly doubled, from 507 lines to 1,110.
From v39 the instruction changed to making it simpler, and every later version removed something. Thirteen rounds of ablation took out the boosted tree stage, a second model class that handled small data separately, half the kernel dictionary, both MAP priors, a binary search over bin counts, and several parameters that turned out to do nothing. The file came back down to 686 lines, below where it started, while keeping the ability to scale that the middle versions bought.
That cost something, and the plot shows it honestly. Mean rank on the small benchmark drifts from 3.71 to 4.60 across the simplification rounds, almost all of it from dropping the separate small-data model in v41. The trade was deliberate: one model instead of two, and still ahead of EBM at 5.00 on the same pool. Each removal was measured first, and the ones that cost real accuracy were put back. Automatic relevance determination, the interaction fitting, and the log target rule all survived that way, having failed to be removable.
The prompts that drove the search, the benchmark harnesses, and the
per-dataset results are all kept alongside the model in
evolve_gam/.
4. Using it
from imodels import GPGamRegressor
from sklearn.datasets import fetch_california_housing
from sklearn.model_selection import train_test_split
X, y = fetch_california_housing(return_X_y=True)
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=42)
model = GPGamRegressor(n_pairs=2).fit(X_train, y_train)
preds = model.predict(X_test)The fitted model is a set of lookup tables, so you can read it directly instead of explaining it after the fact:
grid, values, std = model.shape_function(0, return_std=True) # a curve and its band
model.kernel_weights(0) # {'matern-0.05': ..., 'rbf-0.25': ...}
model.interaction_terms() # [(6, 7), (1, 5)]Reading the fitted model
Below is a GPGam fit to the California housing data: 16,512 training rows, eight features, predicting median house value in hundreds of thousands of dollars. It was asked for two interaction terms, and nothing else was set. On a held-out fifth of the data it reaches an R2 of 0.8247.
This is the whole model. Eight curves, one per feature, each with the Gaussian process posterior drawn two standard deviations wide, and two interaction surfaces, one of them over the state it describes. Nothing else is fitted, so there is nothing else to explain. Click any panel to read it.
- Longitude. Peaks near −122.5°, San Francisco, and falls away eastward. At $452k the widest single curve here.
- Latitude. Falls from south to north and spans $422k, but read it with the map panel: most of the geography lives in the joint term.
- MedInc. Rises across the whole income range, a spread of $203k, steepest through the middle and flattening at both ends.
- AveOccup. Falls from about +0.5 at 1.4 people per household to −0.4 by 4, then flattens. Crowding costs, but only up to a point.
- AveRooms. Almost flat from 2 to 6 rooms, rising only at the top end. Room count separates mansions, not ordinary houses.
- HouseAge. Weak alone, about $19k. Its real contribution is in the interaction panel, where age matters only at low occupancy.
- AveBedrms. Nearly flat, about $18k end to end. Its amplitude was driven near zero rather than removed by any selection step.
- Population. Nearly flat, about $7k, pruned by its own amplitude.
- Latitude × Longitude. The largest term in the model after the two geographic curves. It adds up to +1.14 over the Los Angeles basin and subtracts up to −0.86 in rural northern California, which is what neither marginal curve can express: latitude alone calls the south expensive and misreads the Bay Area, longitude alone calls the coast expensive and misreads the far north. Only the joint term can say coastal and southern. Fewer than half its cells hold data, since the state is a diagonal strip inside its bounding box, and empty cells are left blank rather than drawn as zero.
- HouseAge × AveOccup. Smaller and more specific: house age matters almost entirely among sparsely occupied block groups. Along the age axis at about 1.9 occupants the term swings from −0.37 to +0.28; at a typical 2.8 occupants it barely moves.
The curves are ordered by how far they move. The last two barely move at all:
Population and AveBedrms span under $20k end to end.
No selection step removed them; their kernel amplitudes were driven near zero by
the same likelihood that fit everything else, which is what automatic relevance
determination does.
Every band widens toward the ends of its curve. That is an edge effect rather than a density one, since quantile bins hold roughly equal counts by construction: a bin in the middle of a grid borrows strength from neighbours on both sides, while a bin at the end has neighbours on one side only.
The two surfaces are worth +0.0542 in R2. Without them the same model scores 0.7705, with them 0.8247. Fitting all 24 available pairs instead of two adds only a further +0.014, so nearly all the interaction structure in this data sits in these two terms.