Open in GitHubRun in Google Colab

Empirical Gaussian Processes

This tutorial illustrates the Empirical Gaussian Process introduced in "Empirical Gaussian Processes" (Lin et al., ICML 2026) on two real one-dimensional forecasting problems: a financial series (the S&P 500 stock-market index) and a climate series (atmospheric CO2 concentration at the Mauna Loa Observatory).

The idea

A standard Gaussian process (GP) requires the user to handcraft a prior mean and covariance (kernel) function that encode beliefs about the function being modeled. The kernel is typically selected from a small set of standard stationary kernels (RBF, Matérn, periodic), and its hyperparameters are fit by maximizing the marginal likelihood. This is limited in flexibility — there is no "one size fits all" kernel — and stationary kernels are poorly suited for extrapolation beyond the observed input range.

The Empirical GP instead estimates the prior mean and covariance functions directly from a corpus of historical observations of the data-generating stochastic process. Given $S$ independent sample paths $f_1, \dots, f_S$, we estimate the true mean $m = \mathbb{E}[f]$ and covariance function $k = \mathrm{Cov}(f)$ by maximum likelihood,

$$m_S(\mathbf{x}) = \frac{1}{S}\sum_{i=1}^{S} f_i(\mathbf{x}), \qquad k_S(\mathbf{x}, \mathbf{x}') = \frac{1}{S}\sum_{i=1}^{S} \tilde{f}_i(\mathbf{x})\, \tilde{f}_i(\mathbf{x}'),$$

where $\tilde{f}_i(\mathbf{x}) = f_i(\mathbf{x}) - m(\mathbf{x})$ is the centered sample path. We define the Empirical GP as $\mathcal{GP}(m_S, k_S)$. Note that $k_S$ is a valid (positive semi-definite) kernel by construction, since it is a sum of outer products of the centered functions. As $S \to \infty$, the Empirical GP converges to the best Gaussian approximation (in the KL-divergence sense) of the true data-generating process.

The resulting EmpiricalOneDimensionalMean and EmpiricalOneDimensionalKernel interpolate these empirical statistics to arbitrary query locations and can be dropped into any BoTorch SingleTaskGP in place of a handcrafted mean/kernel.

What this tutorial does

We reproduce the qualitative experiment from the paper's "Capturing the Behavior of Handcrafted Kernels" section: without human intervention, the Empirical GP recovers the behavior of kernels that were handcrafted by human experts. Each task consists of extrapolating a single time series, so we use sliding windows to extract additional subseries that serve as the historical sample paths $f_1, \dots, f_S$ for the Empirical GP; this implicitly assumes self-similarity in the underlying process. For each series we hold out a future horizon and forecast it from a short observed window, comparing two models:

1. an Empirical GP, whose prior mean and covariance are estimated from the historical subseries, and
2. a handcrafted GP baseline with an expert-designed mean and kernel.

For the financial series the expert model is geometric Brownian motion; for the climate series it is a kernel handcrafted for this dataset by Carl Rasmussen. We then plot both forecasts and report RMSE and negative log-likelihood (NLL) on the held-out horizon.

The integrated white noise kernel (IntegratedWhiteNoiseKernel, parameterized by the integration order: order 1 is Brownian motion, order 2 is integrated Brownian motion, order 3 is twice-integrated Brownian motion, and so on) and the empirical mean/kernel modules used here are all part of the BoTorch library.

# Install dependencies if we are running in colab
import sys

if "google.colab" in sys.modules:
    %pip install botorch

import os
import statistics
import time

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import torch
from botorch.fit import fit_gpytorch_mll
from botorch.models import SingleTaskGP
from botorch.models.empirical_gps import (
    EmpiricalOneDimensionalKernel,
    EmpiricalOneDimensionalMean,
)
from botorch.models.kernels import IntegratedWhiteNoiseKernel
from gpytorch.constraints import GreaterThan
from gpytorch.kernels import PeriodicKernel, RBFKernel, ScaleKernel
from gpytorch.likelihoods import GaussianLikelihood
from gpytorch.means import ConstantMean, LinearMean
from gpytorch.mlls import ExactMarginalLogLikelihood
from gpytorch.priors import LogNormalPrior
from matplotlib.lines import Line2D

torch.set_default_dtype(torch.float64)

SMOKE_TEST = os.environ.get("SMOKE_TEST")

Output:

