Gaussian Processes for Modeling and Optimization
# Gaussian Processes for Modeling and Optimization
## Introduction & Motivation
Gaussian processes (GPs) provide non-parametric Bayesian approaches for modeling complex functions with uncertainty quantification. Critical for expensive simulations and experiments, GPs enable efficient exploration through Bayesian optimization, finding optima with minimal evaluations.
Motivation: Use GPs for sample-efficient modeling and optimization.
Applications: Surrogate modeling, Bayesian optimization, hyperparameter tuning, process design.
---
## Core Concepts & Theory
### Kernel Functions
Measuring similarity between inputs.
### Covariance Matrices
Capturing functional relationships.
### Mean Function
Prior expectation of function values.
### Marginal Likelihood
Model evidence for inference.
---
## Mathematical Formulation
Gaussian Process:
$$f(\mathbf{x}) \sim \mathcal{GP}(m(\mathbf{x}), k(\mathbf{x}, \mathbf{x}'))$$
Posterior Predictive:
$$p(f(\mathbf{x}^*)|D) = \mathcal{N}(\mu(\mathbf{x}^*), \sigma^2(\mathbf{x}^*))$$
Acquisition Function (EI):
$$\alpha_{EI}(\mathbf{x}) = (\mu(\mathbf{x}) - f_{best})\Phi(Z) + \sigma(\mathbf{x})\phi(Z)$$
---
## Advanced Theory & Extensions
### Kernel Composition
Building complex kernels.
### Sparse Approximations
Efficient GPs for large datasets.
### Multi-Task Learning
Transfer across related tasks.
---
## Computational Considerations
Covariance Matrix: O(N³) Cholesky decomposition.
Prediction: O(N²) per point.
Hyperparameter Optimization: O(N³·M) for M iterations.
---
## Practical Implementation Strategies
### Kernel Selection
Matern, RBF, and composite kernels.
### Hyperparameter Tuning
Maximum likelihood estimation.
### Bayesian Optimization
Acquisition function optimization.
---
## Benchmark Datasets & Evaluation
Optimization Benchmarks: Standard test functions.
Simulator Data: Physics-based models.
Experimental Records: Laboratory measurements.
---
## Key Challenges & Limitations
### Computational Cost
Scales poorly with dataset size.
### Hyperparameter Sensitivity
Model dependent on kernel choice.
### High Dimensions
Curse of dimensionality.
---
## Hyperparameter Tuning
RBF lengthscale: Problem-dependent.
Noise variance: 0.01-1.0.
Optimization iterations: 10-100.
---
## Real-World Applications & Case Studies
Materials Discovery: Optimal composition search.
Reactor Design: Temperature and pressure optimization.
Drug Development: Efficient screening.
---
## Integration with Other Methods
Gaussian processes + Bayesian optimization; + neural networks; + transfer learning.
---
## Summary & Key Takeaways
GPs enable efficient exploration and uncertainty quantification.
Principles:
1. Kernel: Define similarity measure.
2. Training: Fit GP to data.
3. Prediction: Posterior with uncertainty.
4. Optimization: Bayesian acquisition.
5. Adaptation: Sequential learning.
---
## Appendix: Practical Labs
### Lab 1: Kernel Functions
import numpy as np
class KernelFunctions:
@staticmethod
def rbf_kernel(x1, x2, length_scale=1.0, variance=1.0):
"""Radial basis function kernel"""
dist_sq = np.sum((x1[:, np.newaxis, :] - x2[np.newaxis, :, :]) ** 2, axis=2)
return variance * np.exp(-dist_sq / (2 * length_scale ** 2))
@staticmethod
def matern_kernel(x1, x2, nu=2.5, length_scale=1.0):
"""Matern kernel"""
dist = np.sqrt(np.sum((x1[:, np.newaxis, :] - x2[np.newaxis, :, :]) ** 2, axis=2))
if nu == 1.5:
return (1 + np.sqrt(3)*dist/length_scale) * np.exp(-np.sqrt(3)*dist/length_scale)
else:
# Simplified for nu=2.5
scaled_dist = np.sqrt(5) * dist / length_scale
return (1 + scaled_dist + scaled_dist**2/3) * np.exp(-scaled_dist)
@staticmethod
def linear_kernel(x1, x2, variance=1.0):
"""Linear kernel"""
return variance * np.dot(x1, x2.T)
# Test
X1 = np.array([[0], [1], [2]])
X2 = np.array([[0.5], [1.5]])
K_rbf = KernelFunctions.rbf_kernel(X1, X2, length_scale=0.5)
K_matern = KernelFunctions.matern_kernel(X1, X2)
print(f"✓ RBF kernel matrix shape: {K_rbf.shape}")
print(f" Values: {K_rbf[0]}")
print(f"✓ Matern kernel matrix: {K_matern[0]}")### Lab 2: Gaussian Process Regression
import numpy as np
class GaussianProcessRegression:
def __init__(self, length_scale=1.0, noise_variance=0.01):
self.length_scale = length_scale
self.noise_var = noise_variance
self.X_train = None
self.y_train = None
self.K_inv = None
def rbf_kernel(self, x1, x2):
"""RBF kernel matrix"""
dist_sq = np.sum((x1[:, np.newaxis, :] - x2[np.newaxis, :, :]) ** 2, axis=2)
return np.exp(-dist_sq / (2 * self.length_scale ** 2))
def fit(self, X, y):
"""Fit GP to data"""
self.X_train = X
self.y_train = y
# Compute covariance matrix
K = self.rbf_kernel(X, X)
K += np.eye(len(X)) * self.noise_var
self.K_inv = np.linalg.inv(K)
def predict(self, X_test):
"""Predict on test data"""
K_test = self.rbf_kernel(self.X_train, X_test)
K_test_test = self.rbf_kernel(X_test, X_test)
# Posterior mean
mu = K_test.T @ self.K_inv @ self.y_train
# Posterior variance
sigma_sq = np.diag(K_test_test) - np.sum(K_test * (self.K_inv @ K_test), axis=0)
sigma_sq = np.maximum(sigma_sq, 0)
return mu, np.sqrt(sigma_sq)
# Test
X_train = np.array([[0], [1], [2], [3], [4]])
y_train = np.sin(X_train.flatten()) + np.random.randn(5) * 0.1
gp = GaussianProcessRegression(length_scale=0.5)
gp.fit(X_train, y_train)
X_test = np.array([[0.5], [1.5], [2.5]])
mu, sigma = gp.predict(X_test)
print(f"✓ GP predictions:")
for x, m, s in zip(X_test.flatten(), mu, sigma):
print(f" f({x:.1f}) = {m:.3f} ± {s:.3f}")### Lab 3: Bayesian Optimization
import numpy as np
class BayesianOptimization:
def __init__(self, n_initial=3, n_iterations=10):
self.n_initial = n_initial
self.n_iterations = n_iterations
self.X_obs = None
self.y_obs = None
def acquisition_ei(self, X_test, gp_model, y_best):
"""Expected improvement acquisition"""
mu, sigma = gp_model.predict(X_test)
# Avoid division by zero
sigma = np.maximum(sigma, 1e-9)
# EI
Z = (mu - y_best) / sigma
ei = (mu - y_best) * scipy_norm_cdf(Z) + sigma * scipy_norm_pdf(Z)
return ei
def scipy_norm_cdf(x):
"""Approximate normal CDF"""
return 0.5 * (1 + np.tanh(0.7 * x))
def scipy_norm_pdf(x):
"""Approximate normal PDF"""
return np.exp(-x**2/2) / np.sqrt(2*np.pi)
class SimpleGPForBO:
def __init__(self):
self.X = None
self.y = None
def fit(self, X, y):
self.X = X
self.y = y
def predict(self, X_test):
# Simplified: return mean and constant std
mu = np.mean(self.y)
sigma = np.std(self.y) / (1 + np.linalg.norm(X_test - self.X, axis=1))
return np.ones(len(X_test)) * mu, np.maximum(sigma, 0.1)
def objective_function(x):
"""Function to optimize"""
return -(x[0] - 2)**2 - (x[0] - 3) + 5
# Bayesian optimization
X_initial = np.random.uniform(0, 5, (3, 1))
y_initial = np.array([objective_function(x) for x in X_initial])
X_all = X_initial.copy()
y_all = y_initial.copy()
gp = SimpleGPForBO()
for _ in range(5):
gp.fit(X_all, y_all)
# Find next point by EI
X_candidates = np.random.uniform(0, 5, (50, 1))
mu, sigma = gp.predict(X_candidates)
ei = (mu - y_all.max()) * scipy_norm_cdf((mu - y_all.max()) / (sigma + 1e-6))
next_idx = np.argmax(ei)
x_next = X_candidates[next_idx]
y_next = objective_function(x_next)
X_all = np.vstack([X_all, x_next])
y_all = np.append(y_all, y_next)
print(f"✓ Bayesian optimization:")
print(f" Best x: {X_all[np.argmax(y_all)][0]:.2f}")
print(f" Best y: {y_all.max():.2f}")### Lab 4: GP-Based Process Optimization
import numpy as np
class GPProcessOptimizer:
def __init__(self, process_simulator):
self.simulator = process_simulator
self.X_observed = []
self.y_observed = []
def run_experiment(self, parameters):
"""Run process and record outcome"""
outcome = self.simulator(parameters)
self.X_observed.append(parameters)
self.y_observed.append(outcome)
return outcome
def sequential_optimization(self, n_rounds=10):
"""Optimize process parameters sequentially"""
# Initial experiments
for _ in range(3):
params = np.random.uniform([0, 0], [10, 10], 2)
self.run_experiment(params)
# Iterative optimization
for round in range(n_rounds):
# Find best parameters so far
best_idx = np.argmax(self.y_observed)
best_value = self.y_observed[best_idx]
# Generate candidates
candidates = np.random.uniform([0, 0], [10, 10], (20, 2))
# Simple selection: try nearest to best
distances = np.linalg.norm(candidates - self.X_observed[best_idx], axis=1)
next_candidate_idx = np.argmin(distances)
next_params = candidates[next_candidate_idx]
self.run_experiment(next_params)
return np.array(self.X_observed), np.array(self.y_observed)
def process_simulator(params):
"""Simulate chemical process"""
T, P = params
yield_pred = 50 + 0.1*T + 2*P - 0.01*T*P + np.random.randn() * 2
return np.clip(yield_pred, 0, 100)
optimizer = GPProcessOptimizer(process_simulator)
X_opt, y_opt = optimizer.sequential_optimization(n_rounds=7)
print(f"✓ GP optimization results:")
print(f" Best temperature: {X_opt[np.argmax(y_opt), 0]:.1f}°C")
print(f" Best pressure: {X_opt[np.argmax(y_opt), 1]:.1f} bar")
print(f" Best yield: {y_opt.max():.1f}%")---