GPGam: SOTA generalized additive models via Autoresearch

Chandan Singh · September 2026

A GAM (generalized additive model) predicts by adding up one function per feature, y = f₁(x₁) + f₂(x₂) + ... + fₔ(xₔ). This form allows complete interpretability, as the contribution of each feature can be visualized as a curve. Nevertheless, GAMs maintain fairly strong performance, as f can be modeled nonlinearly, with SOTA GAMs generally modeling them via gradient-boosted trees. Improving the GAM-learning algorithm has been a longstanding open research question. Here, we introduce GPGam, a new GAM-learning algorithm based on Gaussian processes. GPGam is discovered through autoresearch with frontier agents and achieves significant improvents over the SOTA on limited-size held-out tabular regression.

1. Quickstart & example usage

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.204': ..., 'rbf-0.561': ...}, the learned shared scales
model.interaction_terms()      # [(6, 7), (1, 5)]

Reading the fitted model

Below is a GPGam fit to the California housing data. It predicts 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 R² of 0.825.

This is the whole model: eight curves, one per feature, and two interaction heatmaps. No parameters exist besides these curves.

Hover over a panel to read what it shows.

Fig 1. The whole model: a GPGam fit to California housing. Hover over a plot to explain it.

2. Experimental results

GPGam was produced by the autoresearch loop described in Agentic-imodels (see accompanying code on github). Specifically, Fable 5.1 was given the imodels-65 suite of benchmarks from there and instructed to autoresearch a GAM algorithm that maximized predictive performance. The autoresearch loop explored a diverse range of GAM solutions, steadily improving performance and finally culminating in GPGam.

Fig 2. Size and accuracy across the whole search: lines of code (bars) and mean rank on the 65-dataset suite (line, lower is better) at every committed version, from the shallow-tree baseline through spline GAMs and boosted GA2Ms to the Gaussian-process GAM.

Benchmarks

GPGam outperforms interpretable baselines across both visible datasets it saw during the autoresearch loop and held-out datasets where it was only evaluated at the end (see Table 1). During development it saw imodels-65, which contains 65 datasets which were subsampled to have no more than 50 features and 1000 rows. imodels-7 shows 7 key datasets from there evalauted at full size. The other datasets are held out and completely non-overlapping with imodels-65. In all cases, GPGam outperforms baseline interpretable models, including the SOTA GAM algorithm EBM. Each dataset is randomly split into an 80/20 train/test split.

Table 1. Mean rank of GPGam and ten other regression models on four benchmark suites, held-out suites de-duplicated against the development suite.

	imodels-65	imodels-7	TabArena-12	CTR23-23
	Visible		Hidden
	65 sets	7 sets	12 sets	23 sets
Interpretable
GPGam	3.45	2.86	2.50	3.13
EBM	3.94	3.57	3.00	3.48
RuleFit	6.78	7.86	7.42	7.52
Ridge	6.83	9.43	8.92	7.83
Hierarchical shrinkage	7.77	8.43	8.08	7.74
FIGS	7.86	7.29	6.75	6.39
Decision tree	8.11	8.29	8.42	8.48
Not interpretable
TabPFN*	2.51	1.86	3.75	2.87
Gradient boosting	4.78	4.57	4.67	4.57
Random forest	5.14	3.14	3.50	4.87
MLP	8.83	8.71	9.00	9.13

*Due to GPU memory constraints, TabPFN is capped at 2,500 training rows by GPU memory and 500 features, so its numbers may underestimate its performance. The checkpoint is the default regressor of TabPFN-2.5 (tabpfn-v2.5-regressor-v2.5_default.ckpt, tabpfn 7.0.0).

3. Methods: GPGam

Every curve in GPGam is a GP over the quantile bins of its feature, and those curves are fit by maximizing the exact GP marginal likelihood. There is no randomness, so re-fitting on the same data returns the same model.

Binning makes the exact likelihood cheap

Fitting a GAM as a GP normally means an n×n kernel matrix, which is infeasible past a few thousand rows. Every component of a GA²M is a function of one feature, or of two on a grid, so once the features are binned each component holds 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 3 quantities:

C = Z.T @ Z      # bin co-occurrence counts
b = Z.T @ y      # bin sums
y.T @ y

That is P×P numbers plus P more, where P is the total number of bins, which a budget holds to a few thousand. One pass over the data computes all of it, and every optimizer step after that uses those numbers alone. Each step therefore costs O(P³) regardless of the number of rows.

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.

Fig 1 fits one model to a synthetic dataset of 250 rows whose five features differ only in how smooth their effect is: a slow wave, a straight line, a sharp step, a rapid wave, and a feature with no effect.

Fig 3. One GPGam fit to 250 rows whose five features differ only in how smooth their effect is, ordered smoothest to roughest by the kernel balance the likelihood chose. Dots are partial residuals, the dotted line is the truth, and the band is the posterior two standard deviations wide. The same model, untuned, handles all five.

