🤖 Bayesian Optimization Experiment Selection Engine

Seeding the Search — Latin Hypercube and Sobol Sequences Before Any Model Exists

Every Bayesian optimization campaign begins with a problem: the surrogate model needs data to fit, but no data exists yet. Rather than guessing where the optimum might be (introducing human bias) or sampling purely at random (which clusters points and leaves large gaps), self-driving labs use quasi-random space-filling designs to seed the domain with maximally diverse initial points.

• 2(d+1): Typical init budget (points, d = # parameters)
• Sobol / LHS: Common design (low-discrepancy sequences)
• Sobol: BoTorch/Ax default (quasi-Monte Carlo init)
• Latin Hypercube: GPyOpt default (stratified random design)

Why the first points matter more than they look

Latin Hypercube Sampling (LHS, McKay et al. 1979) divides each input dimension into N equal-probability strata and places exactly one sample per stratum per dimension, then randomly pairs strata across dimensions. This guarantees even one-dimensional coverage that pure random sampling cannot: with 10 random points in 1D, expect visible clustering and gaps; with 10 LHS points, every decile of the range is represented exactly once.

Sobol sequences (Sobol 1967) go further, using a deterministic low-discrepancy construction so that even partial prefixes of the sequence (the first 4 of a planned 20 points) remain well-spread — useful because campaigns are often extended on the fly.

Practical defaults used by real frameworks: • GPyOpt: Latin Hypercube by default, optionally random or grid • BoTorch / Ax (Meta): Sobol engine (torch.quasirandom.SobolEngine) for all initial designs, chosen for its scrambling properties and reproducibility • Typical initial batch size: 2×(d+1) to 5×d points for d input dimensions — enough to estimate a first set of GP hyperparameters without wasting the experimental budget before any model-guided search begins

In a real self-driving chemistry lab, each of these "seed" points is a physically executed reaction or synthesis — expensive and slow — so the design must extract maximum information per point. Sobol/LHS designs typically reduce the number of wasted or redundant initial evaluations by 30–50% relative to naive random sampling of the same budget.

Space-filling design is a one-time investment: it is used only until the very first surrogate model can be fit. After that, every subsequent point is chosen adaptively by the acquisition function — the initial design exists purely to give the Gaussian Process enough spread to estimate sensible length-scale and noise hyperparameters.

Gaussian Process Regression — Turning Five Points into a Continuous Belief

The Gaussian Process (GP) is the statistical engine at the heart of Bayesian optimization. Instead of predicting a single value at each candidate point, a GP predicts an entire probability distribution — a mean and a variance — making it possible to reason not just about what the response surface probably looks like, but about how confident that guess actually is.

• Matérn 5/2: Default kernel (BoTorch) (smooth but not infinitely so)
• Type-II MLE: Hyperparameter fit (maximize marginal likelihood)
• O(n³): Fit cost (Cholesky of n×n covariance)
• ~200: Typical n before refit lag (points before sparse GP needed)

Kernel choice and the marginal-likelihood fit

A GP is fully specified by a mean function (usually zero after centering) and a covariance kernel k(x,x′) that encodes how correlated two outputs are expected to be, given the distance between their inputs. The Matérn 5/2 kernel is the standard default in BoTorch and most modern self-driving-lab stacks:

k(x,x′) = σf² (1 + √5r/ℓ + 5r²/3ℓ²) exp(−√5r/ℓ), r = |x−x′|

• ℓ (length-scale): how far apart two points must be before the model treats them as roughly independent — larger ℓ means smoother, more global extrapolation • σf² (signal variance): overall vertical scale of function variation • σn² (noise variance): observation noise floor, added on the diagonal of the covariance matrix

Given training pairs {(xᵢ,yᵢ)}, the posterior mean and variance at a new point x* are:

μ(x*) = k*ᵀ (K + σn²I)⁻¹ y σ²(x*) = k(x*,x*) − k*ᵀ (K + σn²I)⁻¹ k*

where K is the n×n kernel matrix over training points and k* is the vector of kernel values between x* and each training point. Hyperparameters (ℓ, σf², σn²) are not guessed — they are fit by maximizing the log marginal likelihood, an automatic Occam's-razor criterion that trades off data fit against model complexity, exactly as GPyTorch/BoTorch do via L-BFGS on the negative log marginal likelihood.

The defining visual signature of a GP posterior: the uncertainty band pinches to near-zero exactly at observed points (the model has seen ground truth there) and widens smoothly in unexplored regions — this shape is what the acquisition function exploits in the next stage.

Expected Improvement and Upper Confidence Bound — Scoring Every Unmeasured Point

A GP posterior alone does not tell you what to do next — it describes belief, not action. The acquisition function converts that belief into a single scalar score at every candidate point, folding predicted performance and predictive uncertainty into one number that can be maximized to choose the next experiment.

• Jones et al. 1998: EI origin ("Efficient Global Optimization")
• Srinivas et al. 2010: UCB origin (GP-UCB regret bounds)
• Φ, φ: EI closed form (normal CDF / PDF terms)
• exploration margin: ξ role (higher ξ → more exploration)

Expected Improvement — the workhorse acquisition function

Expected Improvement (EI, Jones, Schonlau & Welch 1998) computes the expected amount by which a candidate point x will beat the current best observed value y⁺, under the GP posterior at x:

Z = (μ(x) − y⁺ − ξ) / σ(x) EI(x) = (μ(x) − y⁺ − ξ)·Φ(Z) + σ(x)·φ(Z), for σ(x) > 0 EI(x) = 0, for σ(x) = 0

Φ and φ are the standard normal CDF and PDF. The exploration parameter ξ (xi) sets a margin of required improvement: ξ≈0 makes EI almost purely exploitative (favoring points near the current best mean), while larger ξ (0.01–0.3 in practice) forces EI to reward uncertain, unexplored regions more heavily, since a wide σ(x) can only pay off the higher improvement bar through raw uncertainty.

Upper Confidence Bound (UCB, Srinivas et al. 2010, "GP-UCB") is the other dominant choice, favored for its provable regret bounds:

UCB(x) = μ(x) + κ·σ(x)

κ (kappa) directly and linearly trades off mean against uncertainty — small κ chases the best-known region, large κ pushes toward the least-explored region regardless of predicted value. Both EI and UCB are implemented as standard acquisition classes in GPyOpt, BoTorch (qExpectedImprovement, UpperConfidenceBound), and Ax's Bayesian optimization service API — the choice between them is often empirical, with EI slightly favored for noisy/expensive black-box chemistry and UCB favored when formal exploration guarantees matter.

Choosing the Argmax — From an Abstract Score to One Concrete Experiment

The acquisition surface is itself a continuous, usually multi-modal, non-convex function — finding its maximum is a genuine (if cheap) optimization problem in its own right. Whatever point maximizes it becomes the literal instruction sent to a robotic liquid handler, an autonomous synthesis platform, or a reaction-screening array.

• multi-start L-BFGS: Acquisition optimizer (or CMA-ES for rugged surfaces)
• 10–20: Restarts (BoTorch default) (from Sobol-sampled starts)
• qEI / qUCB: Batch selection (parallel multi-point picks)
• milliseconds–seconds: Selection cost (vs. hours/days per real experiment)

Optimizing the acquisition function itself

Maximizing EI(x) or UCB(x) over the input domain is cheap relative to running a real experiment (milliseconds to seconds vs. hours to days), so it is standard to over-invest in this inner optimization: BoTorch by default draws a large pool of candidate points via Sobol sampling, seeds multiple L-BFGS-B restarts from the most promising subset, and returns the best local optimum found — since acquisition surfaces are frequently multi-modal (several plausible "next best guesses" separated by explored regions).

For discrete or combinatorial parameter spaces (categorical reagents, fixed equipment settings), the acquisition function is instead evaluated exhaustively or via CMA-ES / genetic search over the finite candidate set.

Batch (parallel) selection: real labs often run several reactions simultaneously (multi-well plates, parallel reactors). Naively picking the top-k acquisition points tends to cluster them together (they all look good under the same posterior). Modern frameworks instead use batch-aware acquisition functions — qEI (Monte Carlo Expected Improvement over a joint batch, via BoTorch's q-family) or local penalization — that sequentially "fantasize" each pick's effect on the posterior before choosing the next, ensuring the batch is diverse rather than redundant.

Closing the Loop — Refitting the Model on Every New Measurement

Once the physical experiment at the selected point completes and a real (possibly noisy) measurement comes back, that pair is appended to the training set and the GP is refit — updating not just the local prediction near the new point but, through the shared kernel hyperparameters, the model's beliefs everywhere in the domain.

• seconds: Update latency (GP refit on <200 points)
• every N points: Hyperparameter refit (often not every single one)
• ~10³ points: Scaling limit (before sparse/inducing-point GPs)
• ChemOS, Ada, Emerald Cloud Lab: Real deployments (closed-loop BO platforms)

What actually changes when one point is added

Adding a single (x,y) pair changes three things simultaneously:

1. The n×n kernel matrix K grows to (n+1)×(n+1), requiring a fresh Cholesky decomposition — O(n³) cost that is trivial for tens to low hundreds of points but motivates sparse/inducing-point GP approximations once campaigns exceed roughly 10³ evaluations 2. The posterior mean and variance shift everywhere, not just locally — points near the new sample see the largest correction, but the global length-scale and noise hyperparameters (refit periodically via marginal-likelihood maximization) can shift the whole surface's smoothness assumption 3. The acquisition function is recomputed from scratch on the new posterior, and the entire cycle — fit, score, select, measure, update — repeats

In production self-driving labs this loop is genuinely closed: platforms such as ChemOS (Häse et al.), the Ada system, Emerald Cloud Lab's script-driven workflows, and Chemspeed/Opentrons-integrated BoTorch pipelines feed instrument results directly back into the optimizer with no manual re-entry, typically completing one full fit-acquire-select-measure-update cycle in the time it takes the physical experiment itself to run — meaning the model update step is essentially free relative to the bottleneck of real-world synthesis or measurement time.

Sample Efficiency — Why Bayesian Optimization Beats Grid and Random Search

The entire justification for the added statistical machinery of Bayesian optimization is sample efficiency: reaching a near-optimal answer using dramatically fewer physical experiments than exhaustive or random alternatives. In self-driving labs, where each evaluation can cost real time, reagents, and instrument hours, this efficiency is the whole economic argument for the approach.

• 15–30 evals: BO to within 1% of optimum (typical low-dimensional problems)
• 10²–10³ evals: Grid search, same accuracy (combinatorial explosion in d>2)
• Shields et al. 2021: Reaction yield optimization (Nature; beat expert chemists)
• kᵈ evaluations: Grid cost, d dims, k levels (curse of dimensionality)

Quantifying the efficiency gap

Grid search evaluates every combination of discretized parameter levels: for d independent parameters each discretized into k levels, the required evaluation count is kᵈ — 4 parameters at 10 levels each is already 10,000 evaluations, most of them wasted on clearly suboptimal regions the model never needed to visit. Random search improves slightly on grid search's worst-case redundancy but still ignores everything learned from prior evaluations — every point is chosen with zero regard for what has already been measured.

Bayesian optimization, by contrast, spends its early budget on space-filling exploration (Stage 1) and then rapidly concentrates subsequent evaluations in promising regions once the surrogate is confident enough to make that call. Shields et al. (Nature, 2021, "Bayesian reaction optimization as a tool for chemical synthesis") benchmarked BO against expert human chemists and standard search baselines across real reaction-yield optimization tasks (the Olympus/Summit benchmark suite) and found BO consistently reached near-optimal yields within 30–100 experiments, matching or beating expert intuition while requiring far fewer wet-lab runs than grid-based screening.

The practical rule of thumb widely cited in self-driving-lab literature: BO reaches within 1% of the true global optimum in roughly 15–30 evaluations for problems with up to ~5 continuous parameters, while grid or exhaustive search needs hundreds to thousands of evaluations for comparable coverage and accuracy — an efficiency gain that compounds further as experiment cost rises (multi-day syntheses, expensive reagents, limited instrument time).

Sample efficiency is not a marginal optimization — it is the difference between a discovery campaign that fits inside a PhD student's year and one that does not fit inside a career. A materials or reaction-condition search space with 6 continuous parameters at even coarse 10-level discretization implies 10⁶ grid evaluations; the same space is routinely searched to near-optimality by Bayesian optimization in well under 100 real experiments.