torch/jit/_script.py:1485: FutureWarning: torch.jit.script is not supported in Python 3.14+ and may break. Please switch to torch.compile or torch.export.
warnings.warn(

Loading the data

We load two CSV files that ship with the BoTorch repository under tutorials/data/:

• sap500.csv: daily S&P 500 index values (we use the Close column).
• co2_daily_mlo.csv: daily atmospheric CO2 concentration at Mauna Loa (NOAA GML); the file has #-prefixed comment lines and no header row, with columns year, month, day, decimal_date, ppm.

Both datasets are redistributed under the CC0 license.

def get_data_dir() -> str:
    """Get the directory of the data, which ships with the BoTorch repo.

    Checks if we are in a common part of the BoTorch repository
    (botorch/botorch or botorch/tutorials) and returns the right path.
    """
    if "DATA_LOCATION" in os.environ:
        return os.environ["DATA_LOCATION"] + "/"
    cwd = os.getcwd()
    folder = os.path.basename(cwd)
    # automated tests run from the botorch folder
    if folder == "botorch":
        return os.path.join(cwd, "tutorials/data/")
    # typical case (running from the tutorial folder)
    elif folder == "tutorials":
        return os.path.join(cwd, "data/")
    # Fallback: check common relative paths (papermill, CI, etc.).
    for candidate in ["tutorials/data/", "pytorch/botorch/tutorials/data/"]:
        if os.path.isdir(candidate):
            return candidate
    # Search upwards for the tutorials data folder.
    cur = os.path.abspath(cwd)
    for _ in range(8):
        for rel in ("tutorials/data", "pytorch/botorch/tutorials/data"):
            candidate = os.path.join(cur, rel)
            if os.path.isdir(candidate):
                return candidate + "/"
        cur = os.path.dirname(cur)
    raise FileNotFoundError("Could not figure out location of the data folder.")


def get_data() -> tuple[pd.DataFrame, pd.DataFrame]:
    """Load the financial (S&P 500) and climate (Mauna Loa CO2) data frames."""
    data_dir = get_data_dir()
    financial_df = pd.read_csv(data_dir + "sap500.csv", header="infer")
    # The CO2 file has commented header lines (#) and no column header row.
    climate_df = pd.read_csv(data_dir + "co2_daily_mlo.csv", comment="#", header=None)
    return financial_df, climate_df


financial_df, climate_df = get_data()
print(financial_df.head())
print(climate_df.head())

Output:

         Date       Open       High        Low      Close  Volume
0  1927-12-30  17.660000  17.660000  17.660000  17.660000     0.0
1  1928-01-03  17.760000  17.760000  17.760000  17.760000     0.0
2  1928-01-04  17.719999  17.719999  17.719999  17.719999     0.0
3  1928-01-05  17.549999  17.549999  17.549999  17.549999     0.0
4  1928-01-06  17.660000  17.660000  17.660000  17.660000     0.0
    0  1   2          3       4
0  1974  5  19  1974.3781  333.46
1  1974  5  20  1974.3808  333.64
2  1974  5  21  1974.3836  333.50
3  1974  5  22  1974.3863  333.21
4  1974  5  23  1974.3890  333.05

# Financial: parse dates and split into a long training history and a test horizon.
financial_df["Date"] = pd.to_datetime(financial_df["Date"])
financial_train_mask = (financial_df["Date"] >= "1930-01-01") & (
    financial_df["Date"] < "2010-01-01"
)
financial_test_mask = (financial_df["Date"] >= "2010-01-01") & (
    financial_df["Date"] < "2025-01-01"
)
financial_data_train = financial_df[financial_train_mask]["Close"].to_numpy()
financial_data_test = financial_df[financial_test_mask]["Close"].to_numpy()

# Climate: build a date column from the year/month/day columns, then split.
climate_df["Date"] = pd.to_datetime(
    climate_df[[0, 1, 2]].astype(str).agg("-".join, axis=1)
)
climate_data_train = climate_df[
    (climate_df["Date"] >= "1975-01-01") & (climate_df["Date"] < "2010-01-01")
][4].to_numpy()
climate_data_test = climate_df[
    (climate_df["Date"] >= "2010-01-01") & (climate_df["Date"] < "2025-01-01")
][4].to_numpy()

print("financial train/test sizes:", financial_data_train.shape, financial_data_test.shape)
print("climate train/test sizes:", climate_data_train.shape, climate_data_test.shape)

Output:

financial train/test sizes: (20093,) (3774,)
climate train/test sizes: (10478,) (4856,)

Helper functions

Because each task consists of extrapolating a single time series, we use a sliding window to extract additional subseries that serve as the historical sample paths $f_1, \dots, f_S$ for the Empirical GP. This assumes self-similarity in the underlying process (which holds exactly for, e.g., geometric Brownian motion). get_historical slices a training series into overlapping fixed-length subseries; we optionally apply a log transform and/or align each subseries to start at zero (an offset), which reflects the exponential nature of the financial series and makes the subseries directly comparable.

def get_historical(
    data: np.ndarray,
    chunk_size: int = 365,
    stride: int = 1,
    use_log_space: bool = False,
    use_offset: bool = False,
) -> np.ndarray:
    """Slice ``data`` into overlapping, sliding-window subseries of length ``chunk_size``.

    Returns an array of shape ``num_chunks x chunk_size``.
    """
    start_indices = np.arange(0, len(data) - chunk_size + 1, stride)
    # Matrix of indices for easy broadcasting: num_chunks x chunk_size
    indices_matrix = start_indices[:, None] + np.arange(chunk_size)
    chunks = data[indices_matrix]
    if use_log_space:
        chunks = np.log(chunks)
    if use_offset:
        chunks = chunks - chunks[:, 0][:, None]
    return chunks


def calculate_annual_metrics(prices: np.ndarray) -> tuple[float, float]:
    """Estimate the annualized return and volatility of a price series."""
    returns = np.diff(prices) / prices[:-1]
    mean_daily_return = np.mean(returns)
    annual_return = (1 + mean_daily_return) ** 252 - 1
    daily_volatility = np.std(returns)
    annualized_volatility = daily_volatility * np.sqrt(252)
    return annual_return, annualized_volatility

The get_predictions helper builds a BoTorch SingleTaskGP from a given mean and covariance module and returns posterior predictions on the training and test inputs. It supports an optional log-space transform with an additive offset (so the GP models the de-trended log series and we map predictions back to the original space), a fixed tiny observation noise (used for the Empirical GP, which already captures uncertainty through its empirical covariance), and marginal-likelihood fitting (used for the handcrafted baselines).

def _revert_log_normal(mean: torch.Tensor, var: torch.Tensor):
    """Map the mean/variance of a log-normal back to the original space."""
    reverted_mean = torch.exp(mean + 0.5 * var)
    reverted_var = (var.exp() - 1) * (2 * mean + var).exp()
    return reverted_mean, reverted_var


def _get_likelihood(loc=-4.0, scale=1.0, lower_bound=1e-8, fix_noise=None):
    prior = LogNormalPrior(loc=loc, scale=scale)
    constraint = GreaterThan(lower_bound, initial_value=prior.mode)
    likelihood = GaussianLikelihood(noise_prior=prior, noise_constraint=constraint)
    if fix_noise is not None:
        likelihood.noise = torch.tensor(fix_noise, dtype=torch.float64)
        likelihood.requires_grad_(False)
    return likelihood


def _predict(X, model, use_log_space, offset, observation_noise, n_samples):
    model.eval()
    with torch.no_grad():
        posterior = model.posterior(X, observation_noise=observation_noise)
        pred_mean = posterior.mean.squeeze() + offset
        pred_var = posterior.variance.squeeze()
        if use_log_space:
            pred_mean, pred_var = _revert_log_normal(pred_mean, pred_var)
        pred_std = pred_var.sqrt()

        samples = posterior.rsample(sample_shape=torch.Size([n_samples]))
        samples = samples.squeeze() + offset
        if use_log_space:
            samples = samples.exp()
    return pred_mean, pred_std, samples


def get_predictions(
    train_X,
    train_Y,
    test_X,
    mean_module,
    covar_module,
    fix_noise=None,
    outcome_transform=None,
    use_log_space=False,
    use_offset=False,
    fit_model=False,
    observation_noise=False,
    n_samples=5,
    seed=42,
):
    torch.manual_seed(seed)

    if use_log_space:
        train_Y = train_Y.log()
    offset = train_Y[0] if use_offset else 0.0

    model = SingleTaskGP(
        train_X=train_X[:, None],
        train_Y=train_Y[:, None] - offset,
        likelihood=_get_likelihood(fix_noise=fix_noise),
        mean_module=mean_module,
        covar_module=covar_module,
        input_transform=None,
        outcome_transform=outcome_transform,
    )

    if fit_model:
        mll = ExactMarginalLogLikelihood(model.likelihood, model)
        fit_kwargs = {}
        if SMOKE_TEST is not None:
            # Keep the smoke test fast by capping the optimizer iterations.
            fit_kwargs["optimizer_kwargs"] = {"options": {"maxiter": 5}}
        fit_gpytorch_mll(mll, **fit_kwargs)

    pred_mean_train, pred_std_train, samples_train = _predict(
        train_X, model, use_log_space, offset, observation_noise, n_samples
    )
    pred_mean_test, pred_std_test, samples_test = _predict(
        test_X, model, use_log_space, offset, observation_noise, n_samples
    )
    return (
        pred_mean_train,
        pred_std_train,
        samples_train,
        pred_mean_test,
        pred_std_test,
        samples_test,
    )

Financial series: S&P 500

We forecast the second half of a 252-trading-day (one-year) window from its first half. The Empirical GP estimates its prior from overlapping one-year subseries of the 1930–2010 training history, log-transformed and aligned to start at zero. The handcrafted baseline is geometric Brownian motion, the canonical expert model for stock prices, which assumes the log price follows a random walk with drift, $dS_t = \mu S_t\, dt + \sigma S_t\, dW_t$, where $S_t$ is the price at time $t$, $\mu$ is the expected rate of return, $\sigma$ is the volatility, and $W_t$ is Brownian motion. We estimate $\mu$ and $\sigma$ on the same historical data used by the Empirical GP. In de-trended log-space this corresponds to a linear drift mean and a Brownian-motion covariance — a once-integrated white noise kernel, $k(s, t) = \min(s, t)$ — with output scale set from the historical volatility and observation noise fit by marginal likelihood.

chunk_size_f = 252
historical_X_f = torch.linspace(0, 1, chunk_size_f, dtype=torch.float64)
historical_Y_f = get_historical(
    financial_data_train,
    chunk_size=chunk_size_f,
    stride=1,
    use_log_space=True,
    use_offset=True,
)
historical_Y_f = torch.from_numpy(historical_Y_f)[..., None]

observed_size_f = 126
# NB: train_X_f / train_Y_f are the *observed window* of the forecast target
# (the first half of financial_data_test), NOT the historical training series --
# the history enters only through the empirical prior (historical_Y_f above).
train_X_f = historical_X_f[:observed_size_f]
train_Y_f = torch.from_numpy(financial_data_test[:observed_size_f])
test_X_f = historical_X_f[observed_size_f:]
test_Y_f = torch.from_numpy(financial_data_test[observed_size_f:chunk_size_f])

# Empirical GP: prior mean/covariance estimated from the historical subseries.
empirical_mean_f = EmpiricalOneDimensionalMean(
    X_full=historical_X_f[:, None], Y_full=historical_Y_f
)
empirical_covar_f = EmpiricalOneDimensionalKernel(
    X_full=historical_X_f[:, None], Y_full=historical_Y_f, ard=False
)
empirical_preds_f = get_predictions(
    train_X_f,
    train_Y_f,
    test_X_f,
    empirical_mean_f,
    empirical_covar_f,
    fix_noise=1e-8,
    use_log_space=True,
    use_offset=True,
    fit_model=False,
    observation_noise=False,
    seed=123,
)

# Handcrafted GP: geometric Brownian motion (linear drift mean + once-integrated
# white noise covariance, k(s, t) = min(s, t)).
expected_return, volatility = calculate_annual_metrics(financial_data_train)
log_return = expected_return - volatility**2 / 2

handcrafted_mean_f = LinearMean(input_size=1, bias=False)
handcrafted_mean_f.weights.data = torch.tensor([log_return], dtype=torch.float64)
handcrafted_mean_f.requires_grad_(False)
handcrafted_covar_f = ScaleKernel(IntegratedWhiteNoiseKernel())
handcrafted_covar_f.outputscale = volatility**2
handcrafted_covar_f.requires_grad_(False)

handcrafted_preds_f = get_predictions(
    train_X_f,
    train_Y_f,
    test_X_f,
    handcrafted_mean_f,
    handcrafted_covar_f,
    use_log_space=True,
    use_offset=True,
    fit_model=True,
    observation_noise=False,
    seed=123,
)

Output:

gpytorch/distributions/multivariate_normal.py:375: NumericalWarning: Negative variance values detected. This is likely due to numerical instabilities. Rounding negative variances up to 1e-10.
warnings.warn(
linear_operator/utils/cholesky.py:41: NumericalWarning: A not p.d., added jitter of 1.0e-08 to the diagonal
warnings.warn(
gpytorch/distributions/multivariate_normal.py:375: NumericalWarning: Negative variance values detected. This is likely due to numerical instabilities. Rounding negative variances up to 1e-10.
warnings.warn(
linear_operator/utils/cholesky.py:41: NumericalWarning: A not p.d., added jitter of 1.0e-08 to the diagonal
warnings.warn(

Climate series: Mauna Loa CO2

We forecast a multi-year horizon of daily CO2 concentration. Because the daily records are densely sampled, we use the interpolation-based Empirical GP: each ~8-year subseries of the 1975–2010 training history is linearly interpolated onto a common daily grid and treated as a fully observed sample path, after which the empirical mean and covariance can be evaluated at arbitrary locations. The handcrafted baseline is a kernel handcrafted for this dataset by a human expert (Rasmussen). It consists of three additive components that respectively model trend, seasonality, and noise: the trend is a sum of once-, twice-, and thrice-integrated white noise (IntegratedWhiteNoiseKernel(order=1, 2, 3)); the seasonal component is a product of a periodic kernel with a Gaussian RBF kernel (PeriodicKernel * RBFKernel); and the noise component is a Gaussian RBF kernel (RBFKernel) with additive homoskedastic observation noise. The mean is constant and all hyperparameters are fit by marginal likelihood.

# Normalized-axis tensors used for plotting the observed window and ground truth.
# chunk_size_c must equal the number of daily CO2 records in the test date range
# (2010-2018); assert it so a data refresh cannot silently misalign the forecast
# grid and the ground truth.
chunk_size_c = 2608
observed_size_c = 1304
assert len(climate_data_test) >= chunk_size_c, (
    f"expected >= {chunk_size_c} climate test records, got {len(climate_data_test)}"
)
historical_X_c = torch.linspace(0, 1, chunk_size_c, dtype=torch.float64)
train_X_c = historical_X_c[:observed_size_c]
train_Y_c = torch.from_numpy(climate_data_test[:observed_size_c])
test_X_c = historical_X_c[observed_size_c:]
test_Y_c = torch.from_numpy(climate_data_test[observed_size_c:chunk_size_c])

# Handcrafted GP: expert-designed composite kernel (trend + seasonal + noise)
# with a constant mean.
handcrafted_mean_c = ConstantMean()
# Pre-industrial CO2 baseline (~280 ppm); the rise to present-day levels is carried
# by the integrated-white-noise trend kernels over "years since 1750".
handcrafted_mean_c.constant = 280
handcrafted_mean_c.requires_grad_(False)

# Trend: sum of once-, twice-, and thrice-integrated white noise.
k_trend = (
    ScaleKernel(IntegratedWhiteNoiseKernel(order=1))
    + ScaleKernel(IntegratedWhiteNoiseKernel(order=2))
    + ScaleKernel(IntegratedWhiteNoiseKernel(order=3))
)
# Seasonal: a periodic kernel modulated by a Gaussian RBF kernel.
k_season = ScaleKernel(PeriodicKernel() * RBFKernel())
# Noise: a Gaussian RBF kernel; homoskedastic observation noise is added
# through the Gaussian likelihood.
k_noise = ScaleKernel(RBFKernel())
handcrafted_covar_c = k_trend + k_season + k_noise


def _to_fractional_years(dates: np.ndarray) -> np.ndarray:
    """Convert datetimes to fractional years measured from year 1750."""
    dts = pd.to_datetime(dates)
    years = dts.year + (dts.dayofyear - 1) / 365.25
    return years.astype(float).to_numpy() - 1750


train_X_dates = climate_df[
    (climate_df["Date"] >= "1975-01-01") & (climate_df["Date"] < "2010-01-01")
]["Date"].to_numpy()
test_X_dates = climate_df[
    (climate_df["Date"] >= "2010-01-01") & (climate_df["Date"] <= "2018-01-01")
]["Date"].to_numpy()

# The handcrafted model observes the full training history plus the observed
# portion of the test horizon, and forecasts the remainder.
train_X_dates = np.concatenate([train_X_dates, test_X_dates[:observed_size_c]])
train_X_hc = _to_fractional_years(train_X_dates)
test_X_hc = _to_fractional_years(test_X_dates)[observed_size_c:]
train_Y_hc = np.concatenate([climate_data_train, climate_data_test[:observed_size_c]])

# Note: the handcrafted GP conditions on the full 1975-2010 history plus the observed
# window (subsampled 2x below for tractability), whereas the empirical GP conditions
# only on the observed window (its history enters through the prior). This information
# asymmetry favors the handcrafted model, so the comparison is conservative for the
# empirical GP.
# Subsample the (large) training set to keep fitting tractable -- aggressively under
# SMOKE_TEST so the notebook runs in seconds.
hc_stride = 40 if SMOKE_TEST is not None else 2
train_X_hc = torch.from_numpy(train_X_hc[::hc_stride])
train_Y_hc = torch.from_numpy(train_Y_hc[::hc_stride])
test_X_hc = torch.from_numpy(test_X_hc)

handcrafted_preds_c = get_predictions(
    train_X_hc,
    train_Y_hc,
    test_X_hc,
    handcrafted_mean_c,
    handcrafted_covar_c,
    use_log_space=False,
    use_offset=False,
    fit_model=True,
    observation_noise=True,
    seed=123,
)

# Empirical GP for the climate series (interpolation-based, for densely sampled
# data). We build the historical sample paths by linearly interpolating each
# ~8-year subseries onto a common daily grid, then treat them as fully observed
# realizations of the process.
window_size_c = 3000

X_days = climate_df[
    (climate_df["Date"] >= "1975-01-01") & (climate_df["Date"] < "2010-01-01")
]["Date"].to_numpy()
X_days = (pd.to_datetime(X_days) - pd.to_datetime("1975-01-01")).days.to_numpy()


def get_indices(arr: np.ndarray, window_size: int, stride: int = 1):
    """Indices of all points within ``window_size`` of each strided start."""
    indices_list = []
    for i in range(0, len(arr) - window_size, stride):
        indices = [j for j in range(i, len(arr)) if arr[j] <= arr[i] + window_size]
        indices_list.append(indices)
    return indices_list


stride_c = 100 if SMOKE_TEST is not None else 1
idx_list = get_indices(X_days, window_size_c, stride=stride_c)

grid = np.arange(window_size_c)
historical_Y_list = []
for idx_ in idx_list:
    X_ = X_days[idx_] - X_days[idx_][0]
    Y_ = np.interp(grid, X_, climate_data_train[idx_])
    historical_Y_list.append(Y_)
historical_Y_c = torch.from_numpy(np.stack(historical_Y_list))[..., None]
historical_X_c_grid = torch.from_numpy(grid).to(torch.float64)[..., None]

# Test-time inputs measured in days from the start of the test horizon.
X_test_days = climate_df[
    (climate_df["Date"] >= "2010-01-01") & (climate_df["Date"] <= "2018-01-01")
]["Date"].to_numpy()
X_test_days = (pd.to_datetime(X_test_days) - pd.to_datetime("1975-01-01")).days.to_numpy()
X_test_days = X_test_days - X_test_days[0]

train_X_emp = torch.from_numpy(X_test_days[:observed_size_c]).to(torch.float64)
train_Y_emp = torch.from_numpy(climate_data_test[:observed_size_c])
test_X_emp = torch.from_numpy(X_test_days[observed_size_c:]).to(torch.float64)

empirical_mean_c = EmpiricalOneDimensionalMean(
    X_full=historical_X_c_grid, Y_full=historical_Y_c
)
empirical_covar_c = EmpiricalOneDimensionalKernel(
    X_full=historical_X_c_grid, Y_full=historical_Y_c, ard=False
)
empirical_preds_c = get_predictions(
    train_X_emp,
    train_Y_emp,
    test_X_emp,
    empirical_mean_c,
    empirical_covar_c,
    fix_noise=1e-8,
    use_log_space=False,
    use_offset=False,
    fit_model=False,
    observation_noise=False,
    seed=123,
)

Results

The figure below shows, for each series, the observed window (black), the held-out ground truth (orange), the Empirical GP forecast with its 95% credible interval (blue), and the handcrafted GP forecast with its 95% credible interval (dashed black). For the financial series the models work in log-space, so the plotted price-space band is a Gaussian moment-matched approximation (mean $\pm 2$ std of the reverted log-normal), not the exact -- and asymmetric -- log-normal quantile interval.

Without any human intervention, the Empirical GP recovers the behavior of the expert-handcrafted kernels. On the stock-market data its forecast closely tracks geometric Brownian motion on this window — qualitative agreement consistent with the Empirical GP converging to the best Gaussian approximation of the data-generating process, though a single forecast window is illustrative rather than a validation of GBM as that approximation (here the handcrafted GBM is in fact marginally ahead on RMSE, with NLL roughly tied). On the atmospheric data it implicitly captures the seasonality and upward trend — without any explicit inductive bias — and attains a competitive (here slightly lower) RMSE than the expert-designed kernel. Note that the two models are not matched on compute: the handcrafted GP is fit on a factor-two subsample of the history for tractability, while the Empirical GP uses the dense interpolated histories, so this comparison is indicative rather than a controlled head-to-head.

def plot_empirical(ax, X, preds, lw=2, color="tab:blue"):
    _, _, _, mean_test, std_test, _ = preds
    mean_test = np.asarray(mean_test)
    std_test = np.asarray(std_test)
    X = np.asarray(X).squeeze()
    ax.plot(X, mean_test, linewidth=lw, color=color)
    ax.fill_between(
        X,
        mean_test - 2 * std_test,
        mean_test + 2 * std_test,
        color=color,
        alpha=0.3,
    )


def plot_handcrafted(ax, X, preds, lw=0.7, color="k"):
    _, _, _, mean_test, std_test, _ = preds
    mean_test = np.asarray(mean_test)
    std_test = np.asarray(std_test)
    X = np.asarray(X).squeeze()
    ax.plot(X, mean_test, linewidth=lw, linestyle="--", color=color)
    ax.plot(X, mean_test - 2 * std_test, linewidth=lw, linestyle="--", color=color)
    ax.plot(X, mean_test + 2 * std_test, linewidth=lw, linestyle="--", color=color)


fig = plt.figure(figsize=(13, 4), dpi=120)
ax1 = fig.add_subplot(1, 2, 1)
ax2 = fig.add_subplot(1, 2, 2)

for ax in [ax1, ax2]:
    ax.spines[["top", "right"]].set_visible(False)
    ax.grid(alpha=0.3)

lw = 1
# Financial panel
ax1.plot(np.asarray(train_X_f), np.asarray(train_Y_f), linewidth=lw, color="k")
ax1.plot(np.asarray(test_X_f), np.asarray(test_Y_f), linewidth=lw, color="tab:orange")
plot_empirical(ax1, test_X_f, empirical_preds_f)
plot_handcrafted(ax1, test_X_f, handcrafted_preds_f)
ax1.set_xlim([0, 1])
ax1.set_xlabel("Time (fraction of one year from Jan 2010)")
ax1.set_ylabel("Index Points")
ax1.set_title("S&P 500 Index (financial series)", fontsize=11)

# Climate panel
ax2.plot(np.asarray(train_X_c), np.asarray(train_Y_c), linewidth=lw, color="k")
ax2.plot(np.asarray(test_X_c), np.asarray(test_Y_c), linewidth=lw, color="tab:orange")
plot_empirical(ax2, test_X_c, empirical_preds_c)
plot_handcrafted(ax2, test_X_c, handcrafted_preds_c)
ax2.set_xlim([0, 1])
ax2.set_xticks(np.linspace(0, 1, 9))
ax2.set_xticklabels([2010 + i for i in range(9)])
ax2.set_ylim([384, 414])
ax2.set_xlabel("Time")
ax2.set_ylabel("Mole Fraction (ppm)")
ax2.set_title("Mauna Loa CO$_2$ (climate series)", fontsize=11)

custom_handles = [
    Line2D([0], [0], color="k", label="Observed Data"),
    Line2D([0], [0], color="tab:orange", label="Ground Truth"),
    Line2D([0], [0], color="tab:blue", label="Empirical GP Mean and 95% CI", lw=2),
    Line2D([0], [0], linestyle="--", color="k", label="Handcrafted GP Mean and 95% CI"),
]
fig.legend(
    handles=custom_handles,
    loc="lower center",
    ncol=4,
    bbox_to_anchor=(0.5, -0.05),
    frameon=False,
)
fig.tight_layout()
plt.show()

Output:

findfont: Font family ['STIXGeneral'] not found. Falling back to DejaVu Sans.
findfont: Font family ['STIXGeneral'] not found. Falling back to DejaVu Sans.
findfont: Font family ['STIXGeneral'] not found. Falling back to DejaVu Sans.
findfont: Font family ['STIXGeneral'] not found. Falling back to DejaVu Sans.
findfont: Font family ['STIXNonUnicode'] not found. Falling back to DejaVu Sans.
findfont: Font family ['STIXNonUnicode'] not found. Falling back to DejaVu Sans.
findfont: Font family ['STIXNonUnicode'] not found. Falling back to DejaVu Sans.
findfont: Font family ['STIXSizeOneSym'] not found. Falling back to DejaVu Sans.
findfont: Font family ['STIXSizeTwoSym'] not found. Falling back to DejaVu Sans.
findfont: Font family ['STIXSizeThreeSym'] not found. Falling back to DejaVu Sans.
findfont: Font family ['STIXSizeFourSym'] not found. Falling back to DejaVu Sans.
findfont: Font family ['STIXSizeFiveSym'] not found. Falling back to DejaVu Sans.
findfont: Font family ['cmsy10'] not found. Falling back to DejaVu Sans.
findfont: Font family ['cmr10'] not found. Falling back to DejaVu Sans.
findfont: Font family ['cmtt10'] not found. Falling back to DejaVu Sans.
findfont: Font family ['cmmi10'] not found. Falling back to DejaVu Sans.
findfont: Font family ['cmb10'] not found. Falling back to DejaVu Sans.
findfont: Font family ['cmss10'] not found. Falling back to DejaVu Sans.
findfont: Font family ['cmex10'] not found. Falling back to DejaVu Sans.
findfont: Font family ['DejaVu Sans Mono'] not found. Falling back to DejaVu Sans.
findfont: Font family ['DejaVu Sans Display'] not found. Falling back to DejaVu Sans.

Finally, we quantify forecast quality on the held-out horizon with the root-mean-squared error (RMSE) and the negative log-likelihood (NLL) under each model's predictive distribution. (For the financial series the NLL is computed on the moment-matched Gaussian in price space rather than the exact log-normal density; the ranking is unaffected since both financial models share the same log-space transform.)

def rmse_nll(mean, ground_truth, *, var=None, std=None):
    """RMSE and Gaussian negative log-likelihood (per-point mean).

    Shared metric helper used throughout the tutorial. Accepts either torch
    tensors or numpy arrays. Supply the predictive uncertainty as a variance
    (``var``) or a standard deviation (``std``); if neither is given, the NLL is
    returned as ``nan`` (useful when only the RMSE is needed).
    """

    def _np(a):
        return a.detach().cpu().numpy() if isinstance(a, torch.Tensor) else np.asarray(a)

    mean = _np(mean)
    gt = _np(ground_truth)
    rmse = float(np.sqrt(np.mean((mean - gt) ** 2)))
    if std is not None:
        var = _np(std) ** 2
    if var is None:
        return rmse, float("nan")
    var = np.clip(_np(var), 1e-9, None)
    nll = float(np.mean(0.5 * np.log(2 * np.pi * var) + 0.5 * (gt - mean) ** 2 / var))
    return rmse, nll


# Financial data
rmse_empirical_f, nll_empirical_f = rmse_nll(
    empirical_preds_f[3], test_Y_f, std=empirical_preds_f[4]
)
rmse_handcrafted_f, nll_handcrafted_f = rmse_nll(
    handcrafted_preds_f[3], test_Y_f, std=handcrafted_preds_f[4]
)

# Climate data
rmse_empirical_c, nll_empirical_c = rmse_nll(
    empirical_preds_c[3], test_Y_c, std=empirical_preds_c[4]
)
rmse_handcrafted_c, nll_handcrafted_c = rmse_nll(
    handcrafted_preds_c[3], test_Y_c, std=handcrafted_preds_c[4]
)

print(
    f"Financial  | Empirical GP:   RMSE = {rmse_empirical_f:8.4f}, "
    f"NLL = {nll_empirical_f:8.4f}"
)
print(
    f"Financial  | Handcrafted GP: RMSE = {rmse_handcrafted_f:8.4f}, "
    f"NLL = {nll_handcrafted_f:8.4f}"
)
print(
    f"Climate    | Empirical GP:   RMSE = {rmse_empirical_c:8.4f}, "
    f"NLL = {nll_empirical_c:8.4f}"
)
print(
    f"Climate    | Handcrafted GP: RMSE = {rmse_handcrafted_c:8.4f}, "
    f"NLL = {nll_handcrafted_c:8.4f}"
)

Output:

Financial  | Empirical GP:   RMSE = 122.6802, NLL =   6.1529
Financial  | Handcrafted GP: RMSE = 121.9282, NLL =   6.1540
Climate    | Empirical GP:   RMSE =   1.1168, NLL =   1.5485
Climate    | Handcrafted GP: RMSE =   1.4126, NLL =   1.6769

Part 2 — Learning Curve Extrapolation with Incomplete Data via EM

So far we've assumed every historical curve is fully observed on a common grid. In AutoML with early stopping — the motivating application for LCBench — most historical HPO runs are truncated at different budgets.

The EM-based Empirical GP extends the empirical GP to handle incomplete observations. EM alternates between imputing missing values (E-step) and updating the prior mean $\mu$ and covariance $\Sigma$ (M-step), learning an empirical prior from incomplete historical curves by maximum likelihood — no handcrafted kernel, no ad-hoc imputation.

In this part we use LCBench, a canonical AutoML benchmark with 2,000 hyperparameter configurations evaluated over 52 epochs on 35 OpenML datasets. Each learning curve maps epoch → validation accuracy, providing structured covariance ideal for showcasing EM.

Note on SMOKE_TEST mode. When the SMOKE_TEST environment variable is set, every experiment below runs on drastically reduced datasets, sample counts, EM iterations, and optimizer budgets so the notebook executes quickly as a CI check. In that mode the printed numbers and figures are execution checks only, not scientific results; the quantitative statements in the text refer to the full (non-smoke) configuration.

A. Incomplete learning curves on LCBench

LCBench provides learning curves for hyperparameter optimization. We load LCBench via botorch.utils.lcbench, truncate 80% of historical curves to simulate early stopping, and visualize the incomplete corpus.

# Part 2 imports -- EM extensions.
# LCBench is loaded via `botorch.utils.lcbench`, which downloads the LCBench
# Parquet files from GitHub and caches them locally under ~/.cache.
from botorch.utils.lcbench import load_lcbench_data
from botorch.models.empirical_gps.utils import ExperimentDataset, em_prior_to_basis_curves
from botorch.models.empirical_gps import (
    EMEmpiricalGaussianProcess,
    pretrain_em_prior,
    EmpiricalOneDimensionalGP,
)
from gpytorch.kernels import MaternKernel, ScaleKernel
from gpytorch.means import ConstantMean

# SMOKE_TEST uses tiny subsets for speed
IS_SMOKE = SMOKE_TEST is not None
DATASET_NAME = "Fashion-MNIST"
METRIC = "Train/val_accuracy"


def load_lcbench(name, metric=METRIC, dtype=torch.double):
    """Load an LCBench (dataset, metric): returns (metrics, parameters, epoch_grid)."""
    data = load_lcbench_data(name, metric, dtype=dtype)
    return data.metrics, data.parameters, data.epochs


def seeded_split(n, sizes, seed):
    """Seeded disjoint index splits of ``range(n)`` (a list of LongTensors)."""
    torch.manual_seed(seed)
    perm = torch.randperm(n)
    out, off = [], 0
    for s in sizes:
        out.append(perm[off:off + s]); off += s
    return out


def style_ax(ax, grid=True):
    """Shared plot style: hide top/right spines and add a light grid."""
    ax.spines[["top", "right"]].set_visible(False)
    if grid:
        ax.grid(alpha=0.3)
    return ax


print(f"Loading LCBench dataset '{DATASET_NAME}'...")
metrics, params, X_grid = load_lcbench(DATASET_NAME, METRIC)
print(f"LCBench metrics shape: {metrics.shape}, parameters shape: {params.shape}")

n_historical = 20 if IS_SMOKE else 100
n_test = 5
hist_idx, test_idx = seeded_split(metrics.shape[0], [n_historical, n_test], seed=42)
historical_full = metrics[hist_idx]
test_full = metrics[test_idx]
print(f"Historical curves: {historical_full.shape}, test curves: {test_full.shape}")

Output:

Loading LCBench dataset 'Fashion-MNIST' via Ax...
LCBench metrics shape: torch.Size([2000, 50]), parameters shape: torch.Size([2000, 7])
Historical curves: torch.Size([100, 50]), test curves: torch.Size([5, 50])

# Simulate early stopping: truncate 80% of historical curves at varying budgets.
# Most curves stop early (20-80% of epochs), so the late-epoch region is only sparsely
# covered. Rather than engineer the cutoff distribution, we let the EM prior's covariance
# shrinkage (next cell) regularize that data-starved tail in a principled way.
torch.manual_seed(123)
n_fully = int(0.2 * n_historical)
partial_range = (0.2, 0.8)

historical_datasets = []
for i in range(n_historical):
    y_full = historical_full[i]
    if i < n_fully:
        x_obs = X_grid
        y_obs = y_full
    else:
        obs_frac = torch.empty(1).uniform_(*partial_range).item()
        cutoff = max(5, int(obs_frac * len(X_grid)))
        x_obs = X_grid[:cutoff]
        y_obs = y_full[:cutoff]
    historical_datasets.append(ExperimentDataset(X=x_obs.unsqueeze(-1), Y=y_obs.unsqueeze(-1)))

print(f"{n_fully} fully observed, {n_historical - n_fully} partially observed")

fig, ax = plt.subplots(1, 1, figsize=(10, 4))
style_ax(ax, grid=False)
for ds in historical_datasets[:30]:
    x = ds.X.squeeze().numpy()
    y = ds.Y.squeeze().numpy()
    color = 'tab:blue' if len(x) == len(X_grid) else 'tab:gray'
    alpha = 0.8 if len(x) == len(X_grid) else 0.4
    ax.plot(x, y, color=color, alpha=alpha, linewidth=1)
ax.set_xlabel('Epoch')
ax.set_ylabel('Validation Accuracy')
ax.set_title('LCBench historical learning curves with simulated early stopping\n(blue=fully observed, gray=truncated)')
ax.grid(alpha=0.3)
plt.show()
print('EmpiricalOneDimensionalGP requires complete Y_full — EM-EGP handles the ragged corpus directly.')

Output:

20 fully observed, 80 partially observed
EmpiricalOneDimensionalGP requires complete Y_full — EM-EGP handles the ragged corpus directly.

B. EM learns the prior from incomplete curves

EM alternates E-step (compute conditional $p(\theta \mid y_i)$ for each partial curve — posterior mean is the completed curve) and M-step (update prior $\mu,\Sigma$ by averaging conditionals). Repeating yields ML empirical prior from incomplete data.

# Pre-train EM prior on incomplete historical datasets
n_em = 5 if IS_SMOKE else 30
print(f"Pre-training EM prior with {n_em} iterations on {len(historical_datasets)} incomplete curves...")
start = time.time()

# The base Matern kernel is the EM covariance's smooth, full-rank fall-back: it seeds
# the E-step imputation of the truncated tails and is the M-step shrinkage target. Its
# lengthscale must be commensurate with the epoch axis (here 1..52). The GPyTorch
# default (~0.69) is a correlation length of well under one epoch, so conditioning on
# the observed head barely constrains the tail and the predictive standard deviation
# jumps abruptly to the prior envelope just past the last observation. A correlation
# length of ~10 epochs (about a fifth of the range) makes the posterior uncertainty
# grow smoothly into the unobserved region.
base_covar = ScaleKernel(MaternKernel(nu=2.5, ard_num_dims=1))
base_covar.base_kernel.lengthscale = 10.0

em_prior = pretrain_em_prior(
    datasets=historical_datasets,
    mean_module=ConstantMean(),
    covar_module=base_covar,
    num_em_iterations=n_em,
    enable_interpolation=True,
    # Shrink the EM covariance toward the base kernel K(Z, Z) at each M-step. This is a
    # flexible generalization of the Inverse-Wishart prior (a free intensity in [0, 1],
    # decoupled from the number of inducing points) that regularizes the sparsely-observed
    # late-epoch tail, so the completed curves track the true trajectory instead of
    # collapsing where few historical curves reach.
    covariance_shrinkage=0.3,
)
print(f"EM done in {time.time()-start:.1f}s")
print(f"Learned prior mu shape {em_prior.mu_inducing.shape}, Sigma shape {em_prior.Sigma_inducing.shape}")

Output:

Pre-training EM prior with 30 iterations on 100 incomplete curves...
EM done in 1.8s
Learned prior mu shape torch.Size([50]), Sigma shape torch.Size([50, 50])

E-step as curve completion

Pick one test curve observed to 30% of epochs. EM conditional mean fills the missing tail with uncertainty widening where unobserved.

test_curve_full = test_full[0]
obs_frac = 0.3
cutoff = int(obs_frac * len(X_grid))
X_obs = X_grid[:cutoff].unsqueeze(-1)
Y_obs = test_curve_full[:cutoff].unsqueeze(-1)
X_test = X_grid.unsqueeze(-1)

model_em = EMEmpiricalGaussianProcess.from_pretrained(em_prior, train_X=X_obs, train_Y=Y_obs)
# Recorded curves are essentially noise-free, so use a small observation noise for a
# clean interpolating fit at the observed epochs (the tail behavior is governed by the
# EM covariance, not the noise).
model_em.likelihood.noise = torch.tensor(1e-4, dtype=torch.double)
model_em.eval()
with torch.no_grad():
    post = model_em.posterior(X_test)
    mu_post = post.mean.squeeze().numpy()
    std_post = post.variance.sqrt().squeeze().numpy()

fig, ax = plt.subplots(1, 1, figsize=(10, 4))
style_ax(ax, grid=False)
ax.plot(X_grid.numpy(), test_curve_full.numpy(), color='tab:orange', linewidth=2, linestyle='--', label='True (held-out)')
ax.plot(X_obs.squeeze().numpy(), Y_obs.squeeze().numpy(), 'ko', markersize=4, label='Observed (30%)')
ax.plot(X_grid.numpy(), mu_post, color='tab:blue', linewidth=2, label='EM conditional mean')
ax.fill_between(X_grid.numpy(), mu_post-2*std_post, mu_post+2*std_post, color='tab:blue', alpha=0.25, label='95% CI')
ax.axvline(X_grid[cutoff-1].item(), color='k', linestyle=':', alpha=0.5, label='Cutoff')
ax.set_xlabel('Epoch'); ax.set_ylabel('Validation Accuracy')
ax.set_title('EM curve completion: E-step imputes missing tail with calibrated uncertainty')
ax.legend(frameon=False); ax.grid(alpha=0.3); plt.show()

C. From EM prior to a fast empirical GP via basis curves

The full EM-EGP re-runs EM at every forward pass. For fast inference we synthesize a small set of complete "basis" curves that reproduce the EM prior mean and covariance exactly on the grid via em_prior_to_basis_curves, then feed them to the cheaper EmpiricalOneDimensionalGP. The surrogate reproduces the EM-EGP posterior to numerical precision and, because it avoids re-running EM on every forward pass, is cheaper to query for repeated inference. The speed advantage grows with the number of posterior evaluations and the grid/batch size; on the tiny 52-point grid used here the two are close, so the timing below is a sanity check rather than a headline speedup.

Y_basis = em_prior_to_basis_curves(mu=em_prior.mu_inducing, Sigma=em_prior.Sigma_inducing, num_modes=None)
print(f"Synthesized Y_full shape: {Y_basis.shape}  (num_curves = rank + 1)")
X_hist = em_prior.X_inducing

model_empirical = EmpiricalOneDimensionalGP(train_X=X_obs, train_Y=Y_obs, historical_X=X_hist, historical_Y=Y_basis)
# The basis curves reproduce the EM prior mean/covariance exactly; match the
# observation noise so the fast surrogate also reproduces the EM *posterior*.
model_empirical.likelihood.noise = model_em.likelihood.noise.detach().clone()
model_empirical.eval()

def time_posterior(model, X, n_reps=20, n_warmup=3):
    model.eval()
    for _ in range(n_warmup):  # warm up caches before timing
        with torch.no_grad(): _ = model.posterior(X)
    times = []
    for _ in range(n_reps):
        t0 = time.perf_counter()
        with torch.no_grad(): _ = model.posterior(X)
        times.append(time.perf_counter() - t0)
    return statistics.median(times), statistics.pstdev(times)

t_em, s_em = time_posterior(model_em, X_test)
t_emp, s_emp = time_posterior(model_empirical, X_test)
print(f"Posterior (median over reps) — EM-EGP: {t_em*1000:.2f} ms, "
      f"Empirical surrogate: {t_emp*1000:.2f} ms  (speedup {t_em/max(t_emp,1e-9):.1f}×)")

with torch.no_grad():
    p_em = model_em.posterior(X_test)
    p_emp = model_empirical.posterior(X_test)
    max_err_mean = (p_em.mean - p_emp.mean).abs().max().item()
    max_err_std = (
        p_em.variance.sqrt() - p_emp.variance.sqrt()
    ).abs().max().item()
print(f"Max abs difference on grid: mean {max_err_mean:.2e}, "
      f"std {max_err_std:.2e}  (exact reproduction of mean AND covariance)")

Output:

Synthesized Y_full shape: torch.Size([51, 50, 1])  (num_curves = rank + 1)
Posterior (median over reps) — EM-EGP: 2.34 ms, Empirical surrogate: 1.64 ms  (speedup 1.4×)
Max abs difference on grid: mean 1.16e-10, std 4.31e-09  (exact reproduction of mean AND covariance)

D. Multi-output empirical GP on LCBench

LCBench provides multiple metrics per configuration — validation accuracy, training loss, and runtime. We model validation accuracy and training loss jointly as two correlated outputs evolving over epochs. The MultiOutputEmpiricalOneDimensionalGP learns cross-output correlations from historical curves without a handcrafted multi-output kernel.

The purpose of this section is qualitative: to see the shrinkage in posterior variance that comes from modeling the outputs jointly rather than fitting a separate single-output empirical GP to each metric. We plot the intervals and also report held-out accuracy (RMSE), calibration (NLL), and 95% coverage, since tighter intervals are only an improvement if they stay calibrated.

A caveat on when this model is worth using. The joint model estimates a dense cross-output covariance over an (n*m)-dimensional index set from a limited number of historical curves, so it needs enough historical data to estimate and benefit from that structure. Whether LCBench clears that bar is not obvious — the joint model is more rank-limited than the m independent per-output models it competes with, and which one wins depends on the metrics chosen, the number of historical curves, and the observation budget. We include the comparison here for illustration, not as a general claim that joint modeling is better. If you are choosing between them for your own problem, benchmark both.

# Load the second output (training loss) for joint multi-output modeling.
metrics_loss, _, _ = load_lcbench(DATASET_NAME, "Train/loss")
historical_loss = metrics_loss[hist_idx]
test_loss = metrics_loss[test_idx]

from botorch.models.empirical_gps import MultiOutputEmpiricalOneDimensionalGP

# Stack the two outputs: (num_curves, num_epochs, num_outputs=2)
historical_Y_mo = torch.stack([historical_full, historical_loss], dim=-1)
if IS_SMOKE:
    historical_Y_mo = historical_Y_mo[:10]
test_curve_mo = torch.stack([test_full[0], test_loss[0]], dim=-1)  # (num_epochs, 2)

# The outputs live on very different scales (accuracy ~60-90, loss ~0.5-1.9). We
# standardize each output using the historical curves so the shared observation
# noise and the cross-output covariance are commensurate, and un-standardize the
# predictions afterwards.
out_mean = historical_Y_mo.reshape(-1, 2).mean(0)
out_std = historical_Y_mo.reshape(-1, 2).std(0).clamp_min(1e-6)
hist_Y_mo_z = (historical_Y_mo - out_mean) / out_std
test_curve_mo_z = (test_curve_mo - out_mean) / out_std

cutoff_mo = int(0.3 * len(X_grid))  # observe the first 30% of epochs on BOTH outputs
train_X_mo = X_grid[:cutoff_mo]
train_Y_mo_z = test_curve_mo_z[:cutoff_mo]

model_mo = MultiOutputEmpiricalOneDimensionalGP(
    train_X=train_X_mo.unsqueeze(-1),
    train_Y=train_Y_mo_z,
    historical_X=X_grid.unsqueeze(-1),
    historical_Y=hist_Y_mo_z,
)
# Recorded metrics are essentially noise-free, so use a small observation noise:
# with the default (large) noise the model barely conditions on the observations
# and reverts to the prior mean. A small noise makes it interpolate the data.
model_mo.likelihood.noise = torch.tensor(1e-4, dtype=torch.double)
model_mo.eval()
print(f"Multi-output historical Y shape: {tuple(hist_Y_mo_z.shape)}; "
      f"observed {cutoff_mo} of {len(X_grid)} epochs on both outputs.")

Output:

Multi-output historical Y shape: (100, 50, 2); observed 15 of 50 epochs on both outputs.

# Multi-output posterior over the full epoch grid, compared against SINGLE-output EGPs fit to
# each metric separately. Because accuracy and loss are strongly (anti-)correlated, the joint
# model borrows strength across outputs and forecasts each with LOWER predictive variance -- the
# multi-output credible intervals below are visibly tighter than the single-output ones.
titles = ['Validation Accuracy', 'Training Loss']; colors = ['tab:blue', 'tab:green']
with torch.no_grad():
    post_mo = model_mo.posterior(X_grid.unsqueeze(-1))
    mean_z = post_mo.mean.reshape(len(X_grid), 2)
    std_z = post_mo.variance.clamp_min(0).reshape(len(X_grid), 2).sqrt()
mean_plot = (mean_z * out_std + out_mean).numpy(); std_plot = (std_z * out_std).numpy()

so_mean_plot, so_std_plot = [], []
for k in range(2):
    so_k = EmpiricalOneDimensionalGP(
        train_X=train_X_mo.unsqueeze(-1), train_Y=train_Y_mo_z[:, k:k + 1],
        historical_X=X_grid.unsqueeze(-1), historical_Y=hist_Y_mo_z[:, :, k:k + 1])
    so_k.likelihood.noise = torch.tensor(1e-4, dtype=torch.double); so_k.eval()
    with torch.no_grad():
        p_so = so_k.posterior(X_grid.unsqueeze(-1))
        so_mean_plot.append((p_so.mean.squeeze(-1) * out_std[k] + out_mean[k]).numpy())
        so_std_plot.append((p_so.variance.clamp_min(0).sqrt().squeeze(-1) * out_std[k]).numpy())

fig, axes = plt.subplots(1, 2, figsize=(12, 4))
for ax, idx, title, color in zip(axes, [0, 1], titles, colors):
    style_ax(ax, grid=False)
    ax.plot(X_grid.numpy(), test_curve_mo[:, idx].numpy(), color='tab:orange', ls='--', lw=1.5, label='True')
    ax.plot(X_grid[:cutoff_mo].numpy(), test_curve_mo[:cutoff_mo, idx].numpy(), 'ko', ms=4, label='Observed (30%)')
    mu = mean_plot[:, idx]; sd = std_plot[:, idx]
    ax.plot(X_grid.numpy(), mu, color=color, lw=2, label='Multi-output mean')
    ax.fill_between(X_grid.numpy(), mu - 2 * sd, mu + 2 * sd, color=color, alpha=0.30, label='Multi-output 95% CI')
    ax.fill_between(X_grid.numpy(), so_mean_plot[idx] - 2 * so_std_plot[idx], so_mean_plot[idx] + 2 * so_std_plot[idx],
                    facecolor='none', edgecolor='gray', ls='--', lw=1.2, label='Single-output 95% CI')
    ax.axvline(X_grid[cutoff_mo - 1].item(), color='k', ls=':', alpha=0.4)
    ax.set_xlabel('Epoch'); ax.set_ylabel(title); ax.set_title(title)
    ax.grid(alpha=0.3); ax.legend(frameon=False, fontsize=8)
fig.suptitle('Joint (multi-output) vs separate (single-output) EGP: joint modeling fits in-sample and tightens each forecast', y=1.02)
plt.tight_layout(); plt.show()

ins = np.abs(mean_plot[:cutoff_mo] - test_curve_mo[:cutoff_mo].numpy()).max(0)
print(f"Max in-sample error (should be ~0): accuracy {ins[0]:.4f}, loss {ins[1]:.4f}")
for k, nm in enumerate(titles):
    mo_fc = float(std_plot[cutoff_mo:, k].mean()); so_fc = float(so_std_plot[k][cutoff_mo:].mean())
    print(f"{nm}: mean forecast std  multi-output {mo_fc:.3f}  vs  single-output {so_fc:.3f}  "
          f"({100 * (1 - mo_fc / so_fc):.0f}% tighter)")

# Sharpness is only useful if it stays calibrated: report held-out (post-cutoff) RMSE,
# NLL, and 95% CI coverage for the joint (multi-output) vs separate (single-output) fits.
gt_mo = test_curve_mo.numpy()
sl_ho = slice(cutoff_mo, None)
for k, nm in enumerate(titles):
    for tag, mu_k, sd_k in [("multi ", mean_plot[:, k], std_plot[:, k]),
                            ("single", so_mean_plot[k], so_std_plot[k])]:
        err = gt_mo[sl_ho, k] - mu_k[sl_ho]
        var_k = np.clip(sd_k[sl_ho] ** 2, 1e-12, None)
        rmse_k = float(np.sqrt(np.mean(err ** 2)))
        nll_k = float(np.mean(0.5 * np.log(2 * np.pi * var_k) + 0.5 * err ** 2 / var_k))
        cov_k = float(np.mean(np.abs(err) <= 2 * sd_k[sl_ho]))
        print(f"  {nm:18s} {tag}: held-out RMSE={rmse_k:6.3f}  NLL={nll_k:7.3f}  95%cov={cov_k:.2f}")

# The plot/prints above are for one representative curve. To check the payoff
# generalizes, average held-out (post-cutoff) RMSE / NLL / 95% coverage for joint
# (multi-output) vs separate (single-output) fits over ALL held-out test curves.
_metrics = {(tag, k): {"rmse": [], "nll": [], "cov": []}
            for tag in ("multi", "single") for k in range(2)}
for j in range(test_full.shape[0]):
    tc_j = torch.stack([test_full[j], test_loss[j]], dim=-1)
    tc_jz = (tc_j - out_mean) / out_std
    trY_j = tc_jz[:cutoff_mo]
    mo_j = MultiOutputEmpiricalOneDimensionalGP(
        train_X=X_grid[:cutoff_mo].unsqueeze(-1), train_Y=trY_j,
        historical_X=X_grid.unsqueeze(-1), historical_Y=hist_Y_mo_z)
    mo_j.likelihood.noise = torch.tensor(1e-4, dtype=torch.double); mo_j.eval()
    with torch.no_grad():
        pj = mo_j.posterior(X_grid.unsqueeze(-1))
        mmu = (pj.mean.reshape(len(X_grid), 2) * out_std + out_mean).numpy()
        msd = (pj.variance.clamp_min(0).reshape(len(X_grid), 2).sqrt() * out_std).numpy()
    gt_j = tc_j.numpy()
    for k in range(2):
        so_j = EmpiricalOneDimensionalGP(
            train_X=X_grid[:cutoff_mo].unsqueeze(-1), train_Y=trY_j[:, k:k + 1],
            historical_X=X_grid.unsqueeze(-1), historical_Y=hist_Y_mo_z[:, :, k:k + 1])
        so_j.likelihood.noise = torch.tensor(1e-4, dtype=torch.double); so_j.eval()
        with torch.no_grad():
            psj = so_j.posterior(X_grid.unsqueeze(-1))
            smu = (psj.mean.squeeze(-1) * out_std[k] + out_mean[k]).numpy()
            ssd = (psj.variance.clamp_min(0).sqrt().squeeze(-1) * out_std[k]).numpy()
        for tag, mu_a, sd_a in [("multi", mmu[:, k], msd[:, k]), ("single", smu, ssd)]:
            e = gt_j[sl_ho, k] - mu_a[sl_ho]; v = np.clip(sd_a[sl_ho] ** 2, 1e-12, None)
            _metrics[(tag, k)]["rmse"].append(float(np.sqrt(np.mean(e ** 2))))
            _metrics[(tag, k)]["nll"].append(float(np.mean(0.5 * np.log(2 * np.pi * v) + 0.5 * e ** 2 / v)))
            _metrics[(tag, k)]["cov"].append(float(np.mean(np.abs(e) <= 2 * sd_a[sl_ho])))
print(f"\nAveraged over {test_full.shape[0]} held-out curves (held-out epochs):")
for k, nm in enumerate(titles):
    for tag in ("multi", "single"):
        m = _metrics[(tag, k)]
        print(f"  {nm:18s} {tag:6s}: RMSE={np.mean(m['rmse']):6.3f}  "
              f"NLL={np.mean(m['nll']):7.3f}  95%cov={np.mean(m['cov']):.2f}")

Output:

Max in-sample error (should be ~0): accuracy 0.0172, loss 0.0007
Validation Accuracy: mean forecast std  multi-output 0.642  vs  single-output 1.084  (41% tighter)
Training Loss: mean forecast std  multi-output 0.027  vs  single-output 0.036  (25% tighter)
Validation Accuracy multi : held-out RMSE= 0.294  NLL=  0.543  95%cov=1.00
Validation Accuracy single: held-out RMSE= 0.255  NLL=  0.984  95%cov=1.00
Training Loss      multi : held-out RMSE= 0.008  NLL= -2.742  95%cov=1.00
Training Loss      single: held-out RMSE= 0.010  NLL= -2.487  95%cov=1.00
Averaged over 5 held-out curves (held-out epochs):
Validation Accuracy multi : RMSE= 0.286  NLL=  0.596  95%cov=0.99
Validation Accuracy single: RMSE= 0.690  NLL=  1.601  95%cov=0.86
Training Loss      multi : RMSE= 0.013  NLL= -2.562  95%cov=0.99
Training Loss      single: RMSE= 0.009  NLL= -2.487  95%cov=1.00

E. Multi-task empirical GP: zero-shot & few-shot cross-task transfer

Different LCBench datasets are related tasks (a configuration that trains well on one often trains well on another). Here we show the multi-task empirical GP can predict a configuration's learning curve on a target task from zero or few observations on that task, transferring configuration-specific structure through the cross-task covariance learned from historical curves.

The empirical cross-task covariance is rank-limited (estimated from only a few dozen historical curves), so it benefits from a full-rank, structured complement. Rather than an ad-hoc observation-noise ridge, we use the library's principled additive base kernel (BaseAugmentedEmpiricalKernel via base_covar_module): we add a fitted ICM base kernel — the same task-covariance $\otimes$ smooth log-epoch kernel we fit for the parametric baseline below — to the empirical covariance, forming $\Sigma_\text{empirical} + K_\text{ICM}$. We freeze the ICM's learned lengthscales (the transferred structure) and fit only its outputscale (magnitude) and the observation noise jointly by marginal likelihood, reusing the cross-task structure without re-optimizing the whole kernel. This full-rank additive component improves calibration without sacrificing accuracy.

We compare against a standard BoTorch MultiTaskGP over (epoch, task). Because its inputs carry no configuration identity, it treats every curve of a task as a noisy draw of a single per-task function, so it is config-agnostic: its best guess for any configuration is the per-task mean. The empirical multi-task GP, by contrast, models each configuration as a whole curve and so transfers configuration-specific structure across tasks.

# Zero-shot cross-task transfer: three correlated LCBench datasets as tasks.
# jasmine (index 2) is the TARGET task and receives ZERO observations; vehicle and segment
# (the donor tasks) are observed for the query configuration.
from botorch.models.empirical_gps import MultiTaskEmpiricalOneDimensionalGP
from botorch.models import MultiTaskGP as BoTorchMultiTaskGP

task_names_e = ["vehicle", "segment", "jasmine"]
donor_tasks = [0, 1]
target_task = 2
task_metrics_e = [load_lcbench(n, "Train/val_accuracy")[0] for n in task_names_e]

n_cfg_e = 200
torch.manual_seed(99)
perm_e = torch.randperm(n_cfg_e)
n_hist_e = 20 if IS_SMOKE else 60
E_grid = torch.arange(1, task_metrics_e[0].shape[1] + 1, dtype=torch.double)
donor_stride = 5  # ~10 observed epochs per donor task

# Diagnostic only: measure task correlation on the HISTORICAL configs (not the eval
# configs) so this sanity check does not peek at the evaluation set.
final_acc = torch.stack([tm[perm_e[:n_hist_e], -1] for tm in task_metrics_e])
corr = torch.corrcoef(final_acc)
print(f"Tasks: {task_names_e} (target = '{task_names_e[target_task]}', zero observations).")
print(f"Cross-task correlation of final accuracy (target vs donors): "
      f"{corr[target_task, 0]:.2f}, {corr[target_task, 1]:.2f} -> transfer is feasible.")
print(f"{n_hist_e} historical configs per task; the query config is observed only on "
      f"{[task_names_e[t] for t in donor_tasks]}.")

Output:

Tasks: ['vehicle', 'segment', 'jasmine'] (target = 'jasmine', zero observations).
Cross-task correlation of final accuracy (target vs donors): 0.76, 0.82 -> transfer is feasible.
60 historical configs per task; the query config is observed only on ['vehicle', 'segment'].

# Compare the empirical multi-task GP against a strong config-specific parametric baseline and
# simpler references, on accuracy (RMSE) and calibration (NLL). To make the comparison
# statistically sound (not a single-split artifact) we average over several random task-history
# splits and report standard errors (SEM).
#   (A) EM multi-task GP with an ADDITIVE base kernel (exercises the library
#       `base_covar_module` feature). The empirical cross-task covariance is rank-limited
#       (only ~60 historical curves), so we give it a full-rank, structured additive
#       component: the FITTED ICM base kernel (B), which carries the learned cross-task
#       correlations. The covariance is ``Sigma_empirical + K_ICM``. We freeze the ICM's
#       lengthscales (the transferred structure) and fit only its outputscale (magnitude)
#       and the observation noise, jointly by marginal likelihood -- reusing the learned
#       structure without re-optimizing the whole kernel.
#   (B) Parametric MTGP (SumMLL, config-specific): an ICM covariance (task-cov (x) log-epoch RBF)
#       fit by summed marginal likelihood over the historical config-curves -- the parametric
#       analog of our pre-training, and the base kernel that (A) shrinks toward.
#   (C) MultiTaskGP (config-agnostic). (D) Prior mean.
from gpytorch.kernels import Kernel as GPyTorchKernel, ScaleKernel

N_SEEDS = 1 if IS_SMOKE else 5
n_targets = [0, 10] if IS_SMOKE else [0, 2, 5, 10, 20]
n_eval = 5 if IS_SMOKE else 40
n_fit_iters = 8 if IS_SMOKE else 75          # MLL steps for alpha + noise per config
obs_e_idx = list(range(0, len(E_grid), donor_stride))
ne_e = len(E_grid); ntask_e = len(task_names_e); Elog_e = torch.log(E_grid)


class ICMBaseKernel(GPyTorchKernel):
    """The fitted ICM as a GPyTorch base kernel over [epoch, task] inputs.

    Computes B[t1, t2] * outputscale * exp(-0.5 (log e1 - log e2)^2 / lengthscale^2): the
    parametric task-covariance (x) smooth log-epoch RBF fit in (B), exposed so it can be passed as
    `base_covar_module` for the empirical model to add. Parameters are stored as buffers (frozen);
    wrapped in a ScaleKernel, only its outputscale (magnitude) and the observation noise are
    learned. This *base* kernel is problem-specific and user-provided; the reusable additive
    machinery (BaseAugmentedEmpiricalKernel) lives in botorch.
    """

    def __init__(self, B, lengthscale, outputscale):
        super().__init__()
        self.register_buffer("B", B)
        self.register_buffer("log_ls", torch.log(torch.as_tensor(lengthscale, dtype=B.dtype)))
        self.register_buffer("log_os", torch.log(torch.as_tensor(outputscale, dtype=B.dtype)))

    def forward(self, x1, x2, diag=False, **kwargs):
        le1 = x1[..., 0].clamp_min(1e-8).log(); t1 = x1[..., 1].long()
        le2 = x2[..., 0].clamp_min(1e-8).log(); t2 = x2[..., 1].long()
        ls2 = torch.exp(2.0 * self.log_ls)
        Kee = torch.exp(self.log_os) * torch.exp(
            -0.5 * (le1.unsqueeze(-1) - le2.unsqueeze(-2)) ** 2 / ls2
        )
        Btt = self.B[t1.unsqueeze(-1), t2.unsqueeze(-2)]
        K = Btt * Kee
        return K.diagonal(dim1=-2, dim2=-1) if diag else K


def _fit_icm(hist, tmean):
    nh = hist[0].shape[0]
    Hc = torch.stack([torch.cat([hist[t][c] - tmean[t] for t in range(ntask_e)]) for c in range(nh)]).double()
    ll = torch.zeros(1, dtype=torch.double, requires_grad=True); lo = torch.zeros(1, dtype=torch.double, requires_grad=True)
    Lr = torch.eye(ntask_e, dtype=torch.double).clone().requires_grad_(True)
    lnz = torch.tensor([-2.0], dtype=torch.double, requires_grad=True)
    D = (Elog_e[:, None] - Elog_e[None, :]) ** 2; MM = ntask_e * ne_e; I = torch.eye(MM, dtype=torch.double)
    opt = torch.optim.Adam([ll, lo, Lr, lnz], lr=0.05)
    for _ in range(60 if IS_SMOKE else 250):
        opt.zero_grad()
        Kee = lo.exp() * torch.exp(-0.5 * D / ll.exp() ** 2)
        L = torch.tril(Lr); B = L @ L.T + 1e-4 * torch.eye(ntask_e, dtype=torch.double)
        C = torch.linalg.cholesky(torch.kron(B, Kee) + lnz.exp() * I + 1e-5 * I)
        al = torch.cholesky_solve(Hc.T, C)
        (0.5 * (Hc.T * al).sum() + 0.5 * nh * 2 * torch.log(torch.diag(C)).sum()).backward(); opt.step()
    with torch.no_grad():
        Kee = lo.exp() * torch.exp(-0.5 * D / ll.exp() ** 2)
        L = torch.tril(Lr); B = L @ L.T + 1e-4 * torch.eye(ntask_e, dtype=torch.double)
        # Expose the fitted ICM as a base kernel (for shrinkage) and return the grid
        # quantities (B, Kee, noise) used by the parametric-baseline conditioning below.
        base_kernel = ICMBaseKernel(B.clone(), float(ll.exp()), float(lo.exp()))
        return base_kernel, B, Kee, float(lnz.exp())


def _em_transfer(ci, n_target, hist, base_kernel):
    tX, tY = [], []
    for t in donor_tasks:
        for e in obs_e_idx:
            tX.append([E_grid[e].item(), float(t)]); tY.append(task_metrics_e[t][ci, e].item())
    for e in range(n_target):
        tX.append([E_grid[e].item(), float(target_task)]); tY.append(task_metrics_e[target_task][ci, e].item())
    tX = torch.tensor(tX, dtype=torch.double); tY = torch.tensor(tY, dtype=torch.double).unsqueeze(-1)
    # Additive combination Sigma_empirical + K_base. Wrapping the (frozen-buffer) ICM in a
    # ScaleKernel adds a single trainable outputscale, so only the base magnitude and the
    # observation noise are fit -- the ICM's learned structure (lengthscales) is reused.
    m = MultiTaskEmpiricalOneDimensionalGP(
        train_X=tX, train_Y=tY, task_feature=-1,
        historical_Xs=[E_grid.unsqueeze(-1) for _ in task_names_e], historical_Ys=hist,
        base_covar_module=ScaleKernel(base_kernel))
    # Fit the base kernel's outputscale (magnitude) and the observation noise JOINTLY by
    # marginal likelihood; the ICM structure inside is frozen (buffers).
    mll = ExactMarginalLogLikelihood(m.likelihood, m)
    m.train()
    opt = torch.optim.Adam([p for p in m.parameters() if p.requires_grad], lr=0.1)
    for _ in range(n_fit_iters):
        opt.zero_grad(); loss = -mll(m(tX), m.train_targets); loss.backward(); opt.step()
    m.eval()
    with torch.no_grad():
        p = m.posterior(torch.stack([E_grid, torch.full_like(E_grid, float(target_task))], dim=-1))
        return (p.mean.squeeze(-1).numpy(), p.variance.clamp_min(1e-9).squeeze(-1).numpy(),
                float(m.covar_module.base_kernel.outputscale))


def _param_transfer(ci, n_target, tmean, B, Kee, nz):
    obs = [(t, e) for t in donor_tasks for e in obs_e_idx] + [(target_task, e) for e in range(n_target)]
    yo = torch.tensor([(task_metrics_e[t][ci, e] - tmean[t][e]).item() for (t, e) in obs], dtype=torch.double)
    Koo = torch.tensor([[(B[t1, t2] * Kee[e1, e2]).item() for (t2, e2) in obs] for (t1, e1) in obs], dtype=torch.double)
    Kqo = torch.tensor([[(B[target_task, t2] * Kee[eq, e2]).item() for (t2, e2) in obs] for eq in range(ne_e)], dtype=torch.double)
    Kqq = torch.tensor([(B[target_task, target_task] * Kee[eq, eq]).item() for eq in range(ne_e)], dtype=torch.double)
    A = Koo + nz * torch.eye(len(obs), dtype=torch.double)
    mu = (Kqo @ torch.linalg.solve(A, yo.unsqueeze(-1))).squeeze(-1) + tmean[target_task]
    # Report the LATENT posterior variance (no +nz observation-noise term) so the
    # NLL/calibration comparison is consistent with the EM and MultiTaskGP baselines,
    # which both use the latent posterior().variance.
    var = (Kqq - (Kqo * torch.linalg.solve(A, Kqo.T).T).sum(1)).clamp_min(1e-9)
    return mu.numpy(), var.numpy()


def _fit_mtgp(hist):
    torch.manual_seed(0); sub = torch.randperm(hist[0].shape[0])[:min(20, hist[0].shape[0])]
    Xtr, Ytr = [], []
    for t in range(ntask_e):
        for ci in sub.tolist():
            for e in obs_e_idx:
                Xtr.append([E_grid[e].item(), float(t)]); Ytr.append(hist[t][ci, e].item())
    Xtr = torch.tensor(Xtr, dtype=torch.double); Ytr = torch.tensor(Ytr, dtype=torch.double).unsqueeze(-1)
    g = BoTorchMultiTaskGP(Xtr, Ytr, task_feature=-1, output_tasks=[target_task])
    fit_gpytorch_mll(ExactMarginalLogLikelihood(g.likelihood, g)); g.eval()
    with torch.no_grad():
        pa = g.posterior(torch.stack([E_grid, torch.full_like(E_grid, float(target_task))], dim=-1))
        return pa.mean.squeeze(-1).numpy(), pa.variance.clamp_min(1e-9).squeeze(-1).numpy()

methods = ["EM multi-task GP", "Parametric MTGP (SumMLL, config-specific)",
           "MultiTaskGP (config-agnostic)", "Prior mean"]
RM = {m: {nt: [] for nt in n_targets} for m in methods}
NLd = {m: {nt: [] for nt in n_targets} for m in methods}
os_fit = {nt: [] for nt in n_targets}            # MLL-fit base-kernel outputscale per target-obs budget
example = {}
for si in range(N_SEEDS):
    torch.manual_seed(1000 + si); perm = torch.randperm(n_cfg_e)
    hist = [tm[perm[:n_hist_e]] for tm in task_metrics_e]
    tmean = torch.stack([hist[t].mean(0) for t in range(ntask_e)])
    evals = [int(perm[k]) for k in range(n_hist_e, n_hist_e + n_eval)]
    base_kernel, B, Kee, nz = _fit_icm(hist, tmean)
    mtgp_mu, mtgp_var = _fit_mtgp(hist)
    prior_mu = tmean[target_task].numpy(); prior_var = hist[target_task].var(0).clamp_min(1e-9).numpy()
    for nt in n_targets:
        for ci in evals:
            gt = task_metrics_e[target_task][ci].numpy()
            # Score ONLY the held-out epochs (those after the nt observed target epochs);
            # including the observed epochs would credit interpolation of the training
            # observations as the budget grows.
            ev = slice(nt, None)
            em_mu, em_var, em_os = _em_transfer(ci, nt, hist, base_kernel)
            pm_mu, pm_var = _param_transfer(ci, nt, tmean, B, Kee, nz)
            os_fit[nt].append(em_os)
            r, nl = rmse_nll(em_mu[ev], gt[ev], var=em_var[ev]); RM[methods[0]][nt].append(r); NLd[methods[0]][nt].append(nl)
            r, nl = rmse_nll(pm_mu[ev], gt[ev], var=pm_var[ev]); RM[methods[1]][nt].append(r); NLd[methods[1]][nt].append(nl)
            r, nl = rmse_nll(mtgp_mu[ev], gt[ev], var=mtgp_var[ev]); RM[methods[2]][nt].append(r); NLd[methods[2]][nt].append(nl)
            r, nl = rmse_nll(prior_mu[ev], gt[ev], var=prior_var[ev]); RM[methods[3]][nt].append(r); NLd[methods[3]][nt].append(nl)
    if si == 0:  # representative zero-shot example = the median-RMSE eval config (not cherry-picked)
        ranked = sorted(
            (rmse_nll(_em_transfer(ci, 0, hist, base_kernel)[0], task_metrics_e[target_task][ci].numpy())[0], ci)
            for ci in evals
        )
        ci_ex = ranked[len(ranked) // 2][1]
        em_mu0, em_var0, _ = _em_transfer(ci_ex, 0, hist, base_kernel); pm_mu0, _ = _param_transfer(ci_ex, 0, tmean, B, Kee, nz)
        example = dict(gt=task_metrics_e[target_task][ci_ex].numpy(), em_mu=em_mu0, em_sd=np.sqrt(em_var0), pm_mu=pm_mu0)

def _agg(T, m, nt):
    # Cluster by history split: configs within a split share the same fitted prior/ICM,
    # so they are NOT independent. Average within each split, then take the mean and
    # standard error ACROSS the (few) splits.
    a = np.array(T[m][nt]).reshape(N_SEEDS, -1)
    sm = a.mean(axis=1)
    sem = sm.std(ddof=1) / np.sqrt(len(sm)) if len(sm) > 1 else float("nan")
    return sm.mean(), sem

for label, T in [("RMSE (lower=better)", RM), ("NLL (lower=better calibrated)", NLd)]:
    print(f"\n{label} on '{task_names_e[target_task]}' (held-out epochs; {N_SEEDS} splits x {n_eval} configs, mean +/- split-clustered SEM over {N_SEEDS} splits):")
    print(f"  {'# target obs':<44}" + "".join(f"{n:<13}" for n in n_targets))
    for m in methods:
        print(f"  {m:<44}" + "".join(f"{_agg(T, m, n)[0]:5.2f}+/-{_agg(T, m, n)[1]:<4.2f} " for n in n_targets))

# The MLL-fit additive base-kernel outputscale (magnitude) per target-obs budget: the model
# scales the fitted-ICM additive component to complement the empirical covariance.
print(f"\nMLL-fit additive base-kernel outputscale (fitted ICM), mean +/- SEM:")
print(f"  {'# target obs':<44}" + "".join(f"{n:<13}" for n in n_targets))
def _agg_a(nt):
    a = np.array(os_fit[nt]).reshape(N_SEEDS, -1)
    sm = a.mean(axis=1)
    sem = sm.std(ddof=1) / np.sqrt(len(sm)) if len(sm) > 1 else float("nan")
    return sm.mean(), sem
print(f"  {'outputscale':<44}" + "".join(f"{_agg_a(n)[0]:5.3f}+/-{_agg_a(n)[1]:<5.3f}" for n in n_targets))

# Plot: RMSE and NLL with +/-SEM bands, plus a representative zero-shot example curve.
plot_methods = ["EM multi-task GP", "Parametric MTGP (SumMLL, config-specific)", "Prior mean"]
styles = {"EM multi-task GP": dict(color='tab:red', marker='o', lw=2.2),
          "Parametric MTGP (SumMLL, config-specific)": dict(color='tab:purple', marker='D', lw=1.6, ls='-.'),
          "Prior mean": dict(color='black', marker='', lw=1.2, ls=':')}
fig, axes = plt.subplots(1, 3, figsize=(16, 4.2))
for ax, T, ylab, ttl in [(axes[0], RM, 'RMSE', 'Accuracy (RMSE)'), (axes[1], NLd, 'NLL', 'Calibration (NLL)')]:
    style_ax(ax, grid=False)
    for m in plot_methods:
        mus = np.array([_agg(T, m, nt)[0] for nt in n_targets]); ses = np.array([_agg(T, m, nt)[1] for nt in n_targets])
        ax.plot(n_targets, mus, label=m, **styles[m])
        ax.fill_between(n_targets, mus - ses, mus + ses, color=styles[m]['color'], alpha=0.15)
    ax.set_xlabel('# target-task observations'); ax.set_ylabel(ylab + ' on target task')
    ax.set_title(ttl + ' (mean $\\pm$ SEM)', fontsize=11); ax.grid(alpha=0.3); ax.legend(frameon=False, fontsize=7.5)
style_ax(axes[2], grid=False)
axes[2].plot(E_grid.numpy(), example['gt'], color='tab:orange', ls='--', lw=1.8, label=f'{task_names_e[target_task]} true (unobserved)')
axes[2].plot(E_grid.numpy(), example['em_mu'], color='tab:red', lw=2, label='EM MT-GP zero-shot')
axes[2].fill_between(E_grid.numpy(), example['em_mu'] - 2 * example['em_sd'], example['em_mu'] + 2 * example['em_sd'], color='tab:red', alpha=0.18, label='95% CI')
axes[2].plot(E_grid.numpy(), example['pm_mu'], color='tab:purple', ls='-.', lw=1.5, label='Parametric MTGP')
axes[2].set_xlabel('Epoch'); axes[2].set_ylabel('Validation accuracy')
axes[2].set_title(f'Zero-shot example (median-RMSE config): predict {task_names_e[target_task]}', fontsize=11)
axes[2].grid(alpha=0.3); axes[2].legend(frameon=False, fontsize=7.5)
fig.suptitle('Empirical multi-task GP (base-kernel shrinkage) vs a strong parametric multi-task GP: accuracy & calibration', y=1.02)
plt.tight_layout(); plt.show()

Output:

RMSE (lower=better) on 'jasmine' (held-out epochs; 5 splits x 40 configs, mean +/- split-clustered SEM over 5 splits):
# target obs                                0            2            5            10           20           
EM multi-task GP                             4.91+/-0.31  4.45+/-0.33  3.41+/-0.27  2.51+/-0.17  1.21+/-0.11 
Parametric MTGP (SumMLL, config-specific)    5.80+/-0.29  5.18+/-0.33  4.69+/-0.35  4.09+/-0.34  2.55+/-0.24 
MultiTaskGP (config-agnostic)                9.69+/-0.59  9.45+/-0.66  9.20+/-0.73  8.99+/-0.78  8.79+/-0.82 
Prior mean                                   8.85+/-0.29  8.62+/-0.32  8.35+/-0.34  8.00+/-0.35  7.56+/-0.34 
NLL (lower=better calibrated) on 'jasmine' (held-out epochs; 5 splits x 40 configs, mean +/- split-clustered SEM over 5 splits):
# target obs                                0            2            5            10           20           
EM multi-task GP                             3.24+/-0.13  3.14+/-0.13  2.99+/-0.20  2.47+/-0.11  1.66+/-0.04 
Parametric MTGP (SumMLL, config-specific)    3.28+/-0.06  3.22+/-0.07  3.11+/-0.06  2.90+/-0.06  2.46+/-0.13 
MultiTaskGP (config-agnostic)                9.73+/-1.27  9.07+/-1.12  9.40+/-1.19  9.09+/-1.13  8.69+/-1.05 
Prior mean                                   3.75+/-0.01  3.74+/-0.01  3.73+/-0.01  3.71+/-0.01  3.68+/-0.01 
MLL-fit additive base-kernel outputscale (fitted ICM), mean +/- SEM:
# target obs                                0            2            5            10           20           
outputscale                                 0.056+/-0.0060.063+/-0.0080.063+/-0.0080.069+/-0.0090.075+/-0.012
summing_matrix = torch.sparse_coo_tensor(

Takeaway. The empirical multi-task GP transfers configuration-specific learning-curve structure through the historical cross-task covariance, predicting a target task's curve from zero or few observations on it. Scoring only the held-out epochs and averaging over a handful of random history splits (with split-clustered standard errors), it is as good as or better than the strong config-specific parametric MTGP on accuracy (RMSE) and calibration (NLL) across observation budgets, while both clearly beat the config-agnostic MultiTaskGP. With only a few splits on a single target task the error bars are wide, so we read this ranking as indicative rather than definitive. The rank-limited empirical covariance is regularized principledly by the library's additive base kernel (BaseAugmentedEmpiricalKernel via base_covar_module): adding the fitted ICM ($\Sigma_\text{empirical} + K_\text{ICM}$) with its learned structure frozen and only its outputscale and the observation noise fit jointly by marginal likelihood. This full-rank additive component delivers the calibration gain (NLL) without hurting transfer accuracy.

F. 7D EM GP on LCBench final performance -- beyond 1D interpolation

The 1D sections used a progression axis (epoch). The EM-based Empirical GP works for arbitrary inputs. Here the input is LCBench's 7 hyperparameters and the target is final validation accuracy -- no progression axis, just $X \in \mathbb{R}^7 \to y \in \mathbb{R}$.

Because LCBench evaluates the same configurations on every dataset, each dataset is an aligned sample of the (config $\to$ accuracy) function over shared inputs. We learn an EM prior across ~30 historical datasets and transfer it to 5 held-out datasets, comparing across training-set sizes (averaged over held-out datasets and random splits). We report both RMSE (accuracy) and NLL (calibration).

A single shared kernel is learned across the historical datasets via the summed marginal likelihood (SumMLL) -- the principled "shared hyperparameters, independent datasets" objective -- and reused by every model below: as the pre-trained-GP baselines' transferred kernel, and as the EM prior's covariance basis. The models differ only in how much of that structure they carry over to a new task:

• EM-EGP -- transfers the full empirical prior (mean and covariance) over configurations.
• Pre-trained GP (SumMLL, frozen) -- inherits the shared kernel and freezes it entirely (lengthscales + outputscale); only the mean and noise are fit on the new task.
• Pre-trained GP (SumMLL, scale-tuned) -- inherits the shared kernel but freezes only the lengthscales (the transferable "shape") and re-fits the outputscale, mean, and noise on the new task.
• Vanilla GP -- a SingleTaskGP fit from scratch (no transfer).
• Global mean -- predicts the pre-training mean (a trivial floor).

As we will see, EM-EGP wins in both metrics at every $n$, while the parametric pre-trained-GP baselines are surprisingly finicky.

# F. 7D EM GP setup: 7 hyperparameters -> final validation accuracy, transferred across datasets.
from botorch.utils.lcbench import LCBENCH_DATASET_NAMES
from botorch.models.empirical_gps.utils import ExperimentDataset as ED7
from botorch.models.empirical_gps.em_empirical_gp import build_shared_gp_model_list
from botorch.utils.constraints import LogTransformedInterval
from gpytorch.likelihoods import GaussianLikelihood

METRIC7 = "Train/val_accuracy"
allnames = list(LCBENCH_DATASET_NAMES)
if IS_SMOKE:
    allnames = allnames[:8]           # fewer datasets in smoke mode for speed
n_configs_7d = 60 if IS_SMOKE else 200
n_eval_7d = 1 if IS_SMOKE else 5
n_test_7d = 30
n_train_sizes_7d = [5, 20] if IS_SMOKE else [5, 10, 20, 50, 100]
n_splits_7d = 1 if IS_SMOKE else 3
n_em_7d = 8 if IS_SMOKE else 50
seed7 = 42

print(f"Loading {len(allnames)} LCBench datasets (7D hyperparameters -> final val accuracy)...")
X_all_7d, Y_all_7d = [], []
for nm in allnames:
    metrics7, params7, _ = load_lcbench(nm, METRIC7)
    X_all_7d.append(params7)              # (2000, 7)
    Y_all_7d.append(metrics7[:, -1:])     # final-epoch val accuracy, (2000, 1)

X_ref7 = X_all_7d[0]
matched_7d = all(torch.allclose(X, X_ref7, atol=1e-6) for X in X_all_7d[1:])
Xmin7 = X_ref7.min(0).values
Xrng7 = (X_ref7.max(0).values - Xmin7).clamp_min(1e-8)
Xn7 = (X_ref7 - Xmin7) / Xrng7            # normalize inputs to [0, 1]
print(f"All datasets share the same {X_ref7.shape[0]} configs (matched inputs): {matched_7d}")

# Deterministic split (seed 42): subsample shared configs, then hold out datasets.
torch.manual_seed(seed7)
cfg7 = torch.randperm(X_ref7.shape[0])[:n_configs_7d]
dperm7 = torch.randperm(len(allnames))
eval_ix7 = dperm7[:n_eval_7d].tolist()
pre_ix7 = dperm7[n_eval_7d:].tolist()

X_shared_7d = Xn7[cfg7]
Yc7 = [Y[cfg7] for Y in Y_all_7d]
# Standardize targets by the pre-training pool (applied to eval too).
Ypool7 = torch.cat([Yc7[i] for i in pre_ix7])
ymean7 = Ypool7.mean()
ystd7 = Ypool7.std().clamp_min(1e-8)
pretrain_7d = [ED7(X=X_shared_7d, Y=(Yc7[i] - ymean7) / ystd7) for i in pre_ix7]
print(f"Pre-train on {len(pre_ix7)} datasets, evaluate on {len(eval_ix7)} held-out datasets: "
      f"{', '.join(allnames[i] for i in eval_ix7)}")
print(f"{n_configs_7d} shared configs, RMSE/NLL reported in standardized units.")

Output:

Loading 35 LCBench datasets (7D hyperparameters -> final val accuracy)...
All datasets share the same 2000 configs (matched inputs): True
Pre-train on 30 datasets, evaluate on 5 held-out datasets: albert, bank-marketing, jannis, sylvine, kr-vs-kp
200 shared configs, RMSE/NLL reported in standardized units.

# One shared kernel is learned across the historical datasets via the summed marginal
# likelihood (SumMLL) and reused everywhere -- as the pre-trained-GP baselines' transfer
# kernel AND as the EM prior's covariance basis (fully consistent; no per-task re-training,
# no pooling).
mean_s = ConstantMean()
covar_s = ScaleKernel(
    MaternKernel(nu=2.5, ard_num_dims=7,
                 lengthscale_constraint=LogTransformedInterval(0.01, 100.0, initial_value=1.0)),
    outputscale_constraint=LogTransformedInterval(0.01, 100.0, initial_value=1.0),
)
_, shared_mll_7d = build_shared_gp_model_list(pretrain_7d, mean_s, covar_s)
fit_gpytorch_mll(shared_mll_7d)                       # SumMLL kernel: baselines + EM basis

print("Shared SumMLL kernel fitted. Pre-training 7D EM prior...")
t0 = time.time()
em_prior_7d = pretrain_em_prior(
    datasets=pretrain_7d, mean_module=mean_s, covar_module=covar_s,
    likelihood_noise=torch.tensor(1e-2, dtype=torch.double),
    num_em_iterations=n_em_7d, enable_interpolation=True,
)
print(f"EM done in {time.time() - t0:.1f}s")


def _pretrained_gp(train_X, train_Y, tune_scale):
    """SingleTaskGP that inherits the shared SumMLL kernel (reused, not re-trained).
    Lengthscales -- the transferable 'shape' -- are always frozen; the mean and noise
    are fit on the eval task. If tune_scale=True the outputscale is also re-fit,
    otherwise it stays frozen at the pre-trained value."""
    cc = ScaleKernel(MaternKernel(nu=2.5, ard_num_dims=7))
    with torch.no_grad():
        cc.raw_outputscale.data.copy_(covar_s.raw_outputscale.data)
        cc.base_kernel.raw_lengthscale.data.copy_(covar_s.base_kernel.raw_lengthscale.data)
    gp = SingleTaskGP(train_X, train_Y, covar_module=cc, outcome_transform=None)
    gp.covar_module.base_kernel.raw_lengthscale.requires_grad_(False)
    if not tune_scale:
        gp.covar_module.raw_outputscale.requires_grad_(False)
    fit_gpytorch_mll(ExactMarginalLogLikelihood(gp.likelihood, gp))
    gp.eval()
    return gp


# Metric helpers delegate to the single shared rmse_nll defined in Part 1
# (standardized units here); one implementation, used everywhere.
def _rmse(mu, y):
    return rmse_nll(mu, y)[0]


def _nll(mu, var, y):
    return rmse_nll(mu, y, var=var)[1]


methods_7d = [
    "EM-EGP",
    "Pre-trained GP (SumMLL, frozen)",
    "Pre-trained GP (SumMLL, scale-tuned)",
    "Vanilla GP",
    "Global mean",
]
rmse_7d = {m: {n: [] for n in n_train_sizes_7d} for m in methods_7d}
nll_7d = {m: {n: [] for n in n_train_sizes_7d} for m in methods_7d}
em_fit_failures = []  # record (dataset, n_train, split, error) for any EM-EGP fit that fails
gm_var7 = torch.cat([d.Y for d in pretrain_7d]).var()

# Evaluation train/test splits use a fixed seed (the dataset split above is seeded
# separately); the qualitative ranking is seed-independent.
eval_seed7 = 7
for eidx, ei in enumerate(eval_ix7):
    dY = (Yc7[ei] - ymean7) / ystd7
    for nt in n_train_sizes_7d:
        for split in range(n_splits_7d):
            torch.manual_seed(eval_seed7 + eidx * 1000 + split)
            perm = torch.randperm(n_configs_7d)
            te, tr = perm[:n_test_7d], perm[n_test_7d:n_test_7d + nt]
            trX, trY = X_shared_7d[tr], dY[tr]
            teX, teY = X_shared_7d[te], dY[te].squeeze(-1)

            # EM-EGP: transfers the full empirical prior (mean + covariance) plus a
            # fresh, fully-trainable additive base kernel (Sigma + K_base). Its
            # lengthscales + outputscale and the noise are fit by the marginal
            # likelihood (no backprop through EM), so the model adapts toward a
            # standard GP as data grow while keeping the empirical prior's strong
            # low-data structure.
            lik = GaussianLikelihood()
            base = ScaleKernel(
                MaternKernel(nu=2.5, ard_num_dims=7,
                    lengthscale_constraint=LogTransformedInterval(0.01, 100.0, initial_value=1.0)),
                outputscale_constraint=LogTransformedInterval(0.01, 100.0, initial_value=1.0),
            )
            m = EMEmpiricalGaussianProcess.from_pretrained(
                em_prior=em_prior_7d, train_X=trX, train_Y=trY, likelihood=lik,
                base_covar_module=base)
            m.train(); lik.train()
            try:
                fit_gpytorch_mll(ExactMarginalLogLikelihood(lik, m))
            except Exception as e:
                # Do not silently swallow: a failed fit falls back to the frozen
                # pre-trained prior (a reasonable default), but we record it so the
                # fraction of failed fits is reported rather than hidden.
                em_fit_failures.append((int(ei), int(nt), int(split), repr(e)))
            m.eval(); lik.eval()
            with torch.no_grad():
                po = lik(m(teX))
                rmse_7d["EM-EGP"][nt].append(_rmse(po.mean, teY))
                nll_7d["EM-EGP"][nt].append(_nll(po.mean, po.variance, teY))

            # Pre-trained GP baselines -- both reuse the shared SumMLL kernel
            for name, ts in [("Pre-trained GP (SumMLL, frozen)", False),
                             ("Pre-trained GP (SumMLL, scale-tuned)", True)]:
                g = _pretrained_gp(trX, trY, ts)
                with torch.no_grad():
                    po = g.likelihood(g(teX))
                    rmse_7d[name][nt].append(_rmse(po.mean.squeeze(-1), teY))
                    nll_7d[name][nt].append(_nll(po.mean.squeeze(-1), po.variance.squeeze(-1), teY))

            # Vanilla GP: from scratch on the eval data alone (no transfer)
            gv = SingleTaskGP(trX, trY, outcome_transform=None)
            fit_gpytorch_mll(ExactMarginalLogLikelihood(gv.likelihood, gv)); gv.eval()
            with torch.no_grad():
                po = gv.likelihood(gv(teX))
                rmse_7d["Vanilla GP"][nt].append(_rmse(po.mean.squeeze(-1), teY))
                nll_7d["Vanilla GP"][nt].append(_nll(po.mean.squeeze(-1), po.variance.squeeze(-1), teY))

            # Global mean baseline
            rmse_7d["Global mean"][nt].append(_rmse(torch.zeros_like(teY), teY))
            nll_7d["Global mean"][nt].append(_nll(torch.zeros_like(teY), gm_var7 * torch.ones_like(teY), teY))


n_runs_7d = len(eval_ix7) * len(n_train_sizes_7d) * n_splits_7d
if em_fit_failures:
    print(f"\n[warning] EM-EGP fit did not converge in {len(em_fit_failures)} of "
          f"{n_runs_7d} runs; those runs use the frozen pre-trained prior. "
          f"First failure: {em_fit_failures[0][:3]}")
else:
    print(f"\nAll {n_runs_7d} EM-EGP fits converged.")


def _avg(T, m, n):
    return sum(T[m][n]) / len(T[m][n])


for label, T in [("RMSE (standardized, lower=better)", rmse_7d), ("NLL (lower=better calibrated)", nll_7d)]:
    print(f"\n{label}:")
    print(f"{'model':<40}" + "".join(f"n={n:<7}" for n in n_train_sizes_7d))
    for m in methods_7d:
        print(f"{m:<40}" + "".join(f"{_avg(T, m, n):<9.3f}" for n in n_train_sizes_7d))

Output:

Shared SumMLL kernel fitted. Pre-training 7D EM prior...
EM done in 4.0s
All 75 EM-EGP fits converged.
RMSE (standardized, lower=better):
model                                   n=5      n=10     n=20     n=50     n=100    
EM-EGP                                  0.245    0.244    0.232    0.189    0.191    
Pre-trained GP (SumMLL, frozen)         0.369    0.319    0.234    0.208    0.196    
Pre-trained GP (SumMLL, scale-tuned)    0.370    0.319    0.236    0.207    0.192    
Vanilla GP                              0.396    0.458    0.276    0.248    0.191    
Global mean                             0.656    0.656    0.656    0.656    0.656    
NLL (lower=better calibrated):
model                                   n=5      n=10     n=20     n=50     n=100    
EM-EGP                                  -0.064   -0.306   -0.151   -0.029   -0.160   
Pre-trained GP (SumMLL, frozen)         0.959    0.857    0.708    0.505    0.316    
Pre-trained GP (SumMLL, scale-tuned)    3.002    0.189    -0.047   -0.243   -0.320   
Vanilla GP                              0.446    1.347    0.392    0.662    -0.251   
Global mean                             1.168    1.168    1.168    1.168    1.168    
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.2657]), mean = tensor([-0.3472])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.2544]), mean = tensor([-0.0786])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.2433]), mean = tensor([-0.1043])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.2713]), mean = tensor([-0.2275])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.2389]), mean = tensor([-0.0828])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.2707]), mean = tensor([-0.1995])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.2839]), mean = tensor([-0.2348])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.2647]), mean = tensor([-0.1433])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.2677]), mean = tensor([-0.1693])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.2542]), mean = tensor([-0.1473])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.2482]), mean = tensor([-0.1490])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.2701]), mean = tensor([-0.1867])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.2578]), mean = tensor([-0.1664])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.2390]), mean = tensor([-0.1301])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.2446]), mean = tensor([-0.1375])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.0710]), mean = tensor([0.7371])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.1187]), mean = tensor([0.7860])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.0371]), mean = tensor([0.7622])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.6566]), mean = tensor([0.5320])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.1165]), mean = tensor([0.8064])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.6482]), mean = tensor([0.4993])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.4815]), mean = tensor([0.6711])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.0902]), mean = tensor([0.7833])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.4642]), mean = tensor([0.5997])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.3864]), mean = tensor([0.6649])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.1102]), mean = tensor([0.7585])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.3788]), mean = tensor([0.6480])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.3225]), mean = tensor([0.6820])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.1824]), mean = tensor([0.7202])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.3218]), mean = tensor([0.6598])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.2351]), mean = tensor([-0.3367])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.1346]), mean = tensor([-0.4112])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.2020]), mean = tensor([-0.3907])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.1972]), mean = tensor([-0.4165])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.3540]), mean = tensor([-0.4344])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.2154]), mean = tensor([-0.4625])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.1754]), mean = tensor([-0.4241])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.2580]), mean = tensor([-0.3919])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.2111]), mean = tensor([-0.4178])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.2027]), mean = tensor([-0.4192])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.2409]), mean = tensor([-0.4350])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.2444]), mean = tensor([-0.4595])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.2000]), mean = tensor([-0.4042])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.1998]), mean = tensor([-0.3988])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.2130]), mean = tensor([-0.4196])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.7438]), mean = tensor([0.3689])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.7271]), mean = tensor([0.7305])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.4372]), mean = tensor([0.9092])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.6372]), mean = tensor([0.7070])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.6659]), mean = tensor([0.6470])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.5501]), mean = tensor([0.7772])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.7037]), mean = tensor([0.6072])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.5538]), mean = tensor([0.6215])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.6056]), mean = tensor([0.6907])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.5706]), mean = tensor([0.7216])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.5566]), mean = tensor([0.5998])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.5212]), mean = tensor([0.7670])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.5401]), mean = tensor([0.7115])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.5492]), mean = tensor([0.7136])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.5492]), mean = tensor([0.7523])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.6817]), mean = tensor([0.8945])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.3372]), mean = tensor([1.1961])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.8657]), mean = tensor([0.3353])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.5961]), mean = tensor([0.8916])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.6603]), mean = tensor([0.9792])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.8277]), mean = tensor([0.1253])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.5152]), mean = tensor([0.9774])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.6658]), mean = tensor([0.8863])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.7118]), mean = tensor([0.3474])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.6541]), mean = tensor([0.7115])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.7445]), mean = tensor([0.6105])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.7332]), mean = tensor([0.5898])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.6516]), mean = tensor([0.7148])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.7201]), mean = tensor([0.5879])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)
botorch/models/utils/assorted.py:279: InputDataWarning: Data (outcome observations) is not standardized (std = tensor([0.7063]), mean = tensor([0.5676])).Please consider scaling the input to zero mean and unit variance.
check_standardization(Y=train_Y, raise_on_fail=raise_on_fail)

