Chapter 2: Metaheuristic Optimization

Chapter 1 introduced gradient-based optimization, which is fast and reliable when the objective function is smooth and differentiable. But many real optimization problems in machine learning are not like that: hyperparameter choices are discrete, evaluation results are noisy, and some objectives have no usable gradient at all. This chapter introduces metaheuristics — general-purpose, gradient-free search strategies — and shows you how to implement three classics from scratch: simulated annealing, genetic algorithms, and particle swarm optimization. You will also use SciPy's production-quality differential evolution and learn how to compare these methods fairly.

Learning Objectives

By completing this chapter, you will master the following:

• ✅ Understand when gradient-free methods are needed (non-differentiable, discrete, and noisy objectives)
• ✅ Implement simulated annealing from scratch with the Metropolis criterion and cooling schedules
• ✅ Implement a genetic algorithm (selection, crossover, mutation) for a feature-selection problem
• ✅ Implement particle swarm optimization and visualize how the swarm converges
• ✅ Use scipy.optimize.differential_evolution for practical global optimization
• ✅ Compare metaheuristics against gradient methods fairly, using equal evaluation budgets

2.1 Why Gradient-Free Optimization?

What is a Metaheuristic?

A metaheuristic is a general-purpose search strategy that finds good solutions to an optimization problem using only function evaluations — it never asks for a gradient. The name combines "meta" (above, beyond) and "heuristic" (a practical rule of thumb): a metaheuristic is a higher-level recipe that guides simpler search moves, such as "perturb the current solution" or "combine two good solutions".

This is a real trade-off, not a free upgrade. Gradient methods use derivative information, which tells them the locally best direction to move; metaheuristics must discover good directions by trial and error, so they typically need many more function evaluations. The reason we still need them is that gradient information is often unavailable or misleading.

Three Ways Gradients Fail

The third case is subtle: gradient descent still runs on a multimodal function (a function with many local minima), it just converges to whichever local minimum is closest to the starting point. Metaheuristics are designed to keep exploring instead of committing to the first valley they find.

Seeing the Failure Modes in Code

Let us demonstrate the first two failure modes concretely. We minimize three one-dimensional functions, all with true optimum at $x^* = 2$: a smooth parabola, a "staircase" version of it (piecewise constant), and a noisy version. We compare gradient descent using finite differences against the simplest possible gradient-free method, random search (evaluate random points and keep the best).

# Requirements:
# - Python 3.9+
# - numpy>=1.24.0

"""
Example 1: Where gradients fail

Purpose: Show that finite-difference gradient descent breaks on
         step-shaped and noisy objectives, while random search does not
Target: Intermediate
Execution time: under 5 seconds
"""

import numpy as np

rng = np.random.default_rng(42)

def smooth_objective(x):
    """Differentiable: gradient descent works well here."""
    return (x - 2.0) ** 2

def step_objective(x):
    """Piecewise constant: the gradient is zero almost everywhere."""
    return np.floor(4.0 * (x - 2.0) ** 2) / 4.0

def noisy_objective(x):
    """Noisy evaluation: finite differences amplify the noise."""
    return (x - 2.0) ** 2 + rng.normal(0.0, 0.1)

def finite_diff_descent(f, x0, lr=0.1, eps=1e-4, n_steps=200):
    """Gradient descent using a central finite-difference gradient."""
    x = x0
    for _ in range(n_steps):
        grad = (f(x + eps) - f(x - eps)) / (2.0 * eps)
        x = x - lr * grad
        if abs(x) > 1e6:  # stop if the iterate diverges
            return x, True
    return x, False

def random_search(f, low=-5.0, high=5.0, n_evals=200, seed=0):
    """Baseline gradient-free method: pure random sampling."""
    rs = np.random.default_rng(seed)
    candidates = rs.uniform(low, high, size=n_evals)
    values = np.array([f(c) for c in candidates])
    return candidates[np.argmin(values)]

print("True optimum: x* = 2.0\n")
for name, f in [("smooth", smooth_objective),
                ("step  ", step_objective),
                ("noisy ", noisy_objective)]:
    x_gd, diverged = finite_diff_descent(f, x0=4.0)
    x_rs = random_search(f)
    gd_str = "DIVERGED" if diverged else f"{x_gd:8.4f}"
    print(f"{name} | gradient descent: {gd_str} | random search: {x_rs:8.4f}")

Output:

True optimum: x* = 2.0

smooth | gradient descent:   2.0000 | random search:   2.0265
step   | gradient descent:   4.0000 | random search:   2.2950
noisy  | gradient descent: -41.2292 | random search:   2.2149

Read this table carefully — it summarizes the whole motivation for this chapter:

• Smooth: gradient descent nails the optimum exactly. When gradients work, use them.
• Step: gradient descent never moves. Inside each "stair", the function is flat, so the finite-difference gradient is zero and the iterate stays at its starting point $x = 4$.
• Noisy: gradient descent wanders far from the optimum. The finite difference divides evaluation noise of size $\sim 0.1$ by $2\varepsilon = 2 \times 10^{-4}$, producing gradient estimates of size $\sim 500$ that swamp the true gradient.

Random search, crude as it is, gets within about $0.3$ of the optimum in all three cases, because it only compares function values — an operation that is robust to flatness and moderately robust to noise. Every metaheuristic in this chapter is, at heart, a smarter version of random search.

2.2 Simulated Annealing

The Physical Metaphor

Simulated annealing (SA) is a single-solution metaheuristic inspired by annealing in metallurgy: a metal heated and then cooled slowly settles into a low-energy crystalline state, while rapid cooling freezes in defects. In optimization terms, the "energy" is the objective value $E = f(\mathbf{x})$, and the "temperature" $T$ controls how willing we are to accept moves that make the objective worse.

