DEV Community

Cover image for Bayesian Optimization: The Prospector Who Reads Every Core Sample Before Choosing Where to Drill
Sachin Kr. Rajput
Sachin Kr. Rajput

Posted on

Bayesian Optimization: The Prospector Who Reads Every Core Sample Before Choosing Where to Drill

The One-Line Summary: Bayesian optimization fits a model of your loss surface and drills where that model is both promising and uncertain — and the advantage it bought was not mainly a better answer but a far more reliable one: my thirty-line version landed within 0.0004 across three seeds while random search spanned 0.0128, which is a lottery ticket versus a procedure.


The Parable of the Kestrel Basin Survey

The company has rights to eleven thousand hectares and enough money to drill forty holes. Each hole costs a week and returns one number: how much ore is in the core. There is no map, no seismic, nothing — just eleven thousand hectares and forty chances to find out where the deposit is.

The question the geologist has to answer before the first hole is not where the ore is. It is how she will decide where hole seventeen goes.


The Lottery

The first surveyor's method is to scatter the forty holes at random across the whole basin. It is unbiased, it covers the ground, and it requires no thinking between holes.

THE LOTTERY
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
Forty holes, chosen before any drilling starts:

  hole  1   NW quadrant   0.4 g/t
  hole  2   SE edge       0.1 g/t
  hole  3   NW quadrant   3.1 g/t   <- interesting!
  hole  4   S basin       0.2 g/t
  ...
  hole 40   E ridge       0.3 g/t

  best found: 3.1 g/t at hole 3

Hole 4 was chosen before hole 3 came back. So were
holes 5 through 40. The 3.1 reading changed nothing
about where anyone drilled next.

✗ Thirty-seven holes were sited by a plan that
  could not have known about hole 3.

This is called RANDOM SEARCH.
Enter fullscreen mode Exit fullscreen mode

Every hole is independent of every other. The information arrives and is filed, and the drill goes where the envelope said it would go.


The Prospector

The second surveyor drills eight holes at random to get a feel for the ground, and then stops. She builds a rough model of the whole basin from those eight readings — a guess at the grade everywhere, plus a guess at how unsure she is everywhere. Then she drills hole nine wherever that model says the upside is largest.

"I am not looking for the best hole I have drilled. I am looking for the place my ignorance is most expensive."

THE PROSPECTOR
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
After 8 random holes, she fits a surface:

  location      predicted   uncertainty   score
  near hole 3      2.8          0.4        3.2
  far N ridge      0.9          2.1        3.0   <- !
  near hole 1      0.5          0.3        0.8

  score = predicted + k x uncertainty

She drills the far N ridge, where she predicts
almost nothing, because she barely knows anything
there and a surprise would be worth a lot.

Then she refits the surface with 9 points and asks
again. And again, thirty-one more times.

This is called BAYESIAN OPTIMIZATION.
Enter fullscreen mode Exit fullscreen mode

Two things are being traded against each other, and the constant k decides the exchange rate: exploitation (drill near the good reading) and exploration (drill where you know least). Set k to zero and she drills the same hill forty times. Set it huge and she reinvents the lottery.


What The Survey Actually Showed

Both surveyors were run three times, on three different random starts, and the numbers that came back were not the ones the company expected.

THREE SURVEYS EACH
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
The lottery, three attempts:
  0.8687   0.8747   0.8816     <- spread 0.0128

The prospector, three attempts:
  0.8763   0.8765   0.8768     <- spread 0.0004

On the third attempt the lottery WON, and by more
than the prospector's average margin.

But look at the spreads. The lottery's outcome is
a draw. The prospector's is a result.

✓ The prospector is not reliably better.
  The prospector is reliably.
Enter fullscreen mode Exit fullscreen mode

The company had been asking which surveyor finds more ore. The useful question was which surveyor they could plan around.


Why It Works

THE MATHEMATICS OF ASKING BETTER QUESTIONS
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
RANDOM SEARCH
  x_(n+1) ~ P(x)          drawn from a fixed prior
  → trial n+1 is independent of trials 1..n
  → information accumulates and is never used

