Data ScienceMachine Learning
August 12, 2026 | 26 min read
Gaussian process regression addresses one of the most common problems in applied statistics: building a model that captures a non-linear relationship between variables without committing in advance to a specific functional form. A more traditional approach is to assume a particular non-linear form (sinusoidal, exponential, polynomial, and so on), but unless that form is obvious from the outset, model selection becomes a chore. A non-parametric alternative is to define a set of knots and use splines or kernel regression, though knot placement tends to be ad hoc. A third alternative is a Bayesian non-parametric strategy that models the unknown underlying function directly. For this, we use Gaussian process models.
A walkthrough using scikit-learn 1.7, GPflow 2.11 / TensorFlow 2.21, and PyMC 5.20 / PyTensor, the current releases as of this post’s publication. This post was tested against the following package versions. Installing these exact versions before running the code below will avoid most compatibility issues:
%pip install scikit-learn==1.7 gpflow==2.11 tensorflow==2.21 \
tensorflow-probability==0.25 pymc==5.20 arviz==0.20 "matplotlib<3.11"Calling a Bayesian method "non-parametric" is a bit of a misnomer. It does not mean there are *no* parameters, but rather that the number of parameters grows with the data. Bayesian non-parametric models are infinitely parametric.
Gaussian process regression is a natural fit for time series forecasting because the kernel directly encodes assumptions about how a series varies over time, and the resulting forecast comes with a credible interval that widens appropriately the further it extrapolates beyond the observed data. Composite kernels, built by summing components as demonstrated earlier in this post, let a model capture both a smooth long-term trend and a periodic component, such as daily or seasonal cycles, within a single specification.
A surrogate model stands in for an expensive-to-evaluate function, such as the validation loss of a machine learning model across a grid of hyperparameters, and lets an optimization procedure search efficiently without running every candidate configuration. Gaussian process regression is the surrogate model behind most Bayesian hyperparameter optimization tools. The posterior mean guides the search toward promising regions, while the posterior variance drives exploration of regions that remain uncertain, a strategy generally described as balancing exploitation against exploration.
Beyond machine learning tooling, gaussian process regression is widely used wherever a physical process needs to be modeled from a limited number of expensive measurements, including calibrating simulations in engineering, modeling sensor drift, and interpolating spatial data in geostatistics. The same probabilistic regression framework that produces predictive uncertainty for a synthetic dataset in this post applies directly to these settings, since the underlying assumption, that nearby inputs produce correlated outputs, holds broadly across physical and biological systems.
Adopting a set of Gaussians (a multivariate normal vector) provides several advantages. The marginal distribution of any subset of a multivariate normal is itself normal, and the conditional distribution of a subset given the rest is normal too. A Gaussian process (GP) generalizes the multivariate normal to infinite dimension. It is an infinite collection of random variables, any finite subset of which is jointly Gaussian.
Another way to think of an infinite vector is as a function. So a GP can be described as a distribution over functions. Just as a multivariate normal is fully specified by a mean vector and covariance matrix, a GP is fully specified by a mean function and a covariance function.
One specification of a GP might look like this:
m(x) = 0
k(x, x′) = θ₁ exp(−θ₂/2 (x − x′)²)
Here, the covariance function is a squared exponential. Values of x and x′ that are close together produce values of k closer to one, while values that are far apart produce values closer to zero. It may seem odd to adopt the zero function as the mean, since surely something more elaborate would serve better. Most of the learning in a GP happens through the covariance function and its hyperparameters, so little is gained by specifying a complicated mean function. For a finite set of points, the GP becomes a multivariate normal with mean and covariance given by evaluating those functions at the points.
Standard regression methods, including ordinary least squares and kernel ridge regression, return a single best-fit function and, at most, a confidence interval derived from assumptions about the noise. Gaussian process regression instead returns a full posterior distribution over functions consistent with the observed data. This distinction shows up most clearly in the predictions. At any input, gaussian process regression reports both a point estimate and a variance, so predictive uncertainty grows naturally in regions with little or no data and shrinks near observed points.
This makes gaussian process regression a non-parametric regression method. A polynomial or linear model fits a fixed number of coefficients. A Gaussian process instead grows in effective complexity with the amount of data available, governed by a small number of kernel hyperparameters rather than a hand-chosen basis. The tradeoff is computational. Exact GP regression requires inverting an n-by-n covariance matrix, so it scales cubically with the number of training points, which is one reason kernel choice and hyperparameter optimization, covered later in this post, matter so much in practice.
To make the "distribution over functions" idea concrete, let's draw a realization from a GP prior (before seeing any data). We need a covariance function, here the squared exponential, and a way to evaluate it at a set of points to build a covariance matrix.
import numpy as np
def exponential_cov(x, y, params):
return params[0] * np.exp(-0.5 * params[1] * np.subtract.outer(x, y) ** 2)We'll generate realizations sequentially, point by point, using the conditioning property of the multivariate Gaussian. Here is the function that implements the conditional:
def conditional(x_new, x, y, params):
B = exponential_cov(x_new, x, params)
C = exponential_cov(x, x, params)
A = exponential_cov(x_new, x_new, params)
mu = np.linalg.inv(C).dot(B.T).T.dot(y)
sigma = A - B.dot(np.linalg.inv(C).dot(B.T))
return mu.squeeze(), sigma.squeeze()We will start with a Gaussian process prior with hyperparameters σ_0=1, σ_1=10. We will also assume a zero function as the mean, so we can plot a band that represents one standard deviation from the mean.
import matplotlib.pyplot as plt
np.random.seed(42) # reproducible sequential sampling
θ = [1, 10]
σ_0 = exponential_cov(0, 0, θ)
xpts = np.arange(-3, 3, step=0.01)
plt.errorbar(xpts, np.zeros(len(xpts)), yerr=σ_0, capsize=0)
plt.ylim(-3, 3) # force the y-axis to span -3 to 3
plt.title("GP prior: zero mean ± 1 s.d.");
Let's select an arbitrary starting point to sample, say x=1. Since there are no previous points, we can sample from an unconditional Gaussian:
x = [1.]
y = [np.random.normal(scale=σ_0)]
print(y)[0.4967141530112327]We can now update our confidence band, given the point that we just sampled, using the covariance function to generate new point-wise intervals, conditional on the value [x_0, y_0].
σ_1 = exponential_cov(x, x, θ)
def predict(x, data, kernel, params, sigma, t):
k = [kernel(x, y, params) for y in data]
Sinv = np.linalg.inv(sigma)
y_pred = np.dot(k, Sinv).dot(t)
sigma_new = kernel(x, x, params) - np.dot(k, Sinv).dot(k)
return y_pred, sigma_new
x_pred = np.linspace(-3, 3, 1000)
predictions = [predict(i, x, exponential_cov, θ, σ_1, y) for i in x_pred]y_pred, sigmas = np.transpose(predictions)
plt.errorbar(x_pred, y_pred, yerr=sigmas, capsize=0)
plt.plot(x, y, "ro")
plt.title("Posterior after 1 point");
Conditional on this point and the covariance structure, we have constrained the probable location of additional points. Let's sample another:
m, s = conditional([-0.7], x, y, θ)
y2 = np.random.normal(m, s)
print(y2)-0.1382640378102619This point is added to the realization, and can be used to further update the location of the next point.
x.append(-0.7)
y.append(y2)
σ_2 = exponential_cov(x, x, θ)
predictions = [predict(i, x, exponential_cov, θ, σ_2, y) for i in x_pred]
y_pred, sigmas = np.transpose(predictions)
plt.errorbar(x_pred, y_pred, yerr=sigmas, capsize=0)
plt.plot(x, y, "ro")
plt.title("Posterior after 2 points");
Sequential sampling is just a heuristic to show how the covariance structure works. We can just as easily sample several points at once:
x_more = [-2.1, -1.5, 0.3, 1.8, 2.5]
mu, s = conditional(x_more, x, y, θ)
y_more = np.random.multivariate_normal(mu, s)
print(y_more)[ 0.52217716 -1.52084704 0.32309947 -0.93840169 -1.27434377]x += x_more
y += y_more.tolist()
σ_new = exponential_cov(x, x, θ)
predictions = [predict(i, x, exponential_cov, θ, σ_new, y) for i in x_pred]
y_pred, sigmas = np.transpose(predictions)
plt.errorbar(x_pred, y_pred, yerr=sigmas, capsize=0)
plt.plot(x, y, "ro")
plt.title("Posterior after 7 points");
As the density of points grows, we recover a realization (sample function) from the prior GP.
Though we could extend the code above to introduce data and fit a GP by hand, several libraries specialize in this. We'll demonstrate and compare three:
Generate a reproducible data set including a smooth non-linear function corrupted by Gaussian measurement noise. No particular real-world process motivated the shape of the function below; it was chosen simply because it gives each of the three libraries something non-trivial to fit.
rng = np.random.RandomState(42)
N = 40
x_data = np.sort(rng.uniform(-5, 5, N))
f_true = lambda t: np.sin(t) + 0.5 * np.sin(3 * t) * np.exp(-0.1 * t ** 2)
y_data = f_true(x_data) + rng.normal(0, 0.2, N)
X = x_data.reshape(-1, 1) # 2-D design matrix expected by all three libraries
Y = y_data.reshape(-1, 1) # 2-D targets (GPflow / PyMC tabular form)
plt.plot(x_data, y_data, "ko", ms=4)
plt.xlabel("x"); plt.ylabel("y"); plt.title("Simulated data");
scikit-learn provides a GP module under a consistent API. For regression we use GaussianProcessRegressor specifying a covariance function (kernel). Fitting maximizes the log marginal likelihood, avoiding the cross-validation usually needed to choose hyperparameters. The regressor always assumes a zero mean function.
The radial basis function (RBF) kernel is the simplest starting point. It is a special case of the Matérn family in the limit as roughness approaches infinity, producing infinitely smooth realizations, and is available in scikit-learn as RBF. For data with sharper local structure, the Matérn family is usually a better fit because its roughness parameter can be tuned rather than fixed at infinite smoothness.
A flexible kernel to start with is the Matérn covariance. We build it as a sum of an amplitude (ConstantKernel), a Matern component, and observation noise (WhiteKernel) where Γ is the gamma function and K is a modified Bessel function. Three parameters govern the shape of covariance matrices sampled from this function:
Although all three parameters take non-negative real values in general, when ν = p + 1/2 for an integer p, the Matérn function can be written partly as a polynomial of order p, and it generates realizations that are p-times differentiable. For this reason, ν ∈ {3/2, 5/2} are the values seen most often in practice.
The kernel below is a composite kernel built by multiplying a constant term for amplitude with a Matérn term for the covariance structure, then adding a white-noise term for observation error. scikit-learn supports combining kernels with addition and multiplication directly, which makes it straightforward to build kernels tailored to a specific dataset.
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import Matern, WhiteKernel, ConstantKernel
kernel = ConstantKernel() * Matern(length_scale=2, nu=3/2) + WhiteKernel(noise_level=1)All scikit-learn estimators expect a 2-D array of inputs. We already built X with reshape(-1, 1):
print(X.shape)(40, 1)We instantiate the regressor with our kernel and call fit. We pass normalize_y=True so the zero-mean assumption is reasonable for data not centered at zero.
gp = GaussianProcessRegressor(kernel=kernel, normalize_y=True)
gp.fit(X, y_data)(Note: Some cells below emit convergence or deprecation warnings. These are expected and can be ignored)
Fitted attributes carry a trailing underscore. kernel_ returns the kernel with its optimized hyperparameters:
gp.kernel_1.03**2 * Matern(length_scale=1.11, nu=1.5) + WhiteKernel(noise_level=0.044)predict optionally returns posterior standard deviations alongside the mean, which we use to plot a confidence region around the expected function.
x_pred = np.linspace(-6, 6, 200).reshape(-1, 1)
y_pred, sigma = gp.predict(x_pred, return_std=True)
xx = x_pred.ravel()
plt.figure()
plt.plot(xx, y_pred, label="Expected function")
plt.fill_between(xx, y_pred - 2 * sigma, y_pred + 2 * sigma,
alpha=0.3, label="95% confidence region")
plt.plot(x_data, y_data, "ro", ms=4, label="Observations")
plt.legend(); plt.title("scikit-learn GP fit");
Everything up to this point has treated the outcome as continuous, but plenty of real workflows call for a binary decision instead: a flagged transaction, an adverse event, a document that needs manual review. For cases like these, scikit-learn provides GaussianProcessClassifier, which passes a latent Gaussian process through a link function rather than modeling the outcome directly.
The classifier cannot maximize the marginal likelihood in closed form the way GaussianProcessRegressor does, since a Bernoulli likelihood is not conjugate to a Gaussian prior. Instead, it fits a Laplace approximation to the posterior over the latent function and integrates that approximation to produce class probabilities.
To illustrate the API with data already in hand, we can threshold the same synthetic data set into two classes and use the RBF kernel introduced earlier:
from sklearn.gaussian_process import GaussianProcessClassifier
from sklearn.gaussian_process.kernels import RBF
y_labels = (y_data > 0).astype(int) # binarize the simulated targets for classification
gpc = GaussianProcessClassifier(kernel=RBF())
gpc.fit(X, y_labels)
gpc.predict_proba(x_pred)predict_proba returns the probability that each point in x_pred belongs to each class, rather than the mean and variance that predict returns for the regressor. This distinction matters in practice. A classifier used for flagging borderline cases for review needs predict_proba, not the hard label from predict, because the probability itself is usually the useful signal, not the classification.
Worth noting is that binarizing a continuous simulated series at zero has no real connection to the process that generated it. This example demonstrates the classifier's interface, nothing more. For an actual classification problem, the kernel and the threshold should be chosen based on the outcome being modeled, not retrofitted onto a regression example.
GPflow re-implements the ideas of GPy on top of TensorFlow, enabling modern fitting methods (variational inference and MCMC) for larger models. Its API expects tabular inputs for both predictors and outcomes, so we use the 2-D Y we built earlier.
We again use a Matérn-3/2 covariance. Amplitude is an included parameter (variance), so no separate constant kernel is needed.
import gpflow
from gpflow.utilities import print_summary, set_trainable
k = gpflow.kernels.Matern32(variance=1.0, lengthscales=1.2)(Note: You will see a warning, which is expected behavior. You can continue with the analysis without addressing the warning.)
Since our model has a conjugate Gaussian likelihood, we use the Gaussian process regression (GPR) class. Data is passed as a (X, Y) tuple.
m = gpflow.models.GPR((X, Y), kernel=k)
print_summary(m)
In addition to the Matérn hyperparameters there is a likelihood variance (observation noise). We can assign a starting value directly via the parameter's .assign method.
m.likelihood.variance.assign(0.01)<tf.Variable 'UnreadVariable' shape=() dtype=float64, numpy=-4.600266525158521>The printed value above, roughly −4.6, is not the noise variance itself. GPflow stores each parameter on an unconstrained scale internally and applies a transform, in this case a softplus function, when the parameter is read. The value 0.01 we assigned is recovered correctly whenever the model uses m.likelihood.variance in a computation. Only the raw variable's printed representation looks unfamiliar.
The model is fit by minimizing the negative log marginal likelihood. GPflow delegates to SciPy's L-BFGS-B through gpflow.optimizers.Scipy.
opt = gpflow.optimizers.Scipy()
opt.minimize(m.training_loss, m.trainable_variables)
print_summary(m)
Scipy's L-BFGS-B optimizer used above works well for the modest number of hyperparameters in a Gaussian process. But GPflow's TensorFlow backend also supports gradient-based optimizers such as Adam, which is useful when a GP is embedded inside a larger model with many more parameters, or when mini-batch training is required. The pattern is a standard TensorFlow training loop (add import tensorflow as tf here if it has not already been imported earlier in the notebook):
import tensorflow as tf
optimizer = tf.optimizers.Adam(learning_rate=0.01)
@tf.function
def adam_step():
with tf.GradientTape() as tape:
loss = m.training_loss()
grads = tape.gradient(loss, m.trainable_variables)
optimizer.apply_gradients(zip(grads, m.trainable_variables))
for step in range(500):
adam_step()
print_summary(m)
Adam typically needs more iterations than L-BFGS-B to reach a comparable solution, since it relies only on first-order gradient information and a fixed learning rate rather than an approximate Hessian, but it scales to problems where L-BFGS-B becomes impractical.
The predict_y method returns the predictive mean and variance (including observation noise) on an arbitrary grid.
mean, var = m.predict_y(x_pred)
mean = mean.numpy().ravel(); sd = np.sqrt(var.numpy().ravel())
plt.figure()
plt.plot(xx, mean, label="Expected function")
plt.fill_between(xx, mean - 2 * sd, mean + 2 * sd, alpha=0.3, label="95% CI")
plt.plot(x_data, y_data, "ro", ms=4, label="Observations")
plt.legend(); plt.title("GPflow GPR fit");
scikit-learn's GaussianProcessRegressor is difficult to beat for a quick, zero-mean fit with a standard Gaussian likelihood, but it has no mechanism for placing priors on hyperparameters, and no path to a non-Gaussian likelihood. GPflow's TensorFlow-Probability integration adds both. The next two subsections show a maximum a posteriori fit using priors on the kernel hyperparameters, followed by a fully Bayesian fit using Hamiltonian Monte Carlo and a Student-T likelihood that is more robust to outliers than the Gaussian likelihood scikit-learn assumes.
So far nothing is particularly Bayesian. We just maximized the marginal likelihood. We can assign priors to the hyperparameters. GPflow uses TensorFlow-Probability distributions for priors (the old GPflow.priors.Gamma is replaced by tfp.distributions.Gamma). We can also fix a parameter, e.g. the observation noise, using set_trainable(..., False).
import tensorflow_probability as tfp
k2 = gpflow.kernels.Matern32(variance=1.0, lengthscales=1.2)
m2 = gpflow.models.GPR((X, Y), kernel=k2)
m2.kernel.variance.prior = tfp.distributions.Gamma(np.float64(1.0), np.float64(0.1))
m2.kernel.lengthscales.prior = tfp.distributions.Gamma(np.float64(1.0), np.float64(0.1))
m2.likelihood.variance.assign(0.1)
set_trainable(m2.likelihood.variance, False) # fix the measurement error
opt.minimize(m2.training_loss, m2.trainable_variables)
print_summary(m2)
Adding the log-priors to the objective turns the result into a maximum a posteriori (MAP) estimate.
For a fully Bayesian analysis we use GPMC, which jointly samples the hyperparameters and the latent function. We switch to a Student-T likelihood, which is more robust to outliers, and place priors on the kernel parameters.
import tensorflow as tf
k3 = gpflow.kernels.Matern32(variance=1.0, lengthscales=1.2)
lik = gpflow.likelihoods.StudentT()
m3 = gpflow.models.GPMC((X, Y), kernel=k3, likelihood=lik)
m3.kernel.variance.prior = tfp.distributions.Gamma(np.float64(1.0), np.float64(1.0))
m3.kernel.lengthscales.prior = tfp.distributions.Gamma(np.float64(1.0), np.float64(1.0))
set_trainable(m3.likelihood.scale, False)
# MAP initialization helps the sampler start in a good place
gpflow.optimizers.Scipy().minimize(m3.training_loss, m3.trainable_variables,
options=dict(maxiter=20))(Note: You will see a message, which is expected behavior. You can continue with the analysis without addressing the message.)
hmc_helper = gpflow.optimizers.SamplingHelper(
m3.log_posterior_density, m3.trainable_parameters)
hmc = tfp.mcmc.HamiltonianMonteCarlo(
target_log_prob_fn=hmc_helper.target_log_prob_fn,
num_leapfrog_steps=10, step_size=0.01)
adaptive_hmc = tfp.mcmc.SimpleStepSizeAdaptation(
hmc, num_adaptation_steps=50,
target_accept_prob=tf.cast(0.75, tf.float64), adaptation_rate=0.1)
@tf.function
def run_chain():
return tfp.mcmc.sample_chain(
num_results=300, num_burnin_steps=100,
current_state=hmc_helper.current_state,
kernel=adaptive_hmc, trace_fn=None)
samples = run_chain()
print("HMC complete:", len(samples), "parameter groups,",
int(samples[0].shape[0]), "draws each")(Note: You will see a message, which is expected behavior. You can continue with the analysis without addressing the message.)
We can generate predictions from the posterior by assigning sampled states back into the model and drawing function realizations with predict_f_samples.
realizations = []
for i in range(0, 300, 30):
for var, s in zip(hmc_helper.current_state, samples):
var.assign(s[i])
realizations.append(m3.predict_f_samples(x_pred, 1).numpy().squeeze())
realizations = np.array(realizations)
plt.figure()
plt.plot(xx, realizations.T, "C0", alpha=0.3)
plt.plot(x_data, y_data, "ro", ms=4)
plt.title("GPflow GPMC posterior realizations");
PyMC is a general probabilistic-programming library. Models are declared inside a l Model context.
import pymc as pm
with pm.Model() as gp_fit:
ρ = pm.Gamma("ρ", alpha=1, beta=1) # lengthscale
η = pm.Gamma("η", alpha=1, beta=1) # amplitude
K = η ** 2 * pm.gp.cov.Matern32(input_dim=1, ls=ρ)
M = pm.gp.mean.Zero()
σ = pm.HalfCauchy("σ", beta=2.5) # observation noise
gp = pm.gp.Marginal(mean_func=M, cov_func=K)
y_obs = gp.marginal_likelihood("y_obs", X=X, y=y_data, sigma=σ)The sample function, called inside the model context, fits the model with MCMC. By default PyMC uses the auto-tuning NUTS sampler.
with gp_fit:
trace = pm.sample(400, tune=400, chains=2, cores=1,
target_accept=0.9, random_seed=42, progressbar=False)(Note: You will see a message regarding convergence diagnostics, which is expected behavior. You can continue with the analysis without addressing the message.)
import arviz as az
az.plot_trace(trace, var_names=["ρ", "η", "σ"]);
To generate predictions we sample from the posterior predictive distribution. We add a conditional GP over a grid of new points, then draw from it using the posterior samples.
Z = np.linspace(-6, 6, 200).reshape(-1, 1)
with gp_fit:
f_pred = gp.conditional("f_pred", Z)
ppc = pm.sample_posterior_predictive(trace, var_names=["f_pred"],
random_seed=42, progressbar=False)
samps = ppc.posterior_predictive["f_pred"].stack(s=("chain", "draw")).values # (200, n)
mean_pm = samps.mean(axis=1)
lo, hi = np.percentile(samps, [2.5, 97.5], axis=1)
plt.figure()
plt.plot(Z.ravel(), mean_pm, label="Posterior mean")
plt.fill_between(Z.ravel(), lo, hi, alpha=0.3, label="95% credible region")
plt.plot(x_data, y_data, "ro", ms=4, label="Observations")
plt.legend(); plt.title("PyMC posterior predictive");Running NUTS over every hyperparameter is more expensive than the point estimates scikit-learn and GPflow's Scipy optimizer produce, and for a routine regression task with well-behaved data, that expense often buys little. The fully Bayesian approach earns its cost when the downstream decision depends on getting the uncertainty right rather than just the mean. For example when a surrogate model is used to choose the next experiment in an active learning loop, when the dataset is small enough that hyperparameter uncertainty meaningfully affects the predictions, or when domain knowledge about plausible lengthscales or noise levels is available and worth encoding as a prior rather than discarding.
The three libraries occupy different points on the tradeoff between automation and flexibility, and the right choice depends on how much control the problem demands.
scikit-learn
Fastest path to a fitted gaussian process regressor with a standard workflow
Marginal-likelihood optimization
No
No
GPflow
Larger models, custom likelihoods, or when TensorFlow is already in the stack
L-BFGS-B, Adam, or other TensorFlow optimizers
Yes, via TensorFlow Probability
Yes, via variational inference or MCMC
PyMC5
GP as one component of a larger hierarchical Bayesian model
NUTS (Hamiltonian Monte Carlo)
Yes, on every parameter by default
Yes, through PyMC's broader distribution library
scikit-learn is the right default when the goal is a quick, well-understood fit to continuous data.
GPflow is worth the added setup when the model needs to scale beyond what Scipy's optimizer handles comfortably, or when a non-Gaussian likelihood is required.
PyMC5 is the natural choice when the Gaussian process is embedded in a larger model with its own priors and other Bayesian components, since it treats the GP as just another piece of a shared probabilistic program rather than a standalone fit.
Python users have many options for constructing and fitting Gaussian process models. We fit a GP to continuous data with scikit-learn, then extended to more general forms and more sophisticated fitting algorithms with GPflow and PyMC. This is far from a complete survey. Other actively maintained options today include GPy, GPyTorch (PyTorch-based, scales to large data), Stan, and scikit-learn's classifier variant. Try a few to see which fits your workflow best.
Gaussian process regression in Python refers to fitting a GP, a probabilistic model that defines a distribution over functions, to a set of observed input-output pairs, using one of several available libraries. Rather than assuming a fixed functional form such as a line or a polynomial, a GP is specified by a mean function and a covariance function, commonly called a kernel, and the model infers the most probable functions consistent with the data. In Python, three libraries dominate this space: scikit-learn, which offers a GaussianProcessRegressor class for a fast, standard fit; GPflow, which builds on TensorFlow and supports more advanced fitting methods including variational inference and Markov chain Monte Carlo; and PyMC, a general-purpose probabilistic programming library that treats the GP as one component within a larger Bayesian model. All three return not just a predicted value at a new input but a variance around that prediction, which sets gaussian process regression apart from most other regression methods available in Python.
The unknown function connecting inputs to outputs is treated as a random variable drawn from a distribution over functions, rather than as a fixed quantity to be estimated. Before seeing any data, the model expresses a prior belief about that function through a mean function, typically the zero function, and a covariance function, which encodes an assumption about how similar the outputs of two inputs should be based on how close those inputs are to each other. Once training data arrives, Bayesian inference combines that prior with the observed points to produce a posterior distribution over functions, which remains Gaussian because of a key property of the multivariate normal distribution: any subset of its variables, and any subset conditioned on the rest, is itself normally distributed. Evaluating that posterior at a new input yields both a predicted mean and a variance, which together describe the model's best estimate and its uncertainty about that estimate. The covariance function's hyperparameters, such as the lengthscale that controls how far the influence of one point extends, are typically chosen by maximizing the marginal likelihood of the observed data, a process usually referred to as hyperparameter optimization.
Kernel choice depends on what is already known about the underlying function, since no single kernel serves every case well. The RBF kernel, also called the radial basis function kernel or the squared exponential kernel, is a reasonable default when the function is expected to be very smooth, since it produces infinitely differentiable realizations. The Matérn kernel is a more flexible alternative because its roughness parameter can be tuned to produce realizations that are only as smooth as the data justify, and the two most common settings, corresponding to once- and twice-differentiable functions, cover most practical cases. Periodic kernels are appropriate when the function is known to repeat, such as with seasonal time series, and linear kernels suit functions with an underlying trend. In practice, composite kernels built by adding or multiplying simpler kernels together, for example a Matérn term for local structure summed with a periodic term for seasonality, often outperform any single kernel, since real-world functions rarely conform to one assumption alone. The kernel's hyperparameters should always be fit to the data rather than fixed by hand, since even the right kernel family performs poorly with poorly chosen lengthscale or amplitude values.
Predictive uncertainty is the main reason to reach for Gaussian process regression over other methods, provided the training set stays small enough for the computational cost to remain manageable. Because exact gaussian process regression requires inverting a covariance matrix whose size grows with the number of training points, it scales cubically and becomes impractical for large datasets without approximations such as sparse or inducing-point methods. Within that range, it is a non-parametric regression method, so its effective complexity grows with the data rather than being fixed in advance, and it naturally reports wider uncertainty in regions with little or no data, which linear models and most tree-based methods do not do without additional calibration. It is particularly well suited to problems such as Bayesian hyperparameter optimization, where a surrogate model needs to identify not just where the function is likely to be best but where it remains most uncertain, and to scientific or engineering applications where each observation is expensive and every prediction should come with an honest error bar, the kind of rigor MLOps practices expect in production.
The most significant limitation of gaussian process regression is computational cost. Exact inference requires inverting an n-by-n covariance matrix, so both training time and memory use scale poorly as the number of observations grows into the thousands or beyond. Standard implementations in scikit-learn and GPflow's exact GPR class are best suited to datasets of a few hundred to a few thousand points. Larger problems generally require sparse approximations or inducing-point methods, which trade some accuracy for scalability. A second limitation is sensitivity to kernel choice. An inappropriate kernel, or one with poorly optimized hyperparameters, can produce either overconfident or badly underfit predictions, and diagnosing the difference requires some familiarity with what each kernel implies about the function being modeled. Gaussian process regression also assumes, by default, that observation noise is Gaussian and constant across the input space, an assumption that a Student-T or other custom likelihood, available in GPflow and PyMC, can relax but only at the cost of exact inference, requiring approximate methods such as variational inference or Markov chain Monte Carlo instead.
Chris Fonnesbeck is a professor of biostatistics at Vanderbilt University and, as of recent, Principal Quantitative Analyst at the Philadelphia Phillies.