Why would we ever accept a worse solution? Because that is exactly how you escape a local minimum: sometimes you must climb out of the current valley to reach a deeper one. A method that only accepts improvements (a greedy search) gets trapped in the first valley it enters.

The Metropolis Criterion

At each step, SA proposes a random perturbation of the current solution and computes the change in objective value $\Delta E = f(\mathbf{x}_{\text{new}}) - f(\mathbf{x}_{\text{current}})$. The proposal is accepted with probability given by the Metropolis criterion:

$$ P(\text{accept}) = \begin{cases} 1 & \text{if } \Delta E \le 0 \\[4pt] \exp\!\left(-\dfrac{\Delta E}{T}\right) & \text{if } \Delta E > 0 \end{cases} $$

Improvements are always accepted. Deteriorations are accepted with a probability that shrinks exponentially with how bad they are ($\Delta E$) and how cold the system is ($T$). At high temperature, SA behaves like a random walk that explores broadly; as $T \to 0$, it becomes a pure greedy search that only accepts improvements.

Cooling Schedules

A cooling schedule specifies how the temperature decreases over iterations $k$. Three common choices:

The logarithmic schedule illustrates an honest lesson about theory versus practice: it is the only schedule with a proof of convergence to the global optimum, but the proof requires so many iterations that nobody uses it. Practitioners use exponential schedules and tune $\alpha$ empirically.

Test Problem: The Rastrigin Function

We test SA on the Rastrigin function, a standard multimodal benchmark:

$$ f(\mathbf{x}) = 10n + \sum_{i=1}^{n} \left[ x_i^2 - 10\cos(2\pi x_i) \right] $$

It is a parabola with a cosine ripple superimposed, creating a regular grid of local minima (roughly one near every integer coordinate). The global minimum is $f(\mathbf{0}) = 0$. In 2-D on $[-5.12, 5.12]^2$ there are over a hundred local minima — a genuine trap for greedy methods.

From-Scratch Implementation

# Requirements:
# - Python 3.9+
# - numpy>=1.24.0

"""
Example 2: Simulated annealing on the 2-D Rastrigin function

Purpose: Implement SA with the Metropolis criterion and an exponential
         cooling schedule; compare against greedy search
Target: Intermediate
Execution time: 10-30 seconds
"""

import numpy as np

def rastrigin(x):
    """Rastrigin function: global minimum f(0,...,0) = 0."""
    x = np.asarray(x, dtype=float)
    return 10.0 * x.size + np.sum(x**2 - 10.0 * np.cos(2.0 * np.pi * x))

def simulated_annealing(f, x0, T0=10.0, alpha=0.999, step_size=0.5,
                        n_iter=20000, seed=0):
    """Minimize f with simulated annealing.

    T0:        initial temperature
    alpha:     exponential cooling factor, T <- T * alpha each iteration
    step_size: std of the Gaussian proposal perturbation
    """
    rng = np.random.default_rng(seed)
    x = np.array(x0, dtype=float)
    fx = f(x)
    best_x, best_f = x.copy(), fx
    T = T0
    history = np.empty(n_iter)
    for k in range(n_iter):
        candidate = x + rng.normal(0.0, step_size, size=x.shape)
        f_cand = f(candidate)
        delta = f_cand - fx
        # Metropolis criterion: always accept improvements, sometimes
        # accept deteriorations with probability exp(-delta / T)
        if delta <= 0.0 or rng.random() < np.exp(-delta / T):
            x, fx = candidate, f_cand
            if fx < best_f:
                best_x, best_f = x.copy(), fx
        T = max(T * alpha, 1e-12)
        history[k] = best_f
    return best_x, best_f, history

def greedy_search(f, x0, step_size=0.2, n_iter=20000, seed=0):
    """Same proposal mechanism, but only accepts improvements (T = 0)."""
    rng = np.random.default_rng(seed)
    x = np.array(x0, dtype=float)
    fx = f(x)
    for _ in range(n_iter):
        candidate = x + rng.normal(0.0, step_size, size=x.shape)
        f_cand = f(candidate)
        if f_cand < fx:
            x, fx = candidate, f_cand
    return x, fx

# Single SA run from a distant start
x0 = np.array([4.5, 4.5])
best_x, best_f, history = simulated_annealing(rastrigin, x0, seed=3)
print(f"Start:          x = [4.5, 4.5],  f = {rastrigin(x0):.3f}")
print(f"SA best found:  x = [{best_x[0]:.4f}, {best_x[1]:.4f}],  f = {best_f:.6f}")
print(f"Global optimum: x = [0, 0],  f = 0")

# Fair SA vs greedy comparison: same proposals, same budget, 10 seeds
sa_finals = [simulated_annealing(rastrigin, x0, step_size=0.2, seed=s)[1]
             for s in range(10)]
greedy_finals = [greedy_search(rastrigin, x0, step_size=0.2, seed=s)[1]
                 for s in range(10)]
print(f"\nWith step_size = 0.2 (10 random seeds):")
print(f"SA:     mean best f = {np.mean(sa_finals):7.3f} +/- {np.std(sa_finals):.3f}")
print(f"Greedy: mean best f = {np.mean(greedy_finals):7.3f} +/- {np.std(greedy_finals):.3f}")

Output:

Start:          x = [4.5, 4.5],  f = 80.500
SA best found:  x = [0.0004, -0.0027],  f = 0.001511
Global optimum: x = [0, 0],  f = 0

