Chapter 1: Fundamentals of Bayesian Optimization

1.1 Black-Box Optimization Problems

In process industries, there are many "black-box" problems where the relationship between inputs and outputs is complex and cannot be described analytically. For example, when searching for optimal operating conditions for a chemical reactor, the relationship between parameters such as temperature, pressure, and reaction time and the yield can only be evaluated through experiments or simulations.

Example 1: Formulation of Black-Box Problem (Chemical Reactor)

We define a problem that optimizes three parameters: temperature, pressure, and catalyst concentration, assuming a polymerization reaction process.

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm

# ===================================
# Example 1: Black-box optimization problem for chemical reactor
# ===================================

class ChemicalReactor:
    """Chemical reactor black-box simulator

    Parameters:
        - Temperature (T): 300-400 K
        - Pressure (P): 1-5 bar
        - Catalyst concentration (C): 0.1-1.0 mol/L

    Objective: Maximize yield
    """

    def __init__(self, noise_level=0.02):
        self.noise_level = noise_level
        self.bounds = np.array([[300, 400], [1, 5], [0.1, 1.0]])
        self.dim_names = ['Temperature (K)', 'Pressure (bar)', 'Catalyst (mol/L)']
        self.optimal_x = np.array([350, 3.0, 0.5])  # True optimal solution

    def evaluate(self, x):
        """Calculate yield (equivalent to actual experiment/simulation)

        Args:
            x: [temperature, pressure, catalyst concentration]

        Returns:
            yield: Yield [0-1] + noise
        """
        T, P, C = x

        # Temperature dependency (Arrhenius type)
        T_opt = 350
        temp_factor = np.exp(-((T - T_opt) / 30)**2)

        # Pressure dependency (parabolic type)
        P_opt = 3.0
        pressure_factor = 1 - 0.3 * ((P - P_opt) / 2)**2

        # Catalyst concentration dependency (Langmuir type)
        catalyst_factor = C / (0.2 + C)

        # Interaction term (synergistic effect of temperature and pressure)
        interaction = 0.1 * np.sin((T - 300) / 50 * np.pi) * (P - 1) / 4

        # Yield calculation
        yield_val = 0.85 * temp_factor * pressure_factor * catalyst_factor + interaction

        # Add noise (measurement error)
        noise = np.random.normal(0, self.noise_level)

        return float(np.clip(yield_val + noise, 0, 1))

    def plot_landscape(self, fixed_catalyst=0.5):
        """Visualize objective function (catalyst concentration fixed)"""
        T_range = np.linspace(300, 400, 50)
        P_range = np.linspace(1, 5, 50)
        T_grid, P_grid = np.meshgrid(T_range, P_range)

        Y_grid = np.zeros_like(T_grid)
        for i in range(len(T_range)):
            for j in range(len(P_range)):
                Y_grid[j, i] = self.evaluate([T_grid[j, i], P_grid[j, i], fixed_catalyst])

        fig = plt.figure(figsize=(10, 4))

        # 3D surface
        ax1 = fig.add_subplot(121, projection='3d')
        surf = ax1.plot_surface(T_grid, P_grid, Y_grid, cmap=cm.viridis, alpha=0.8)
        ax1.set_xlabel('Temperature (K)')
        ax1.set_ylabel('Pressure (bar)')
        ax1.set_zlabel('Yield')
        ax1.set_title('Chemical Reactor Response Surface')

        # Contour
        ax2 = fig.add_subplot(122)
        contour = ax2.contourf(T_grid, P_grid, Y_grid, levels=20, cmap='viridis')
        ax2.set_xlabel('Temperature (K)')
        ax2.set_ylabel('Pressure (bar)')
        ax2.set_title('Yield Contour (Catalyst=0.5 mol/L)')
        plt.colorbar(contour, ax=ax2, label='Yield')

        plt.tight_layout()
        return fig

# Usage example
reactor = ChemicalReactor(noise_level=0.02)

# Evaluate at initial conditions
x_initial = np.array([320, 2.0, 0.3])
yield_initial = reactor.evaluate(x_initial)
print(f"Initial conditions: T={x_initial[0]}K, P={x_initial[1]}bar, C={x_initial[2]}mol/L")
print(f"Initial yield: {yield_initial:.3f}")

# Evaluate near optimal conditions
x_optimal = np.array([350, 3.0, 0.5])
yield_optimal = reactor.evaluate(x_optimal)
print(f"\nOptimal conditions: T={x_optimal[0]}K, P={x_optimal[1]}bar, C={x_optimal[2]}mol/L")
print(f"Optimal yield: {yield_optimal:.3f}")

# Visualization
fig = reactor.plot_landscape()
plt.show()

1.2 Sequential Design Strategy

The core of Bayesian optimization lies in "Sequential Design". By leveraging previous observation results, it enables efficient optimization by intelligently selecting the next point to evaluate.

Example 2: Sequential Design vs Random Sampling

We demonstrate how optimization performance greatly differs depending on strategy, even with the same number of evaluations.

import numpy as np
import matplotlib.pyplot as plt

# ===================================
# Example 2: Demonstration of sequential design strategy
# ===================================

def simple_objective(x):
    """1D test function (analytical solution available)"""
    return -(x - 3)**2 * np.sin(5 * x) + 2

