Introduction

Imagine you are training a neural network and you need to pick the right learning rate, the number of hidden layers, and the dropout probability. Each combination requires a full training run that can take hours — sometimes days — on expensive hardware. Trying every combination is simply impossible.

This is the setting Bayesian optimization was built for: find the best input to an expensive black-box function using as few evaluations as possible. The word black-box means you can measure the output for any input you choose, but you cannot see a formula or gradient — only the result.

The key idea is to be smart about where you look next. Instead of searching at random (which wastes evaluations) or following a gradient (which doesn't exist), Bayesian optimization builds a probabilistic model of the unknown function, uses that model to decide the most promising next point, evaluates the function there, then updates the model and repeats. Every evaluation makes the model sharper.

Introduced in its modern form by Jonas Mockus in the 1970s and popularized for machine learning by Snoek, Larochelle, and Adams (2012), Bayesian optimization is now the backbone of hyperparameter tuning systems like Google Vizier and Optuna.

Try It: Minimize with a Surrogate

The demo below runs Bayesian optimization on a noisy one-dimensional function. The blue curve is the true (hidden) function; the orange curve is the Gaussian process surrogate the algorithm has built from the evaluations so far. The green band shows uncertainty — wider means less explored.

<p class="hint">The <b>blue curve</b> is the hidden true function. The <b>orange curve</b>
is the Gaussian process surrogate. The <b>green band</b> is uncertainty (±2σ).
Dots are evaluations so far. Click <b>Next evaluation</b> to let the algorithm
pick the point with highest Expected Improvement.</p>
<canvas id="chart" width="560" height="220"></canvas>
<div class="info" id="info">Press "Next evaluation" to start.</div>
<div class="btns">
  <button id="btn-next" type="button">Next evaluation</button>
  <button id="btn-reset" type="button" class="ghost">Reset</button>
</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 { display: block; width: 100%; max-width: 560px; border: 1px solid #dde3ea; border-radius: 8px; background: #f8fafc; }
.info { font-size: .92rem; font-weight: 600; min-height: 1.4em; margin: .5rem 0; color: #1d3557; }
.btns { display: flex; gap: .5rem; flex-wrap: wrap; }
button { font: 600 14px system-ui; padding: .45rem .9rem; 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] ?? ''));
// Bayesian Optimization demo — 1-D Gaussian Process surrogate with
// Expected Improvement acquisition function.
//
// True function: f(x) = sin(3x) + x*cos(x) + small Gaussian noise
// We want to find the MINIMUM.
// GP uses a squared-exponential (RBF) kernel.
// Expected Improvement (EI) picks the next x.

const W = 560, H = 220;
const XMIN = 0, XMAX = 6.5;
const NOISE = 0.12;        // observation noise std-dev
const LENGTHSCALE = 0.8;   // RBF kernel length-scale
const SIGNAL_VAR = 1.5;    // RBF kernel signal variance
const NOISE_VAR = NOISE * NOISE;
const NX = 300;            // grid resolution for plotting

// True objective (unknown to the algorithm, known to us for display)
function trueF(x) {
  return Math.sin(3 * x) + x * Math.cos(x);
}

// Noisy observation
function observe(x) {
  return trueF(x) + NOISE * randn();
}

// Box-Muller normal sample
function randn() {
  let u = 0, v = 0;
  while (u === 0) u = Math.random();
  while (v === 0) v = Math.random();
  return Math.sqrt(-2 * Math.log(u)) * Math.cos(2 * Math.PI * v);
}

// RBF kernel
function kernel(x1, x2) {
  const d = x1 - x2;
  return SIGNAL_VAR * Math.exp(-0.5 * d * d / (LENGTHSCALE * LENGTHSCALE));
}

// Build covariance matrix K(X,X) + noise*I
function buildK(xs) {
  const n = xs.length;
  const K = [];
  for (let i = 0; i < n; i++) {
    K[i] = [];
    for (let j = 0; j < n; j++) {
      K[i][j] = kernel(xs[i], xs[j]) + (i === j ? NOISE_VAR : 0);
    }
  }
  return K;
}

// Cholesky decomposition L s.t. K = L * L^T
function cholesky(K) {
  const n = K.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 = K[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-10)) : s / L[j][j];
    }
  }
  return L;
}