With step_size = 0.2 (10 random seeds):
SA:     mean best f =   2.529 +/- 2.783
Greedy: mean best f =  35.820 +/- 6.748

Two observations:

• From a starting point with $f = 80.5$, SA reaches the global basin and finishes at $f \approx 0.0015$ — essentially the global optimum of a function with over a hundred local minima.
• The multi-seed comparison shows why: with small steps (0.2), greedy search gets stuck in a local minimum near the start (mean $f \approx 36$), while SA's occasional uphill moves let it hop between basins (mean $f \approx 2.5$).

One honest caveat: if you rerun the comparison with step_size=0.5, greedy search performs almost as well as SA on this particular problem, because a large Gaussian step can jump directly between Rastrigin's closely spaced basins. The advantage of accepting uphill moves is largest when the step size is small relative to the spacing of local minima — which, in high-dimensional real problems where you cannot see the landscape, is the common situation.

Comparing Cooling Schedules Experimentally

Which schedule should you use? Instead of trusting folklore, let us measure. We drive the same SA core with four precomputed temperature sequences and average over 10 seeds.

# Requirements:
# - Python 3.9+
# - numpy>=1.24.0

"""
Example 3: Cooling schedule comparison on the 2-D Rastrigin function

Purpose: Compare exponential, linear, logarithmic, and constant
         temperature schedules under an equal iteration budget
Target: Intermediate
Execution time: 1-2 minutes
"""

import numpy as np

def rastrigin(x):
    x = np.asarray(x, dtype=float)
    return 10.0 * x.size + np.sum(x**2 - 10.0 * np.cos(2.0 * np.pi * x))

def sa_with_schedule(f, x0, temps, step_size=0.5, seed=0):
    """Simulated annealing driven by a precomputed temperature array."""
    rng = np.random.default_rng(seed)
    x = np.array(x0, dtype=float)
    fx = f(x)
    best_f = fx
    for T in temps:
        candidate = x + rng.normal(0.0, step_size, size=x.shape)
        f_cand = f(candidate)
        delta = f_cand - fx
        if delta <= 0.0 or rng.random() < np.exp(-delta / max(T, 1e-12)):
            x, fx = candidate, f_cand
            best_f = min(best_f, fx)
    return best_f

K = 20000
T0 = 10.0
k = np.arange(K)
schedules = {
    "exponential (alpha=0.999)": T0 * 0.999**k,
    "linear": T0 * (1.0 - k / K),
    "logarithmic": T0 / np.log(k + 2.0),
    "constant T=1": np.full(K, 1.0),
}

n_seeds = 10
x0 = np.array([4.5, 4.5])
print(f"{'schedule':<28}{'mean best f':>14}{'std':>10}")
for name, temps in schedules.items():
    finals = [sa_with_schedule(rastrigin, x0, temps, seed=s)
              for s in range(n_seeds)]
    print(f"{name:<28}{np.mean(finals):>14.4f}{np.std(finals):>10.4f}")

Output:

schedule                       mean best f       std
exponential (alpha=0.999)           0.0069    0.0058
linear                              0.0290    0.0157
logarithmic                         0.0261    0.0276
constant T=1                        0.0164    0.0155

The exponential schedule wins here, consistent with common practice. Notice, though, that even a constant temperature performed respectably — better than linear and logarithmic on this problem. This is a useful reality check: on an easy 2-D benchmark, the choice of schedule matters less than getting the step size and temperature scale roughly right. Do not over-tune the schedule before checking the basics.

2.3 Genetic Algorithms

From One Solution to a Population

A genetic algorithm (GA) is a population-based metaheuristic inspired by natural selection. Instead of improving a single solution, a GA maintains a population of candidate solutions and evolves it over generations. The vocabulary comes from biology:

• Chromosome (individual): the encoding of one candidate solution — here, a binary vector.
• Fitness: a score measuring how good an individual is (higher is better, by convention).
• Selection: choosing parents, biased toward fitter individuals. We use tournament selection: pick $k$ individuals at random and keep the fittest one.
• Crossover (recombination): combining two parents into offspring. One-point crossover cuts both parents at a random position and swaps the tails.
• Mutation: randomly flipping bits with small probability, which maintains diversity and lets the population reach solutions that crossover alone cannot construct.
• Elitism: copying the best individual unchanged into the next generation, so the best-so-far fitness never decreases.

Why GAs Shine on Discrete Problems: Feature Selection

Feature selection — choosing which input columns a model should use — is a perfect GA showcase. With 20 features there are $2^{20} \approx 10^6$ possible subsets. The objective (validation error of a model trained on the chosen subset) is defined on binary vectors, so it has no gradient, and exhaustive search quickly becomes impossible as the feature count grows. A binary chromosome maps directly onto the problem: bit $i$ says whether feature $i$ is used.

We build a synthetic regression dataset where only 5 of 20 features carry signal, and check whether the GA finds exactly those 5.

# Requirements:
# - Python 3.9+
# - numpy>=1.24.0

"""
Example 4: Genetic algorithm for feature selection

Purpose: Implement selection, crossover, and mutation from scratch and
         recover the informative features of a synthetic dataset
Target: Intermediate
Execution time: 10-30 seconds
"""

import numpy as np

rng = np.random.default_rng(0)

# Synthetic regression data: 20 features, only 5 carry signal
n_samples, n_features = 300, 20
informative = [0, 3, 7, 12, 18]
X = rng.normal(size=(n_samples, n_features))
true_coef = np.zeros(n_features)
true_coef[informative] = [3.0, -2.0, 1.5, 2.5, -1.0]
y = X @ true_coef + rng.normal(0.0, 0.5, size=n_samples)