class SequentialDesigner:
    """Optimization by sequential design"""

    def __init__(self, objective_func, bounds, n_initial=3):
        self.objective = objective_func
        self.bounds = bounds
        self.X_observed = []
        self.Y_observed = []

        # Initial points (Latin Hypercube Sampling)
        np.random.seed(42)
        for _ in range(n_initial):
            x = np.random.uniform(bounds[0], bounds[1])
            y = objective_func(x)
            self.X_observed.append(x)
            self.Y_observed.append(y)

    def select_next_point(self):
        """Select next evaluation point (simplified version: balance of exploration and exploitation)"""
        # Generate candidate points
        candidates = np.linspace(self.bounds[0], self.bounds[1], 100)

        # Distance from observed points (exploration)
        min_distances = []
        for c in candidates:
            distances = [abs(c - x) for x in self.X_observed]
            min_distances.append(min(distances))

        # Expected improvement from current best value (exploitation)
        best_y = max(self.Y_observed)
        improvements = [max(0, self.objective(c) - best_y) for c in candidates]

        # Score calculation (60% exploration + 40% exploitation)
        scores = 0.6 * np.array(min_distances) + 0.4 * np.array(improvements)

        return candidates[np.argmax(scores)]

    def optimize(self, n_iterations=10):
        """Execute optimization"""
        for i in range(n_iterations):
            x_next = self.select_next_point()
            y_next = self.objective(x_next)
            self.X_observed.append(x_next)
            self.Y_observed.append(y_next)

            current_best = max(self.Y_observed)
            print(f"Iteration {i+1}: x={x_next:.2f}, y={y_next:.3f}, best={current_best:.3f}")

        return self.X_observed, self.Y_observed

# Comparison experiment
bounds = [0, 5]
n_total = 13  # 3 initial points + 10 additional points

# 1. Sequential design
print("=== Sequential Design ===")
seq_designer = SequentialDesigner(simple_objective, bounds, n_initial=3)
X_seq, Y_seq = seq_designer.optimize(n_iterations=10)

# 2. Random sampling
print("\n=== Random Sampling ===")
np.random.seed(42)
X_random = np.random.uniform(bounds[0], bounds[1], n_total)
Y_random = [simple_objective(x) for x in X_random]
for i, (x, y) in enumerate(zip(X_random, Y_random), 1):
    current_best = max(Y_random[:i])
    print(f"Sample {i}: x={x:.2f}, y={y:.3f}, best={current_best:.3f}")

# Visualization
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# True function
x_true = np.linspace(0, 5, 200)
y_true = simple_objective(x_true)

# Sequential design
axes[0].plot(x_true, y_true, 'k-', linewidth=2, label='True function', alpha=0.7)
axes[0].scatter(X_seq, Y_seq, c=range(len(X_seq)), cmap='viridis',
                s=100, edgecolor='black', linewidth=1.5, label='Sequential samples', zorder=5)
axes[0].scatter(X_seq[np.argmax(Y_seq)], max(Y_seq), color='red', s=300,
                marker='*', edgecolor='black', linewidth=2, label='Best found', zorder=6)
axes[0].set_xlabel('x')
axes[0].set_ylabel('f(x)')
axes[0].set_title('Sequential Design (Best: {:.3f})'.format(max(Y_seq)))
axes[0].legend()
axes[0].grid(alpha=0.3)

# Random sampling
axes[1].plot(x_true, y_true, 'k-', linewidth=2, label='True function', alpha=0.7)
axes[1].scatter(X_random, Y_random, c=range(len(X_random)), cmap='plasma',
                s=100, edgecolor='black', linewidth=1.5, label='Random samples', zorder=5)
axes[1].scatter(X_random[np.argmax(Y_random)], max(Y_random), color='red', s=300,
                marker='*', edgecolor='black', linewidth=2, label='Best found', zorder=6)
axes[1].set_xlabel('x')
axes[1].set_ylabel('f(x)')
axes[1].set_title('Random Sampling (Best: {:.3f})'.format(max(Y_random)))
axes[1].legend()
axes[1].grid(alpha=0.3)

plt.tight_layout()
plt.show()

# Convergence comparison
best_seq = [max(Y_seq[:i+1]) for i in range(len(Y_seq))]
best_random = [max(Y_random[:i+1]) for i in range(len(Y_random))]

plt.figure(figsize=(8, 5))
plt.plot(best_seq, 'o-', linewidth=2, markersize=8, label='Sequential Design')
plt.plot(best_random, 's-', linewidth=2, markersize=8, label='Random Sampling')
plt.xlabel('Number of Evaluations')
plt.ylabel('Best Objective Value')
plt.title('Convergence Comparison')
plt.legend()
plt.grid(alpha=0.3)
plt.show()

1.3 Exploration-Exploitation Trade-off

The most important concept in Bayesian optimization is the balance between "Exploration" and "Exploitation".

• Exploitation: Intensively investigate near the current best point to pursue local improvements
• Exploration: Investigate unknown regions to discover better global optimal solutions

Example 3: Visualization of Exploration vs Exploitation

We visually understand the impact of different balances on optimization.

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

# ===================================
# Example 3: Exploration-Exploitation Trade-off
# ===================================