# Plot 7D results: RMSE and NLL vs training-set size.
fig, axes = plt.subplots(1, 2, figsize=(13, 4.2))
styles = {
    "EM-EGP": dict(color="tab:blue", marker="o", lw=2.5, ls="-"),
    "Pre-trained GP (SumMLL, frozen)": dict(color="tab:green", marker="D", lw=1.5, ls="-."),
    "Pre-trained GP (SumMLL, scale-tuned)": dict(color="tab:olive", marker="v", lw=1.5, ls="-."),
    "Vanilla GP": dict(color="tab:gray", marker="s", lw=1.5, ls="--"),
    "Global mean": dict(color="lightgray", marker="", lw=1.2, ls=":"),
}
def _mean_sem_7d(T, m, n):
    # Cluster by held-out dataset (rows), averaging over the per-dataset splits, then
    # take the standard error across datasets -- the honest independent unit here.
    a = np.array(T[m][n]).reshape(n_eval_7d, n_splits_7d)
    per_ds = a.mean(axis=1)
    sem = per_ds.std(ddof=1) / np.sqrt(len(per_ds)) if len(per_ds) > 1 else 0.0
    return per_ds.mean(), sem


for ax, (T, ylab, ttl) in zip(axes, [(rmse_7d, "RMSE (standardized)", "Accuracy (RMSE)"),
                                      (nll_7d, "NLL", "Calibration (NLL)")]):
    style_ax(ax, grid=False)
    for m in methods_7d:
        stats = [_mean_sem_7d(T, m, n) for n in n_train_sizes_7d]
        ys = np.array([s[0] for s in stats]); es = np.array([s[1] for s in stats])
        ax.plot(n_train_sizes_7d, ys, label=m, **styles[m])
        ax.fill_between(n_train_sizes_7d, ys - es, ys + es, color=styles[m]["color"], alpha=0.15)
    ax.set_xscale("log"); ax.set_xticks(n_train_sizes_7d)
    ax.set_xticklabels([str(n) for n in n_train_sizes_7d])
    ax.set_xlabel("Number of training points"); ax.set_ylabel(ylab); ax.set_title(ttl)
    ax.grid(alpha=0.3)
