Bayesian Optimization from scratch: a GP surrogate plus an acquisition function finds a black-box peak in a dozen evals

Yesterday I built a Gaussian Process that gives a calibrated uncertainty band. Today I turn that band into a decision. When every evaluation of f(x) is expensive — a training run, a lab experiment, an A/B test — you can't afford a grid search. Bayesian Optimization fits the GP as a cheap surrogate, then an acquisition function reads its mean and band to score every candidate by how promising it is, spends one real evaluation at the argmax, refits, and repeats. In a dozen steps it homes in on the peak while random search is still wandering. Every curve on the demo page is computed live with the same Cholesky-based GP from yesterday. Here's the whole loop.

The expensive black-box and a few seed points

BO is for functions that are expensive to evaluate and whose formula you don't know. You can only afford a handful of calls, so start with two or three seed evaluations to give the surrogate something to fit.

import numpy as np
def objective(x):     # EXPENSIVE black-box; formula unknown to us
    return 1.7*np.exp(-(x-7)**2/1.6) + 1.15*np.exp(-(x-2.6)**2/2.5) - 0.15
bounds = (0.0, 10.0)
X = np.array([1.6, 8.4])                    # a few seed evaluations
y = np.array([objective(x) for x in X])     # the only truth we have

The surrogate is just the Gaussian Process

The surrogate is a cheap stand-in we can query for free. A GP is perfect: at every input it returns a posterior mean and a calibrated ±σ band. This is exactly the previous day's engine — kernel, three matrices, one Cholesky solve.

def gp_posterior(X, y, Xs, noise=1e-3):
    K   = rbf(X, X) + noise*np.eye(len(X))       # train-train
    Ks  = rbf(Xs, X)                             # test-train
    Kss = rbf(Xs, Xs)                            # test-test
    L     = np.linalg.cholesky(K)                # K = L Lᵀ
    alpha = np.linalg.solve(L.T, np.linalg.solve(L, y))
    mean  = Ks @ alpha                           # posterior mean
    v     = np.linalg.solve(L, Ks.T)
    std   = np.sqrt(np.clip(np.diag(Kss) - np.sum(v*v, 0), 0, None))
    return mean, std                             # the surrogate's belief

Expected Improvement — the acquisition function

An acquisition function scores every candidate x by how worth sampling it is, folding mean and uncertainty into one number. Expected Improvement asks: how much do we expect to beat the best value so far, f⁺? It rewards a high mean (exploit) and a wide band (explore) at once, and drops to zero at points already measured — so it never wastes an evaluation.

from scipy.stats import norm
def expected_improvement(mean, std, f_best, xi=0.01):
    std = np.maximum(std, 1e-9)
    imp = mean - f_best - xi          # xi nudges toward more exploration
    z   = imp / std
    return imp*norm.cdf(z) + std*norm.pdf(z)   # E[max(0, f - f_best)]

UCB is the same idea with one explicit dial: μ + κ·σ, an optimist's score where small κ exploits the current best and large κ chases the widest bands. Both are cheap, so you can maximise them over a dense grid.

The loop: fit → maximise acquisition → evaluate → repeat

This is the entire method. Fit the surrogate to everything seen so far, score a dense grid of candidates, pick the acquisition's argmax, spend one real (expensive) evaluation there, append it, go again. A dozen iterations is usually plenty — exactly the loop the demo runs on every click.

Xs = np.linspace(*bounds, 400)                 # a fine candidate grid
for t in range(12):                            # a handful of expensive evals
    mean, std = gp_posterior(X, y, Xs)         # 1. fit surrogate to data so far
    acq       = expected_improvement(mean, std, y.max())   # 2. score candidates
    x_next    = Xs[np.argmax(acq)]             # 3. pick the acquisition's argmax
    y_next    = objective(x_next)              # 4. spend ONE real evaluation
    X = np.append(X, x_next); y = np.append(y, y_next)     # 5. repeat

Why it beats grid and random search

The demo runs BO against a random search with the exact same budget and the exact same seed points, and plots the best value found per evaluation. Both climb toward the true maximum, but BO gets there in far fewer evaluations because it spends every one where its uncertainty-aware model says the payoff is highest — random just guesses. That gap is the value of a surrogate: after each new point, the band collapses around it and the acquisition re-aims.

from skopt import gp_minimize          # scikit-optimize; BoTorch/Ax for torch/GPU
res = gp_minimize(lambda x: -objective(x[0]), [bounds], n_calls=15, acq_func="EI")
print(-res.fun, res.x)                 # best value, best x — in ~15 evals

scikit-optimize's gp_minimize (or BoTorch/Ax on the GPU) does all of this and tunes the kernel too. Reach for BO when each evaluation is genuinely expensive and the budget is small: hyperparameter tuning, materials and drug screening, robotics, A/B tests. The cost is the GP's O(N³) Cholesky, so for high-dimensional or batched settings you move to Thompson sampling, TuRBO, or q-EI. But the core is small: a surrogate for the band, an acquisition to turn the band into a decision, and that explore-vs-exploit balance is the whole game.

Press "Sample next" and watch the hidden objective get found:
https://dev48v.infy.uk/ml/day38-bayesian-optimization.html