// Solve L * y = b (forward substitution)
function solveL(L, b) {
  const n = b.length;
  const y = new Array(n);
  for (let i = 0; i < n; i++) {
    let s = b[i];
    for (let j = 0; j < i; j++) s -= L[i][j] * y[j];
    y[i] = s / L[i][i];
  }
  return y;
}

// Solve L^T * x = y (back substitution)
function solveLT(L, y) {
  const n = y.length;
  const x = new Array(n);
  for (let i = n - 1; i >= 0; i--) {
    let s = y[i];
    for (let j = i + 1; j < n; j++) s -= L[j][i] * x[j];
    x[i] = s / L[i][i];
  }
  return x;
}

// GP posterior mean and variance at a set of test points
function gpPredict(trainX, trainY, testX) {
  const n = trainX.length;
  if (n === 0) {
    // Prior: mean 0, var = SIGNAL_VAR
    return testX.map(() => ({ mean: 0, std: Math.sqrt(SIGNAL_VAR) }));
  }
  const K = buildK(trainX);
  const L = cholesky(K);
  const alpha = solveLT(L, solveL(L, trainY));

  return testX.map(xs => {
    // k_* = [k(xs, x1), ..., k(xs, xn)]
    const kstar = trainX.map(x => kernel(xs, x));
    // mean = kstar^T * alpha
    let mean = 0;
    for (let i = 0; i < n; i++) mean += kstar[i] * alpha[i];
    // variance = k(xs,xs) - kstar^T * (K^{-1} kstar)
    const v = solveL(L, kstar);
    let var_ = kernel(xs, xs);
    for (let i = 0; i < n; i++) var_ -= v[i] * v[i];
    var_ = Math.max(var_, 0);
    return { mean, std: Math.sqrt(var_) };
  });
}

// Standard normal CDF (approximation)
function normCDF(z) {
  const t = 1 / (1 + 0.2316419 * Math.abs(z));
  const poly = t * (0.319381530 + t * (-0.356563782 + t * (1.781477937 + t * (-1.821255978 + t * 1.330274429))));
  const p = 1 - (1 / Math.sqrt(2 * Math.PI)) * Math.exp(-0.5 * z * z) * poly;
  return z >= 0 ? p : 1 - p;
}

// Standard normal PDF
function normPDF(z) {
  return (1 / Math.sqrt(2 * Math.PI)) * Math.exp(-0.5 * z * z);
}

// Expected Improvement (minimization)
function ei(mean, std, fBest) {
  if (std < 1e-10) return 0;
  const z = (fBest - mean) / std;
  return (fBest - mean) * normCDF(z) + std * normPDF(z);
}

// ── State ─────────────────────────────────────────────────────────────────
let obsX = [], obsY = [];
let evalCount = 0;

// Build x-grid for plotting/acquisition
const grid = [];
for (let i = 0; i < NX; i++) grid.push(XMIN + (XMAX - XMIN) * i / (NX - 1));

// Pre-compute true function on grid (for display)
const trueGrid = grid.map(trueF);
const trueMin = Math.min(...trueGrid);
const trueMinX = grid[trueGrid.indexOf(trueMin)];

// ── Drawing ───────────────────────────────────────────────────────────────
const canvas = document.getElementById('chart');
const ctx = canvas.getContext('2d');
const DPR = window.devicePixelRatio || 1;
canvas.width = W * DPR; canvas.height = H * DPR;
canvas.style.width = W + 'px'; canvas.style.height = H + 'px';
ctx.scale(DPR, DPR);

const PAD = { top: 12, right: 14, bottom: 32, left: 42 };
const PW = W - PAD.left - PAD.right;
const PH = H - PAD.top - PAD.bottom;

function toCanvasX(x) { return PAD.left + (x - XMIN) / (XMAX - XMIN) * PW; }
function toCanvasY(y, ymin, ymax) { return PAD.top + PH - (y - ymin) / (ymax - ymin) * PH; }

