What it actually does
Instead of fitting one function, put a probability distribution over all the functions that could have produced the data, and let the observations narrow it. The kernel is what makes this concrete: it states how correlated the function values at two inputs should be, which is the same as stating how smooth, how wiggly, and how periodic you believe the function is. Conditioning that prior on the data gives a posterior that is again Gaussian, so the answer is closed-form — a mean and a variance at every input. The variance is the part that earns its keep: it is small where data is dense, and it widens automatically wherever you leave the data behind, which no point-estimate regressor will tell you.
- Function values at nearby inputs are correlated in the way the kernel says — the kernel is the model, and choosing RBF asserts infinite differentiability whether or not that was intended
- Observation noise is Gaussian and, in the standard formulation, the same everywhere; heteroscedastic noise needs an explicit extension
- The prior mean is correct where there is no data, since that is exactly what the posterior returns there
- The process is stationary unless the kernel says otherwise — a stationary kernel cannot represent a function that is smooth in one region and rough in another
- The sample is small enough for an n-by-n factorization, which is a constraint on the model rather than on the hardware
Likelihood
How probable is the observed data under this prior over functions, with the latent function integrated out entirely rather than estimated. The two terms are a fit-versus-complexity trade made explicit: the first rewards explaining the data, the second penalizes a prior loose enough to explain any data. That balance is the reason a GP can select its own hyperparameters without a validation set — the marginal likelihood already contains an Occam term. It is not a free lunch, though: the surface is non-convex, and with many hyperparameters it can and does overfit itself.
- the n-by-n kernel matrix, K_ij = k(x_i, x_j) — the prior covariance between function values
- observation noise, added to the diagonal; also the numerical regularizer that keeps K_y invertible
- kernel hyperparameters — lengthscale, amplitude, period, noise — fitted by maximizing this quantity
- the complexity penalty: a kernel flexible enough to explain anything has a large determinant, and pays for it here
How it is minimized
Closed-form posterior via Cholesky factorization; gradient ascent on the log marginal likelihood for the kernel hyperparameters
- covariance between the query point and every training point — how much each observation is allowed to say about it
- prior variance at the query point, before any data is taken into account
- posterior mean: a weighted average of the training targets, with kernel-derived weights
- posterior variance — the prior variance minus whatever the data explained, so it can only shrink
There is nothing iterative in the prediction: conditioning a joint Gaussian gives another Gaussian, so both the mean and the variance are exact algebra. The whole implementation question is how the inverse is handled, and the answer is that it never is. Factorize K_y once with a Cholesky and every subsequent quantity follows from triangular solves — the posterior mean, the posterior variance, and the log-determinant, which is just twice the sum of the log diagonal. Forming the inverse explicitly is both slower and numerically worse, and on a near-singular kernel matrix it is the difference between a usable model and noise. The iteration that does exist is over hyperparameters: gradient ascent on the marginal likelihood, which is non-convex and needs restarts.
kernel family
The modelling assumption, and the only one that matters. RBF asserts an infinitely differentiable function; Matérn 3/2 or 5/2 asserts something rougher and is usually the more honest default for real data
lengthscale
How far the function has to move before values decorrelate — short means wiggly, long means smooth. Learned from the marginal likelihood, per dimension if you want automatic relevance determination (log-scale search around the median pairwise distance)
signal variance
Amplitude of the function; sets the prior spread and therefore how far the posterior can move from the mean
noise variance (sigma_n^2)
Observation noise. Doing double duty as the numerical jitter that keeps the kernel matrix positive definite (1e-6 to the target variance; never exactly zero)
mean function
Usually zero, which matters more than it sounds: predictions revert to it away from data, so a zero-mean GP forecasts zero once extrapolating
The posterior is exact — no convergence question at all, given hyperparameters. Everything that can go wrong is numerical or is about the hyperparameters. Numerically, the kernel matrix becomes ill-conditioned as soon as two training inputs are close together, because their rows are nearly identical; the Cholesky then fails outright, and the standard remedy is jitter on the diagonal, which is a small deliberate lie about the noise level. On hyperparameters, the marginal likelihood is non-convex with real local optima that correspond to genuinely different interpretations — one explaining the data as a smooth trend with high noise, another as a wiggly function with low noise — and the optimizer picks whichever basin it started in. Restarts are not optional.
Training: O(n^3) for the Cholesky, O(n^2) memory for the kernel matrix. That is the defining fact — worse than the SVM next door — and it puts exact GP regression out of reach somewhere around 10,000 points. Prediction, with the factorization cached: O(n) per query for the mean and O(n^2) for the variance, which is why the variance is the expensive half. Sparse and inducing-point approximations reduce training to O(n·m^2) for m inducing points and are what any GP at scale is actually doing.
How, where, and why it is used
Time-Series Forecasting
ViableTreat time as the input and compose the kernel to match the structure: a long-lengthscale Matérn for trend, a periodic kernel for seasonality, a product of periodic and slowly-varying for seasonality that changes shape, plus a white-noise term. The composition IS the model specification — it is written down rather than searched for — and the posterior then supplies a mean forecast with genuinely calibrated intervals.
- Small-sample forecasting where a calibrated interval matters more than the point estimate — clinical, environmental, and sensor-calibration series
- Irregularly sampled and gappy series, which a GP handles natively because it never assumed a fixed sampling grid
- Interpolation and smoothing between observations, where the posterior mean is the optimal smoother under the assumed kernel
- Series where the kernel structure is a deliverable in itself: showing that a periodic component with a learned period explains the data is an argument, not just a fit
Its advantage over every other forecaster here is the interval — not a bootstrap approximation but the actual posterior, and one that widens correctly as the horizon extends. It also handles irregular sampling without interpolation and encodes structure explicitly. The two reasons not to use it are hard ones: O(n^3) rules out long series outright, and the zero-mean prior means the forecast reverts to the mean beyond the data, so a GP with an RBF kernel will not extrapolate a trend — it will flatten. Encoding the trend in the kernel or the mean function is mandatory, not optional.
- Compose the kernel deliberately: trend times periodic plus noise, rather than an RBF and hope
- Prefer Matérn over RBF unless the process really is infinitely smooth — the RBF prior is far stronger than most people intend it to be
- Standardize the target, since the signal variance and noise are both learned on its scale
- Use an explicit mean function or a linear kernel component when the series trends, or the forecast reverts to zero
Rolling-origin backtesting scored on a proper scoring rule that rewards the whole distribution — CRPS or the negative log predictive density — not RMSE, which ignores the model’s main output. Check interval coverage explicitly: an 80% interval should contain about 80% of the held-out points, and a GP whose coverage is wrong has a mis-specified kernel rather than a tuning problem.
- Extrapolating with a stationary kernel and a zero mean, which flattens toward zero and reports growing uncertainty around a forecast nobody believes
- Optimizing hyperparameters on the whole series before splitting — the marginal likelihood is a fit statistic, and it leaks
- Hitting the O(n^3) wall silently by resampling to a finer grid
- A local optimum that explains all structure as noise, producing a flat mean with wide bands that looks like humility and is actually a failed fit
Anomaly Detection
AdaptedFit on clean data, then score a new observation by its standardized residual against the predictive distribution: how many predictive standard deviations away it is. The difference from residual scoring on any other regressor is that the denominator varies with the input — a deviation in a well-sampled region is significant, and the same deviation where the model has little data is not.
- Sensor calibration and drift monitoring, where the expected value is a smooth function of operating conditions
- Monitoring a metric against a smooth, seasonally-structured baseline with irregular observation times
- Screening for inputs outside the training distribution — the predictive variance alone is a usable novelty score, with no residual required
The input-dependent variance is the entire argument: a fixed threshold on raw residuals over-alarms in sparse regions and under-alarms in dense ones, and the GP gives the right denominator for free. It is also the honest option when data is scarce, since it degrades to wide intervals rather than to confident nonsense. Against it: O(n^3) makes the reference window small, so a GP detector monitors a well-understood low-dimensional process rather than a high-volume stream, where an isolation forest is the correct tool.
- Fit on a confirmed-clean window, and keep it small enough that a refit stays affordable — the cubic cost is the real design constraint here
- Include the operating-condition covariates that legitimately move the baseline, or their effect is scored as anomaly
- Learn the noise level rather than fixing it; the threshold is expressed in units of predictive standard deviation, so an underestimated noise term makes everything an anomaly
Precision@k against confirmed incidents, plus a calibration check on clean data: the standardized residuals should look standard normal, and if they do not, the kernel is wrong and any threshold on them is arbitrary.
- Contamination absorbed into the noise term, which widens the intervals until nothing is anomalous — a quieter failure than the boundary being dragged
- Reading the predictive variance as a confidence in the model rather than in the function value; it is conditional on the kernel being right, and says nothing about the kernel being wrong
- A reference window that grows over time until the refit stops fitting in the maintenance budget
Optimization
Primary fitThis is the surrogate in Bayesian optimization. When each evaluation of an objective is expensive — a training run, a physical experiment, a simulation — fit a GP to the evaluations you have, then use its posterior to choose where to evaluate next: an acquisition function such as expected improvement trades exploiting the region with the best predicted mean against exploring the region with the largest variance. The optimizer is then over the cheap acquisition surface, not the expensive objective.
- Hyperparameter optimization, where a single evaluation is a full training run
- Experimental design and A/B allocation with a limited budget of expensive trials
- Materials, chemical, and process parameter search — the domain where this originated, under the name kriging
- Simulation-based design where each run costs hours and the search space is continuous and low-dimensional
It is the default surrogate for one reason: the acquisition function needs calibrated uncertainty, and a GP supplies it in closed form. A point-estimate model cannot say where it is ignorant, and without that there is no principled way to explore. The conditions are equally specific — the objective must be expensive relative to the surrogate, the dimension low (roughly under 20, since the lengthscale becomes unidentifiable above that), and the budget small enough that O(n^3) stays trivial. When evaluations are cheap, random search or an evolutionary method beats this outright, because the surrogate overhead buys nothing.
- Scale every search dimension to a common range, since a shared lengthscale otherwise encodes the arbitrary units of the search space
- Search log-scale parameters in log space — a learning rate is not smooth in linear coordinates
- Use Matérn 5/2 rather than RBF; the standard recommendation in this domain precisely because real objective surfaces are not infinitely smooth
- Model integer and categorical dimensions explicitly rather than rounding a continuous relaxation, which makes the surrogate wrong in a way it cannot detect
Best-value-found against evaluation count, averaged over repeated runs with different seeds — a single run of a stochastic search says nothing. Compare against random search on the same budget, which is a stronger baseline than most reports admit and wins outright in high dimension.
- Over-exploiting: an acquisition function that collapses onto the incumbent optimum and stops exploring, which is what a badly-fitted noise term causes
- Duplicate or near-duplicate evaluations making the kernel matrix singular, so the surrogate fails numerically at exactly the point it has become confident
- Applying it above roughly twenty dimensions, where the surrogate learns nothing and the whole apparatus becomes an expensive random search
- Ignoring evaluation noise: a deterministic-noise assumption on a stochastic objective produces a surrogate that chases its own measurement error
Control & Operations
ViableTwo roles. As a learned dynamics model inside model-predictive control, where the posterior variance lets the controller be cautious in states it has not visited — the basis of the data-efficient model-based methods. And as a response-surface model in process operations, mapping controllable settings to an outcome so that the setting can be chosen with a stated confidence rather than a point prediction.
- Data-efficient model-based control, where a GP dynamics model needs far fewer trials than a neural one
- Process and yield optimization in manufacturing, where each trial run is expensive and uncertainty must be reported
- Robotics and physical systems with low state dimension and costly real-world interaction
Sample efficiency and honest uncertainty are exactly the properties this domain needs: physical trials cost money, and a controller that does not know what it does not know will confidently drive a system somewhere it has never been. The GP supplies both. The limits are as usual structural — state dimension must be small, the cubic cost bounds how much interaction history can be retained, and the smoothness assumption fails at contacts and discontinuities, which is precisely where physical systems are hardest.
- Model the state delta rather than the next state, so the prior mean of zero corresponds to "nothing changes" instead of "the system returns to the origin"
- Scale each state and action dimension, since a shared lengthscale over mixed units is meaningless
- Keep a bounded, subsampled history or move to a sparse GP; an unbounded interaction buffer meets the cubic wall quickly
Multi-step rollout error rather than one-step prediction error, since compounding is what actually breaks a learned dynamics model. Check that the predictive variance grows along the rollout — a model whose uncertainty stays flat over a long horizon is mis-specified.
- Smoothness assumptions failing at contact events and saturation limits, exactly where the controller needs accuracy most
- Compounding error over a horizon while the reported uncertainty stays optimistic because each step was conditioned on a point estimate
- The training buffer growing until the refit no longer fits inside the control loop
What it costs to run
O(n^3) time and O(n^2) memory, which at n = 10,000 is an 800 MB kernel matrix and a factorization measured in minutes. Hyperparameter fitting repeats that factorization at every gradient step, so the practical cost is a multiple of it. This is the number that decides whether a GP is a candidate.
With the Cholesky factor and the weight vector cached at fit time, the posterior mean is O(n) per query — a dot product against the training set — and the variance is O(n^2), which makes uncertainty the expensive half. The model also carries the entire training set, as the kernel methods next door do.
Infrequent, and forced to be: a refit is cubic and there is no exact incremental update, though a rank-one Cholesky update handles single-point additions cheaply and is what online GP implementations use. Where data arrives continuously, a sparse GP with fixed inducing points is the version that can keep up.
- Track empirical coverage of the predictive intervals — an 80% interval that covers 60% of new points means the kernel is wrong, and it surfaces before the mean error moves
- Monitor the condition number of the kernel matrix across refits; a rising one means training inputs are clustering and the factorization is about to need more jitter
- Watch the fitted noise variance over time, since a rising estimate is usually the model absorbing structure it can no longer explain
- Alert on queries whose predictive variance approaches the prior variance — those are inputs the model has no information about, and its mean there is the prior, not a prediction
- Jitter must be persisted with the model and applied identically at scoring time; a model factorized with jitter and scored without it is a different model
- The posterior mean reverts to the prior mean away from the data, so a zero-mean GP silently predicts zero in unfamiliar regions — the variance is the only thing that says so
- Never form the inverse: an explicit K^-1 is slower and loses roughly half the available precision relative to solving with the Cholesky factor
- Cache the factorization and the weight vector at fit time; re-solving per query turns an O(n) prediction into an O(n^3) one and is a common accidental regression
- The whole training set ships inside the model, which makes deletion requests and data licensing model problems, the same way they are for kNN and SVM
Pros and cons, in context
Calibrated uncertainty at every input, in closed form
This is the reason to choose it, and it is what makes Bayesian optimization and cautious control possible at all. Worth nothing if only the point prediction is consumed, where cheaper models match it.
Works well with very little data
A dozen points and a well-chosen kernel produce a usable model with honest error bars. Decisive when each observation is expensive; irrelevant when data is abundant, which is also when the cubic cost forbids it.
Assumptions are stated in the kernel rather than buried
Smoothness, periodicity, and additive structure are written down and can be argued about — and the marginal likelihood scores them. That is a genuine scientific advantage over a model whose inductive bias is implicit in its architecture.
Hyperparameters can be fitted without a validation set
The marginal likelihood contains its own complexity penalty, so there is no fold-splitting step. The caveat is real: it is non-convex, and with many kernel hyperparameters it will overfit itself.
O(n^3) time and O(n^2) memory
A hard ceiling near 10,000 points, and the reason nearly every GP in production is a sparse approximation. It is the first thing to check and it settles most arguments before accuracy is discussed.
Does not extrapolate — the posterior reverts to the prior mean
With a stationary kernel and a zero mean, the forecast decays to zero beyond the data. Correct behaviour for a model that knows it is ignorant, and a trap for anyone who wanted a trend extended.
Everything depends on the kernel, and the default is usually wrong
The RBF kernel asserts an infinitely differentiable function, which almost no real process is; a Matérn is the more defensible default. Choosing the kernel is modelling work, not configuration, and there is no way to avoid doing it.
Standard formulation assumes Gaussian, homoscedastic noise
Non-Gaussian likelihoods — counts, classification, censored data — lose the closed form and need approximate inference, at which point most of the simplicity that motivated the model is gone.
Intuitive, idiomatic, optimized
The same algorithm in three languages, each progressing through Kent Beck's three stages. Every optimization names the lever it pulls and the trade-off it accepts — cost profiles are illustrative, not measured benchmarks.
The math, transcribed literally. No dependencies.
"""Gaussian process regression - the posterior equations, transcribed.
mu_* = k_*^T (K + sigma_n^2 I)^-1 y
sigma_*^2 = k_** - k_*^T (K + sigma_n^2 I)^-1 k_*
The inverse is never formed. K_y is symmetric positive definite, so a Cholesky
factorization plus two triangular solves gives the same answer at half the
cost and far better conditioning - and the log-determinant, needed for the
marginal likelihood, is just twice the sum of the log diagonal.
"""
import math
def rbf(a, b, lengthscale, amplitude):
"""k(a, b) = amplitude * exp(-||a - b||^2 / (2 * lengthscale^2))"""
squared = 0.0
for j in range(len(a)):
difference = a[j] - b[j]
squared += difference * difference
return amplitude * math.exp(-0.5 * squared / (lengthscale * lengthscale))
def cholesky(matrix):
"""Lower-triangular L with L L^T = matrix. Fails loudly if not SPD."""
n = len(matrix)
lower = [[0.0] * n for _ in range(n)]
for i in range(n):
for j in range(i + 1):
total = matrix[i][j]
for k in range(j):
total -= lower[i][k] * lower[j][k]
if i == j:
if total <= 0.0:
raise ValueError("kernel matrix is not positive definite; add jitter")
lower[i][j] = math.sqrt(total)
else:
lower[i][j] = total / lower[j][j]
return lower
def forward_substitution(lower, b):
"""Solve L z = b."""
n = len(b)
z = [0.0] * n
for i in range(n):
total = b[i]
for k in range(i):
total -= lower[i][k] * z[k]
z[i] = total / lower[i][i]
return z
def back_substitution(lower, z):
"""Solve L^T w = z."""
n = len(z)
w = [0.0] * n
for step in range(n):
i = n - 1 - step
total = z[i]
for k in range(i + 1, n):
total -= lower[k][i] * w[k]
w[i] = total / lower[i][i]
return w
def fit_predict(X_train, y_train, X_test, lengthscale=1.0, amplitude=1.0, noise=1e-4):
n = len(X_train)
# K_y = K + sigma_n^2 I. The noise term is also what keeps this invertible
# when two training inputs are close together.
k_y = [[0.0] * n for _ in range(n)]
for i in range(n):
for j in range(n):
k_y[i][j] = rbf(X_train[i], X_train[j], lengthscale, amplitude)
k_y[i][i] += noise
lower = cholesky(k_y)
# alpha = K_y^-1 y, computed once and reused for every query point.
alpha = back_substitution(lower, forward_substitution(lower, y_train))
means = []
variances = []
for query in X_test:
k_star = [rbf(query, X_train[i], lengthscale, amplitude) for i in range(n)]
mean = 0.0
for i in range(n):
mean += k_star[i] * alpha[i]
# v = L^-1 k_*, so k_*^T K_y^-1 k_* = ||v||^2 - one solve, no inverse.
v = forward_substitution(lower, k_star)
explained = 0.0
for i in range(n):
explained += v[i] * v[i]
means.append(mean)
variances.append(amplitude - explained)
return means, variancesO(n^3) for the factorization and O(n^2) per query for the variance, in interpreter loops over lists of lists.
How this language actually wants it written.
"""Gaussian process regression - typed, factorization cached, kernel injected."""
from collections.abc import Callable
from dataclasses import dataclass
import numpy as np
from numpy.typing import NDArray
from scipy.linalg import cho_factor, cho_solve, solve_triangular
Vector = NDArray[np.float64]
Matrix = NDArray[np.float64]
Kernel = Callable[[Matrix, Matrix], Matrix]
def matern52(lengthscale: float, amplitude: float = 1.0) -> Kernel:
"""Matern 5/2 - twice differentiable, and the honest default.
The RBF kernel asserts an infinitely differentiable function, which almost
no real process is. Matern 5/2 assumes far less and is the standard choice
in Bayesian optimization for exactly that reason.
"""
def kernel(a: Matrix, b: Matrix) -> Matrix:
squared = (
np.einsum("ij,ij->i", a, a)[:, None]
- 2.0 * (a @ b.T)
+ np.einsum("ij,ij->i", b, b)[None, :]
)
distance = np.sqrt(np.maximum(squared, 0.0)) / lengthscale
scaled = np.sqrt(5.0) * distance
return amplitude * (1.0 + scaled + 5.0 / 3.0 * distance**2) * np.exp(-scaled)
return kernel
@dataclass(frozen=True)
class GaussianProcess:
"""A fitted GP. The Cholesky factor and the weights are the model.
Both are computed once at fit time: re-solving per query would turn an O(n)
prediction into an O(n^3) one, which is the most common accidental
regression in a GP implementation.
"""
X_train: Matrix
factor: tuple[Matrix, bool] # cho_factor output
weights: Vector # alpha = K_y^-1 y
kernel: Kernel
prior_variance: float
log_marginal_likelihood: float
def predict(self, X: Matrix, return_std: bool = False) -> tuple[Vector, Vector | None]:
if X.ndim != 2 or X.shape[1] != self.X_train.shape[1]:
raise ValueError(
f"expected (m, {self.X_train.shape[1]}) inputs, got {X.shape}"
)
cross = self.kernel(X, self.X_train) # (m, n)
mean = cross @ self.weights
if not return_std:
return mean, None
# v = L^-1 k_*, so the explained variance is the squared row norm.
lower, _ = self.factor
solved = solve_triangular(lower, cross.T, lower=True)
explained = np.einsum("ij,ij->j", solved, solved)
variance = np.maximum(self.prior_variance - explained, 0.0)
return mean, np.sqrt(variance)
def fit(
X: Matrix,
y: Vector,
kernel: Kernel | None = None,
noise: float = 1e-4,
jitter: float = 1e-10,
) -> GaussianProcess:
"""Condition the prior on the data. Raises ValueError on malformed input."""
if X.ndim != 2:
raise ValueError(f"X must be 2-D, got shape {X.shape}")
if X.shape[0] != y.shape[0]:
raise ValueError(f"X has {X.shape[0]} rows but y has {y.shape[0]}")
if noise <= 0.0:
raise ValueError("noise must be positive; it is also the numerical regularizer")
# None rather than a default kernel instance: a mutable default would be
# shared across every call to this function.
active_kernel = matern52(lengthscale=1.0) if kernel is None else kernel
n = X.shape[0]
k_y = active_kernel(X, X) + (noise + jitter) * np.eye(n)
try:
factor = cho_factor(k_y, lower=True)
except np.linalg.LinAlgError as error:
raise ValueError(
"kernel matrix is not positive definite - duplicate or near-duplicate "
"training inputs; increase jitter"
) from error
weights = cho_solve(factor, y)
# log|K_y| = 2 * sum(log diag(L)) - free, given the factor already computed.
lower, _ = factor
log_det = 2.0 * np.log(np.diag(lower)).sum()
log_marginal = -0.5 * float(y @ weights) - 0.5 * log_det - 0.5 * n * np.log(2.0 * np.pi)
return GaussianProcess(
X_train=X,
factor=factor,
weights=weights,
kernel=active_kernel,
prior_variance=float(active_kernel(X[:1], X[:1])[0, 0]),
log_marginal_likelihood=log_marginal,
)What changed: The structural change is that fitting and predicting become separate: the Cholesky factor and the weight vector are computed once and stored on the model, so a prediction is a dot product rather than a re-factorization — the difference between O(n) and O(n^3) per query, and the most common accidental regression in a GP implementation. The kernel becomes an injected callable rather than a hard-coded RBF, and the default moves to Matern 5/2, because RBF asserts infinite differentiability that almost no real process has. The hand-rolled triangular solves become LAPACK calls, the log marginal likelihood is computed from the factor that already exists, and the non-positive-definite case — which the previous stage raised as a bare arithmetic failure — becomes a specific error naming its actual cause.
O(n^3) once at fit; O(n) per query for the mean and O(n^2) for the variance, in LAPACK.
Faster, and the justification for each change.
"""Gaussian process regression - blocked prediction, fused variance, in place."""
import numpy as np
from numpy.typing import NDArray
from scipy.linalg import cho_factor, cho_solve, solve_triangular
Vector = NDArray[np.float64]
Matrix = NDArray[np.float64]
class BlockedGaussianProcess:
"""Exact GP prediction over large query batches.
The variance is the expensive half - O(n^2) per query - and the naive
implementation computes it one query at a time, re-entering LAPACK for each.
Here the whole block is solved at once: L^-1 K_*^T is a single triangular
solve over an (n, m) right-hand side, and the explained variance is the
squared column norms of the result.
Nothing here makes exact GP regression sub-cubic. Past roughly ten thousand
training points the answer is an inducing-point approximation, not a better
kernel evaluation.
"""
def __init__(self, X: Matrix, y: Vector, lengthscale: float, amplitude: float,
noise: float) -> None:
self._X = np.ascontiguousarray(X, dtype=np.float64)
self._lengthscale = lengthscale
self._amplitude = amplitude
# Squared norms of the training inputs, precomputed once: every kernel
# evaluation for every future query reuses them.
self._train_norms = np.einsum("ij,ij->i", self._X, self._X)
k_y = self._kernel(self._X, self._train_norms)
k_y[np.diag_indices_from(k_y)] += noise
self._factor = cho_factor(k_y, lower=True)
self._weights = cho_solve(self._factor, np.ascontiguousarray(y, dtype=np.float64))
def _kernel(self, B: Matrix, b_norms: Vector) -> Matrix:
"""RBF against the training set, via the squared-norm expansion.
||a - b||^2 = ||a||^2 - 2 a . b + ||b||^2, so the only term touching
both operands is a GEMM.
"""
squared = self._X @ B.T
squared *= -2.0 # in place
squared += self._train_norms[:, None]
squared += b_norms[None, :]
np.maximum(squared, 0.0, out=squared) # cancellation guard
squared *= -0.5 / (self._lengthscale * self._lengthscale)
np.exp(squared, out=squared)
squared *= self._amplitude
return squared # (n, m)
def predict(self, X: Matrix, batch: int = 4096) -> tuple[Vector, Vector]:
queries = np.ascontiguousarray(X, dtype=np.float64)
m = queries.shape[0]
mean = np.empty(m, dtype=np.float64) # allocated once
std = np.empty(m, dtype=np.float64)
lower, _ = self._factor
# Batching bounds the (n, batch) cross-covariance block, which is what
# actually limits how many queries can be scored in one call.
for start in range(0, m, batch):
stop = min(start + batch, m)
block = queries[start:stop]
block_norms = np.einsum("ij,ij->i", block, block)
cross = self._kernel(block, block_norms) # (n, b)
mean[start:stop] = cross.T @ self._weights
# One triangular solve for the whole block, not one per query.
solved = solve_triangular(lower, cross, lower=True, overwrite_b=True)
explained = np.einsum("ij,ij->j", solved, solved)
np.subtract(self._amplitude, explained, out=explained)
np.maximum(explained, 0.0, out=explained)
np.sqrt(explained, out=std[start:stop])
return mean, stdWhat changed: The kernel evaluation stops being a broadcast subtraction and becomes a matrix product, by expanding the squared distance so the only term involving both operands is a GEMM — with the training-set squared norms precomputed once and reused by every future query. The predictive variance, which is the expensive half, stops being a per-query triangular solve and becomes one solve over the whole query block, so LAPACK is entered once per batch rather than once per point. Batching bounds the (n, m) cross-covariance block, and the kernel arithmetic runs entirely through in-place operations because that block is the largest array in the computation. What is deliberately not claimed: none of this changes the cubic fit cost.
Vectorize to NumPy/BLAS (@, np.dot, einsum)
The cross-covariance becomes one GEMM through the squared-norm expansion, and the explained variance becomes an einsum reduction over the solved block with no intermediate.
Trade-off: The expansion is less numerically stable than a direct subtraction — cancellation between large squared norms yields small negative distances, which is why the clip before the exponential is load-bearing rather than defensive.
Batch work to amortize interpreter overhead
One triangular solve serves a whole block of queries instead of one per query, and the (n, batch) cross-covariance bounds peak memory independently of how many points are scored.
Trade-off: The batch size is an untunable-by-default parameter: too small and the LAPACK call is dominated by dispatch, too large and the cross-covariance block no longer fits in cache.
Pre-allocate output arrays and use in-place operations
Mean and standard-deviation outputs are allocated once for the batch, and the kernel is built through in-place multiply, add, and exp rather than six full-size temporaries.
Trade-off: The in-place chain reuses one buffer for the squared distance, then the scaled distance, then the kernel — fast, materially harder to read, and silently wrong if the lines are reordered.
O(n^3) once at fit, then O(n*m) for the means and O(n^2*m) for the variances in BLAS. Illustrative, not a measured benchmark.