Kernel Autoresearch for Open-Ended Model Discovery

Kernels encode the inductive bias of a wide range of machine learning (ML) models, from Gaussian processes (GPs) and support vector machines to kernel ridge regression and kernel two-sample tests. A kernel defines which inputs a model treats as similar, so it determines what the model can learn from data.

The choice of a kernel matters most when data are scarce. In Bayesian optimization (BO) and experimental design, where each evaluation is costly or time-consuming, we need models that predict well from only a few observations. With so little data, a poorly matched kernel cannot be corrected by more observations, and every later decision inherits its mistakes.

Automated kernel design usually searches a fixed grammar. It starts from base kernels such as RBF and Matérn and combines them with sums and products, guided by greedy selection, Bayesian optimization, neural composition, or language models. A grammar guarantees that every candidate is valid, but it also limits the search to structures that the base kernels and operators can express.

A natural alternative is to let coding agents write kernels as programs, which removes that limit. It also removes the guarantee. An arbitrary program need not be symmetric and positive semidefinite, and a numerical check on sampled inputs cannot prove that it is. A kernel that fails on other scales or dimensions silently breaks the model built on it. The problem is to search an open space of kernel programs while keeping every accepted kernel valid.

Kernaut (Kernel Autoresearch) addresses this with construction contracts. The agent writes only a component, and trusted backend code assembles the kernel with rules that preserve validity. The rest of this page describes the method, tests whether the discovered kernels transfer to tasks withheld from search, and examines what one of them encodes.

Kernaut treats that similarity rule as a short Python program. The agent writes a component, and trusted backend code assembles the kernel through a construction contract.

The four contracts let agents write new feature maps and input transformations, as well as combine existing kernels:

• Feature map. Map each input to a vector. The interpreter takes the inner product of two vectors.
• Spectral. Propose frequencies and raw weights. The interpreter squares the weights and combines cosine components.
• Pullback. Transform each input, then apply a trusted library kernel to the transformed inputs.
• Closure. Combine trusted or accepted kernels through sums, products, nonnegative scaling, and input transformations.

These constructions preserve positive semidefiniteness (PSD), the matrix property required of a kernel. Kernel validity also requires finite outputs, correct dimensions, and pure components. A pure component depends only on its input and parameters, without batch dependence, randomness, or state from earlier calls.

Each discovery campaign runs a sequence of fresh agent sessions. A session either creates a program or revises an archived program. A MAP-Elites archive retains strong candidates with distinct kernel behavior ¹. Archived examples and structured rejection reasons guide later proposals. Numerical fingerprints discourage candidates that behave like known kernels.

Search and parameter tuning use meta-training data only. Validation then selects one frozen program, which transfers to held-out tasks without structural changes or further LLM calls. Each new GP fits only covariance amplitude and observation noise.

Contracts are needed because numerical checks alone cannot certify an unrestricted program. A covariance matrix can pass a PSD test on one set of inputs and fail on another. In three unrestricted-program cohorts, additional checks falsified 22–58% of the programs that passed the initial screen.

Kernaut checks execution and output shapes at Tier 0, sampled PSD at Tier 1, and construction plus consistency at Tier 2. The validity argument comes from trusted assembly under the contract assumptions. Source inspection and behavioral checks test conformance, but do not prove arbitrary Python programs pure over every input.

The experiments ask three questions. Does one frozen kernel structure transfer beyond the tasks used during search? Does open-ended synthesis matter? What do the discovered kernels encode? We test transfer on black-box optimization, greenhouse-gas forecasting, enzyme kinetics, and glucose dynamics. We report the continuous ranked probability score (CRPS), which evaluates a predictive distribution against the observed value ². Lower CRPS is better. Optimization uses regret AUC, the mean normalized simple regret over BO steps, where lower is also better.

We start with black-box optimization (BBO), the setting that motivated this work. Search and validation use five function families. Testing uses 60 tasks from six other families, with dimensions up to eight. Each predictive episode fits on 8 + 2d observations and scores 32 held-out points, where d is the input dimension. Each discovery campaign has 12 runs.

Validation selected the dual warp–fold (DWF) kernel in the first ensemble campaign. DWF combines a mild coordinate warp with a triangular fold. Its test CRPS was 0.601, compared with 0.612 for fixed Matérn-5/2. Its regret AUC was 0.671, compared with 0.693.

DWF is a representative discovery. It is not a sum or product of standard kernels. The agent wrote only the pointwise transform below, which warps each coordinate slightly and folds it around a mirror point, so values at mirrored inputs are encouraged to agree. The trusted interpreter applies a Matérn-5/2 kernel to the output, so validity follows from the pullback contract and not from the agent’s code. The section on human refinement analyzes its geometry.

def transform_point(x, parameters):
    x = np.asarray(x, dtype=float).ravel()
    warp_rate = float(parameters["warp_rate"])
    fold_frequency = float(parameters["fold_frequency"])
    fold_phase = float(parameters["fold_phase"])

    # mild monotone warp of each coordinate
    warp = (1.0 - np.exp(-warp_rate * x)) / (1.0 - np.exp(-warp_rate))

    # triangular fold: mirrored inputs get similar values
    fold = (1.0 / np.pi) * np.arcsin(
        np.sin(np.pi * fold_frequency * x)
    ) * np.sin(fold_phase)

    return np.concatenate([warp, fold])