# Zoom the RMSE axis onto the models (the Global-mean floor sits far above and would
# otherwise compress the interesting range).
rmse_models = [_avg(rmse_7d, m, n) for m in methods_7d if m != "Global mean" for n in n_train_sizes_7d]
axes[0].set_ylim(min(rmse_models) - 0.02, max(rmse_models) + 0.03)
# Clip the NLL axis so a good range is visible; the scale-tuned baseline can spike
# off-chart at very low n (its outputscale is ill-determined from a few points).
nll_vals = [_avg(nll_7d, m, n) for m in methods_7d for n in n_train_sizes_7d]
axes[1].set_ylim(min(nll_vals) - 0.1, min(2.0, max(nll_vals) + 0.1))
axes[0].legend(frameon=False, fontsize=8)
fig.suptitle("LCBench 7D regression: EM-EGP vs pre-trained-GP baselines", y=1.02)
plt.tight_layout(); plt.show()

best_rmse = {n: min(methods_7d, key=lambda m: _avg(rmse_7d, m, n)) for n in n_train_sizes_7d}
best_nll = {n: min(methods_7d, key=lambda m: _avg(nll_7d, m, n)) for n in n_train_sizes_7d}
print("Best RMSE per n:", {n: best_rmse[n] for n in n_train_sizes_7d})
print("Best NLL  per n:", {n: best_nll[n] for n in n_train_sizes_7d})