n_train = 200
X_train, X_val = X[:n_train], X[n_train:]
y_train, y_val = y[:n_train], y[n_train:]

def fitness(mask):
    """Negative validation MSE, with a small penalty per feature."""
    if mask.sum() == 0:
        return -np.inf
    cols = mask.astype(bool)
    coef, *_ = np.linalg.lstsq(X_train[:, cols], y_train, rcond=None)
    mse = np.mean((X_val[:, cols] @ coef - y_val) ** 2)
    return -mse - 0.02 * mask.sum()

def tournament_select(pop, scores, rng, k=3):
    """Pick k random individuals, return a copy of the fittest."""
    idx = rng.choice(len(pop), size=k, replace=False)
    return pop[idx[np.argmax(scores[idx])]].copy()

def one_point_crossover(p1, p2, rng):
    """Cut both parents at one point and swap the tails."""
    cut = rng.integers(1, len(p1))
    c1 = np.concatenate([p1[:cut], p2[cut:]])
    c2 = np.concatenate([p2[:cut], p1[cut:]])
    return c1, c2

def mutate(mask, rng, rate=0.05):
    """Flip each bit independently with probability `rate`."""
    flip = rng.random(len(mask)) < rate
    mask = mask.copy()
    mask[flip] = 1 - mask[flip]
    return mask

def genetic_algorithm(pop_size=40, n_generations=30, seed=1):
    rng = np.random.default_rng(seed)
    pop = (rng.random((pop_size, n_features)) < 0.5).astype(int)
    for gen in range(n_generations):
        scores = np.array([fitness(ind) for ind in pop])
        if gen % 10 == 0 or gen == n_generations - 1:
            print(f"Generation {gen:2d}: best fitness = {scores.max():.4f}, "
                  f"mean = {scores[np.isfinite(scores)].mean():.4f}")
        new_pop = [pop[np.argmax(scores)].copy()]  # elitism: keep the best
        while len(new_pop) < pop_size:
            p1 = tournament_select(pop, scores, rng)
            p2 = tournament_select(pop, scores, rng)
            c1, c2 = one_point_crossover(p1, p2, rng)
            new_pop.append(mutate(c1, rng))
            if len(new_pop) < pop_size:
                new_pop.append(mutate(c2, rng))
        pop = np.array(new_pop)
    scores = np.array([fitness(ind) for ind in pop])
    return pop[np.argmax(scores)], scores.max()

best_mask, best_score = genetic_algorithm()
selected = np.where(best_mask == 1)[0].tolist()
print(f"\nSelected features:    {selected}")
print(f"Informative features: {informative}")
print(f"All-features fitness: {fitness(np.ones(n_features, dtype=int)):.4f}")
print(f"GA best fitness:      {best_score:.4f}")

Output:

Generation  0: best fitness = -3.3562, mean = -13.1745
Generation 10: best fitness = -0.3652, mean = -1.8847
Generation 20: best fitness = -0.3652, mean = -1.2705
Generation 29: best fitness = -0.3652, mean = -1.6963

Selected features:    [0, 3, 7, 12, 18]
Informative features: [0, 3, 7, 12, 18]
All-features fitness: -0.6799
GA best fitness:      -0.3652

The GA recovers exactly the 5 informative features out of $2^{20}$ possible subsets, and its solution scores better than using all 20 features (which overfits the noise columns). Note the small per-feature penalty in the fitness function: without it, the GA has no incentive to drop harmless-but-useless features. Designing the fitness function is where most of your problem knowledge enters a GA.

2.4 Particle Swarm Optimization

Learning from Neighbors

Particle swarm optimization (PSO), introduced by Kennedy and Eberhart in 1995, is a population-based method inspired by bird flocking. Each particle $i$ has a position $\mathbf{x}_i$ (a candidate solution) and a velocity $\mathbf{v}_i$, and remembers the best position it has personally visited, $\mathbf{p}_i$. The swarm as a whole shares the best position anyone has found, $\mathbf{g}$. At each iteration:

$$ \mathbf{v}_i^{t+1} = \omega\, \mathbf{v}_i^{t} + c_1 r_1 \left(\mathbf{p}_i - \mathbf{x}_i^{t}\right) + c_2 r_2 \left(\mathbf{g} - \mathbf{x}_i^{t}\right) $$

$$ \mathbf{x}_i^{t+1} = \mathbf{x}_i^{t} + \mathbf{v}_i^{t+1} $$

where $r_1, r_2 \sim \mathcal{U}(0, 1)$ are fresh random numbers each step (drawn per dimension). The three terms have intuitive names:

• Inertia $\omega\, \mathbf{v}_i^t$: keep moving in the current direction. Typical $\omega \approx 0.4$–$0.9$; larger values favor exploration.
• Cognitive term $c_1 r_1 (\mathbf{p}_i - \mathbf{x}_i)$: pull back toward the particle's own best memory.
• Social term $c_2 r_2 (\mathbf{g} - \mathbf{x}_i)$: pull toward the swarm's best discovery. Typical $c_1 = c_2 \approx 1.5$–$2.0$.

Implementation with Trajectory Visualization

# Requirements:
# - Python 3.9+
# - numpy>=1.24.0
# - matplotlib>=3.7.0

"""
Example 5: Particle swarm optimization with trajectory visualization

Purpose: Implement PSO on the 2-D Rastrigin function and visualize how
         the swarm contracts onto the global optimum
Target: Intermediate
Execution time: 10-30 seconds
"""

import numpy as np
import matplotlib.pyplot as plt

