bayesian-optimization

name	bayesian-optimization
description	Bayesian optimization — Gaussian process surrogate, acquisition functions (EI/UCB/PI), kernel selection, hyperparameter tuning, black-box optimization, engineering design applications.
metadata	{"priority":7,"promptSignals":{"phrases":["Bayesian optimization","Gaussian process surrogate","acquisition function","expected improvement","black box optimization","surrogate model optimization","design space exploration"],"minScore":3}}

Bayesian Optimization — Complete Skill

Problem Statement

Optimize an expensive-to-evaluate black-box function f(x): x* = argmax f(x), x ∈ D ⊆ R^d

Key characteristics:

• f(x) is expensive to evaluate (e.g., CFD, FEA, physical experiments)
• No gradient information available (or unreliable)
• Function may be noisy: y_i = f(x_i) + ε_i, ε_i ~ N(0, σ_n²)
• Goal: find global optimum with minimum evaluations (typically 50–200 total)

Gaussian Process Surrogate Model

GP Definition: f(x) ~ GP(μ(x), k(x, x'))

μ(x) = prior mean function (often μ = 0 after normalization) k(x, x') = covariance kernel function

GP Posterior after observations D = {(x_i, y_i)}: μ_n(x) = k(x, X) [K + σ_n² I]⁻¹ y σ²_n(x) = k(x, x) - k(x, X) [K + σ_n² I]⁻¹ k(X, x)

K = [k(x_i, x_j)] = n×n kernel matrix k(x, X) = [k(x, x₁), ..., k(x, x_n)] = 1×n vector

Prediction: GP returns mean μ_n(x) and variance σ²_n(x) at any point

Kernel Functions (Covariance)

Squared Exponential (RBF/Gaussian) — most common

k(x, x') = σ_f² × exp(-||x - x'||²/(2ℓ²)) ℓ = length scale (controls smoothness); σ_f² = signal variance Infinitely differentiable → very smooth functions; may over-smooth engineering responses

Matérn Kernel (better for engineering)

Matérn 5/2: k(x,x') = σ_f²(1 + √5 r/ℓ + 5r²/(3ℓ²)) exp(-√5 r/ℓ) [r = ||x-x'||] Matérn 3/2: k(x,x') = σ_f²(1 + √3 r/ℓ) exp(-√3 r/ℓ) Twice differentiable (5/2) or once differentiable (3/2) → better for physical responses Recommendation: Matérn 5/2 for most engineering applications

Automatic Relevance Determination (ARD)

Different length scale per dimension: r² = Σ(x_d - x'_d)²/ℓ_d² Allows feature importance estimation; longer ℓ_d → dimension d less important

Composite Kernels

k = k₁ + k₂ (additive structure) k = k₁ × k₂ (product structure) Combine to model different behavior at different scales

Hyperparameter Learning

Optimize kernel hyperparameters θ = {σ_f, ℓ, σ_n} by maximizing log marginal likelihood: log p(y|X, θ) = -½ y^T (K + σ_n² I)⁻¹ y - ½ log|K + σ_n² I| - n/2 log(2π)

Optimize with L-BFGS: multiple restarts to avoid local optima (3–10 starts typical) Full Bayesian (MCMC): sample θ instead of point estimate; better uncertainty quantification

Acquisition Functions

Expected Improvement (EI) — most widely used

α_EI(x) = E[max(f(x) - f*, 0)] = (μ_n(x) - f*) × Φ(Z) + σ_n(x) × φ(Z) Z = (μ_n(x) - f*) / σ_n(x) f* = current best observed value

Properties: balances exploitation (high μ_n) and exploration (high σ_n) Handles noise: EI computed using GP uncertainty naturally

Upper Confidence Bound (UCB)

α_UCB(x) = μ_n(x) + κ × σ_n(x) κ = exploration parameter (typically 2–5) Higher κ → more exploration; lower → more exploitation Theoretical guarantees: sublinear regret bounds (Srinivas et al.)

Probability of Improvement (PI)

α_PI(x) = Φ((μ_n(x) - f* - ξ) / σ_n(x)) ξ = improvement threshold (usually 0.01–0.1) More exploitative than EI; may miss global optimum

Thompson Sampling

Sample f̃ ~ GP posterior; maximize f̃ Inherently stochastic; naturally parallel (sample multiple and evaluate simultaneously)

Knowledge Gradient (KG)

Optimal single-observation expected improvement for finite horizon Computationally expensive; best for noisy objectives

Optimization Loop

Initialize: x_1,...,x_n0 via Latin hypercube sampling
For i = n0+1, 2, ... until budget:
    1. Fit GP to current data D = {(x_j, y_j)}
    2. Optimize acquisition: x_i = argmax α(x; GP)
    3. Evaluate: y_i = f(x_i)  [expensive step]
    4. Update D = D ∪ (x_i, y_i)
Return x* = argmax_{x_j∈D} y_j

Inner optimization (acquisition function): Multi-start gradient (L-BFGS) or evolutionary (CMA-ES) for non-convex α(x) Acquisition is cheap to evaluate (closed form) → exhaustive inner optimization OK

Handling Constraints

Hard constraints (known feasibility function): constrained EI: multiply EI by P(feasible) P(feasible) = GP posterior probability that constraint satisfied

Soft constraints (expensive constraint evaluations): Multi-objective formulation: maximize f while satisfying g(x) ≤ 0 Infill criterion: constrained EI; PoF × EI

Known constraint regions: simply restrict search domain; re-parameterize if possible

Multi-Objective Bayesian Optimization

Optimize f₁(x), ..., f_k(x) simultaneously Hypervolume Expected Improvement (EHVI): expected increase in dominated hypervolume ParEGO: scalarize objectives with random weights; multiple Bayesian loops

Parallel / Batch Bayesian Optimization

Evaluate q points simultaneously (batch size q > 1) q-EI (joint expected improvement): q-EI = E[max(max{f(x₁),...,f(x_q)} - f*, 0)]

Kriging Believer / Constant Liar: sequential greedy; assign fantasy values to pending evaluations Local Penalization: penalize acquisition near already selected points

Engineering Applications

Hyperparameter tuning (ML): 50–200 evaluations vs. grid search (thousands) Material property optimization: alloy composition, heat treatment parameters CFD shape optimization: airfoil profile, turbine blade, nozzle Manufacturing process: welding parameters, heat treat, casting conditions Structural optimization: wall thicknesses, material selection Experimental design: lab conditions, reagent concentrations

Software Tools

BoTorch (Meta): PyTorch-based; GPU-accelerated; q-EI, multi-objective GPyOpt: scikit-learn compatible; simple API Spearmint: original Bayesian optimization package (Adams group) Ax (Meta): high-level platform; ABN test, ML hyperparameters, experiments SMAC3: Sequential Model-based Algorithm Configuration; handles categorical Dragonfly: high-dimensional; neural architecture search

Scalability Limits

Standard GP: O(n³) training, O(n²) prediction n > 1000 observations: use sparse GP (inducing points; GPy, GPflow) or random features High-dimensional (d > 20): dimensionality reduction (REMBO, ALEBO); additive models; neural network kernels (Deep GP)

Output

Provide: kernel selection (with rationale), initial DoE size n₀ (typically max(5d, 10)), acquisition function, convergence criterion (budget N_max or improvement threshold), expected reduction in f (from initial to optimum), recommended batch size q, constraint handling method, and software library suggestion for implementation.