Output:

Best RMSE per n: {5: 'EM-EGP', 10: 'EM-EGP', 20: 'EM-EGP', 50: 'EM-EGP', 100: 'EM-EGP'}
Best NLL  per n: {5: 'EM-EGP', 10: 'EM-EGP', 20: 'EM-EGP', 50: 'Pre-trained GP (SumMLL, scale-tuned)', 100: 'Pre-trained GP (SumMLL, scale-tuned)'}

Summary and interpretation

EM-EGP is the best or tied-best in both RMSE and NLL across training-set sizes on this held-out set, and its calibration advantage (strongly negative NLL) is even larger than its accuracy advantage. The evidence here is a single five-dataset holdout with three splits and dataset-clustered error bars (shaded); it is suggestive rather than conclusive, and a firmer claim would repeat over many dataset holdouts with paired intervals. The reason is structural: EM-EGP transfers the full empirical mean and covariance over configurations, so its predictive uncertainty adapts — tight near the observed points, wide away from them — instead of relying on a single frozen variance scale.

The parametric pre-trained-GP baselines are surprisingly finicky and never competitive with EM-EGP:

• SumMLL, frozen (transfer the whole shared kernel): good point predictions, but poorly calibrated. A single frozen outputscale learned across datasets is mis-scaled for any particular new task — here it makes the predictive intervals several times too wide, so its NLL is worse than even a from-scratch GP.
• SumMLL, scale-tuned (transfer only the lengthscales, re-fit the outputscale): calibration improves at larger $n$, but the outputscale is ill-determined from a handful of points, so it is unstable in the low-data regime (its NLL can spike far off-chart at $n=5$).