def rastrigin(x):
    x = np.asarray(x, dtype=float)
    return 10.0 * x.size + np.sum(x**2 - 10.0 * np.cos(2.0 * np.pi * x))

def pso(f, bounds, n_particles=30, n_iter=100, omega=0.7,
        c1=1.5, c2=1.5, seed=2, record=False):
    """Minimize f with particle swarm optimization.

    bounds: list of (low, high) per dimension
    omega:  inertia weight
    c1, c2: cognitive and social acceleration coefficients
    """
    rng = np.random.default_rng(seed)
    dim = len(bounds)
    low = np.array([b[0] for b in bounds])
    high = np.array([b[1] for b in bounds])

    x = rng.uniform(low, high, size=(n_particles, dim))
    v = rng.uniform(-1.0, 1.0, size=(n_particles, dim))
    p_best = x.copy()
    p_best_f = np.array([f(xi) for xi in x])
    g_idx = np.argmin(p_best_f)
    g_best, g_best_f = p_best[g_idx].copy(), p_best_f[g_idx]

    trajectory = [x.copy()] if record else None
    for _ in range(n_iter):
        r1 = rng.random((n_particles, dim))
        r2 = rng.random((n_particles, dim))
        v = (omega * v
             + c1 * r1 * (p_best - x)
             + c2 * r2 * (g_best - x))
        x = np.clip(x + v, low, high)
        fx = np.array([f(xi) for xi in x])
        improved = fx < p_best_f
        p_best[improved] = x[improved]
        p_best_f[improved] = fx[improved]
        g_idx = np.argmin(p_best_f)
        if p_best_f[g_idx] < g_best_f:
            g_best, g_best_f = p_best[g_idx].copy(), p_best_f[g_idx]
        if record:
            trajectory.append(x.copy())
    return g_best, g_best_f, trajectory

bounds = [(-5.12, 5.12), (-5.12, 5.12)]
g_best, g_best_f, traj = pso(rastrigin, bounds, record=True)
print(f"PSO best: x = [{g_best[0]:.4f}, {g_best[1]:.4f}], f = {g_best_f:.6f}")

# Visualize the swarm at three moments
traj = np.array(traj)
xg = np.linspace(-5.12, 5.12, 200)
Xg, Yg = np.meshgrid(xg, xg)
Z = 20 + Xg**2 + Yg**2 - 10 * (np.cos(2*np.pi*Xg) + np.cos(2*np.pi*Yg))

fig, axes = plt.subplots(1, 3, figsize=(15, 4.5))
for ax, it in zip(axes, [0, 10, 100]):
    ax.contourf(Xg, Yg, Z, levels=30, cmap="viridis")
    ax.scatter(traj[it][:, 0], traj[it][:, 1],
               c="white", edgecolors="black", s=30)
    ax.set_title(f"Iteration {it}")
    ax.set_xlabel("$x_1$")
    ax.set_ylabel("$x_2$")
plt.tight_layout()
plt.savefig("pso_trajectory.png", dpi=120)
plt.show()

Output:

PSO best: x = [-0.0000, 0.0000], f = 0.000000

The three panels of the saved figure tell the story: at iteration 0 the particles are scattered uniformly; by iteration 10 they cluster around a few promising basins; by iteration 100 the swarm has contracted onto the global optimum at the origin. This collapse is PSO's strength and its weakness — once the swarm has contracted, it has little ability to explore elsewhere. If the global best found early is in the wrong basin, the whole swarm can converge there prematurely. Larger $\omega$, more particles, or restarting the swarm are the standard remedies.

2.5 Differential Evolution with SciPy

The Algorithm

Differential evolution (DE), proposed by Storn and Price in 1997, is a population-based method for continuous spaces that is remarkably effective for its simplicity. For each individual $\mathbf{x}_i$ in the population, DE builds a mutant vector from three other randomly chosen individuals:

$$ \mathbf{v}_i = \mathbf{x}_{r_1} + F \cdot \left(\mathbf{x}_{r_2} - \mathbf{x}_{r_3}\right) $$

where $F \in (0, 2]$ is the mutation factor (also called differential weight). The clever part: the perturbation $\mathbf{x}_{r_2} - \mathbf{x}_{r_3}$ is a difference of population members, so its scale and direction automatically adapt to the current spread of the population — large exploratory steps early, small refining steps as the population converges. The mutant is then mixed with $\mathbf{x}_i$ coordinate-by-coordinate with crossover probability $CR$ (binomial crossover), and the trial vector replaces $\mathbf{x}_i$ only if it is at least as good.

Using scipy.optimize.differential_evolution

You do not need to implement DE yourself — SciPy ships a well-tested implementation. Here is a complete example on the 5-dimensional Rastrigin function, together with the fairest gradient-based competitor: L-BFGS-B restarted from 20 random points.

# Requirements:
# - Python 3.9+
# - numpy>=1.24.0
# - scipy>=1.11.0

"""
Example 6: Differential evolution with SciPy vs multi-start L-BFGS-B

Purpose: Solve the 5-D Rastrigin problem with
         scipy.optimize.differential_evolution and compare cost and
         quality against a restarted gradient method
Target: Intermediate
Execution time: 10-30 seconds
"""

import numpy as np
from scipy.optimize import differential_evolution, minimize

def rastrigin(x):
    x = np.asarray(x, dtype=float)
    return 10.0 * x.size + np.sum(x**2 - 10.0 * np.cos(2.0 * np.pi * x))

dim = 5
bounds5 = [(-5.12, 5.12)] * dim