class ExplorationExploitationDemo:
    """Visualize exploration-exploitation balance"""

    def __init__(self):
        # Sample function: has multiple local optima
        self.x_range = np.linspace(0, 10, 200)
        self.true_func = lambda x: np.sin(x) + 0.3 * np.sin(3*x) + 0.5 * np.cos(5*x)

        # Observed points (5 points)
        self.X_obs = np.array([1.0, 3.0, 4.5, 7.0, 9.0])
        self.Y_obs = self.true_func(self.X_obs) + np.random.normal(0, 0.1, len(self.X_obs))

    def predict_with_uncertainty(self, x):
        """Simple prediction mean and uncertainty (Gaussian process approximation)"""
        # Distance-based weighted average
        weights = np.exp(-((self.X_obs[:, None] - x) / 1.0)**2)
        weights = weights / (weights.sum(axis=0) + 1e-10)

        mean = (weights.T @ self.Y_obs)

        # Uncertainty (depends on distance from observed points)
        min_dist = np.min(np.abs(self.X_obs[:, None] - x), axis=0)
        uncertainty = 0.5 * (1 - np.exp(-min_dist / 2.0))

        return mean, uncertainty

    def exploitation_strategy(self):
        """Exploitation strategy: Select point with maximum predicted mean"""
        mean, _ = self.predict_with_uncertainty(self.x_range)
        return self.x_range[np.argmax(mean)]

    def exploration_strategy(self):
        """Exploration strategy: Select point with maximum uncertainty"""
        _, uncertainty = self.predict_with_uncertainty(self.x_range)
        return self.x_range[np.argmax(uncertainty)]

    def balanced_strategy(self, alpha=0.5):
        """Balanced strategy: UCB (Upper Confidence Bound)"""
        mean, uncertainty = self.predict_with_uncertainty(self.x_range)
        ucb = mean + alpha * uncertainty
        return self.x_range[np.argmax(ucb)]

    def visualize(self):
        """Visualization"""
        mean, uncertainty = self.predict_with_uncertainty(self.x_range)

        fig, axes = plt.subplots(2, 2, figsize=(14, 10))

        # (1) True function and prediction
        ax = axes[0, 0]
        ax.plot(self.x_range, self.true_func(self.x_range), 'k--',
                linewidth=2, label='True function', alpha=0.7)
        ax.plot(self.x_range, mean, 'b-', linewidth=2, label='Predicted mean')
        ax.fill_between(self.x_range, mean - uncertainty, mean + uncertainty,
                        alpha=0.3, label='Uncertainty (±1Ã)')
        ax.scatter(self.X_obs, self.Y_obs, c='red', s=100, zorder=5,
                  edgecolor='black', linewidth=1.5, label='Observations')
        ax.set_xlabel('x')
        ax.set_ylabel('f(x)')
        ax.set_title('Model Prediction with Uncertainty')
        ax.legend()
        ax.grid(alpha=0.3)

        # (2) Exploitation strategy
        ax = axes[0, 1]
        x_exploit = self.exploitation_strategy()
        ax.plot(self.x_range, mean, 'b-', linewidth=2, label='Predicted mean')
        ax.scatter(self.X_obs, self.Y_obs, c='gray', s=80, zorder=4, alpha=0.5)
        ax.axvline(x_exploit, color='red', linestyle='--', linewidth=2,
                  label=f'Next point (Exploitation)\nx={x_exploit:.2f}')
        ax.set_xlabel('x')
        ax.set_ylabel('Predicted mean')
        ax.set_title('Exploitation Strategy: Maximize Mean')
        ax.legend()
        ax.grid(alpha=0.3)

        # (3) Exploration strategy
        ax = axes[1, 0]
        x_explore = self.exploration_strategy()
        ax.plot(self.x_range, uncertainty, 'g-', linewidth=2, label='Uncertainty')
        ax.scatter(self.X_obs, np.zeros_like(self.X_obs), c='gray', s=80,
                  zorder=4, alpha=0.5, label='Observations')
        ax.axvline(x_explore, color='blue', linestyle='--', linewidth=2,
                  label=f'Next point (Exploration)\nx={x_explore:.2f}')
        ax.set_xlabel('x')
        ax.set_ylabel('Uncertainty')
        ax.set_title('Exploration Strategy: Maximize Uncertainty')
        ax.legend()
        ax.grid(alpha=0.3)

        # (4) Balanced strategy (UCB)
        ax = axes[1, 1]
        alpha = 1.5
        x_balanced = self.balanced_strategy(alpha=alpha)
        ucb = mean + alpha * uncertainty
        ax.plot(self.x_range, mean, 'b-', linewidth=1.5, label='Mean', alpha=0.7)
        ax.plot(self.x_range, ucb, 'purple', linewidth=2, label=f'UCB (±={alpha})')
        ax.fill_between(self.x_range, mean, ucb, alpha=0.2, color='purple')
        ax.scatter(self.X_obs, self.Y_obs, c='gray', s=80, zorder=4, alpha=0.5)
        ax.axvline(x_balanced, color='purple', linestyle='--', linewidth=2,
                  label=f'Next point (Balanced)\nx={x_balanced:.2f}')
        ax.set_xlabel('x')
        ax.set_ylabel('Value')
        ax.set_title('Balanced Strategy: UCB = Mean + ± × Uncertainty')
        ax.legend()
        ax.grid(alpha=0.3)

        plt.tight_layout()
        return fig

# Execution
demo = ExplorationExploitationDemo()

print("=== Exploration vs Exploitation ===")
print(f"Exploitation (best predicted point): x = {demo.exploitation_strategy():.2f}")
print(f"Exploration (maximum uncertainty): x = {demo.exploration_strategy():.2f}")
print(f"Balanced (UCB, ±=0.5): x = {demo.balanced_strategy(alpha=0.5):.2f}")
print(f"Balanced (UCB, ±=1.5): x = {demo.balanced_strategy(alpha=1.5):.2f}")

fig = demo.visualize()
plt.show()

1.4 Basic Loop of Bayesian Optimization

The Bayesian optimization algorithm repeats the following four steps:

1. Build surrogate model: Approximate the objective function from observed data
2. Calculate acquisition function: Quantify the promise of the next point to evaluate
3. Select next point: Choose the point that maximizes the acquisition function
4. Evaluate and update: Actually evaluate and add to observed data

Example 4: Simple Bayesian Optimization Implementation

We understand the overall flow with a minimal implementation.

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
from scipy.optimize import minimize

# ===================================
# Example 4: Simple Bayesian Optimization Loop
# ===================================

