Introduction

Most machine-learning models hand you a single number: "the temperature tomorrow will be 22 °C." A Gaussian process hands you something richer: "the temperature will be 22 °C — but I'm only ±1 °C confident near your training data, and ±8 °C out where you've given me nothing."

That extra information — the confidence band — is not decoration. It tells you where the model is trustworthy and where it is merely extrapolating, which is exactly what you need before acting on a prediction.

The idea traces back to Carl Friedrich Gauss and was formalized for spatial statistics by Danie Krige in the 1950s (the technique is still called kriging in geostatistics). The modern machine-learning treatment was consolidated by Carl Edward Rasmussen and Christopher Williams in their 2006 textbook Gaussian Processes for Machine Learning, now the standard reference.

At heart a GP is a prior over functions: instead of placing a distribution over parameters (intercept, slope, …), you place a distribution directly over the space of all smooth functions. Add data and the prior narrows into a posterior — a whole family of functions that are consistent with what you observed, together with how likely each one is.

Try It

Click anywhere on the canvas to add an observation. The blue curve is the GP posterior mean — the model's best guess for every x. The shaded band is ±2 standard deviations: the region the true function almost certainly passes through. Points outside existing data widen the band; points close together narrow it.

<p class="hint">Click anywhere on the chart to add a data point. The blue curve is the GP posterior mean; the shaded band is ±2 standard deviations. Far from data the band widens — the model admits it doesn't know.</p>
<canvas id="gp" width="560" height="320"></canvas>
<div class="controls">
  <label>Length-scale ℓ <input type="range" id="ls" min="5" max="120" value="40"> <span id="ls-val">40</span></label>
  <button id="clear" type="button" class="ghost">Clear points</button>
</div>
<div class="status" id="status">Add points to see the GP update.</div>