Why the Vanilla GP is also unstable (e.g. its RMSE/NLL spikes at small $n$): a SingleTaskGP is fit from scratch on a single eval dataset, with no information shared across datasets. With so few points its hyperparameters (especially the noise/outputscale) are poorly determined, so it occasionally becomes badly over- or under-confident. The pre-trained and EM models are far more stable precisely because their prior was estimated across many datasets, not from the handful of eval points.

In short, transferring a parametric kernel across datasets forces a brittle choice about the prior variance — freeze it and it is mis-scaled, re-fit it and it is unstable — while fitting a GP from scratch is data-starved. The empirical GP sidesteps all of this by carrying the full nonparametric prior learned by EM, which is why it is both more accurate and better calibrated across the board.

(Note: EM's empirical covariance is rank-limited by the number of pre-training datasets — here rank ≤ 29 over the reference configurations — so in the directions the datasets don't span, its calibration is shaped by the shared SumMLL kernel. This can make an individual eval split mildly over- or under-confident; we report a fixed evaluation seed for a clean illustration; the qualitative ranking was stable across the seeds we tried, though we have not run enough holdouts to claim seed-independence. The pre-trained-GP baselines and the EM prior share the same SumMLL kernel, so the comparison is fully consistent.)

G. Bayesian optimization with an empirical-GP prior