class SimpleBayesianOptimization:
    """Simple Bayesian optimization (1D)"""

    def __init__(self, objective_func, bounds, noise_std=0.01):
        self.objective = objective_func
        self.bounds = bounds
        self.noise_std = noise_std

        self.X_obs = []
        self.Y_obs = []

        # Initial sampling (3 points)
        np.random.seed(42)
        for _ in range(3):
            x = np.random.uniform(bounds[0], bounds[1])
            y = self.evaluate(x)
            self.X_obs.append(x)
            self.Y_obs.append(y)

    def evaluate(self, x):
        """Evaluate objective function (with noise)"""
        return self.objective(x) + np.random.normal(0, self.noise_std)

    def gaussian_kernel(self, x1, x2, length_scale=0.5):
        """RBF kernel"""
        return np.exp(-0.5 * ((x1 - x2) / length_scale)**2)

    def predict(self, x_test):
        """Prediction by Gaussian process (simplified version)"""
        X_obs_array = np.array(self.X_obs).reshape(-1, 1)
        Y_obs_array = np.array(self.Y_obs).reshape(-1, 1)
        x_test_array = np.array(x_test).reshape(-1, 1)

        # Kernel matrix
        K = np.zeros((len(self.X_obs), len(self.X_obs)))
        for i in range(len(self.X_obs)):
            for j in range(len(self.X_obs)):
                K[i, j] = self.gaussian_kernel(self.X_obs[i], self.X_obs[j])

        # Add noise term
        K += self.noise_std**2 * np.eye(len(self.X_obs))

        # Kernel with test point
        k_star = np.array([self.gaussian_kernel(self.X_obs[i], x_test)
                          for i in range(len(self.X_obs))])

        # Predicted mean
        K_inv = np.linalg.inv(K)
        mean = k_star.T @ K_inv @ Y_obs_array

        # Predicted variance
        k_star_star = self.gaussian_kernel(x_test, x_test)
        variance = k_star_star - k_star.T @ K_inv @ k_star
        std = np.sqrt(np.maximum(variance, 0))

        return mean.flatten(), std.flatten()

    def expected_improvement(self, x):
        """Expected Improvement acquisition function"""
        mean, std = self.predict(x)
        best_y = max(self.Y_obs)

        # EI calculation
        with np.errstate(divide='ignore', invalid='ignore'):
            z = (mean - best_y) / (std + 1e-9)
            ei = (mean - best_y) * norm.cdf(z) + std * norm.pdf(z)
            ei[std == 0] = 0.0

        return -ei  # Convert to minimization problem

    def select_next_point(self):
        """Select next evaluation point (maximize EI)"""
        result = minimize(
            lambda x: self.expected_improvement(x),
            x0=np.random.uniform(self.bounds[0], self.bounds[1]),
            bounds=[self.bounds],
            method='L-BFGS-B'
        )
        return result.x[0]

    def optimize(self, n_iterations=10, verbose=True):
        """Execute optimization"""
        for i in range(n_iterations):
            # Select next point
            x_next = self.select_next_point()
            y_next = self.evaluate(x_next)

            # Update data
            self.X_obs.append(x_next)
            self.Y_obs.append(y_next)

            current_best = max(self.Y_obs)
            if verbose:
                print(f"Iter {i+1}: x={x_next:.3f}, y={y_next:.3f}, best={current_best:.3f}")

        best_idx = np.argmax(self.Y_obs)
        return self.X_obs[best_idx], self.Y_obs[best_idx]

    def plot_progress(self):
        """Visualize optimization progress"""
        x_plot = np.linspace(self.bounds[0], self.bounds[1], 200)
        y_true = [self.objective(x) for x in x_plot]
        mean, std = self.predict(x_plot)

        fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8))

        # Surrogate model
        ax1.plot(x_plot, y_true, 'k--', linewidth=2, label='True function', alpha=0.7)
        ax1.plot(x_plot, mean, 'b-', linewidth=2, label='GP mean')
        ax1.fill_between(x_plot, mean - 2*std, mean + 2*std, alpha=0.3, label='95% CI')
        ax1.scatter(self.X_obs, self.Y_obs, c='red', s=100, zorder=5,
                   edgecolor='black', linewidth=1.5, label='Observations')
        best_idx = np.argmax(self.Y_obs)
        ax1.scatter(self.X_obs[best_idx], self.Y_obs[best_idx],
                   c='gold', s=300, marker='*', zorder=6,
                   edgecolor='black', linewidth=2, label='Best')
        ax1.set_xlabel('x')
        ax1.set_ylabel('f(x)')
        ax1.set_title('Gaussian Process Surrogate Model')
        ax1.legend()
        ax1.grid(alpha=0.3)

        # Acquisition function
        ei_values = [-self.expected_improvement(x) for x in x_plot]
        ax2.plot(x_plot, ei_values, 'g-', linewidth=2, label='Expected Improvement')
        ax2.fill_between(x_plot, 0, ei_values, alpha=0.3, color='green')
        ax2.axvline(self.X_obs[-1], color='red', linestyle='--',
                   linewidth=2, label=f'Last selected: x={self.X_obs[-1]:.3f}')
        ax2.set_xlabel('x')
        ax2.set_ylabel('EI(x)')
        ax2.set_title('Acquisition Function (Expected Improvement)')
        ax2.legend()
        ax2.grid(alpha=0.3)

        plt.tight_layout()
        return fig

# Test function
def test_function(x):
    return -(x - 2.5)**2 * np.sin(10 * x) + 3

# Execution
print("=== Simple Bayesian Optimization ===\n")
bo = SimpleBayesianOptimization(test_function, bounds=[0, 5], noise_std=0.05)
x_best, y_best = bo.optimize(n_iterations=12, verbose=True)

print(f"\n=== Final Result ===")
print(f"Best x: {x_best:.4f}")
print(f"Best y: {y_best:.4f}")

fig = bo.plot_progress()
plt.show()

1.5 Comparison: BO vs Grid Search vs Random Search

We quantitatively evaluate the superiority of Bayesian optimization.