Slow wave and straight line. Almost all the weight goes to the smooth kernel: 621 times as much for the wave, 195 times for the line. There is nothing sharp to fit, so the rough kernel is left unused.
Rapid wave. The balance tips the other way: the rough kernel carries 1.4 times the weight of the smooth one, and the curve turns ten times across the range instead of smoothing the oscillation away. The shared lengthscales, learned from all five features at once, sit between what this feature and the slow wave would each choose alone, and the amplitude split does the rest.
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 43 times more weight here than it does for the slow wave, and that is what buys the jump.
No effect. The fitted curve spans 0.09 against about 2 for the others. Nothing removed this feature; the same likelihood that fit the rest found no reason to give it amplitude, and the hierarchical prior keeps it at a floor rather than exactly zero.

Three more choices that a GAM usually leaves to its user are made by that likelihood too. The two lengthscales are not fixed: one Matérn and one squared exponential lengthscale are learned and shared by every feature, with a weak prior around their starting values so that a dataset of forty rows cannot pull them anywhere. The likelihood consistently prefers smoother curves than the constants a search had settled on. Each kernel's log-amplitudes are also shrunk toward their centre across features, a hierarchical prior. Without it, which of two near-duplicate features carries an effect is decided by the train/test split, and the model's error on a redrawn split can change by half; with it the choice is stable.

Interactions. Interactions are screened on the residual of the main effects and the top candidates are fit as two-dimensional surfaces over grids of quantile bins. Up to 48 of them are fit jointly in blocks that share shrinkage, with the grid resolution of each block chosen by comparing marginal likelihoods. That joint fit is cubic in the number of cells, which put a hard cap on the number of interactions, and the cap turned out to be where the model lost to EBM on wide data: EBM's default fits five interactions per feature, and on a dataset with 48 features that is 240 surfaces against our 48. Above a thousand rows, GPGam now fits the remaining interactions by backfitting: each further surface is an exact 2-D Gaussian process fit to the residual of the joint model, its resolution picked by marginal likelihood, swept three times with a shared noise level. The cost is linear in the number of surfaces, so five per feature is affordable, and below a thousand rows nothing changes.

GPGam extends a line of work on GP-based additive models:

Work	How it relates
What the additive components are
Accuracy versus interpretability with GPs (Plate 1999)	Fits a series of GPs constrained to be progressively more additive, trading accuracy for readability. GPGam fixes that constraint at order two.
Bayesian functional ANOVA with GP priors (Kaufman et al. 2010)	GP prior on each term of an ANOVA decomposition, posterior by MCMC. Same decomposition as GPGam, which gets it in closed form from the sufficient statistics.
Additive GPs (Duvenaud et al. 2011)	One amplitude per interaction order, shared across features. GPGam gives each feature its own amplitudes and stops at order two.
Orthogonal additive kernels (Lu et al. 2022)	Constrains components to be orthogonal, which identifies the decomposition. GPGam does not constrain them, and instead reports each curve's variance about its own mean.
Hierarchical additive interaction GPs on grid data (Ishida et al. 2023)	Centred hierarchical ANOVA kernel with main and pairwise terms, made fast by Kronecker structure on gridded data. GPGam bins to a grid too, but reduces to `Z'Z`, `Z'y` and `y'y`, and does not centre.
How the GP is made affordable
KISS-GP (Wilson et al. 2015)	Places inducing points on a grid and interpolates onto it. GPGam assigns each row to a bin exactly, so its statistics carry no interpolation error.
Scalable GAMs with sparse variational GPs (Adam et al. 2018)	The same model: every component of the GAM is a GP. Reaches scale with sparse variational inducing points, where GPGam gets the exact marginal likelihood from binning.
GP regression for binned data (Smith et al. 2018)	Exact GP inference when observations arrive already binned, at a cost set by the bin count. GPGam bins deliberately to get the same property.
Hilbert space methods for reduced-rank GPs (Solin et al. 2020)	Approximates a kernel by a truncated eigenfunction basis. GPGam's bin grid is also a finite basis, but exact for the piecewise-constant function it defines.
GP Neural Additive Models (Zhang et al. 2024)	GP shape functions approximated by random Fourier features, fit by convex optimization. GPGam keeps the exact GP posterior over bins and fits amplitudes by exact marginal likelihood.
Additive GPs in applications
Additive GP regression for longitudinal data (Cheng et al. 2019)	Uses the additive GP decomposition itself as the interpretation, on longitudinal covariates.
Geographical GP GAM (Comber et al. 2023)	GP smooths parameterised by location, which makes the coefficients vary over space. Same construction as GPGam, applied to spatial rather than tabular features.

Citation

GPGam came out of the autoresearch loop described in Agentic-imodels. If you use either, please cite:

@misc{singh2026agenticimodels,
      title={Agentic-imodels: Evolving agentic interpretability tools via autoresearch},
      author={Chandan Singh and Yan Shuo Tan and Weijia Xu and Zelalem Gero and Weiwei Yang and Michel Galley and Jianfeng Gao},
      year={2026},
      eprint={2605.03808},
      archivePrefix={arXiv},
      primaryClass={cs.AI},
      url={https://arxiv.org/abs/2605.03808},
}