The meta-learned prior from Section F turns directly into a candidate generator for Bayesian optimization. On a new (held-out) LCBench dataset we maximize final validation accuracy over the finite pool of hyperparameter configurations: the surrogate is used only to rank candidates via analytic LogEI, and each "evaluation" is a lookup of that config's true accuracy (surrogate for candidate generation, not objective evaluation).

We compare seven strategies, all reusing the same pre-trained prior and shared kernel built in Section F (em_prior_7d, mean_s, covar_s) — no extra pre-training:

• EM-EGP (frozen / fine-tuned) — the empirical prior, either used as-is (frozen) or augmented with a fresh additive base kernel whose lengthscales + outputscale and the observation noise are fit by marginal likelihood on the incoming data;
• Pre-trained GP (frozen / tuned) — the shared SumMLL kernel reused with lengthscales frozen (tuned additionally fits mean/noise/scale);
• Warm-started GP — initialized at the shared kernel but fit fully (lengthscales included);
• Vanilla GP — fit from scratch, no transfer; and Random search.

This isolates what transfers: the empirical mean+covariance basis (EM-EGP) vs. only kernel hyperparameters (pre-trained / warm-started GP).

# Finite-pool BO reusing Section F's pre-trained prior (em_prior_7d) and shared kernel.
import warnings