Example 5: Performance Comparison of Three Methods

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
from scipy.stats import norm

# ===================================
# Example 5: BO vs Grid Search vs Random Search
# ===================================

# Test function (2D)
def branin_function(x):
    """Branin function (global optimization benchmark)"""
    x1, x2 = x[0], x[1]
    a, b, c = 1, 5.1/(4*np.pi**2), 5/np.pi
    r, s, t = 6, 10, 1/(8*np.pi)

    term1 = a * (x2 - b*x1**2 + c*x1 - r)**2
    term2 = s * (1 - t) * np.cos(x1)
    term3 = s

    return -(term1 + term2 + term3)  # Convert to maximization problem

class OptimizationComparison:
    """Comparison experiment of three methods"""

    def __init__(self, objective, bounds, budget=30):
        self.objective = objective
        self.bounds = np.array(bounds)  # [[x1_min, x1_max], [x2_min, x2_max]]
        self.budget = budget
        self.dim = len(bounds)

    def grid_search(self):
        """Grid search"""
        n_per_dim = int(np.ceil(self.budget ** (1/self.dim)))

        grid_1d = [np.linspace(b[0], b[1], n_per_dim) for b in self.bounds]
        grid = np.meshgrid(*grid_1d)

        X_grid = np.column_stack([g.ravel() for g in grid])[:self.budget]
        Y_grid = [self.objective(x) for x in X_grid]

        return X_grid, Y_grid

    def random_search(self, seed=42):
        """Random search"""
        np.random.seed(seed)
        X_random = np.random.uniform(
            self.bounds[:, 0], self.bounds[:, 1],
            size=(self.budget, self.dim)
        )
        Y_random = [self.objective(x) for x in X_random]

        return X_random, Y_random

    def bayesian_optimization(self, seed=42):
        """Bayesian optimization (simplified version)"""
        np.random.seed(seed)

        # Initial samples (5 points)
        n_init = 5
        X = np.random.uniform(self.bounds[:, 0], self.bounds[:, 1],
                             size=(n_init, self.dim))
        Y = [self.objective(x) for x in X]

        # Sequential optimization
        for _ in range(self.budget - n_init):
            # Simple GP prediction
            def gp_mean_std(x_test):
                distances = np.linalg.norm(X - x_test, axis=1)
                weights = np.exp(-distances**2 / 2.0)
                weights = weights / (weights.sum() + 1e-10)

                mean = weights @ Y
                std = 1.0 * np.exp(-np.min(distances) / 1.5)

                return mean, std

            # EI acquisition function
            def neg_ei(x):
                mean, std = gp_mean_std(x)
                best_y = max(Y)
                z = (mean - best_y) / (std + 1e-9)
                ei = (mean - best_y) * norm.cdf(z) + std * norm.pdf(z)
                return -ei

            # Next point selection
            result = minimize(
                neg_ei,
                x0=np.random.uniform(self.bounds[:, 0], self.bounds[:, 1]),
                bounds=self.bounds,
                method='L-BFGS-B'
            )

            x_next = result.x
            y_next = self.objective(x_next)

            X = np.vstack([X, x_next])
            Y.append(y_next)

        return X, Y

    def compare(self, n_trials=5):
        """Compare average performance over multiple trials"""
        results = {
            'Grid Search': [],
            'Random Search': [],
            'Bayesian Optimization': []
        }

        for trial in range(n_trials):
            print(f"\n=== Trial {trial + 1}/{n_trials} ===")

            # Grid Search
            X_grid, Y_grid = self.grid_search()
            best_grid = [max(Y_grid[:i+1]) for i in range(len(Y_grid))]
            results['Grid Search'].append(best_grid)
            print(f"Grid Search best: {max(Y_grid):.4f}")

            # Random Search
            X_rand, Y_rand = self.random_search(seed=trial)
            best_rand = [max(Y_rand[:i+1]) for i in range(len(Y_rand))]
            results['Random Search'].append(best_rand)
            print(f"Random Search best: {max(Y_rand):.4f}")

            # Bayesian Optimization
            X_bo, Y_bo = self.bayesian_optimization(seed=trial)
            best_bo = [max(Y_bo[:i+1]) for i in range(len(Y_bo))]
            results['Bayesian Optimization'].append(best_bo)
            print(f"Bayesian Optimization best: {max(Y_bo):.4f}")

        return results

    def plot_comparison(self, results):
        """Visualize comparison results"""
        fig, ax = plt.subplots(figsize=(10, 6))

        colors = {'Grid Search': 'blue', 'Random Search': 'orange',
                 'Bayesian Optimization': 'green'}

        for method, trials in results.items():
            trials_array = np.array(trials)
            mean_curve = trials_array.mean(axis=0)
            std_curve = trials_array.std(axis=0)

            x_axis = np.arange(1, len(mean_curve) + 1)
            ax.plot(x_axis, mean_curve, linewidth=2.5, label=method,
                   color=colors[method], marker='o', markersize=4)
            ax.fill_between(x_axis, mean_curve - std_curve, mean_curve + std_curve,
                           alpha=0.2, color=colors[method])

        ax.set_xlabel('Number of Evaluations', fontsize=12)
        ax.set_ylabel('Best Objective Value Found', fontsize=12)
        ax.set_title('Optimization Performance Comparison (Mean ± Std over 5 trials)',
                    fontsize=13, fontweight='bold')
        ax.legend(fontsize=11)
        ax.grid(alpha=0.3)

        return fig

# Run experiment
bounds = [[-5, 10], [0, 15]]  # Domain of Branin function
comparison = OptimizationComparison(branin_function, bounds, budget=30)

print("Running optimization comparison...")
results = comparison.compare(n_trials=5)