BAYESIAN OPTIMIZATION
  1. fit a SURROGATE model  f-hat(x)  to all
     trials so far, returning mean AND spread
  2. maximise an ACQUISITION function over it,
     e.g. UCB:   mu(x) + k * sigma(x)
  3. evaluate the true objective there, once
  4. refit and repeat

  REQUIREMENTS
  1. The objective must be EXPENSIVE relative to
     fitting the surrogate. Otherwise just sample.
  2. The surrogate must report uncertainty, not
     only a prediction.
  3. Trials must be SEQUENTIAL. Parallelism costs
     you the thing that makes it work.
  4. Roughly 8-15 random trials first, or the
     surrogate is fitting noise.
Enter fullscreen mode Exit fullscreen mode

What Is Bayesian Optimization?

It is a loop with two models in it. The objective is the thing you actually care about and cannot afford to evaluate often — a cross-validated model score, a simulation, a physical experiment. The surrogate is a cheap statistical model of the objective, fitted to every evaluation so far, which must return both a prediction and an uncertainty.

The acquisition function turns those two numbers into a single score that says how much you want to try a candidate point. Upper Confidence Bound is the simplest honest one:

aUCB(x)=μ(x)+kσ(x) a_{\text{UCB}}(x) = \mu(x) + k \cdot \sigma(x)

Expected Improvement is the more common choice and it does the same job with a different shape — it integrates how much better than the current best a point might be, rather than adding a fixed multiple of the standard deviation.

The classical surrogate is a Gaussian process, which is where "Bayesian" comes from: the GP is a prior over functions, updated by observations into a posterior [1]. In practice the surrogate is often a tree ensemble or, in Optuna's default, a Tree-structured Parzen Estimator — which does not model p(yx)p(y|x) at all but instead models p(xy)p(x|y) , building separate densities over the good trials and the bad ones and sampling where their ratio is favourable [2].

Step by Step

  1. Evaluate ninitn_{\text{init}} random points. Eight to fifteen.
  2. Fit the surrogate on everything observed so far.
  3. Maximise the acquisition function over a large set of cheap candidates.
  4. Evaluate the true objective at the argmax. Once.
  5. Add it to the observations and go to 2.

Step 3 is the part that surprises people: you are running an optimisation inside an optimisation. That is fine, because the inner one is over a model that costs microseconds, and the outer one is over a model fit that costs seconds or hours.


TPE Against Random, Three Budgets

Optuna's TPE sampler against Optuna's random sampler, identical objective, five seeds each, so the spread is visible.

import warnings, numpy as np; warnings.filterwarnings("ignore")
import optuna; optuna.logging.set_verbosity(optuna.logging.WARNING)
from sklearn.datasets import make_friedman1
from sklearn.model_selection import train_test_split, KFold, cross_val_score
from sklearn.ensemble import GradientBoostingRegressor

X, y = make_friedman1(n_samples=800, noise=1.0, random_state=42)
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.3, random_state=42)
cv = KFold(3, shuffle=True, random_state=0)

def score(lr, depth, sub, mf):
    m = GradientBoostingRegressor(n_estimators=60, learning_rate=lr, max_depth=depth,
                                  subsample=sub, max_features=mf, random_state=0)
    return cross_val_score(m, Xtr, ytr, cv=cv, scoring='r2', n_jobs=-1).mean()

def objective(t):
    return score(t.suggest_float("lr", 0.02, 0.25, log=True),
                 t.suggest_int("depth", 2, 5),
                 t.suggest_float("sub", 0.6, 1.0),
                 t.suggest_float("mf", 0.5, 1.0))

def run(sampler_name, budget, seed):
    s = optuna.samplers.TPESampler(seed=seed) if sampler_name == "tpe" \
        else optuna.samplers.RandomSampler(seed=seed)
    st = optuna.create_study(direction="maximize", sampler=s)
    st.optimize(objective, n_trials=budget, show_progress_bar=False)
    return st.best_value