# --- Differential evolution ---
result = differential_evolution(
    rastrigin,
    bounds5,
    strategy="best1bin",     # mutation strategy: best + 1 difference, binomial crossover
    maxiter=300,             # generations
    popsize=20,              # population = popsize * dim individuals
    mutation=(0.5, 1.0),     # F drawn from this range each generation (dithering)
    recombination=0.7,       # crossover probability CR
    tol=1e-8,
    seed=7,
    polish=True,             # finish with a local L-BFGS-B refinement
)
print("=== Differential Evolution (5-D Rastrigin) ===")
print(f"Best f:  {result.fun:.8f}")
print(f"Best x:  {np.round(result.x, 5)}")
print(f"Function evaluations: {result.nfev}")

# --- Multi-start L-BFGS-B: the fair gradient-based baseline ---
rng = np.random.default_rng(0)
finals, total_nfev = [], 0
for _ in range(20):
    x0 = rng.uniform(-5.12, 5.12, size=dim)
    res = minimize(rastrigin, x0=x0, method="L-BFGS-B", bounds=bounds5)
    finals.append(res.fun)
    total_nfev += res.nfev
finals = np.array(finals)
print("\n=== L-BFGS-B, 20 random restarts ===")
print(f"Best of 20 restarts: {finals.min():.4f}")
print(f"Median result:       {np.median(finals):.4f}")
print(f"Runs reaching f < 1: {(finals < 1).sum()} / 20")
print(f"Total evaluations:   {total_nfev}")

Output:

=== Differential Evolution (5-D Rastrigin) ===
Best f:  0.00000000
Best x:  [ 0.  0. -0.  0.  0.]
Function evaluations: 17906

=== L-BFGS-B, 20 random restarts ===
Best of 20 restarts: 10.9445
Median result:       37.3107
Runs reaching f < 1: 0 / 20
Total evaluations:   1980

This is the fair comparison the chapter promised, and both sides of it matter:

• DE finds the exact global optimum; none of the 20 gradient restarts even gets close (best: $f \approx 10.9$). On highly multimodal functions, restarting a local method is not enough.
• But DE spent about 9× more function evaluations (17,906 vs 1,980). If each evaluation were a 30-minute model training run, that difference is days of compute. Metaheuristics buy global search ability with evaluations.

Method Comparison Table

2.6 Practical Guidance

The Optimizers Have Hyperparameters Too

There is an irony at the heart of this chapter: we often reach for metaheuristics to tune hyperparameters, yet metaheuristics have hyperparameters of their own ($T_0$ and $\alpha$ for SA; population size, crossover and mutation rates for GA; $\omega$, $c_1$, $c_2$ for PSO; $F$ and $CR$ for DE). These settings genuinely matter. Let us measure how much, by running SciPy's DE under a fixed budget of roughly 15,000 evaluations with different settings:

# Requirements:
# - Python 3.9+
# - numpy>=1.24.0
# - scipy>=1.11.0

"""
Example 7: Sensitivity of DE to its own hyperparameters

Purpose: Show that optimizer settings change results substantially even
         under an identical evaluation budget
Target: Intermediate
Execution time: 1-3 minutes
"""

import numpy as np
from scipy.optimize import differential_evolution

def rastrigin(x):
    x = np.asarray(x, dtype=float)
    return 10.0 * x.size + np.sum(x**2 - 10.0 * np.cos(2.0 * np.pi * x))

dim = 5
bounds5 = [(-5.12, 5.12)] * dim
budget = 15000  # approximate function-evaluation budget

def run_de_budget(mutation, recombination, popsize, seed):
    # population per generation = popsize * dim, so cap generations
    maxiter = budget // (popsize * dim) - 1
    res = differential_evolution(
        rastrigin, bounds5,
        mutation=mutation, recombination=recombination,
        popsize=popsize, maxiter=maxiter, tol=0, seed=seed,
        polish=False,   # no local refinement: measure DE alone
    )
    return res.fun

settings = [
    ("default   (F=0.5-1.0, CR=0.7, pop=15)", (0.5, 1.0), 0.7, 15),
    ("high CR   (F=0.5-1.0, CR=0.95, pop=15)", (0.5, 1.0), 0.95, 15),
    ("low F     (F=0.2, CR=0.7, pop=15)", 0.2, 0.7, 15),
    ("large pop (F=0.5-1.0, CR=0.7, pop=40)", (0.5, 1.0), 0.7, 40),
]

print(f"Budget: about {budget} evaluations, 5-D Rastrigin, 10 seeds\n")
print(f"{'setting':<42}{'mean f':>10}{'std':>8}")
for name, mut, cr, pop in settings:
    finals = [run_de_budget(mut, cr, pop, seed=s) for s in range(10)]
    print(f"{name:<42}{np.mean(finals):>10.4f}{np.std(finals):>8.4f}")

Output:

Budget: about 15000 evaluations, 5-D Rastrigin, 10 seeds

setting                                       mean f     std
default   (F=0.5-1.0, CR=0.7, pop=15)         0.0995  0.2985
high CR   (F=0.5-1.0, CR=0.95, pop=15)        1.4924  1.7937
low F     (F=0.2, CR=0.7, pop=15)             5.7708  3.9247
large pop (F=0.5-1.0, CR=0.7, pop=40)         0.9058  0.5534

Under the identical budget, results span nearly two orders of magnitude. Three practical lessons:

• Start from library defaults. SciPy's defaults (dithered $F$, $CR = 0.7$) were the best setting here — they encode a lot of accumulated experience.
• Too-small mutation factors are the classic failure. $F = 0.2$ makes steps too small to escape local structure; the population converges prematurely.
• Bigger populations are not free. A larger population explores more per generation but runs fewer generations under a fixed budget. Population size trades breadth for depth.