# Final performance summary
print("\n=== Final Performance Summary ===")
for method, trials in results.items():
    final_values = [trial[-1] for trial in trials]
    print(f"{method:25s}: {np.mean(final_values):.4f} ± {np.std(final_values):.4f}")

fig = comparison.plot_comparison(results)
plt.show()

1.6 Convergence Analysis and Iteration Tracking

We quantitatively evaluate optimization progress and learn methods for convergence determination.

Example 6: Convergence Diagnostic Tool

import numpy as np
import matplotlib.pyplot as plt

# ===================================
# Example 6: Convergence analysis and iteration tracking
# ===================================

class ConvergenceAnalyzer:
    """Analyze convergence of Bayesian optimization"""

    def __init__(self, X_history, Y_history, true_optimum=None):
        self.X_history = np.array(X_history)
        self.Y_history = np.array(Y_history)
        self.true_optimum = true_optimum
        self.n_iter = len(Y_history)

    def compute_metrics(self):
        """Calculate convergence metrics"""
        # Cumulative best value
        cumulative_best = [max(self.Y_history[:i+1]) for i in range(self.n_iter)]

        # Improvement (progress at each iteration)
        improvements = [0]
        for i in range(1, self.n_iter):
            improvements.append(max(0, cumulative_best[i] - cumulative_best[i-1]))

        # Optimality gap (if true optimum is known)
        if self.true_optimum is not None:
            optimality_gap = [self.true_optimum - cb for cb in cumulative_best]
        else:
            optimality_gap = None

        # Convergence rate (standard deviation of recent 5 improvements)
        convergence_rate = []
        window = 5
        for i in range(self.n_iter):
            if i < window:
                convergence_rate.append(np.nan)
            else:
                recent_improvements = improvements[i-window+1:i+1]
                convergence_rate.append(np.std(recent_improvements))

        return {
            'cumulative_best': cumulative_best,
            'improvements': improvements,
            'optimality_gap': optimality_gap,
            'convergence_rate': convergence_rate
        }

    def is_converged(self, tolerance=1e-3, patience=5):
        """Convergence determination"""
        metrics = self.compute_metrics()
        improvements = metrics['improvements']

        # Improvement is less than tolerance for recent patience iterations
        if len(improvements) < patience:
            return False

        recent_improvements = improvements[-patience:]
        return all(imp < tolerance for imp in recent_improvements)

    def plot_diagnostics(self):
        """Diagnostic plots"""
        metrics = self.compute_metrics()

        fig, axes = plt.subplots(2, 2, figsize=(14, 10))
        iterations = np.arange(1, self.n_iter + 1)

        # (1) Transition of cumulative best value
        ax = axes[0, 0]
        ax.plot(iterations, metrics['cumulative_best'], 'b-', linewidth=2, marker='o')
        if self.true_optimum is not None:
            ax.axhline(self.true_optimum, color='red', linestyle='--',
                      linewidth=2, label=f'True optimum: {self.true_optimum:.3f}')
            ax.legend()
        ax.set_xlabel('Iteration')
        ax.set_ylabel('Best Value Found')
        ax.set_title('Convergence: Cumulative Best')
        ax.grid(alpha=0.3)

        # (2) Improvement at each iteration
        ax = axes[0, 1]
        ax.bar(iterations, metrics['improvements'], color='green', alpha=0.7)
        ax.set_xlabel('Iteration')
        ax.set_ylabel('Improvement')
        ax.set_title('Improvement per Iteration')
        ax.grid(alpha=0.3, axis='y')

        # (3) Optimality gap (log scale)
        ax = axes[1, 0]
        if metrics['optimality_gap'] is not None:
            gap = np.array(metrics['optimality_gap'])
            gap[gap <= 0] = 1e-10  # Handle negative values
            ax.semilogy(iterations, gap, 'r-', linewidth=2, marker='s')
            ax.set_xlabel('Iteration')
            ax.set_ylabel('Optimality Gap (log scale)')
            ax.set_title('Distance to True Optimum')
            ax.grid(alpha=0.3, which='both')
        else:
            ax.text(0.5, 0.5, 'True optimum unknown',
                   ha='center', va='center', fontsize=14)
            ax.set_xticks([])
            ax.set_yticks([])

        # (4) Convergence rate (variability of improvements)
        ax = axes[1, 1]
        valid_idx = ~np.isnan(metrics['convergence_rate'])
        ax.plot(iterations[valid_idx], np.array(metrics['convergence_rate'])[valid_idx],
               'purple', linewidth=2, marker='d')
        ax.axhline(1e-3, color='orange', linestyle='--',
                  linewidth=2, label='Convergence threshold')
        ax.set_xlabel('Iteration')
        ax.set_ylabel('Convergence Rate (Std of recent improvements)')
        ax.set_title('Convergence Rate Indicator')
        ax.legend()
        ax.grid(alpha=0.3)

        plt.tight_layout()
        return fig

    def print_summary(self):
        """Summary report"""
        metrics = self.compute_metrics()

        print("=== Convergence Analysis Summary ===")
        print(f"Total iterations: {self.n_iter}")
        print(f"Best value found: {max(self.Y_history):.6f}")
        print(f"Final improvement: {metrics['improvements'][-1]:.6f}")

        if self.true_optimum is not None:
            final_gap = self.true_optimum - max(self.Y_history)
            print(f"True optimum: {self.true_optimum:.6f}")
            print(f"Optimality gap: {final_gap:.6f} ({final_gap/self.true_optimum*100:.2f}%)")

        converged = self.is_converged()
        print(f"Converged: {'Yes' if converged else 'No'}")

        # Iteration with maximum improvement
        max_imp_iter = np.argmax(metrics['improvements'])
        print(f"Largest improvement at iteration: {max_imp_iter + 1} "
              f"(”y = {metrics['improvements'][max_imp_iter]:.6f})")