print("TPE vs RANDOM, equal budget, 5 seeds each")
print("-" * 62)
print(f"{'budget':>7}{'TPE mean':>11}{'TPE spread':>12}{'rand mean':>11}"
      f"{'rand spread':>13}{'diff':>9}")
for budget in (20, 50, 120):
    tpe = [run("tpe", budget, s) for s in range(5)]
    rnd = [run("rand", budget, s) for s in range(5)]
    print(f"{budget:>7}{np.mean(tpe):>11.4f}{max(tpe)-min(tpe):>12.4f}"
          f"{np.mean(rnd):>11.4f}{max(rnd)-min(rnd):>13.4f}"
          f"{np.mean(tpe)-np.mean(rnd):>+9.4f}")
Enter fullscreen mode Exit fullscreen mode
TPE vs RANDOM, equal budget, 5 seeds each
--------------------------------------------------------------
 budget   TPE mean  TPE spread  rand mean  rand spread     diff
     20     0.8737      0.0146     0.8704       0.0077  +0.0032
     50     0.8781      0.0074     0.8730       0.0125  +0.0051
    120     0.8803      0.0090     0.8741       0.0107  +0.0063

Enter fullscreen mode Exit fullscreen mode

TPE wins at all three budgets, and — this is the part that matters — the margin grows with budget: +0.0032, +0.0051, +0.0063.

That growth is the signature of a method that is actually learning. Yesterday's grid-versus-random comparison produced gaps that flipped sign with dimensionality and vanished into seed noise, because neither method used its own results. Here the gap widens monotonically, because every extra trial makes the surrogate better, which makes the next trial better chosen. Random search at 120 trials (0.8741) has barely improved on random search at 50 (0.8730), while TPE gained 0.0022 over the same stretch.

Being honest about the caveat: at budget 20 the gap is +0.0032 while TPE's own spread is 0.0146. One run at a small budget tells you nothing. The direction is consistent across all three budgets, which is the evidence — not the size of any single comparison.


From Scratch: A Surrogate and a Rule

Thirty lines. Fit a random forest to the trials so far, use the disagreement between its trees as the uncertainty, and go wherever mean plus 1.5 standard deviations is highest.

import warnings, numpy as np; warnings.filterwarnings("ignore")
from sklearn.datasets import make_friedman1
from sklearn.model_selection import train_test_split, KFold, cross_val_score
from sklearn.ensemble import GradientBoostingRegressor, RandomForestRegressor

X, y = make_friedman1(n_samples=800, noise=1.0, random_state=42)
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.3, random_state=42)
cv = KFold(3, shuffle=True, random_state=0)
BOUNDS = np.array([[np.log(0.02), np.log(0.25)], [2, 5], [0.6, 1.0], [0.5, 1.0]])

def score(v):
    return cross_val_score(GradientBoostingRegressor(
        n_estimators=60, learning_rate=float(np.exp(v[0])), max_depth=int(round(v[1])),
        subsample=float(v[2]), max_features=float(v[3]), random_state=0),
        Xtr, ytr, cv=cv, scoring='r2', n_jobs=-1).mean()

def sample(rng, n):
    return rng.uniform(BOUNDS[:, 0], BOUNDS[:, 1], size=(n, 4))

def bayes_opt(budget, seed, n_init=8, kappa=1.5, n_cand=400):
    """Fit a forest on the trials so far; go where mean + kappa*sd is highest."""
    rng = np.random.RandomState(seed)
    Xs = sample(rng, n_init); ys = np.array([score(v) for v in Xs])
    for _ in range(budget - n_init):
        sur = RandomForestRegressor(n_estimators=60, min_samples_leaf=2,
                                    random_state=0, n_jobs=-1).fit(Xs, ys)
        cand = sample(rng, n_cand)
        per_tree = np.stack([t.predict(cand) for t in sur.estimators_])
        ucb = per_tree.mean(axis=0) + kappa * per_tree.std(axis=0)
        nxt = cand[int(np.argmax(ucb))]
        Xs = np.vstack([Xs, nxt]); ys = np.append(ys, score(nxt))
    return ys.max()