function drawAll() {
  ctx.clearRect(0, 0, W, H);

  // GP prediction on grid
  const preds = gpPredict(obsX, obsY, grid);
  const means = preds.map(p => p.mean);
  const stds  = preds.map(p => p.std);

  // y axis range: union of true + GP ±2σ
  const allY = [...trueGrid, ...means.map((m, i) => m + 2 * stds[i]),
                             ...means.map((m, i) => m - 2 * stds[i]), ...obsY];
  const yMin = Math.min(...allY) - 0.3;
  const yMax = Math.max(...allY) + 0.3;

  // Uncertainty band ±2σ
  ctx.beginPath();
  ctx.moveTo(toCanvasX(grid[0]), toCanvasY(means[0] + 2 * stds[0], yMin, yMax));
  for (let i = 1; i < NX; i++)
    ctx.lineTo(toCanvasX(grid[i]), toCanvasY(means[i] + 2 * stds[i], yMin, yMax));
  for (let i = NX - 1; i >= 0; i--)
    ctx.lineTo(toCanvasX(grid[i]), toCanvasY(means[i] - 2 * stds[i], yMin, yMax));
  ctx.closePath();
  ctx.fillStyle = 'rgba(34,197,94,0.15)';
  ctx.fill();

  // True function (blue)
  ctx.beginPath();
  ctx.strokeStyle = '#2563eb'; ctx.lineWidth = 1.8;
  ctx.setLineDash([6, 4]);
  grid.forEach((x, i) => {
    const cx = toCanvasX(x), cy = toCanvasY(trueGrid[i], yMin, yMax);
    i === 0 ? ctx.moveTo(cx, cy) : ctx.lineTo(cx, cy);
  });
  ctx.stroke();
  ctx.setLineDash([]);

  // GP mean (orange)
  ctx.beginPath();
  ctx.strokeStyle = '#f97316'; ctx.lineWidth = 2;
  grid.forEach((x, i) => {
    const cx = toCanvasX(x), cy = toCanvasY(means[i], yMin, yMax);
    i === 0 ? ctx.moveTo(cx, cy) : ctx.lineTo(cx, cy);
  });
  ctx.stroke();

  // Axes
  ctx.strokeStyle = '#94a3b8'; ctx.lineWidth = 1;
  ctx.beginPath();
  ctx.moveTo(PAD.left, PAD.top); ctx.lineTo(PAD.left, PAD.top + PH);
  ctx.lineTo(PAD.left + PW, PAD.top + PH);
  ctx.stroke();

  // y=0 line
  if (yMin < 0 && yMax > 0) {
    const y0 = toCanvasY(0, yMin, yMax);
    ctx.strokeStyle = '#cbd5e1'; ctx.lineWidth = 0.8;
    ctx.setLineDash([3, 3]);
    ctx.beginPath(); ctx.moveTo(PAD.left, y0); ctx.lineTo(PAD.left + PW, y0); ctx.stroke();
    ctx.setLineDash([]);
  }

  // Axis labels
  ctx.fillStyle = '#64748b'; ctx.font = '11px system-ui';
  ctx.textAlign = 'center';
  [0, 1, 2, 3, 4, 5, 6].forEach(v => {
    if (v <= XMAX) {
      const cx = toCanvasX(v);
      ctx.fillText(v, cx, PAD.top + PH + 18);
    }
  });
  ctx.textAlign = 'right';
  [-3, -2, -1, 0, 1, 2, 3, 4].forEach(v => {
    if (v >= yMin - 0.1 && v <= yMax + 0.1) {
      ctx.fillText(v, PAD.left - 5, toCanvasY(v, yMin, yMax) + 4);
    }
  });

  // Observed points
  obsX.forEach((x, i) => {
    const cx = toCanvasX(x), cy = toCanvasY(obsY[i], yMin, yMax);
    ctx.beginPath();
    ctx.arc(cx, cy, 4, 0, 2 * Math.PI);
    ctx.fillStyle = '#1d3557'; ctx.fill();
    ctx.strokeStyle = '#fff'; ctx.lineWidth = 1.5; ctx.stroke();
  });

  // EI acquisition: find next point (highlight)
  if (obsX.length > 0) {
    const fBest = Math.min(...obsY);
    let bestEI = -Infinity, bestIdx = 0;
    for (let i = 0; i < NX; i++) {
      const v = ei(means[i], stds[i], fBest);
      if (v > bestEI) { bestEI = v; bestIdx = i; }
    }
    const nx = grid[bestIdx];
    const cx = toCanvasX(nx);
    ctx.strokeStyle = '#16a34a'; ctx.lineWidth = 1.5;
    ctx.setLineDash([4, 3]);
    ctx.beginPath(); ctx.moveTo(cx, PAD.top); ctx.lineTo(cx, PAD.top + PH); ctx.stroke();
    ctx.setLineDash([]);
    ctx.fillStyle = '#16a34a'; ctx.font = '10px system-ui'; ctx.textAlign = 'left';
    ctx.fillText('next →', cx + 3, PAD.top + 12);
  }

  // Legend
  ctx.font = '10px system-ui'; ctx.textAlign = 'left';
  const lx = PAD.left + 8, ly = PAD.top + 14;
  ctx.strokeStyle = '#2563eb'; ctx.lineWidth = 1.8; ctx.setLineDash([5, 3]);
  ctx.beginPath(); ctx.moveTo(lx, ly); ctx.lineTo(lx + 20, ly); ctx.stroke();
  ctx.setLineDash([]);
  ctx.fillStyle = '#2563eb'; ctx.fillText('true f(x)', lx + 24, ly + 4);
  ctx.strokeStyle = '#f97316'; ctx.lineWidth = 2;
  ctx.beginPath(); ctx.moveTo(lx, ly + 14); ctx.lineTo(lx + 20, ly + 14); ctx.stroke();
  ctx.fillStyle = '#f97316'; ctx.fillText('GP surrogate', lx + 24, ly + 18);
}