# Demo execution (using results from Example 4)
np.random.seed(42)

def test_func(x):
    return -(x - 2.5)**2 * np.sin(10 * x) + 3

# Run BO and obtain history
from scipy.stats import norm
from scipy.optimize import minimize

X_hist, Y_hist = [], []
bounds = [0, 5]

# Initial sampling
for _ in range(3):
    x = np.random.uniform(bounds[0], bounds[1])
    X_hist.append(x)
    Y_hist.append(test_func(x) + np.random.normal(0, 0.02))

# Sequential optimization (15 iterations)
for iteration in range(15):
    # Simple GP prediction
    def gp_predict(x_test):
        dists = np.abs(np.array(X_hist) - x_test)
        weights = np.exp(-dists**2 / 0.5)
        mean = weights @ Y_hist / (weights.sum() + 1e-10)
        std = 0.5 * (1 - np.exp(-np.min(dists) / 1.0))
        return mean, std

    # EI acquisition function
    def neg_ei(x):
        mean, std = gp_predict(x)
        z = (mean - max(Y_hist)) / (std + 1e-9)
        ei = (mean - max(Y_hist)) * norm.cdf(z) + std * norm.pdf(z)
        return -ei

    # Next point selection
    res = minimize(neg_ei, x0=np.random.uniform(bounds[0], bounds[1]),
                  bounds=[bounds], method='L-BFGS-B')

    x_next = res.x[0]
    y_next = test_func(x_next) + np.random.normal(0, 0.02)

    X_hist.append(x_next)
    Y_hist.append(y_next)

# Convergence analysis
analyzer = ConvergenceAnalyzer(X_hist, Y_hist, true_optimum=3.62)
analyzer.print_summary()
fig = analyzer.plot_diagnostics()
plt.show()

1.7 Practical Example of Hyperparameter Tuning

We implement hyperparameter tuning of machine learning models as a representative application of Bayesian optimization.

Example 7: Hyperparameter Optimization of Random Forest

import numpy as np
import matplotlib.pyplot as plt
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import cross_val_score
from sklearn.datasets import make_regression
from scipy.stats import norm
from scipy.optimize import minimize

# ===================================
# Example 7: Hyperparameter tuning
# ===================================

# Generate sample dataset (assuming sensor data from chemical process)
np.random.seed(42)
X_data, y_data = make_regression(n_samples=200, n_features=10, noise=10, random_state=42)

