Optimal Experimental Design

Optimal experimental design (OED) provides a principled framework for allocating experimental resources to maximize the information gained about model parameters. Rather than relying on classical designs (factorial, Latin square), OED uses mathematical optimization to construct designs tailored to specific models, objectives, and constraints. This lesson develops the theory of optimal design criteria, information-theoretic foundations, and computational algorithms.

Information Matrix

Optimal Design Criteria

The three classical optimality criteria are matrix orderings of , each optimizing a different scalar function of the parameter covariance.

Classical Optimal Designs

Sequential Design

Bayesian Optimal Design

Coordinate Exchange Algorithm

Design Diagnostics

Python Implementation

import numpy as np
from scipy.optimize import minimize
from itertools import product
import matplotlib.pyplot as plt

np.random.seed(42)

# --- Design Space ---
# 2-factor design space: x1, x2 in [-1, 1]
candidate_set = np.array(list(product(np.linspace(-1, 1, 20), repeat=2)))
n_factors = 2

# --- Model: Quadratic Regression ---
# y = b0 + b1*x1 + b2*x2 + b11*x1^2 + b22*x2^2 + b12*x1*x2 (p=6 parameters)
def design_matrix(X):
    """Build quadratic regression design matrix."""
    n = X.shape[0]
    return np.column_stack([
        np.ones(n), X[:, 0], X[:, 1],
        X[:, 0]**2, X[:, 1]**2, X[:, 0]*X[:, 1]
    ])

def d_criterion(X, n_factors=2):
    """Negative D-criterion (for minimization)."""
    X_design = design_matrix(X)
    M = X_design.T @ X_design
    sign, logdet = np.linalg.slogdet(M)
    if sign <= 0:
        return 1e10
    return -logdet

# --- Coordinate Exchange Algorithm ---
def coordinate_exchange(n_runs, candidate_set, n_iter=100, n_restarts=10):
    """Find D-optimal design via coordinate exchange."""
    best_design = None
    best_det = -np.inf

    for restart in range(n_restarts):
        # Random initialization
        idx = np.random.choice(len(candidate_set), n_runs, replace=False)
        design = candidate_set[idx].copy()

        for iteration in range(n_iter):
            improved = False
            for i in range(n_runs):
                # Remove point i
                X_rest = np.delete(design, i, axis=0)
                M_rest = design_matrix(X_rest).T @ design_matrix(X_rest)

                # Find best replacement
                best_x = design[i]
                best_det_local = np.linalg.slogdet(M_rest + design_matrix(design[i:i+1]).T @ design_matrix(design[i:i+1]))[1]

                for j, x_cand in enumerate(candidate_set):
                    X_new = design_matrix(x_cand.reshape(1, -1))
                    det_new = np.linalg.slogdet(M_rest + X_new.T @ X_new)[1]
                    if det_new > best_det_local:
                        best_det_local = det_new
                        best_x = x_cand
                        improved = True

                design[i] = best_x

            if not improved:
                break

        # Evaluate final design
        X_final = design_matrix(design)
        det_final = np.linalg.slogdet(X_final.T @ X_final)[1]

        if det_final > best_det:
            best_det = det_final
            best_design = design.copy()

    return best_design, best_det

# --- Generate D-Optimal Design ---
n_runs = 12  # More than p=6 for lack-of-fit test
optimal_design, d_crit = coordinate_exchange(n_runs, candidate_set)

print("=== D-Optimal Design (Quadratic Regression) ===")
print(f"Number of runs: {n_runs}")
print(f"Number of parameters: 6")
print(f"D-criterion (log det M): {d_crit:.4f}")
print(f"\nDesign points:")
X_opt = design_matrix(optimal_design)
M_opt = X_opt.T @ X_opt
print(np.array2string(optimal_design, precision=3, suppress_small=True))

# --- D-Efficiency ---
def d_efficiency(design, reference_det, p=6):
    """Compute D-efficiency relative to reference."""
    X = design_matrix(design)
    det_current = np.linalg.slogdet(X.T @ X)[1]
    return np.exp((det_current - reference_det) / p)

print(f"\nD-efficiency: {d_efficiency(optimal_design, d_crit):.4f}")

# --- Prediction Variance Profile ---
x1_grid = np.linspace(-1, 1, 50)
x2_grid = np.linspace(-1, 1, 50)
X1, X2 = np.meshgrid(x1_grid, x2_grid)
spv = np.zeros_like(X1)

M_inv = np.linalg.inv(M_opt)
for i in range(50):
    for j in range(50):
        x_test = np.array([1, X1[i,j], X2[i,j], X1[i,j]**2, X2[i,j]**2, X1[i,j]*X2[i,j]])
        spv[i, j] = n_runs * x_test @ M_inv @ x_test

print(f"\n=== Prediction Variance ===")
print(f"Min SPV: {spv.min():.2f}")
print(f"Max SPV: {spv.max():.2f}")
print(f"Mean SPV: {spv.mean():.2f}")
print(f"Theoretical max (G-opt): p = 6")

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

# Design points
axes[0].scatter(optimal_design[:, 0], optimal_design[:, 1], s=100, c='red', 
                edgecolors='black', zorder=5, label='Design points')
axes[0].scatter(candidate_set[:, 0], candidate_set[:, 1], s=10, c='gray', alpha=0.3, label='Candidates')
axes[0].set_xlabel('x₁')
axes[0].set_ylabel('x₂')
axes[0].set_title(f'D-Optimal Design (n={n_runs}, log det = {d_crit:.2f})')
axes[0].legend()
axes[0].set_xlim([-1.1, 1.1])
axes[0].set_ylim([-1.1, 1.1])
axes[0].grid(True, alpha=0.3)

# Prediction variance profile
im = axes[1].contourf(X1, X2, spv, levels=20, cmap='viridis')
plt.colorbar(im, ax=axes[1], label='SPV(x)')
axes[1].scatter(optimal_design[:, 0], optimal_design[:, 1], s=60, c='red', 
                edgecolors='white', zorder=5)
axes[1].set_xlabel('x₁')
axes[1].set_ylabel('x₂')
axes[1].set_title('Scaled Prediction Variance Profile')

plt.tight_layout()
plt.savefig("optimal_design.png", dpi=150, bbox_inches="tight")
plt.show()

# --- Comparison: Random vs Optimal ---
n_trials = 100
random_dets = []
for _ in range(n_trials):
    idx = np.random.choice(len(candidate_set), n_runs, replace=False)
    X_rand = design_matrix(candidate_set[idx])
    det_val = np.linalg.slogdet(X_rand.T @ X_rand)[1]
    random_dets.append(det_val)

print(f"\n=== Design Comparison ===")
print(f"D-optimal log det: {d_crit:.4f}")
print(f"Random designs: mean={np.mean(random_dets):.4f}, std={np.std(random_dets):.4f}")
print(f"D-efficiency of random vs optimal: {np.exp((np.mean(random_dets) - d_crit)/6)*100:.1f}%")

Summary