Budget Considerations

The single most important number in gradient-free optimization is your evaluation budget: how many times can you afford to evaluate the objective? Plan your method around it:

Two budget rules that prevent most self-deception:

1. Compare methods at equal evaluation counts, never at equal iteration counts. One GA "generation" costs a whole population of evaluations; one SA "iteration" costs one. Example 6 counted nfev for exactly this reason.
2. Report multiple seeds. Metaheuristics are stochastic; a single lucky run proves nothing. The mean-and-standard-deviation tables in Examples 3 and 7 are the honest format.

The No-Free-Lunch Theorem

The no-free-lunch (NFL) theorem for optimization (Wolpert & Macready, 1997) states that, averaged over all possible objective functions, every optimization algorithm — including pure random search — performs exactly the same. No optimizer is universally best.

It is worth being precise about what this does and does not mean:

• What it means: any method that beats random search on some class of problems must lose to it on another. When a paper or vendor claims an optimizer is "the best", the only meaningful question is: best on which problem class?
• What it does not mean: that the choice of method is hopeless in practice. Real objective functions are not drawn uniformly from all possible functions — they have structure (smoothness, low effective dimensionality, clustered good regions). An optimizer succeeds precisely by exploiting structure that your problem actually has: DE exploits continuity, GAs exploit the tendency of good partial solutions to combine well, SA exploits the fact that good solutions are near other good solutions.

The practical corollary: match the method's assumptions to your problem's structure, and always keep random search as a baseline. If your carefully tuned metaheuristic does not clearly beat random search at the same budget, it is not exploiting any structure — use the simpler method.

2.7 Chapter Summary

What We Learned

1. When gradient-free methods are needed

   • Discrete/piecewise-constant objectives: the gradient is zero or undefined
   • Noisy evaluations: finite differences amplify the noise catastrophically
   • Multimodal landscapes: gradients lead to the nearest local minimum, not the best one

2. Simulated annealing

   • Metropolis criterion: accept worse solutions with probability $\exp(-\Delta E / T)$
   • Cooling schedules: exponential is the practical default; the theoretically guaranteed logarithmic schedule is too slow to use
   • Accepting uphill moves is what lets SA escape local minima that trap greedy search

3. Genetic algorithms

   • Population + tournament selection + crossover + mutation + elitism
   • Natural fit for discrete encodings — our GA recovered the exact informative feature subset out of $2^{20}$ candidates
   • The fitness function is where your problem knowledge lives

4. Particle swarm optimization

   • Velocity update = inertia + cognitive pull + social pull
   • Fast convergence, but the swarm can contract prematurely onto the wrong basin

5. Differential evolution and fair comparison

   • Difference-vector mutation self-adapts step sizes; scipy.optimize.differential_evolution is a strong default
   • DE solved 5-D Rastrigin exactly where 20 gradient restarts all failed — but used 9× the evaluations
   • Compare at equal evaluation budgets, across multiple seeds, against a random-search baseline
   • No free lunch: every optimizer wins only by exploiting problem structure

Next Chapter

Metaheuristics assume evaluations are cheap enough to spend by the thousand. When each evaluation is a model training run or a laboratory experiment, we need a method that squeezes maximum information out of every single evaluation. In Chapter 3, you will learn Bayesian optimization:

• Gaussian process surrogate models
• Acquisition functions (expected improvement, upper confidence bound)
• The exploration-exploitation trade-off, made explicit
• Hyperparameter tuning with realistic evaluation budgets

Exercises

Problem 1 (Difficulty: easy)

The proposal step size is as important as the cooling schedule in simulated annealing. Run the simulated_annealing function from Example 2 on the 2-D Rastrigin function with step_size values 0.05, 0.5, and 3.0 (keeping all other settings at their defaults, seed=1). Report the best objective value for each, and explain the pattern you observe.

# Requirements:
# - Python 3.9+
# - numpy>=1.24.0

import numpy as np

# (reuse rastrigin and simulated_annealing from Example 2)

def rastrigin(x):
    x = np.asarray(x, dtype=float)
    return 10.0 * x.size + np.sum(x**2 - 10.0 * np.cos(2.0 * np.pi * x))

def simulated_annealing(f, x0, T0=10.0, alpha=0.999, step_size=0.5,
                        n_iter=20000, seed=0):
    rng = np.random.default_rng(seed)
    x = np.array(x0, dtype=float)
    fx = f(x)
    best_x, best_f = x.copy(), fx
    T = T0
    for k in range(n_iter):
        candidate = x + rng.normal(0.0, step_size, size=x.shape)
        f_cand = f(candidate)
        delta = f_cand - fx
        if delta <= 0.0 or rng.random() < np.exp(-delta / T):
            x, fx = candidate, f_cand
            if fx < best_f:
                best_x, best_f = x.copy(), fx
        T = max(T * alpha, 1e-12)
    return best_x, best_f

for ss in [0.05, 0.5, 3.0]:
    _, bf = simulated_annealing(rastrigin, np.array([4.5, 4.5]),
                                step_size=ss, seed=1)
    print(f"SA step_size={ss}: best f = {bf:.4f}")

SA step_size=0.05: best f = 24.8739
SA step_size=0.5: best f = 0.0008
SA step_size=3.0: best f = 0.0632

Problem 2 (Difficulty: medium)

Explain the role of the temperature $T$ in the Metropolis criterion $P = \exp(-\Delta E / T)$. What does simulated annealing become in the two limits $T \to \infty$ and $T \to 0$? Using this, explain why a good cooling schedule must be neither too fast nor too slow.