function pickNext() {
  const preds = gpPredict(obsX, obsY, grid);
  const means = preds.map(p => p.mean);
  const stds  = preds.map(p => p.std);

  let xNext;
  if (obsX.length === 0) {
    // First point: random
    xNext = XMIN + Math.random() * (XMAX - XMIN);
  } else {
    const fBest = Math.min(...obsY);
    let bestEI = -Infinity, bestIdx = 0;
    for (let i = 0; i < NX; i++) {
      const v = ei(means[i], stds[i], fBest);
      if (v > bestEI) { bestEI = v; bestIdx = i; }
    }
    xNext = grid[bestIdx];
  }
  const yNext = observe(xNext);
  obsX.push(xNext);
  obsY.push(yNext);
  evalCount++;
  return { xNext, yNext };
}

function updateInfo() {
  const info = document.getElementById('info');
  if (obsX.length === 0) {
    info.textContent = 'Press "Next evaluation" to start.';
    return;
  }
  const bestIdx = obsY.indexOf(Math.min(...obsY));
  const bestX = obsX[bestIdx].toFixed(3);
  const bestY = obsY[bestIdx].toFixed(3);
  const trueMinStr = trueMin.toFixed(3);
  info.textContent = fmt('Eval #%evalCount% — best so far: f(%bestX%) = %bestY%  |  true min ≈ %trueMinStr% at x ≈ %trueMinX%', { evalCount: evalCount, bestX: bestX, bestY: bestY, trueMinStr: trueMinStr, trueMinX: trueMinX.toFixed(3) });
}

document.getElementById('btn-next').addEventListener('click', () => {
  pickNext();
  drawAll();
  updateInfo();
});

document.getElementById('btn-reset').addEventListener('click', () => {
  obsX = []; obsY = []; evalCount = 0;
  drawAll();
  updateInfo();
});

drawAll();
updateInfo();

// Self-test (runs at load, not shown to user): verify GP and EI are working.
(function selfTest() {
  // With a known observation, GP mean at that x should be close to observed y.
  const tx = [2.0], ty = [trueF(2.0)];
  const pred = gpPredict(tx, ty, [2.0]);
  const meanAt2 = pred[0].mean;
  const stdAt2  = pred[0].std;
  // mean at observed point should be very close to observed value
  if (Math.abs(meanAt2 - ty[0]) > 0.5) {
    console.error('GP self-test FAILED: mean at obs point =', meanAt2, 'expected ~', ty[0]);
  } else {
    console.log('GP self-test passed: mean at obs =', meanAt2.toFixed(4), '(expected ~', ty[0].toFixed(4) + ')');
  }
  // Std at observed point should be low (noise-conditioned)
  if (stdAt2 > 0.5) {
    console.error('GP self-test FAILED: std at obs point too large =', stdAt2);
  } else {
    console.log('GP self-test passed: std at obs =', stdAt2.toFixed(4), '(expected low)');
  }
  // EI at a new point should be > 0 when we can improve
  const testEI = ei(-5, 0.5, ty[0]); // mean << fBest => big EI
  if (testEI <= 0) {
    console.error('EI self-test FAILED: expected positive EI, got', testEI);
  } else {
    console.log('EI self-test passed: EI =', testEI.toFixed(4));
  }
})();
Result
Console

