Why Gaussian Processes Beat Neural Networks for Small‑Sample Regression with Uncertainty

Hello, I'm cos大壮. Today we discuss an important algorithmic model: Gaussian Processes (GP).

Why consider Gaussian Processes?

When modeling an unknown continuous function, practitioners often pick a function family—polynomials, neural networks, or kernel methods—and fit the data. This requires manual model selection and hyper‑parameter tuning, and it rarely yields predictive uncertainty.

Gaussian Process Regression (GPR) offers an attractive alternative:

Represent the function as a distribution: define a GP prior over functions.

Apply Bayesian inference to obtain a posterior, giving both predictive mean and variance.

Regularization is built‑in: the kernel controls smoothness and complexity, while marginal likelihood automatically selects hyper‑parameters.

The output is a probability distribution, providing natural confidence estimates.

GPR excels on small‑ to medium‑scale data (up to a few thousand points) where uncertainty matters, such as experimental design, Bayesian optimization, spatial statistics (kriging), and risk calibration.

Core Theory

A Gaussian Process is a stochastic process where any finite collection of function values follows a multivariate Gaussian distribution. We denote the random function as f(·); for any input set X, the vector f(X) is Gaussian with mean function m(X) and covariance (kernel) k(X, X'). A common zero‑mean prior is m(x)=0.

The kernel must be symmetric and positive‑semi‑definite to ensure the covariance matrix is PSD.

Common Kernels

RBF (Gaussian) kernel : length‑scale ℓ and signal variance σ².

Matérn kernel : provides flexible smoothness (e.g., ν=3/2 or ν=5/2).

Rational Quadratic kernel : a scale mixture of RBFs.

Periodic kernel : captures repeating patterns with period p.

Kernels can be added (modeling sums of functions) or multiplied (modulation/gating effects) to express richer priors.

Observation Noise Model

We usually assume i.i.d. Gaussian noise: y = f(x) + ε, with ε ∼ N(0, σ_n²). The training input matrix is denoted X, and the observation vector y.

From Prior to Posterior: Conditional Gaussian

Given training data (X, y) and test inputs X_* , define:

K_tt: training‑training kernel matrix.

K_t*: training‑test kernel matrix.

K_**: test‑test kernel matrix.

With a zero‑mean prior, the joint distribution of training and test outputs is Gaussian. Conditioning yields the posterior mean μ_* = K_t*ᵀ K_tt⁻¹ y and covariance Σ_* = K_** – K_t*ᵀ K_tt⁻¹ K_t* . Adding observation noise to the diagonal of K_tt improves numerical stability.

For noisy predictions, add σ_n² to the predictive variance.

Numerically, we avoid explicit matrix inversion and use Cholesky decomposition: K_tt = L Lᵀ, then solve triangular systems for efficiency and stability.

Marginal Likelihood and Hyper‑parameter Learning

Let θ denote kernel and noise hyper‑parameters (e.g., length‑scale, variance, noise level). The log marginal likelihood (LML) is:

log p(y|θ) = -0.5 * yᵀ K_tt⁻¹ y - 0.5 * log|K_tt| - (n/2) * log 2π

Maximizing LML (type‑II maximum likelihood) yields optimal hyper‑parameters. Its gradient with respect to θ can be derived analytically, and optimization typically uses gradient‑based methods such as L‑BFGS‑B with multiple random restarts to avoid local optima.

LML balances data fit (first term) against model complexity (second term).

Relation to Kernel Ridge Regression (KRR)

KRR predicts: μ_* = K_t* (K_tt + λI)⁻¹ y. If we identify λ with the noise variance σ_n², the predictive means of KRR and GPR coincide, but GPR additionally provides posterior variance and automatically tunes λ via marginal likelihood.

Computational Complexity and Numerical Issues

Training costs O(n³) due to Cholesky decomposition and O(n²) memory. Predictive mean costs O(n n_* ) and variance O(n_*²). GPs work well up to a few thousand points; larger datasets require sparse approximations.

Use Cholesky instead of matrix inverse.

Add a small jitter to the kernel diagonal for numerical stability.

Standardize inputs (zero mean, unit variance) and optionally outputs (normalize_y=True).

Optimize hyper‑parameters in log‑space with sensible bounds.

Full Example

We construct a 1‑D regression task with a complex underlying function (linear trend, low‑ and high‑frequency components, local Gaussian peaks) and heteroscedastic noise. The code uses scikit‑learn's GaussianProcessRegressor with a rich kernel composed of Constant * (RBF + RationalQuadratic + ExpSineSquared) + WhiteKernel.

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
import seaborn as sns
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, RationalQuadratic, ExpSineSquared, WhiteKernel, ConstantKernel
from sklearn.metrics import mean_squared_error
from sklearn.utils import check_random_state

# 1) Random seed
rng = check_random_state(42)

# 2) True function and heteroscedastic noise
def f_true(x):
    return 0.5*x + np.sin(1.5*x) + 0.3*np.exp(-0.5*(x-2.0)**2) + 0.5*np.sin(4*x)

def noise_std(x):
    return 0.20 + 0.10*((np.sin(x) + 1.0)**2)

# 3) Generate data
n_train_total = 1000
X_all = rng.uniform(-6.0, 6.0, size=n_train_total)
X_all = np.sort(X_all)
y_all = f_true(X_all) + rng.normal(0.0, noise_std(X_all))
X_all_2d = X_all.reshape(-1, 1)

X_test = np.linspace(-6.0, 6.0, 400)
X_test_2d = X_test.reshape(-1, 1)
y_test_true = f_true(X_test)

# Subsample for base training set
n_train_base = 160
idx_train = rng.choice(np.arange(n_train_total), size=n_train_base, replace=False)
idx_train.sort()
X_train = X_all[idx_train]
y_train = y_all[idx_train]
X_train_2d = X_train.reshape(-1, 1)

# Rich kernel
kernel_rich = (ConstantKernel(1.0, (1e-3, 1e3)) *
               (RBF(length_scale=1.5, length_scale_bounds=(1e-2, 1e2)) +
                RationalQuadratic(alpha=1.0, length_scale=1.0,
                                   alpha_bounds=(1e-3, 1e3),
                                   length_scale_bounds=(1e-2, 1e2)) +
                ExpSineSquared(length_scale=1.0, periodicity=3.0,
                               length_scale_bounds=(1e-2, 1e2),
                               periodicity_bounds=(0.5, 10.0))) +
               WhiteKernel(noise_level=0.1, noise_level_bounds=(1e-5, 1e1)))

gpr_rich = GaussianProcessRegressor(kernel=kernel_rich,
                                    alpha=0.0,
                                    normalize_y=True,
                                    n_restarts_optimizer=8,
                                    random_state=42)

# Fit and predict
gpr_rich.fit(X_train_2d, y_train)
y_pred_mean, y_pred_std = gpr_rich.predict(X_test_2d, return_std=True)

# Evaluation
rmse_test = np.sqrt(mean_squared_error(y_test_true, y_pred_mean))
var_pred = np.maximum(y_pred_std**2, 1e-9)
nlpd = 0.5*np.log(2*np.pi*var_pred) + 0.5*((y_test_true - y_pred_mean)**2)/var_pred
nlpd_mean = np.mean(nlpd)

print("==== Metrics ====")
print("Rich kernel GPR optimal kernel:", gpr_rich.kernel_)
print(f"Rich kernel GPR test RMSE: {rmse_test:.4f}")
print(f"Rich kernel GPR test NLPD (approx): {nlpd_mean:.4f}")

The code prints the learned kernel, test RMSE, and approximate negative log‑predictive density (NLPD).

Visualization Analysis

Posterior Prediction shows the mean fit and a 95 % confidence band. Uncertainty is low near training points and higher in sparse regions, illustrating Bayesian non‑parametric behavior. Because the data contain heteroscedastic noise but the standard GP assumes homoscedasticity, the confidence band can be mis‑calibrated, motivating heteroscedastic extensions.

Simplified Kernel NLML Landscape visualizes the negative log‑marginal likelihood over log(length‑scale) and log(noise) for a simple Constant × RBF + White kernel. The global basin indicates the optimal hyper‑parameters, while local minima reveal potential traps for gradient‑based optimizers.

Learning Curve plots test RMSE (mean ± 1 σ) versus training sample size, revealing whether the model is data‑hungry and whether performance approaches the irreducible noise floor.

Standardized Residual Histogram vs. Standard Normal PDF checks calibration: well‑calibrated uncertainties produce residuals that follow N(0, 1). Deviations indicate over‑ or under‑confidence, especially in heteroscedastic regions.

Model Training Process and Key Details

1. Data preprocessing and normalization : Standardize inputs to zero mean and unit variance; optionally normalize outputs (normalize_y=True) for stable optimization.

2. Kernel selection and composition :

Start with a simple baseline kernel (RBF + White).

Inspect residuals; if multiple scales or periodicities appear, add RationalQuadratic or ExpSineSquared kernels.

Use additive kernels for sum of effects, multiplicative kernels for gating or modulation.

Set reasonable bounds for hyper‑parameters to avoid edge‑effects.

3. Hyper‑parameter optimization :

Maximize marginal likelihood with several random restarts (n_restarts_optimizer ≥ 5).

For very small datasets, perform a coarse grid search followed by local refinement.

Be aware of multiple local optima in complex kernel landscapes.

4. Uncertainty and numerical stability :

Use Cholesky decomposition instead of explicit matrix inversion.

Add a small jitter to the kernel diagonal (scikit‑learn does this internally, but you can increase the WhiteKernel lower bound).

Predictive variance corresponds to the latent function; add observation noise if you need variance of y.

5. Evaluation and diagnostics :

Point‑wise metrics: RMSE, MAE.

Probabilistic metric: NLPD (negative log‑predictive density).

Calibration tools: reliability curves, standardized residual histograms.

Possible Optimizations

Standardize inputs (z‑score) and center outputs.

Use many optimizer restarts (≥ 5) to increase chance of finding a good solution.

Set sensible bounds for kernel hyper‑parameters.

Incorporate physical priors (e.g., known periodicity) early in the kernel.

Optimize in log‑space (default in scikit‑learn) for robustness.

For larger datasets, consider sparse approximations (inducing points, SVGP).

Conclusion

Theoretically, GPR combines a Gaussian prior with conditional Gaussian inference to yield closed‑form predictive means and variances, while marginal likelihood provides automatic regularization and hyper‑parameter selection.

Our analysis shows that standard homoscedastic GP can be mis‑calibrated under heteroscedastic noise, suggesting extensions such as heteroscedastic GP or more realistic noise models.

For small‑ to medium‑scale problems, GPR delivers excellent point predictions together with well‑calibrated uncertainty, offering strong interpretability and probabilistic consistency.