def random_search(budget, seed):
    rng = np.random.RandomState(seed)
    return max(score(v) for v in sample(rng, budget))

print("FROM SCRATCH: FOREST SURROGATE + UCB vs RANDOM   (3 seeds, budget 50)")
print("-" * 68)
bo = [bayes_opt(50, s) for s in range(3)]
rs = [random_search(50, s) for s in range(3)]
for s, (b, r) in enumerate(zip(bo, rs)):
    print(f"  seed {s}   mine {b:.4f}   random {r:.4f}   diff {b-r:+.4f}")
print()
print(f"  mine   mean {np.mean(bo):.4f}   spread {max(bo)-min(bo):.4f}")
print(f"  random mean {np.mean(rs):.4f}   spread {max(rs)-min(rs):.4f}")
print(f"  mine - random = {np.mean(bo)-np.mean(rs):+.4f}")
print(f"  (optuna TPE at budget 50 scored 0.8781 mean over 5 seeds)")
Enter fullscreen mode Exit fullscreen mode
FROM SCRATCH: FOREST SURROGATE + UCB vs RANDOM   (3 seeds, budget 50)
--------------------------------------------------------------------
  seed 0   mine 0.8763   random 0.8687   diff +0.0076
  seed 1   mine 0.8765   random 0.8747   diff +0.0019
  seed 2   mine 0.8768   random 0.8816   diff -0.0048

  mine   mean 0.8765   spread 0.0004
  random mean 0.8750   spread 0.0128
Enter fullscreen mode Exit fullscreen mode

Read the seed rows first, because they are embarrassing and they are the point. On seed 2, random search beat my optimizer by 0.0048 — more than my average margin over it. If I had run one seed and it happened to be seed 2, the honest headline would have been that Bayesian optimization lost.

Now read the spreads. My three runs landed at 0.8763, 0.8765 and 0.8768 — a spread of 0.0004. Random search's three runs spanned 0.0128, thirty-two times wider.

That reframes the entire method. The mean advantage was +0.0015, which is nothing. The variance reduction was a factor of thirty-two, which is everything. Random search occasionally finds something excellent and occasionally finds something mediocre, and you cannot tell which happened without a baseline. The surrogate search converges to nearly the same answer every time, because after the initial random phase it is following the data rather than the dice.

That is a different value proposition from the one usually advertised. Bayesian optimization is not primarily a way to get a better model. It is a way to stop your tuning result from being a function of your seed.

And the honest scoreboard: Optuna's TPE reached 0.8781 at the same budget, against my 0.8765. Thirty lines gets you most of the mechanism and none of the years of tuning that went into a real sampler. Use the library.


When Does It Help Most?

Random search Bayesian optimization
Uses previous trials no yes
Parallelises perfectly poorly, it is sequential
Gain grows with budget no (0.8730 → 0.8741) yes (+0.0032 → +0.0063)
Run-to-run spread 0.0128 0.0004
Overhead per trial none a surrogate fit
Needs warm-up trials no 8–15
Good below ~20 trials yes not really

Reach for Bayesian optimization when each evaluation is genuinely expensive — minutes or more — when you have 30 or more trials of budget, when trials run sequentially anyway, and when you need the result to be reproducible rather than lucky.

Reach for random search when evaluations are cheap, when you have hundreds of cores and want them all busy, or when your budget is under about twenty trials and the surrogate will not have enough data to be worth its overhead.

The failure mode worth naming: if your objective is fast, the surrogate fit and the candidate scan can cost more wall clock than the evaluations they are meant to save. Bayesian optimization is an answer to expensive objectives, and it is a tax on cheap ones.


Quick Reference Card

BAYESIAN OPTIMIZATION: CHEAT SHEET
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
THE LOOP
  1. n_init random trials  (8-15)
  2. fit surrogate -> mean + uncertainty
  3. maximise acquisition, e.g. mu + k*sigma
  4. evaluate objective ONCE there
  5. refit, repeat