Click Next evaluation to let the algorithm pick the most promising point (highest Expected Improvement) and sample the true function there. Watch how the orange surrogate gets closer to the blue truth as evaluations accumulate — and how the algorithm stops wasting effort in regions it has already mapped.

The Real Complexity

Bayesian optimization is not a hard open problem — it is a well-understood, practical algorithm. Its status: a solved heuristic with provable regret bounds (Srinivas et al., 2010).

Two costs dominate:

  • Surrogate fitting: the standard Gaussian process (GP) stores an n × n covariance matrix and inverts it at every step. That inversion is O(n3)O(n^{3}) in the number of past evaluations. In practice this caps usage at a few thousand points; beyond that, sparse or deep-learning surrogates are substituted.
  • Acquisition function optimization: choosing the next point requires maximizing a cheap closed-form function (like Expected Improvement or Upper Confidence Bound). Because this inner loop is cheap and smooth, standard gradient methods handle it efficiently.

The total evaluation budget n is usually tiny (tens to a few hundred) because each real evaluation is expensive. So the O(n3)O(n^{3}) GP cost is usually irrelevant compared to the cost of running the actual experiment.

Regret bounds give theoretical guarantees: with a squared-exponential kernel, the cumulative gap from the true optimum grows only as O(√(n log n)), meaning the algorithm converges, just not at an exponential rate.

Compare to alternatives: grid search wastes exponentially many evaluations in high dimensions; random search is often surprisingly competitive but has no memory; gradient descent needs a differentiable function. Bayesian optimization sits in the sweet spot when evaluations are costly and the budget is small.

Related ideas appear in non-convex optimization and the bandit framework used by recommendation systems.

Where It Matters

Any domain where each experiment is expensive and the search space is moderate is a natural home for Bayesian optimization:

  • Machine learning hyperparameters: learning rate, architecture depth, regularization strength. Systems like Google Vizier, Optuna, and Ax are built on it. A single well-placed run can beat a dozen randomly chosen ones.
  • Drug discovery and molecular design: each lab synthesis is costly. BO guides which molecule to test next, steering toward desirable properties like binding affinity or solubility.
  • Materials science and chemistry: finding alloy compositions or battery electrolyte formulas. Each physical experiment can take days; BO reduces the trial count dramatically.
  • Robotics policy search: tuning the control parameters of a physical robot without breaking it. Each rollout is a real-world experiment.
  • A/B testing and product optimization: when running a new variant is costly or risky, BO selects the most informative variant to test rather than a random one.
  • Simulation calibration: fitting simulator parameters to match real-world observations, where each simulation run is slow.

The field sits at the intersection of statistics (Gaussian processes, Bayesian inference) and optimization (acquisition maximization). Understanding it unlocks a family of "smart experiment design" methods used across science and engineering.

Conclusion

Bayesian optimization is the answer to a deceptively simple question: how do you find the best configuration when you can only afford a handful of experiments?

The answer is elegant: keep a running model of what you believe the function looks like, express uncertainty honestly, and always spend your next experiment where the payoff is highest given what you already know. The Gaussian process surrogate does the bookkeeping; the acquisition function makes the decision.

The algorithm is practical today in everything from AutoML pipelines to drug discovery labs. Its costs are real — the O(n3)O(n^{3}) GP and the need for a sensible prior — but for the regimes it was designed for, nothing beats it.

Next time you see a well-tuned neural network that took only twenty training runs to configure, there is a good chance Bayesian optimization was the quiet strategist behind the scenes.

Share this article

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

Comments

Loading comments...

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