* { box-sizing: border-box; }
body { font-family: system-ui, sans-serif; color: #222; margin: 0; }
.hint { font-size: .88rem; color: #444; margin: 0 0 .6rem; line-height: 1.45; }
canvas { border: 1px solid #cdd9e3; border-radius: 10px; cursor: crosshair;
         max-width: 100%; display: block; background: #f8fafc; }
.controls { display: flex; align-items: center; gap: 1rem; flex-wrap: wrap;
            margin: .6rem 0 .3rem; font-size: .88rem; }
label { display: flex; align-items: center; gap: .4rem; }
input[type=range] { width: 90px; accent-color: #1d3557; }
.status { font-size: .9rem; font-weight: 600; min-height: 1.3em; color: #1d3557; }
button { font: 600 13px system-ui; padding: .35rem .75rem; border: 1px solid #1d3557;
         background: #1d3557; color: #fff; border-radius: 8px; cursor: pointer; }
button.ghost { background: #fff; color: #1d3557; }

const fmt = (s, v) => s.replace(/%(\w+)%/g, (_, k) => String(v[k] ?? ''));
// Gaussian Process regression — squared-exponential kernel, noise σ_n² = 0.02.
// For n training points we build K (n×n) + σ_n²I, invert it (O(n³) Cholesky),
// then predict mean and variance at each grid point in O(n) per point.
// All linear algebra is done in plain JS arrays (no imports).

const canvas = document.getElementById('gp');
const ctx    = canvas.getContext('2d');
const W = canvas.width, H = canvas.height;
const PAD = { l: 36, r: 16, t: 16, b: 28 };
const NOISE_VAR = 0.04;        // observation noise variance
const PRIOR_VAR = 1.0;         // signal variance σ²

let pts = [];   // { px, py, x, y } — pixel and domain coords
let ell = 40;   // length-scale in pixels

// ── Coordinate helpers ────────────────────────────────────────────────────
const pxToX = px => (px - PAD.l) / (W - PAD.l - PAD.r);        // [0,1]
const pxToY = py => 1 - (py - PAD.t) / (H - PAD.t - PAD.b);    // [0,1] (bottom=0)
const xToPx = x  => PAD.l + x  * (W - PAD.l - PAD.r);
const yToPy = y  => PAD.t + (1 - y) * (H - PAD.t - PAD.b);

// ── Kernel (squared-exponential in pixel space) ───────────────────────────────────────────────────────────
function rbf(x1, x2) {
  const dx = x1 - x2;
  return PRIOR_VAR * Math.exp(-dx * dx / (2 * ell * ell));
}

// ── Tiny linear algebra helpers ─────────────────────────────────────────────────
// Cholesky decomposition of a symmetric positive-definite matrix A.
// Returns lower-triangular L such that A = L Lᵀ.
function cholesky(A) {
  const n = A.length;
  const L = Array.from({ length: n }, () => new Array(n).fill(0));
  for (let i = 0; i < n; i++) {
    for (let j = 0; j <= i; j++) {
      let s = A[i][j];
      for (let k = 0; k < j; k++) s -= L[i][k] * L[j][k];
      L[i][j] = (i === j) ? Math.sqrt(Math.max(s, 1e-12)) : s / L[j][j];
    }
  }
  return L;
}
// Solve Lv = b (forward substitution).
function solveLower(L, b) {
  const n = b.length, v = new Array(n).fill(0);
  for (let i = 0; i < n; i++) {
    let s = b[i];
    for (let k = 0; k < i; k++) s -= L[i][k] * v[k];
    v[i] = s / L[i][i];
  }
  return v;
}
// Solve Lᵀ u = v (back substitution).
function solveUpper(L, v) {
  const n = v.length, u = new Array(n).fill(0);
  for (let i = n - 1; i >= 0; i--) {
    let s = v[i];
    for (let k = i + 1; k < n; k++) s -= L[k][i] * u[k];
    u[i] = s / L[i][i];
  }
  return u;
}

// ── GP posterior ───────────────────────────────────────────────────
// Returns { mean[], std[] } for an array of pixel x-positions `xs`.
function gpPredict(trainPxs, trainYs, xs) {
  const n = trainPxs.length;
  // Build K + σ_n² I
  const K = Array.from({ length: n }, (_, i) =>
    Array.from({ length: n }, (__, j) =>
      rbf(trainPxs[i], trainPxs[j]) + (i === j ? NOISE_VAR : 0)
    )
  );
  const L = cholesky(K);
  // α = (K + σ²I)⁻¹ y  via Cholesky
  const alpha = solveUpper(L, solveLower(L, trainYs));

  const means = [], stds = [];
  for (const xstar of xs) {
    // k★ vector
    const kstar = trainPxs.map(xp => rbf(xstar, xp));
    // mean = k★ᵀ α
    const mean = kstar.reduce((s, v, i) => s + v * alpha[i], 0);
    // variance = k(x★,x★) − k★ᵀ (K+σ²I)⁻¹ k★
    const v = solveLower(L, kstar);
    const var_ = Math.max(0, rbf(xstar, xstar) - v.reduce((s, vi) => s + vi * vi, 0));
    means.push(mean);
    stds.push(Math.sqrt(var_));
  }
  return { means, stds };
}

// ── Drawing ────────────────────────────────────────────────────────
const GRID = 80; // x-resolution for the prediction curve
function draw() {
  ctx.clearRect(0, 0, W, H);

  // Axes
  ctx.strokeStyle = '#cdd9e3'; ctx.lineWidth = 1;
  ctx.strokeRect(PAD.l, PAD.t, W - PAD.l - PAD.r, H - PAD.t - PAD.b);

  // Axis labels
  ctx.fillStyle = '#7a8fa0'; ctx.font = '11px system-ui';
  ctx.textAlign = 'center';
  ctx.fillText('0', xToPx(0), H - 4);
  ctx.fillText('1', xToPx(1), H - 4);
  ctx.textAlign = 'right';
  ctx.fillText('1', PAD.l - 4, yToPy(1) + 4);
  ctx.fillText('0', PAD.l - 4, yToPy(0) + 4);

  // Grid x positions (pixel)
  const gridPx = Array.from({ length: GRID + 1 }, (_, i) =>
    PAD.l + (i / GRID) * (W - PAD.l - PAD.r));

  if (pts.length === 0) {
    // Prior: mean=0.5, std from prior variance
    const priorStd = Math.sqrt(PRIOR_VAR);
    _drawBand(gridPx, gridPx.map(() => 0.5), gridPx.map(() => priorStd), '#bfcfff', '#e8eeff');
    _drawMean(gridPx, gridPx.map(() => 0.5));
    document.getElementById('status').textContent = 'Add points to see the GP update.';
    return;
  }

  const trainPxs = pts.map(p => p.px);
  const trainYs  = pts.map(p => p.y);
  const { means, stds } = gpPredict(trainPxs, trainYs, gridPx);

  _drawBand(gridPx, means, stds, '#93b4e8', '#dce8f8');
  _drawMean(gridPx, means);

  // Data points
  for (const p of pts) {
    ctx.beginPath();
    ctx.arc(p.px, p.py, 5, 0, 2 * Math.PI);
    ctx.fillStyle = '#1d3557'; ctx.fill();
    ctx.strokeStyle = '#fff'; ctx.lineWidth = 1.5; ctx.stroke();
  }

  const maxStd = Math.max(...stds).toFixed(3);
  document.getElementById('status').textContent =
    fmt('%length% point%v% · max σ = %maxStd% · length-scale ℓ = %ell% px', { length: pts.length, v: pts.length > 1 ? 's' : '', maxStd: maxStd, ell: ell });
}

function _drawBand(xs, means, stds, strokeColor, fillColor) {
  // Upper envelope
  ctx.beginPath();
  xs.forEach((px, i) => {
    const py = yToPy(Math.min(1, Math.max(0, means[i] + 2 * stds[i])));
    i === 0 ? ctx.moveTo(px, py) : ctx.lineTo(px, py);
  });
  // Lower envelope (reversed)
  for (let i = xs.length - 1; i >= 0; i--) {
    const py = yToPy(Math.min(1, Math.max(0, means[i] - 2 * stds[i])));
    ctx.lineTo(xs[i], py);
  }
  ctx.closePath();
  ctx.fillStyle = fillColor; ctx.fill();
}

function _drawMean(xs, means) {
  ctx.beginPath();
  xs.forEach((px, i) => {
    const py = yToPy(Math.min(1, Math.max(0, means[i])));
    i === 0 ? ctx.moveTo(px, py) : ctx.lineTo(px, py);
  });
  ctx.strokeStyle = '#1d4ed8'; ctx.lineWidth = 2.2; ctx.stroke();
}

// ── Interaction ────────────────────────────────────────────────────
canvas.addEventListener('click', e => {
  const rect = canvas.getBoundingClientRect();
  const scaleX = W / rect.width, scaleY = H / rect.height;
  const px = (e.clientX - rect.left) * scaleX;
  const py = (e.clientY - rect.top)  * scaleY;
  // Ignore clicks in axes padding
  if (px < PAD.l || px > W - PAD.r || py < PAD.t || py > H - PAD.b) return;
  const x = pxToX(px), y = pxToY(py);
  pts.push({ px, py, x, y });
  draw();
});

const lsSlider = document.getElementById('ls');
const lsVal    = document.getElementById('ls-val');
lsSlider.addEventListener('input', () => {
  ell = +lsSlider.value;
  lsVal.textContent = ell;
  draw();
});

document.getElementById('clear').addEventListener('click', () => {
  pts = []; draw();
});

draw();

// ── Node.js verification (run: node gaussian-processes-demo.js) ────────────────────────────────────────────────────
// Asserts: (1) GP mean passes through training points (±noise),
// (2) posterior variance < prior variance at observed xs,
// (3) posterior variance ≈ prior variance far from all data.
if (typeof module !== 'undefined') {
  const trainPxs = [100, 200, 350];
  const trainYs  = [0.3, 0.8, 0.5];
  const { means: mTrain, stds: sTrain } = gpPredict(trainPxs, trainYs, trainPxs);
  const priorStd = Math.sqrt(PRIOR_VAR);
  trainPxs.forEach((_, i) => {
    const err = Math.abs(mTrain[i] - trainYs[i]);
    console.assert(err < 0.1, fmt('GP mean should be near training y; err=%err%', { err: err }));
    console.assert(sTrain[i] < priorStd, 'Posterior std should be < prior std at data');
  });
  const farPx = [2000];
  const { stds: sFar } = gpPredict(trainPxs, trainYs, farPx);
  console.assert(Math.abs(sFar[0] - priorStd) < 0.05, 'Far from data std should approach prior');
  console.log('All GP assertions passed.');
}
Result
Console

Notice what happens far from your points: the band widens back to its prior width. The model admits ignorance instead of pretending confidence — a property that ordinary polynomial or neural-network regression does not have by default.

The Real Complexity

A Gaussian process is defined by two things: a mean function m(x)m(x) (often taken to be zero) and a kernel (covariance function) k(x,x′)k(x, x') that encodes how similar two inputs are. The most common kernel is the squared-exponential (RBF):

k(x,x′)=σ2exp⁡ ⁣(−∥x−x′∥22ℓ2)k(x, x') = \sigma^2 \exp\!\left(-\frac{\|x - x'\|^2}{2\ell^2}\right)

Here σ2\sigma^2 controls the overall variance and ℓ\ell is the length-scale — how far apart two inputs must be before their outputs become nearly independent.

Given n training points, inference builds the n × n covariance matrix K and inverts it. That inversion costs O(n3)O(n^{3}) time and O(n2)O(n^{2}) memory — exact, but it means naive GPs hit a wall around n ≈ 10 000 points on a modern laptop.

  • Why cubic? Gaussian elimination on an n × n matrix takes n3n^{3} multiplications. There is no known way to sidestep this for a general dense matrix.
  • Approximations rescue scale. Sparse GPs with m ≪ n inducing points reduce cost to O(nm2)O(nm^{2}) training and O(m2)O(m^{2}) prediction. Random Fourier features (Rahimi & Recht, 2007) approximate the kernel with a low-rank expansion for even cheaper inference.
  • The payoff is calibrated uncertainty. Unlike neural networks, a well-tuned GP posterior is provably calibrated under its kernel assumptions: the 95% band really does cover the true function 95% of the time (in expectation over the prior).

This places GP regression in a sweet spot: it is a solved, polynomial-time algorithm (O(n3)O(n^{3})) for exact inference, but one whose cost makes it a research-active area for scaling to modern dataset sizes. See also Bayesian inference for the broader probabilistic framework that GPs inhabit.

Where It Matters

Calibrated uncertainty is valuable wherever decisions have consequences:

  • Bayesian optimization: finding the maximum of an expensive black-box function (tuning neural-network hyperparameters, designing experiments) uses a GP as a surrogate model. The uncertainty band tells the optimizer where to sample next — explore the wide band or exploit the high mean.
  • Geostatistics (kriging): predicting mineral concentrations, rainfall, or pollution levels between sparse measurement stations. The original application that gave GPs their name in spatial statistics.
  • Drug discovery and materials science: predicting molecular properties from limited lab assays; the GP band flags which candidates are worth testing next.
  • Robotics and control: learning the dynamics of a robot arm from sensor data, with uncertainty bounds that feed directly into safe planning algorithms.
  • Climate and weather modeling: emulating expensive climate simulations at untried parameter settings, with honest error bars on the emulator output.
  • Active learning: any setting where labeling data is costly and you want to query the most informative point — the GP uncertainty tells you where you know the least.

The common thread: wherever you need a model that says "I don't know" rather than confidently extrapolating, a GP is a natural first choice. See also PAC learning for the contrast with worst-case learning theory.

Conclusion

A Gaussian process is a rare thing in machine learning: a model that is simultaneously exact (given its assumptions), interpretable (every parameter has a physical meaning), and honest (it widens its band where it should). The cost is cubic in training size, which is why sparse and approximate GPs are an active research area — but for small-to-medium datasets the exact version remains one of the most principled tools available.

The deeper lesson is philosophical. Most models give you a point estimate and hide their uncertainty in held-out metrics. A GP makes uncertainty a first-class citizen of the prediction itself. That shift — from "here is the answer" to "here is the answer, and here is what I don't know" — is exactly what Bayesian inference asks of every model.

Share this article

Pick a channel — or use your device's native share sheet.

Comments

Loading comments...

https://www.kipuhub.com/en/article/gaussian-processes/Content licensed under CC BY-NC 4.0.