MEASURED (GBM on Friedman #1, 4 params)
  budget 20   TPE 0.8737   random 0.8704   +0.0032
  budget 50   TPE 0.8781   random 0.8730   +0.0051
  budget 120  TPE 0.8803   random 0.8741   +0.0063
     -> the margin GROWS with budget

  from scratch, budget 50, 3 seeds
     mine   0.8765 mean, spread 0.0004
     random 0.8750 mean, spread 0.0128
     -> 32x tighter, not 32x better

KNOBS
  n_init      8-15; fewer and you fit noise
  kappa       0 = greedy, huge = random search
  surrogate   GP (smooth, small d) / TPE / forest

WHEN TO USE
   objective costs minutes+, budget 30+
   sequential trials, want reproducibility
WHEN NOT
   cheap objective (surrogate overhead dominates)
   budget < 20
   hundreds of idle cores

OPTUNA
  optuna.create_study(sampler=TPESampler(seed=0))
Enter fullscreen mode Exit fullscreen mode

Key Takeaways

  1. It beat random at every budget tested — +0.0032, +0.0051 and +0.0063 at 20, 50 and 120 trials.

  2. The margin grows with budget, which is the signature of learning. Random search went 0.8730 to 0.8741 from 50 to 120 trials; TPE went 0.8781 to 0.8803.

  3. The real prize is variance, not mean — my from-scratch version spread 0.0004 across three seeds where random search spread 0.0128.

  4. A single seed can reverse the conclusion — on seed 2 random search beat my optimizer by 0.0048.

  5. Thirty lines gets you most of it — a forest surrogate plus UCB reached 0.8765 against Optuna TPE's 0.8781.

  6. But use the library — that remaining 0.0016 is years of sampler engineering you do not need to repeat.

  7. It is sequential by design, so it trades away the one thing random search is perfect at. With idle cores, random search may finish sooner in wall clock.

  8. kappa is the whole personality — zero makes it greedy and it re-drills one hill; enormous makes it random search with extra steps.


The One-Sentence Summary

Bayesian optimization replaces "where should I look next?" with "where is my ignorance most expensive?" by fitting a model of the loss surface and drilling where predicted value plus uncertainty is highest — and measured over three budgets its advantage over random search grew monotonically from +0.0032 to +0.0063 because every trial improves the model that chooses the next one, but the finding that actually changes how you should think about it is that my thirty-line version landed within 0.0004 across three seeds while random search spanned 0.0128, which means its real product is not a better answer but an answer that does not depend on luck.


What's Next?

  1. The 5-minute Optuna setup — tomorrow. Replacing your grid search tonight, with pruning.
  2. Nested cross-validation — tuning and reporting without burning the same data twice.
  3. Early stopping inside a search — killing hopeless trials at iteration 10 instead of 200.
  4. Multi-objective tuning — when you need accuracy and latency.

Follow me for the next article in the Hyperparameter Tuning series!


Let's Connect!

If the Kestrel Basin survey made the exploration-exploitation trade click, drop a heart!

Questions? Ask in the comments — I read and respond to every one.

What's your n_init? Mine is 10 and I arrived at it badly — I ran with 3 for a year, kept getting searches that fixated on whatever the first lucky trial found, and never connected the two until I plotted where the trials went. 🪨


References

[1] J. Snoek, H. Larochelle and R. P. Adams, Practical Bayesian Optimization of Machine Learning Algorithms (2012), NeurIPS 2012

[2] J. Bergstra, R. Bardenet, Y. Bengio and B. Kégl, Algorithms for Hyper-Parameter Optimization (2011), NeurIPS 2011


The thing I did not expect from these measurements is that the interesting quantity was the spread rather than the mean. We are trained to compare methods by average performance, and on average these two are close enough that you could argue either way with one seed each — which is exactly what most blog comparisons do. The variance is where the difference lives, and variance is invisible unless you deliberately run the same thing several times and resist the urge to report only the best one.


Send this to whoever is still running a 512-point grid overnight.

Top comments (0)