Problem 3 (Difficulty: hard)

Perform a fair comparison of simulated annealing, PSO, differential evolution, and random search on the 2-D Ackley function, another standard multimodal benchmark with global minimum $f(\mathbf{0}) = 0$:

$$ f(\mathbf{x}) = -20 \exp\!\left(-0.2\sqrt{\tfrac{1}{n}\sum_i x_i^2}\right) - \exp\!\left(\tfrac{1}{n}\sum_i \cos(2\pi x_i)\right) + 20 + e $$

Give every method the same budget of about 6,000 function evaluations, run 5 seeds each, and report mean and standard deviation of the best objective value. Which methods solve the problem, and what do you conclude?

# Requirements:
# - Python 3.9+
# - numpy>=1.24.0
# - scipy>=1.11.0

import numpy as np
from scipy.optimize import differential_evolution

def ackley(x):
    x = np.asarray(x, dtype=float)
    n = x.size
    return (-20.0 * np.exp(-0.2 * np.sqrt(np.sum(x**2) / n))
            - np.exp(np.sum(np.cos(2.0 * np.pi * x)) / n)
            + 20.0 + np.e)

def simulated_annealing(f, x0, T0=10.0, alpha=0.999, step_size=0.5,
                        n_iter=20000, seed=0):
    rng = np.random.default_rng(seed)
    x = np.array(x0, dtype=float)
    fx = f(x)
    best_x, best_f = x.copy(), fx
    T = T0
    for k in range(n_iter):
        candidate = x + rng.normal(0.0, step_size, size=x.shape)
        f_cand = f(candidate)
        delta = f_cand - fx
        if delta <= 0.0 or rng.random() < np.exp(-delta / T):
            x, fx = candidate, f_cand
            if fx < best_f:
                best_x, best_f = x.copy(), fx
        T = max(T * alpha, 1e-12)
    return best_x, best_f

def pso(f, bounds, n_particles=30, n_iter=100, omega=0.7,
        c1=1.5, c2=1.5, seed=2):
    rng = np.random.default_rng(seed)
    dim = len(bounds)
    low = np.array([b[0] for b in bounds])
    high = np.array([b[1] for b in bounds])
    x = rng.uniform(low, high, size=(n_particles, dim))
    v = rng.uniform(-1.0, 1.0, size=(n_particles, dim))
    p_best = x.copy()
    p_best_f = np.array([f(xi) for xi in x])
    g_idx = np.argmin(p_best_f)
    g_best, g_best_f = p_best[g_idx].copy(), p_best_f[g_idx]
    for _ in range(n_iter):
        r1 = rng.random((n_particles, dim))
        r2 = rng.random((n_particles, dim))
        v = omega * v + c1 * r1 * (p_best - x) + c2 * r2 * (g_best - x)
        x = np.clip(x + v, low, high)
        fx = np.array([f(xi) for xi in x])
        improved = fx < p_best_f
        p_best[improved] = x[improved]
        p_best_f[improved] = fx[improved]
        g_idx = np.argmin(p_best_f)
        if p_best_f[g_idx] < g_best_f:
            g_best, g_best_f = p_best[g_idx].copy(), p_best_f[g_idx]
    return g_best, g_best_f

bounds = [(-5.0, 5.0)] * 2
n_seeds = 5
budget = 6000

sa_res = [simulated_annealing(ackley, np.array([4.0, -4.0]),
                              n_iter=budget, seed=s)[1]
          for s in range(n_seeds)]
pso_res = [pso(ackley, bounds, n_particles=30,
               n_iter=budget // 30, seed=s)[1]
           for s in range(n_seeds)]
de_res = []
for s in range(n_seeds):
    r = differential_evolution(ackley, bounds, popsize=15,
                               maxiter=budget // (15 * 2) - 1,
                               tol=0, seed=s, polish=False)
    de_res.append(r.fun)
rs_res = []
for s in range(n_seeds):
    rg = np.random.default_rng(s)
    cands = rg.uniform(-5, 5, size=(budget, 2))
    rs_res.append(min(ackley(c) for c in cands))

for name, res in [("SA", sa_res), ("PSO", pso_res),
                  ("DE", de_res), ("random", rs_res)]:
    print(f"{name:<8} mean best f = {np.mean(res):.4f} +/- {np.std(res):.4f}")

SA       mean best f = 0.0362 +/- 0.0176
PSO      mean best f = 0.0000 +/- 0.0000
DE       mean best f = 0.0000 +/- 0.0000
random   mean best f = 0.3056 +/- 0.1684

References

1. Kirkpatrick, S., Gelatt, C. D., & Vecchi, M. P. (1983). Optimization by Simulated Annealing. Science, 220(4598), 671-680.
2. Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H., & Teller, E. (1953). Equation of State Calculations by Fast Computing Machines. Journal of Chemical Physics, 21(6), 1087-1092.
3. Holland, J. H. (1975). Adaptation in Natural and Artificial Systems. University of Michigan Press.
4. Kennedy, J., & Eberhart, R. (1995). Particle Swarm Optimization. Proceedings of ICNN'95, 1942-1948.
5. Storn, R., & Price, K. (1997). Differential Evolution - A Simple and Efficient Heuristic for Global Optimization over Continuous Spaces. Journal of Global Optimization, 11, 341-359.
6. Wolpert, D. H., & Macready, W. G. (1997). No Free Lunch Theorems for Optimization. IEEE Transactions on Evolutionary Computation, 1(1), 67-82.
7. Geman, S., & Geman, D. (1984). Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images. IEEE Transactions on Pattern Analysis and Machine Intelligence, 6(6), 721-741.