Bayesian Optimization
Bayesian optimization (BO) is a sample-efficient strategy for finding the optimum of expensive black-box functions. It maintains a probabilistic surrogate model of the objective and uses an acquisition function to decide where to evaluate next.
1. Problem Formulation
1.1 Black-Box Optimization
Consider optimizing a function where:
- is expensive to evaluate (e.g., training a neural network)
- may be noisy
- We want to find with few evaluations
1.2 Bayesian Optimization Loop
At each iteration :
- Surrogate model: Update posterior given data
- Acquisition function: Compute that balances exploration and exploitation
- Evaluate:
- Observe:
- Update:
2. Surrogate Models
2.1 Gaussian Processes
The most common surrogate model. A GP places a prior over functions:
Posterior: After observing :
where:
2.2 Random Forests (SMAC)
The Sequential Model-based Algorithm Configuration (SMAC) uses random forests:
where are individual regression trees. Variance is estimated by:
2.3 Tree-structured Parzen Estimator (TPE)
TPE models and separately:
- : modeled with Parzen estimators using kernel density estimation
- : modeled similarly
- : split at quantile
The acquisition function is derived from Bayes' rule:
3. Acquisition Functions
3.1 Expected Improvement (EI)
The most popular acquisition function. Let and be the GP posterior mean and standard deviation, and be the current best.
Improvement:
Expected Improvement:
where , and are the standard normal CDF and PDF.
Properties:
- Balances exploration (high ) and exploitation (low )
- Computationally efficient
- Can be optimized analytically in some cases
3.2 Upper Confidence Bound (UCB)
where is a parameter controlling exploration. For minimization, use the negative.
Theoretical guarantees: With , UCB achieves regret:
where is the maximum information gain.
3.3 Probability of Improvement (PI)
Limitation: Tends to get stuck in local optima because it only considers probability, not magnitude of improvement.
3.4 Knowledge Gradient (KG)
The expected increase in the best observed value:
Advantage: Accounts for the value of information, not just immediate improvement. Disadvantage: Computationally expensive (requires optimization over in expectation).
3.5 Entropy Search (ES)
Maximizes the expected reduction in entropy of the location of the minimizer :
3.6 Predictive Entropy Search (PES)
Simpler alternative to ES that maximizes:
4. Theoretical Analysis
4.1 Regret Bounds
Define simple regret as and cumulative regret as .
Theorem (Srinivas et al., 2010): For GP-UCB with :
with probability at least , where:
is the maximum information gain.
4.2 Information Gain Bounds
For common kernels:
- RBF:
- Matérn:
- Periodic:
4.3 No-Free-Lunch
Theorem: No optimization algorithm can be universally better than random search for all objective functions.
Implication: BO works well when the objective has structure (smoothness, etc.) that matches the GP prior.
5. Multi-Fidelity Optimization
5.1 Setup
Consider objectives with multiple fidelities where fidelity has cost and evaluates a cheaper approximation .
5.2 Successive Halving
- Start with configurations at lowest fidelity
- Evaluate all, keep top half
- Increase fidelity, repeat
- Continue until one configuration remains
Complexity: evaluations at lowest fidelity
5.3 Hyperband
Combines successive halving with random search:
where , , and is the halving rate.
5.4 BOHB (Bayesian Optimization + Hyperband)
Combines BO with Hyperband:
- Use BO to suggest configurations
- Use Hyperband to evaluate them efficiently
- Use completed evaluations to update the BO model
6. High-Dimensional BO
6.1 Linear Embedding (REMBO)
Assume the objective depends on a low-dimensional subspace:
where with .
6.2 Additive Models
Decompose the objective:
where are subsets of dimensions.
6.3 Trust Regions
Restrict BO to a region around the current best:
Expand/contract based on success/failure.
7. Parallel and Distributed BO
7.1 Batch Acquisition Functions
To evaluate points in parallel:
This is the sequential greedy approach.
7.2 Kriging Believer
Use the posterior mean as a pseudo-observation:
then update the GP and select the next point.
7.3 Constant Liar
Assume a constant value for pending evaluations:
8. Code: Bayesian Optimization
import numpy as np
from scipy.optimize import minimize
from scipy.stats import norm
from scipy.linalg import cho_solve, cho_factor
class GaussianProcessForBO:
"""Gaussian Process for Bayesian Optimization."""
def __init__(self, length_scale=1.0, signal_variance=1.0, noise=1e-3):
self.length_scale = length_scale
self.signal_variance = signal_variance
self.noise = noise
self.X_train = None
self.y_train = None
self.L = None
self.alpha = None
def _kernel(self, X1, X2):
"""RBF kernel."""
sq_dists = np.sum(X1**2, axis=1, keepdims=True) + \
np.sum(X2**2, axis=1) - 2 * X1 @ X2.T
return self.signal_variance * np.exp(-0.5 * sq_dists / self.length_scale**2)
def fit(self, X, y):
"""Fit GP to data."""
self.X_train = X
self.y_train = y
K = self._kernel(X, X) + self.noise * np.eye(len(X))
self.L = np.linalg.cholesky(K)
self.alpha = cho_solve((self.L, True), y)
def predict(self, X, return_std=True):
"""Predict at new points."""
K_s = self._kernel(X, self.X_train)
K_ss = self._kernel(X, X)
mu = K_s @ self.alpha
if return_std:
v = np.linalg.solve(self.L, K_s.T)
var = np.diag(K_ss) - np.sum(v**2, axis=0)
var = np.maximum(var, 1e-10)
return mu, np.sqrt(var)
return mu
def log_marginal_likelihood(self):
"""Compute log marginal likelihood."""
y = self.y_train
n = len(y)
lml = -0.5 * y @ self.alpha - \
np.sum(np.log(np.diag(self.L))) - \
n / 2 * np.log(2 * np.pi)
return lml
class ExpectedImprovement:
"""Expected Improvement acquisition function."""
def __init__(self, gp, xi=0.01):
self.gp = gp
self.xi = xi # Exploration-exploitation trade-off
def __call__(self, X):
"""Compute EI at points X."""
mu, sigma = self.gp.predict(X, return_std=True)
# Current best
f_best = np.min(self.gp.y_train)
# Standardized improvement
z = (f_best - mu - self.xi) / sigma
sigma = np.maximum(sigma, 1e-10)
# EI formula
ei = sigma * (z * norm.cdf(z) + norm.pdf(z))
return ei
def negative(self, X):
"""Negative EI for minimization."""
return -self.__call__(X.reshape(1, -1))[0]
class UpperConfidenceBound:
"""Upper Confidence Bound acquisition function."""
def __init__(self, gp, beta=2.0):
self.gp = gp
self.beta = beta
def __call__(self, X):
"""Compute UCB at points X."""
mu, sigma = self.gp.predict(X, return_std=True)
# For minimization, use lower confidence bound
ucb = mu - self.beta * sigma
return ucb
def negative(self, X):
return -self.__call__(X.reshape(1, -1))[0]
class ProbabilityOfImprovement:
"""Probability of Improvement acquisition function."""
def __init__(self, gp, xi=0.01):
self.gp = gp
self.xi = xi
def __call__(self, X):
"""Compute PI at points X."""
mu, sigma = self.gp.predict(X, return_std=True)
f_best = np.min(self.gp.y_train)
z = (f_best - mu - self.xi) / sigma
sigma = np.maximum(sigma, 1e-10)
pi = norm.cdf(z)
return pi
def negative(self, X):
return -self.__call__(X.reshape(1, -1))[0]
class KnowledgeGradient:
"""Knowledge Gradient acquisition function."""
def __init__(self, gp):
self.gp = gp
def __call__(self, X):
"""Compute KG at points X."""
mu, sigma = self.gp.predict(X, return_std=True)
# Current best
f_best = np.min(self.gp.y_train)
# Expected improvement in best value
# Approximation: E[max(f_best, mu - sigma*z)] - f_best
# where z ~ N(0,1)
z = (f_best - mu) / sigma
sigma = np.maximum(sigma, 1e-10)
kg = sigma * (z * norm.cdf(z) + norm.pdf(z))
return kg
def negative(self, X):
return -self.__call__(X.reshape(1, -1))[0]
class BayesianOptimizer:
"""Bayesian Optimization main class."""
def __init__(self, objective, bounds, gp_params=None, acq_type='EI', n_initial=5):
"""
Args:
objective: Function to minimize
bounds: List of (min, max) for each dimension
gp_params: Dict of GP parameters
acq_type: 'EI', 'UCB', 'PI', or 'KG'
n_initial: Number of initial random points
"""
self.objective = objective
self.bounds = np.array(bounds)
self.d = len(bounds)
self.n_initial = n_initial
self.X = []
self.y = []
self.gp = GaussianProcessForBO(**(gp_params or {}))
self.acq_type = acq_type
self.history = []
def _random_sample(self, n=1):
"""Sample random points within bounds."""
X = np.random.uniform(self.bounds[:, 0], self.bounds[:, 1],
size=(n, self.d))
return X
def _initial_sample(self):
"""Generate initial samples."""
X = self._random_sample(self.n_initial)
for x in X:
y = self.objective(x)
self.X.append(x)
self.y.append(y)
self.history.append((x.copy(), y))
self.X = np.array(self.X)
self.y = np.array(self.y)
def _optimize_acquisition(self, n_restarts=10):
"""Optimize acquisition function."""
if self.acq_type == 'EI':
acq = ExpectedImprovement(self.gp, xi=0.01)
elif self.acq_type == 'UCB':
acq = UpperConfidenceBound(self.gp, beta=2.0)
elif self.acq_type == 'PI':
acq = ProbabilityOfImprovement(self.gp, xi=0.01)
elif self.acq_type == 'KG':
acq = KnowledgeGradient(self.gp)
else:
raise ValueError(f"Unknown acquisition type: {self.acq_type}")
best_x = None
best_acq = -np.inf
for _ in range(n_restarts):
# Random restart
x0 = self._random_sample(1)[0]
# Optimize
result = minimize(acq.negative, x0, method='L-BFGS-B',
bounds=self.bounds)
if -result.fun > best_acq:
best_acq = -result.fun
best_x = result.x
return best_x
def optimize(self, n_iterations=50):
"""Run Bayesian optimization."""
# Initial sampling
self._initial_sample()
for i in range(n_iterations):
# Fit GP
self.gp.fit(self.X, self.y)
# Optimize acquisition
x_next = self._optimize_acquisition()
# Evaluate objective
y_next = self.objective(x_next)
# Update data
self.X = np.vstack([self.X, x_next.reshape(1, -1)])
self.y = np.append(self.y, y_next)
self.history.append((x_next.copy(), y_next))
print(f"Iteration {i+1}/{n_iterations}: "
f"f = {y_next:.4f}, best = {np.min(self.y):.4f}")
# Final fit
self.gp.fit(self.X, self.y)
best_idx = np.argmin(self.y)
return self.X[best_idx], self.y[best_idx]
class ParallelBayesianOptimizer:
"""Parallel Bayesian Optimization with batch acquisition."""
def __init__(self, objective, bounds, batch_size=4, acq_type='EI'):
self.objective = objective
self.bounds = np.array(bounds)
self.batch_size = batch_size
self.acq_type = acq_type
self.X = []
self.y = []
self.gp = GaussianProcessForBO()
def _sequential_greedy_batch(self):
"""Select batch using sequential greedy approach."""
batch = []
X_pending = self.X.copy()
y_pending = self.y.copy()
for _ in range(self.batch_size):
# Fit GP with pending points
self.gp.fit(X_pending, y_pending)
# Get acquisition function
if self.acq_type == 'EI':
acq = ExpectedImprovement(self.gp, xi=0.01)
else:
acq = UpperConfidenceBound(self.gp, beta=2.0)
# Optimize acquisition
best_x = None
best_acq = -np.inf
for _ in range(5):
x0 = np.random.uniform(self.bounds[:, 0], self.bounds[:, 1])
result = minimize(acq.negative, x0, method='L-BFGS-B',
bounds=self.bounds)
if -result.fun > best_acq:
best_acq = -result.fun
best_x = result.x
batch.append(best_x)
# Add to pending (with pseudo-observation)
y_pseudo = self.gp.predict(best_x.reshape(1, -1))[0]
X_pending = np.vstack([X_pending, best_x.reshape(1, -1)])
y_pending = np.append(y_pending, y_pseudo)
return np.array(batch)
def optimize(self, n_iterations=10):
"""Run parallel BO."""
# Initial random sampling
X_init = np.random.uniform(self.bounds[:, 0], self.bounds[:, 1],
size=(self.batch_size, self.d))
for x in X_init:
y = self.objective(x)
self.X.append(x)
self.y.append(y)
self.X = np.array(self.X)
self.y = np.array(self.y)
for i in range(n_iterations):
# Select batch
batch = self._sequential_greedy_batch()
# Evaluate batch
for x in batch:
y = self.objective(x)
self.X = np.vstack([self.X, x.reshape(1, -1)])
self.y = np.append(self.y, y)
print(f"Iteration {i+1}: best = {np.min(self.y):.4f}")
best_idx = np.argmin(self.y)
return self.X[best_idx], self.y[best_idx]
class HyperbandOptimizer:
"""Hyperband optimizer for multi-fidelity optimization."""
def __init__(self, objective, config_space, max_resource=81, eta=3):
"""
Args:
objective: Function(config, resource) -> loss
config_space: Dict of parameter ranges
max_resource: Maximum resource (e.g., epochs)
eta: Halving rate
"""
self.objective = objective
self.config_space = config_space
self.max_resource = max_resource
self.eta = eta
self.history = []
def _sample_config(self):
"""Sample random configuration."""
config = {}
for param, (low, high) in self.config_space.items():
if isinstance(low, int) and isinstance(high, int):
config[param] = np.random.randint(low, high + 1)
else:
config[param] = np.random.uniform(low, high)
return config
def _successive_halving(self, n_configs, resource, keep_ratio=1/self.eta):
"""Run successive halving."""
configs = [self._sample_config() for _ in range(n_configs)]
while len(configs) > 1 and resource <= self.max_resource:
# Evaluate all configs
losses = []
for config in configs:
loss = self.objective(config, resource)
losses.append(loss)
# Keep top fraction
n_keep = max(1, int(len(configs) * keep_ratio))
indices = np.argsort(losses)[:n_keep]
configs = [configs[i] for i in indices]
# Increase resource
resource = int(resource * self.eta)
return configs[0], min(losses)
def optimize(self, n_brackets=4):
"""Run Hyperband."""
best_config = None
best_loss = np.inf
for s in range(n_brackets, -1, -1):
n_configs = int(np.ceil(self.max_resource /
(self.max_resource * self.eta**(-s)) *
self.eta**s))
resource = int(self.max_resource * self.eta**(-s))
config, loss = self._successive_halving(n_configs, resource)
if loss < best_loss:
best_loss = loss
best_config = config
print(f"Bracket s={s}: n_configs={n_configs}, "
f"resource={resource}, best_loss={best_loss:.4f}")
return best_config, best_loss
class BOHB:
"""Bayesian Optimization + Hyperband."""
def __init__(self, objective, config_space, max_resource=81, eta=3):
self.objective = objective
self.config_space = config_space
self.max_resource = max_resource
self.eta = eta
self.X = [] # Configurations
self.y = [] # Losses
self.resources = [] # Resources used
self.gp = GaussianProcessForBO()
def _config_to_vector(self, config):
"""Convert config dict to vector."""
vec = []
for param in sorted(self.config_space.keys()):
low, high = self.config_space[param]
vec.append((config[param] - low) / (high - low))
return np.array(vec)
def _vector_to_config(self, vec):
"""Convert vector to config dict."""
config = {}
for i, param in enumerate(sorted(self.config_space.keys())):
low, high = self.config_space[param]
config[param] = low + vec[i] * (high - low)
return config
def _suggest_config(self):
"""Suggest configuration using BO."""
if len(self.X) < 10:
# Random exploration
config = {}
for param, (low, high) in self.config_space.items():
if isinstance(low, int) and isinstance(high, int):
config[param] = np.random.randint(low, high + 1)
else:
config[param] = np.random.uniform(low, high)
return config
# Fit GP
X = np.array([self._config_to_vector(x) for x in self.X])
y = np.array(self.y)
self.gp.fit(X, y)
# Optimize EI
acq = ExpectedImprovement(self.gp, xi=0.01)
best_x = None
best_acq = -np.inf
for _ in range(10):
x0 = np.random.uniform(0, 1, size=len(self.config_space))
result = minimize(acq.negative, x0, method='L-BFGS-B',
bounds=[(0, 1)] * len(self.config_space))
if -result.fun > best_acq:
best_acq = -result.fun
best_x = result.x
return self._vector_to_config(best_x)
def optimize(self, n_iterations=50):
"""Run BOHB."""
for i in range(n_iterations):
# Suggest configuration
config = self._suggest_config()
# Run with early stopping (simplified Hyperband)
resource = self.max_resource
keep_ratio = 1 / self.eta
while resource > 1:
loss = self.objective(config, resource)
if loss > np.percentile(self.y + [loss], 100 * keep_ratio):
break # Early stop
resource = max(1, int(resource / self.eta))
# Store result
self.X.append(config)
self.y.append(loss)
self.resources.append(resource)
print(f"Iteration {i+1}: loss={loss:.4f}, "
f"resource={resource}, best={np.min(self.y):.4f}")
best_idx = np.argmin(self.y)
return self.X[best_idx], self.y[best_idx]
# Example usage
if __name__ == "__main__":
# Define objective function
def objective(x):
return (x[0] - 0.5)**2 + (x[1] - 0.3)**2 + 0.1 * np.sin(10 * x[0])
bounds = [(0, 1), (0, 1)]
# Run Bayesian Optimization
bo = BayesianOptimizer(objective, bounds, acq_type='EI', n_initial=5)
best_x, best_y = bo.optimize(n_iterations=20)
print(f"\nBest solution: x = {best_x}, f(x) = {best_y}")
9. Summary
Bayesian optimization provides a principled framework for sample-efficient optimization:
- Surrogate models (GP, random forests, TPE) capture beliefs about the objective
- Acquisition functions (EI, UCB, PI, KG) balance exploration and exploitation
- Multi-fidelity methods (Hyperband, BOHB) leverage cheap approximations
- Theoretical guarantees (regret bounds) provide worst-case assurances
- Parallel extensions enable distributed optimization
The key insight is that by maintaining a probabilistic model and using information-theoretic acquisition functions, BO can find global optima with far fewer evaluations than gradient-based or random methods.