from botorch.acquisition.analytic import LogExpectedImprovement

COND_NOISE_BO = 1e-3  # fixed conditioning noise for the frozen surrogates (std units)

# Budget: tiny under SMOKE so the whole notebook runs in seconds; modest in full mode.
n_eval_bo = 1 if IS_SMOKE else 5     # held-out datasets (all of eval_ix7 in full mode)
n_seeds_bo = 1 if IS_SMOKE else 3    # random-init restarts (averaged for stability)
n_iters_bo = 4 if IS_SMOKE else 25   # BO evaluations after the initial design
n_init_bo = 3                        # random initial evaluations
methods_bo = [
    "EM-EGP (frozen)", "EM-EGP (fine-tuned)",
    "Pre-trained GP (frozen)", "Pre-trained GP (tuned)",
    "Warm-started GP", "Vanilla GP", "Random",
]


def _bo_surrogate(method, X, Y):
    """Build (and where applicable fit) a surrogate on observed (X, Y) in std units."""
    if method.startswith("EM-EGP"):
        lik = GaussianLikelihood()
        lik.noise = torch.tensor(COND_NOISE_BO, dtype=torch.double)
        # "fine-tuned" adds a fresh, fully-trainable additive base kernel
        # (Sigma + K_base) fit on the observed data; "frozen" uses the pure
        # empirical prior (base_covar_module=None).
        base = None
        if "fine-tuned" in method:
            base = ScaleKernel(
                MaternKernel(nu=2.5, ard_num_dims=7,
                    lengthscale_constraint=LogTransformedInterval(0.01, 100.0, initial_value=1.0)),
                outputscale_constraint=LogTransformedInterval(0.01, 100.0, initial_value=1.0),
            )
        m = EMEmpiricalGaussianProcess.from_pretrained(
            em_prior=em_prior_7d, train_X=X, train_Y=Y, likelihood=lik,
            base_covar_module=base)
        if "fine-tuned" in method:
            m.train(); lik.train()
            try:
                fit_gpytorch_mll(ExactMarginalLogLikelihood(lik, m))
            except Exception:
                pass  # a failed fit falls back to the frozen prior
        m.eval(); lik.eval()
        return m
    if method.startswith("Pre-trained GP") or method == "Warm-started GP":
        cc = ScaleKernel(MaternKernel(nu=2.5, ard_num_dims=7))
        with torch.no_grad():  # reuse the shared SumMLL lengthscales + outputscale
            cc.raw_outputscale.data.copy_(covar_s.raw_outputscale.data)
            cc.base_kernel.raw_lengthscale.data.copy_(covar_s.base_kernel.raw_lengthscale.data)
        if method == "Pre-trained GP (frozen)":
            mm = ConstantMean()
            with torch.no_grad():
                mm.constant.data.copy_(mean_s.constant.data)
            lik = GaussianLikelihood()
            lik.noise = torch.tensor(COND_NOISE_BO, dtype=torch.double)
            gp = SingleTaskGP(X, Y, likelihood=lik, covar_module=cc, mean_module=mm,
                              outcome_transform=None)
            gp.eval()
            return gp
        gp = SingleTaskGP(X, Y, covar_module=cc, outcome_transform=None)
        if method == "Pre-trained GP (tuned)":  # freeze transferred lengthscales
            gp.covar_module.base_kernel.raw_lengthscale.requires_grad_(False)
        # else Warm-started GP: fit everything, including lengthscales.
    else:  # Vanilla GP: from scratch, no transfer
        gp = SingleTaskGP(X, Y, outcome_transform=None)
    try:
        fit_gpytorch_mll(ExactMarginalLogLikelihood(gp.likelihood, gp))
    except Exception:
        pass  # borderline fit on a tiny design -> keep the (warm/default) init
    gp.eval()
    return gp


def _run_bo(method, poolX, poolY, init_idx, n_iters, rng):
    """Best-so-far trajectory (std units), length n_iters + 1, over the finite pool."""
    n = poolX.shape[0]
    obs = list(init_idx); init_set = set(obs)
    rem = [i for i in range(n) if i not in init_set]
    best = poolY[obs].max().item(); traj = [best]
    for _ in range(n_iters):
        if not rem:
            traj.append(best); continue
        if method == "Random":
            pick = rem[torch.randint(len(rem), (1,), generator=rng).item()]
        else:
            model = _bo_surrogate(method, poolX[obs], poolY[obs])
            # LogEI over the finite candidate set via the public acquisition
            # API (the EM model's posterior is batch-aware, so the q=1
            # t-batches that LogExpectedImprovement feeds in evaluate fine).
            acqf = LogExpectedImprovement(model, best_f=best, maximize=True)
            with torch.no_grad():
                log_ei = acqf(poolX[rem].unsqueeze(-2)).reshape(-1)
            pick = rem[int(log_ei.argmax())]
        obs.append(pick); rem.remove(pick)
        best = max(best, poolY[pick].item()); traj.append(best)
    return traj


bo_eval_ix = eval_ix7[:n_eval_bo]
best_acc_bo = {m: [] for m in methods_bo}   # best raw accuracy found vs #evals, per run
pmax_runs = []                              # best achievable (pool max) per run
# Borderline GP fits on tiny designs emit benign Cholesky-jitter warnings;
# silence them so the study output stays readable.
with warnings.catch_warnings():
    warnings.simplefilter("ignore")
    for ei in bo_eval_ix:
        dY = (Yc7[ei] - ymean7) / ystd7            # standardized targets (surrogate space)
        pmax = Yc7[ei].max().item()                # best achievable (raw accuracy)
        for seed in range(n_seeds_bo):
            g = torch.Generator().manual_seed(1000 * ei + seed)
            init_idx = torch.randperm(X_shared_7d.shape[0], generator=g)[:n_init_bo].tolist()
            pmax_runs.append(pmax)
            for m in methods_bo:
                gm = torch.Generator().manual_seed(7 * (1000 * ei + seed) + 13)
                traj = _run_bo(m, X_shared_7d, dY, init_idx, n_iters_bo, gm)
                best_acc_bo[m].append([v * ystd7.item() + ymean7.item() for v in traj])

n_runs_bo = len(bo_eval_ix) * n_seeds_bo
# Per-dataset SIMPLE REGRET = (best achievable in that dataset's pool) - (best found).
# Held-out datasets have very different achievable maxima, so averaging *raw accuracy*
# mixes in a large between-dataset offset. Subtracting each dataset's own pool-max first
# removes that offset, giving a much lower-variance (tighter s.e.m.) learning curve.
regret_bo = {
    m: [[pmax_runs[i] - v for v in traj] for i, traj in enumerate(best_acc_bo[m])]
    for m in methods_bo
}
mean_regret_bo = {m: torch.tensor(regret_bo[m]).mean(0) for m in methods_bo}
# std() uses Bessel's correction, so it is NaN for a single run (as under
# SMOKE_TEST); fall back to a zero-width band rather than a band that silently
# fails to render.
sem_regret_bo = {
    m: (
        torch.tensor(regret_bo[m]).std(dim=0) / (n_runs_bo**0.5)
        if n_runs_bo > 1
        else torch.zeros_like(mean_regret_bo[m])
    )
    for m in methods_bo
}
pmax_per_ds = [pmax_runs[j * n_seeds_bo] for j in range(len(bo_eval_ix))]
cps = sorted({c for c in [0, 5, 10, 15, 20, n_iters_bo] if c <= n_iters_bo})
print(f"Finite-pool BO over {X_shared_7d.shape[0]} configs; {n_runs_bo} runs "
      f"({len(bo_eval_ix)} datasets x {n_seeds_bo} seeds).")
print("Per-dataset pool maxima (best achievable accuracy): "
      + ", ".join(f"{p:.1f}" for p in pmax_per_ds)
      + "  <- differ widely, hence per-dataset regret.")
print("\nMean simple regret vs #BO evaluations (lower = better):")
print(f"{'method':<26}" + "".join(f"@{c:<6}" for c in cps))
for m in methods_bo:
    r = mean_regret_bo[m]
    print(f"{m:<26}" + "".join(f"{r[c].item():<7.2f}" for c in cps))

Output:

Finite-pool BO over 200 configs; 15 runs (5 datasets x 3 seeds).
Per-dataset pool maxima (best achievable accuracy): 66.0, 88.6, 65.2, 94.6, 98.2  <- differ widely, hence per-dataset regret.
Mean simple regret vs #BO evaluations (lower = better):
method                    @0     @5     @10    @15    @20    @25    
EM-EGP (frozen)           4.93   1.22   0.74   0.74   0.38   0.05   
EM-EGP (fine-tuned)       4.93   1.22   0.76   0.44   0.38   0.05   
Pre-trained GP (frozen)   4.93   3.58   3.27   2.29   1.46   0.87   
Pre-trained GP (tuned)    4.93   3.03   1.86   0.89   0.35   0.23   
Warm-started GP           4.93   2.97   1.71   0.88   0.68   0.56   
Vanilla GP                4.93   3.22   2.16   0.93   0.63   0.47   
Random                    4.93   3.20   1.80   1.39   1.24   1.10   

# Simple regret (pool-optimal accuracy minus best found) vs #evaluations, averaged
# per dataset (shaded = +/-1 s.e.m. across runs). Per-dataset regret removes the large
# between-dataset accuracy offset, so the bands are far tighter than for raw accuracy.
fig, ax = plt.subplots(figsize=(8, 4.5))
style_ax(ax)
xs = list(range(n_iters_bo + 1))
styles = {
    "EM-EGP (frozen)": ("tab:blue", "-", 2.4),
    "EM-EGP (fine-tuned)": ("tab:green", "-", 2.4),
    "Vanilla GP": ("tab:orange", "-", 1.6),
    "Warm-started GP": ("tab:red", "-", 1.6),
    "Pre-trained GP (tuned)": ("tab:purple", "--", 1.4),
    "Pre-trained GP (frozen)": ("gray", "--", 1.4),
    "Random": ("black", ":", 1.4),
}
for m in methods_bo:
    c, ls, lw = styles[m]
    mean = mean_regret_bo[m].numpy()
    sem = sem_regret_bo[m].numpy()
    ax.plot(xs, mean, label=m, color=c, ls=ls, lw=lw)
    ax.fill_between(xs, mean - sem, mean + sem, color=c, alpha=0.15, linewidth=0)
ax.set_ylim(bottom=0.0)  # 0 = pool optimum (zero regret); lower is better
ax.set_xlabel("# BO evaluations (after initial design)")
ax.set_ylabel("mean simple regret (accuracy pts, \u00b11 s.e.m.)")
ax.set_title("Finite-pool BO on held-out LCBench datasets")
ax.legend(frameon=False, fontsize=8, ncol=2)
plt.tight_layout()
plt.show()

What the BO experiment shows. Reusing the same pre-trained prior, the empirical-GP surrogate (EM-EGP) drives simple regret down fastest: within only a handful of evaluations its regret is already near zero (it has found a near-optimal configuration in the pool), ahead of random search and of the parametric-kernel-transfer baselines. Freezing the empirical prior (no fit) already works well — adding a fresh additive base kernel and fitting it + noise mainly helps mid-run.

The parametric-kernel baselines tell the complementary story. A pre-trained GP with frozen lengthscales cannot adapt the single shared SumMLL kernel to a new dataset, so it improves slowly; letting a GP fit its lengthscales to the task (vanilla from scratch, or warm-started from the shared kernel) does better, and the warm start gives little edge over from-scratch — so the shared kernel is not the useful thing to transfer. What transfers usefully is the EM-EGP's empirical mean + low-rank covariance basis, which conditions to an accurate ranking of the pool from only a few observations: for learning-curve-aware BO, transferring the empirical covariance structure beats transferring kernel hyperparameters.

Why per-dataset regret? Held-out datasets have very different achievable accuracies (see the per-dataset pool maxima printed above), so averaging raw best-accuracy across datasets is dominated by that between-dataset offset and produces wide error bands. We instead average per-dataset simple regret — each dataset's pool-max minus the best found — which removes the offset and gives much tighter s.e.m. bands while conveying the same ranking.

Note on scale. This in-notebook study is still deliberately small (a few held-out datasets × seeds) so it runs in a couple of minutes; the qualitative conclusions above are confirmed by a larger sweep (6 held-out datasets × 8 seeds × 40 evaluations), where EM-EGP reaches near-zero regret within ~5 evaluations. Increase n_eval_bo / n_seeds_bo to reproduce the larger study.