class HyperparameterOptimizer:
    """Bayesian optimization of Random Forest hyperparameters"""

    def __init__(self, X_train, y_train):
        self.X_train = X_train
        self.y_train = y_train

        # Three hyperparameters to optimize
        self.param_bounds = {
            'n_estimators': [10, 200],      # Number of decision trees
            'max_depth': [3, 20],            # Maximum depth
            'min_samples_split': [2, 20]     # Minimum samples for split
        }

        self.X_obs = []
        self.Y_obs = []

        # Initial random sampling (5 points)
        for _ in range(5):
            params = self._sample_random_params()
            score = self._evaluate(params)
            self.X_obs.append(params)
            self.Y_obs.append(score)

    def _sample_random_params(self):
        """Randomly sample hyperparameters"""
        return [
            np.random.randint(self.param_bounds['n_estimators'][0],
                            self.param_bounds['n_estimators'][1] + 1),
            np.random.randint(self.param_bounds['max_depth'][0],
                            self.param_bounds['max_depth'][1] + 1),
            np.random.randint(self.param_bounds['min_samples_split'][0],
                            self.param_bounds['min_samples_split'][1] + 1)
        ]

    def _evaluate(self, params):
        """Evaluate hyperparameter performance (5-fold CV)"""
        n_est, max_d, min_split = [int(p) for p in params]

        model = RandomForestRegressor(
            n_estimators=n_est,
            max_depth=max_d,
            min_samples_split=min_split,
            random_state=42,
            n_jobs=-1
        )

        # Cross-validation R² score
        scores = cross_val_score(model, self.X_train, self.y_train,
                                cv=5, scoring='r2')
        return scores.mean()

    def _gp_predict(self, params_test):
        """Simple Gaussian process prediction"""
        X_obs_array = np.array(self.X_obs)
        params_array = np.array(params_test).reshape(1, -1)

        # Normalization
        X_obs_norm = (X_obs_array - X_obs_array.mean(axis=0)) / (X_obs_array.std(axis=0) + 1e-10)
        params_norm = (params_array - X_obs_array.mean(axis=0)) / (X_obs_array.std(axis=0) + 1e-10)

        # Distance-based weights
        dists = np.linalg.norm(X_obs_norm - params_norm, axis=1)
        weights = np.exp(-dists**2 / 2.0)

        mean = weights @ self.Y_obs / (weights.sum() + 1e-10)
        std = 0.2 * (1 - np.exp(-np.min(dists) / 1.5))

        return mean, std

    def _expected_improvement(self, params):
        """EI acquisition function"""
        mean, std = self._gp_predict(params)
        best_y = max(self.Y_obs)

        z = (mean - best_y) / (std + 1e-9)
        ei = (mean - best_y) * norm.cdf(z) + std * norm.pdf(z)

        return -ei  # Minimization problem

    def optimize(self, n_iterations=15):
        """Execute Bayesian optimization"""
        print("=== Hyperparameter Optimization ===\n")

        for i in range(n_iterations):
            # Select next hyperparameters
            bounds_array = [
                self.param_bounds['n_estimators'],
                self.param_bounds['max_depth'],
                self.param_bounds['min_samples_split']
            ]

            result = minimize(
                self._expected_improvement,
                x0=self._sample_random_params(),
                bounds=bounds_array,
                method='L-BFGS-B'
            )

            params_next = [int(p) for p in result.x]
            score_next = self._evaluate(params_next)

            self.X_obs.append(params_next)
            self.Y_obs.append(score_next)

            current_best = max(self.Y_obs)
            print(f"Iter {i+1}: n_est={params_next[0]}, max_depth={params_next[1]}, "
                  f"min_split={params_next[2]} ’ R²={score_next:.4f} (best={current_best:.4f})")

        # Best parameters
        best_idx = np.argmax(self.Y_obs)
        best_params = self.X_obs[best_idx]
        best_score = self.Y_obs[best_idx]

        print(f"\n=== Best Hyperparameters ===")
        print(f"n_estimators: {best_params[0]}")
        print(f"max_depth: {best_params[1]}")
        print(f"min_samples_split: {best_params[2]}")
        print(f"Best R² score: {best_score:.4f}")

        return best_params, best_score

    def plot_optimization_history(self):
        """Visualize optimization history"""
        fig, axes = plt.subplots(2, 2, figsize=(14, 10))

        iterations = np.arange(1, len(self.Y_obs) + 1)
        cumulative_best = [max(self.Y_obs[:i+1]) for i in range(len(self.Y_obs))]

        # (1) Transition of R² score
        ax = axes[0, 0]
        ax.plot(iterations, self.Y_obs, 'o-', linewidth=2, markersize=8,
               label='Observed R²', alpha=0.7)
        ax.plot(iterations, cumulative_best, 's-', linewidth=2.5, markersize=8,
               color='red', label='Best R²')
        ax.set_xlabel('Iteration')
        ax.set_ylabel('R² Score')
        ax.set_title('Optimization Progress')
        ax.legend()
        ax.grid(alpha=0.3)

        # (2) Exploration trajectory of n_estimators
        ax = axes[0, 1]
        n_estimators = [x[0] for x in self.X_obs]
        scatter = ax.scatter(iterations, n_estimators, c=self.Y_obs,
                           cmap='viridis', s=100, edgecolor='black', linewidth=1.5)
        ax.set_xlabel('Iteration')
        ax.set_ylabel('n_estimators')
        ax.set_title('Parameter Exploration: n_estimators')
        plt.colorbar(scatter, ax=ax, label='R² Score')
        ax.grid(alpha=0.3)

        # (3) Exploration trajectory of max_depth
        ax = axes[1, 0]
        max_depth = [x[1] for x in self.X_obs]
        scatter = ax.scatter(iterations, max_depth, c=self.Y_obs,
                           cmap='viridis', s=100, edgecolor='black', linewidth=1.5)
        ax.set_xlabel('Iteration')
        ax.set_ylabel('max_depth')
        ax.set_title('Parameter Exploration: max_depth')
        plt.colorbar(scatter, ax=ax, label='R² Score')
        ax.grid(alpha=0.3)

        # (4) Parameter space exploration (2D projection: n_estimators vs max_depth)
        ax = axes[1, 1]
        scatter = ax.scatter(n_estimators, max_depth, c=self.Y_obs,
                           s=150, cmap='viridis', edgecolor='black', linewidth=1.5)
        best_idx = np.argmax(self.Y_obs)
        ax.scatter(n_estimators[best_idx], max_depth[best_idx],
                  s=500, marker='*', color='red', edgecolor='black', linewidth=2,
                  label='Best', zorder=5)
        ax.set_xlabel('n_estimators')
        ax.set_ylabel('max_depth')
        ax.set_title('Parameter Space Exploration')
        plt.colorbar(scatter, ax=ax, label='R² Score')
        ax.legend()
        ax.grid(alpha=0.3)

        plt.tight_layout()
        return fig

# Execution
optimizer = HyperparameterOptimizer(X_data, y_data)
best_params, best_score = optimizer.optimize(n_iterations=15)

fig = optimizer.plot_optimization_history()
plt.show()

# Comparison with baseline
print("\n=== Comparison with Default Parameters ===")
default_model = RandomForestRegressor(random_state=42, n_jobs=-1)
default_score = cross_val_score(default_model, X_data, y_data, cv=5, scoring='r2').mean()
print(f"Default R² score: {default_score:.4f}")
print(f"Optimized R² score: {best_score:.4f}")
print(f"Improvement: {(best_score - default_score) / default_score * 100:.2f}%")

Learning Objectives Review

Upon completing this chapter, you will be able to explain and implement the following:

Basic Understanding

•  Can explain the characteristics and challenges of black-box optimization problems
•  Can describe the advantages of sequential design strategy compared to random search
•  Understand the concept of exploration-exploitation trade-off

Practical Skills

•  Can formulate black-box problems in chemical processes
•  Can implement simple Bayesian optimization loops
•  Can compare and evaluate the performance of three methods (BO/Grid/Random)
•  Can analyze optimization progress using convergence diagnostic tools

Application Skills

•  Can apply to hyperparameter tuning of machine learning models
•  Can select appropriate optimization methods for practical problems
•  Can evaluate the reliability of optimization results

Next Steps

In Chapter 1, we learned the basic concepts and implementation of Bayesian optimization. In the next chapter, we will learn in detail about "Gaussian Processes", the core technology of Bayesian optimization.

References

1. Montgomery, D. C. (2019). Design and Analysis of Experiments (9th ed.). Wiley.
2. Box, G. E. P., Hunter, J. S., & Hunter, W. G. (2005). Statistics for Experimenters: Design, Innovation, and Discovery (2nd ed.). Wiley.
3. Seborg, D. E., Edgar, T. F., Mellichamp, D. A., & Doyle III, F. J. (2016). Process Dynamics and Control (4th ed.). Wiley.
4. McKay, M. D., Beckman, R. J., & Conover, W. J. (2000). "A Comparison of Three Methods for Selecting Values of Input Variables in the Analysis of Output from a Computer Code." Technometrics, 42(1), 55-61.