Across seven ensemble campaigns, six selected kernels improved test CRPS over Matérn-5/2. All seven GPT-6 Astra campaigns improved CRPS, with gains from 5.8% to 11.2%. These repeated campaigns measure discovery variability, whereas Figure 3 compares individual frozen kernels.

We also compared DWF with Few-Shot Bayesian Optimization (FSBO), which meta-learns a deep kernel ³. At the matched budget of 518 meta-training points, DWF had lower CRPS and regret AUC in both FSBO evaluation modes. At 64,000 points, neither FSBO mode differed significantly from DWF on either metric.

Transfer alone does not show that open-ended programs are necessary. The matched ablation keeps the proposer ensemble, run budget, and evaluation protocol fixed. Restricting agents to compositions of library kernels raises test CRPS from 0.571 to 0.618. Full search wins 11 of 12 paired campaigns against this closure-only arm.

Independent proposals also perform worse, but that arm changes feedback, novelty screening across runs, and candidate retention together. The comparison measures their combined contribution, rather than isolating archive sharing.

To move beyond optimization benchmarks, we forecast real greenhouse-gas records. Search uses monthly CO₂, CH₄, and N₂O records. Testing uses SF₆, CFC-12, and CFC-11, with 30 windows per record and a 48-month horizon. Some windows reverse time, so the benchmark includes both forecasting and backcasting.

Validation selected a residual period-bank kernel. It combines a Matérn component, a smooth residual component, and low-frequency sinusoids that act as slow trends. The selected kernel ranks first among 58 frozen candidates on both CFC records.

Across all three records, the selected kernel has mean CRPS 0.091, compared with 0.120 for AutoGP, which infers a structure for each window ⁴. The overall gain comes mainly from reversed windows. On forward windows, its CRPS is 0.034, behind CAKE at 0.022. Windows overlap within each record, so these are not independent record-level replications.

Enzyme kinetics. ChemBench maps seven experimental variables to a reaction rate. Search uses ten rate-law families. Testing uses 75 episodes from five unseen mechanisms, with 20 observations for fitting and 20 for scoring in each episode.

The validation-selected frontier kernel reaches CRPS 0.312, compared with 0.567 for optimized automatic relevance determination (ARD). ARD learns a separate lengthscale for each input. The selected kernel improves over all six learned references in paired tests with Holm correction. The selected ensemble kernel is competitive with the strongest learned references.

The frontier program combines substrate-response shapes inspired by the training mechanisms, local sensitivities, and modulation by enzyme loading. These features provide an inspectable biochemical prior. Predictive transfer does not establish recovery of the true mechanism.

Glucose dynamics. GlucoseBench uses simulated patients. Search uses ten children, validation uses ten adolescents, and testing uses ten adults. Each forecast uses only training episodes from its target patient, so the experiment tests transfer of kernel structure.

With two training episodes, the five discovered kernels reduce geometric-mean normalized MSE by 38–67% against relative-time ARD-Matérn. With six episodes and a 45-minute prefix, the reductions narrow to 0–23%. The selected ensemble kernel encodes a delayed meal response and a fading insulin effect.

These glucose results describe ten simulated adults. Across 35 paired patient-level comparisons, no adjusted p-value falls below 0.05 after Holm correction. The reported percentage reductions therefore describe the observed benchmark means.

Across the four domains, frozen discovered kernels transfer to tasks withheld from search, although the glucose gains are descriptive and forward forecasting windows favor CAKE. Because the kernels are short programs, we can also read and refine them. We examine DWF, introduced in the black-box optimization results. Its warp stays near the identity, while its fold bends the input geometry around a peak at 0.513. The original program applies a Matérn kernel jointly to the warped and folded coordinates.

A human refinement gives the warp and fold separate kernels, then averages their covariances. This additive construction preserves PSD and strengthens the reflection preference.

The CRPS gain has a paired 95% bootstrap interval excluding zero. The regret reduction remains descriptive because the multiplicity-adjusted test does not cross 0.05. Smoothing the fold preserves the CRPS gain, so the sharp corner is unnecessary for this improvement.

The experiments test small-data Gaussian processes. The BBO benchmark has at most eight observed dimensions, and the scientific benchmarks use simulated mechanisms or patients. Natural next benchmarks include higher-dimensional, constrained, and mixed categorical and continuous optimization. Larger datasets are also within reach, because feature-map and spectral contracts produce explicit finite features that reduce exact inference from cubic to linear cost in the number of observations.

Each construction contract pairs a component interface with an assembly rule and a proof of validity, so new contracts widen the program space without changing the search loop. Candidates include multi-output kernels through intrinsic coregionalization, feature maps for sequences and graphs, continuous spectral measures, state-space kernels from linear ODEs, and recurrent feature maps over event histories such as meals and insulin doses. The last two would let the agent propose mechanistic structure for glucose and enzyme data directly.

Two further tests would sharpen what the agent contributes. Anonymizing scientific variable names measures how much the proposer relies on domain knowledge, and an adversarial proposer probes whether search-time fitness and held-out performance can diverge. Discovered kernels also apply beyond Gaussian processes, for example to kernel ridge regression or to semantic kernels for language model uncertainty. The documentation describes the implementation and workflow.