Black-box optimization, visualized · metaheuristics from the Delta competition

00 The black box we were handed

Minimize f : ℝ^d → ℝ, given nothing but the right to query it.

All of this came out of Delta 2026, a competition. Task D: retime the traffic signals of a real city so a fleet of vehicles loses as little time as possible over a four-hour simulation. There was no analytic objective. There was a mesoscopic traffic simulator. You submit a full signal plan, it runs the city, and it returns a single scalar (total delay). Lower is better.

341

signals (the decision vector)

5,215

vehicles simulated

14,400 s

simulated per evaluation

~5 s

wall-clock per query

What "black box" means here

You hand the objective a full decision vector and it hands back one number. You do not get a formula, a derivative, or any look inside. The only way to learn anything is to ask, and every question costs about five seconds of simulation:

fig. The interface in full: a 341-number plan goes in, one scalar comes out. Nothing inside is exposed, so the classical tools that need $\nabla f$, such as gradient descent, Newton's method, or linear programming, do not apply here.

That makes the objective black-box, high-dimensional, and expensive all at once. You cannot enumerate (the configuration space is astronomically large), you cannot follow a gradient (there is none to compute), and you cannot afford many evaluations. This is the exact niche of metaheuristics: stochastic search strategies that trade optimality guarantees for a very good solution within a fixed budget. Each one is a different policy for the same dilemma, exploration versus exploitation: sampling unknown regions against refining the best point found so far. The sections below give each policy as an update rule, a schematic, and a demo you can drive.

01 Why local descent is not enough

Start with a picture of what goes wrong.

Greedy descent, and where it stops

Suppose, just for a moment, that you could differentiate. The obvious move is gradient descent: from wherever you are, take a small step straight downhill, and repeat until the slope flattens to zero.

$$ x \;\leftarrow\; x - \eta\, f'(x), \qquad\text{stop when } f'(x) = 0 $$

$\eta$ is the step size. The update only ever moves toward lower $f$, so it can stop at any point where the slope vanishes. That might be the global minimum, but it might just as easily be a local one.

Below is a 1-D objective with several dips of different depths. Each gray point starts at a fixed spot and runs exactly that update. The coloured trail traces where it ends up. Every point reaches a minimum, but only the nearest one, with no knowledge that a deeper basin exists elsewhere.

fig. Thirteen greedy descents on one bumpy objective. Each settles into the basin it happened to start in; only the few that began inside the global basin (green) reach the true optimum, the rest stop short (red).

The basin of attraction each point lands in is decided entirely by where it started. For a convex objective there is only one basin and this is fine. For a multimodal one, local descent is just a lucky dip into whatever basin you happened to initialise inside.

This is the design problem. Every method below adds a mechanism to escape or avoid local minima: a population that samples many basins at once, stochastic acceptance of uphill moves, a self-adapting step, or a probabilistic model of where to look next. Keep the explore/exploit tension in mind; it is the axis every one of them is tuned along.

02 Particle Swarm Optimization

A population of velocity-driven particles, coupled through a shared best.

The idea

Scatter a handful of particles across the search space and let each one move under its own velocity. A particle remembers the best point it has visited, and the swarm as a whole keeps the best point any particle has found. On each step a particle is pulled toward both of those points while keeping some of the velocity it already had. That is the whole method: three pulls on a moving point, and no gradient anywhere.

One particle's update

Particle $i$ carries a position $\mathbf{x}_i$ and a velocity $\mathbf{v}_i$. Each iteration the velocity is rebuilt from three terms: the inertia it already had, a cognitive pull toward its own best $\mathbf{p}_i$, and a social pull toward the global best $\hat{\mathbf{g}}$. The position then follows the velocity:

$$ \mathbf{v}_i \;\leftarrow\; \underbrace{\textcolor{#f0b72f}{w\,\mathbf{v}_i}}_{\text{inertia}} \;+\; \underbrace{\textcolor{#56d364}{c_1\, r_1 \odot (\mathbf{p}_i - \mathbf{x}_i)}}_{\text{cognitive}} \;+\; \underbrace{\textcolor{#5ec8ff}{c_2\, r_2 \odot (\hat{\mathbf{g}} - \mathbf{x}_i)}}_{\text{social}} $$ $$ \mathbf{x}_i \;\leftarrow\; \mathbf{x}_i + \mathbf{v}_i $$

$r_1, r_2 \sim U(0,1)^d$ are drawn fresh per coordinate and $\odot$ is elementwise, so the two pulls jitter independently along each axis. Stable defaults: $c_1 = c_2 \approx 2$, with $w$ on a decreasing schedule from $\sim\!0.9$ down to $\sim\!0.4$.

How the velocity is assembled

Geometrically the update is nothing but vector addition. Lay the three contributions tip-to-tail, and the new velocity is the single arrow that closes the path from the old position to the new one:

fig. 1. The new velocity $\mathbf{v}_i$ is the tip-to-tail sum of carried momentum (gold), the pull toward the particle's own best $\mathbf{p}_i$ (green), and the pull toward the swarm best $\hat{\mathbf{g}}$ (blue). Their resultant carries the particle from $\mathbf{x}_i$ to $\mathbf{x}_i + \mathbf{v}_i$. The random weights $r_1, r_2$ rescale the two coloured pulls every step.

What each knob does

The dial that sets exploration against exploitation is the inertia $w$. High $w$ keeps velocities large and the swarm roams; low $w$ lets the two pulls win and the swarm collapses onto $\hat{\mathbf{g}}$. Making $c_2$ large relative to $c_1$ turns every particle into a follower of the global best. That converges fast, but it invites premature convergence: the swarm agrees before it has explored.

live. The two bold arrows are drawn on one highlighted particle, so you can watch its velocity being assembled; the faint blue arrows are the velocity field $\mathbf{v}_i$ of the whole swarm. The side panel tracks the swarm spread, which collapses as the social term takes over.

Drive it Push inertia $w$ above 1 on Rastrigin: velocities grow and the swarm refuses to settle. Drop it to 0.3 and it collapses into the nearest basin. Zero out $c_1$ and max $c_2$ to watch the spread collapse as every particle chases $\hat{\mathbf{g}}$. Use step to inspect a single iteration.

03 Genetic Algorithm

Selection, recombination, mutation, elitism. A real-coded evolutionary loop.

The loop

A population of candidate solutions evolves generation by generation. Each new generation is built by three operators in order: pick promising parents, blend them into children, then perturb the children a little. A few of the very best survive untouched (elitism), so progress is never lost.

one generationfor each offspring slot (after copying the elite): p₁, p₂ ← tournament(k), tournament(k) # selection child ← BLX-α(p₁, p₂) # recombination child ← child + 𝒟(0, σ²) with prob p_mut # mutation next generation ← elite ∪ offspring

Selection: a tournament

To pick one parent, draw $k$ individuals at random and keep the fittest of them. Nothing more. The tournament size $k$ is the single dial of selection pressure: $k=1$ is a random pick (no pressure), while large $k$ almost always returns a top individual (harsh pressure, fast convergence, more risk of collapsing onto a local optimum).

fig. 1. Three individuals are sampled (gold rings); the greenest wins and becomes a parent. Raising $k$ samples more each time, so the winner is almost always near the top, with more pressure and less diversity.

Recombination: BLX-α crossover

Given two parents, each child coordinate is drawn uniformly from the interval the parents span on that axis, widened by a fraction $\alpha$ on both sides. The widening is what lets a child fall a little outside the parents, so the population can still reach values neither parent had:

$$ \text{child}_j \;\sim\; U\big[\, m_j - \alpha d_j,\ \ M_j + \alpha d_j \,\big], \qquad m_j=\min(p_{1j},p_{2j}),\ \ M_j=\max(p_{1j},p_{2j}),\ \ d_j = M_j - m_j $$

fig. 2. The green band is the parent interval $[m_j, M_j]$; the blue band extends it by $\alpha d_j$ on each side. The child is drawn uniformly across the whole blue band, so it can land beyond either parent (here, to the left). That mild extrapolation keeps the search from shrinking too fast.

Mutation & elitism

Finally each child coordinate is, with probability $p_{\text{mut}}$, nudged by Gaussian noise of scale $\sigma$. This is the GA's exploration valve: too small and the population stalls, too large and it never settles.

$$ \text{child} \;\leftarrow\; \text{child} + \mathcal{N}(0,\sigma^2) \quad\text{per coordinate, with probability } p_{\text{mut}} $$

live. The cloud is colour-coded by fitness; one mating per generation is highlighted (two gold parents give a green child). The panel's diversity reading is the spread that lets a GA cover several basins at once.

Drive it Send mutation $\sigma$ toward 0: the population converges quickly, then stalls, losing the diversity it needs to improve. Raise it and it never settles. Increase tournament $k$ to intensify selection pressure (faster convergence, higher risk of a local optimum). On Rastrigin the GA's diversity routinely finds the global basin where a greedier method settles.

04 Differential Evolution

Mutation by the scaled difference of two random members. Self-scaling, near parameter-free.

Mutation: a difference vector

DE's signature is that its step size is not a parameter at all. It is read off the population itself. To mutate, pick three other members at random and add a $F$-scaled copy of the difference between two of them to the third:

$$ \mathbf{v} \;=\; \mathbf{x}_a \;+\; \textcolor{#f0b72f}{F}\,\textcolor{#d2a8ff}{(\mathbf{x}_b - \mathbf{x}_c)} $$

fig. 1. The donor $\mathbf{v}$ is the base point $\mathbf{x}_a$ plus the difference $\mathbf{x}_b - \mathbf{x}_c$ scaled by $F$. That difference is large when the population is spread out and small once it has clustered, so the step self-scales: big exploratory moves early, fine refinement late, with no schedule to tune.

Crossover: binomial mixing

The donor does not replace the target wholesale. Coordinate by coordinate, a trial vector $\mathbf{u}$ takes the donor's value with probability $CR$ and otherwise keeps the target's. One coordinate $j_{\text{rand}}$ is forced to come from the donor, so the trial always differs from the target in at least one place:

$$ u_j \;=\; \begin{cases} \textcolor{#56d364}{v_j} & \text{if } \text{rand}_j \le CR \ \text{ or }\ j = j_{\text{rand}} \\[2pt] x_{ij} & \text{otherwise} \end{cases} $$

fig. 2. Per coordinate: $\checkmark$ means $\text{rand}\le CR$, so the trial copies the donor; $\times$ means it keeps the target. The gold $\star$ is the forced coordinate $j_{\text{rand}}$. High $CR$ takes mostly donor (good when coordinates are separable); low $CR$ stays close to the target (good when they are coupled).

Selection: greedy and elitist

The trial is evaluated once and replaces the target only if it is no worse. Nothing ever regresses, so the population's best is monotone:

$$ \mathbf{x}_i \;\leftarrow\; \mathbf{u} \quad\text{if}\quad f(\mathbf{u}) \le f(\mathbf{x}_i) $$

live. The donor construction is drawn on one target each step: the purple difference vector and the green donor. The panel's mean $|\mathbf{x}_b - \mathbf{x}_c|$ is the step size, and you can watch it decay on its own as the cloud contracts.

Drive it Switch to Rosenbrock: DE threads the ill-conditioned valley well, because its difference vectors line up with the population's local geometry. $F$ sets step aggressiveness; $CR$ the mixing rate (high $CR$ suits separable problems, low $CR$ suits coupled coordinates).

05 Simulated Annealing

A single Markov chain with a temperature-controlled Metropolis acceptance.

The Metropolis rule

One walker proposes a random neighbour each step. If the move improves the objective ($\Delta E \le 0$) it is always taken. A worsening move is not rejected outright. It is accepted with a probability that shrinks as the move gets worse and grows with the temperature $T$:

$$ P(\text{accept}) = \begin{cases} 1 & \Delta E \le 0 \\[2pt] \textcolor{#5ec8ff}{\exp(-\Delta E / T)} & \Delta E > 0 \end{cases}, \qquad \Delta E = f(x_{\text{new}}) - f(x_{\text{cur}}) $$

That exponential is what lets annealing climb out of a local minimum: while the temperature is high it walks uphill out of a basin fairly often, and once it is low it almost never does. The shape of the curve says the rest:

fig. Acceptance probability for a worsening move. Improving moves ($\Delta E \le 0$, green shelf) are always taken. For uphill moves, a hot chain still says yes fairly often; a cold one almost always says no. Cooling is just sliding from the top curve down to the bottom one over the run.

Cooling

The temperature is lowered geometrically, so the chain explores freely early and commits late. This schedule is the one knob that really matters:

$$ T \;\leftarrow\; \alpha\, T, \qquad 0 \ll \alpha < 1 $$

Too fast (small $\alpha$) and the chain quenches into the first basin it finds; too slow and it wastes the evaluation budget wandering. Annealing buys the ability to escape local minima with a single point, and carries no population.

live. The panel shows $\Delta E$ and $P(\text{accept})$ each step. The yellow trail is the chain; red flickers are rejected proposals. Early on, while hot, you will catch it accepting clearly uphill moves.

Drive it Set cooling $\alpha$ to a fast 0.90: the chain freezes early, often in a local minimum. Push it to 0.999 and it explores far longer before committing. Watch the live $P(\text{accept})$ collapse toward 0 as $T$ falls.

06 CMA-ES

Adapt a sampling distribution instead of moving points around. Usually the best choice on smooth problems.

Sample, rank, re-fit

Instead of tracking individual points, CMA-ES tracks a distribution: a multivariate Gaussian $\mathcal{N}(\mathbf{m}, \sigma^2 \mathbf{C})$ over the search space. Each generation it draws $\lambda$ candidates from that Gaussian, ranks them by objective value, and slides the mean toward a weighted average of the best $\mu$:

$$ \mathbf{x}_k = \mathbf{m} + \sigma\,\mathcal{N}(\mathbf{0}, \mathbf{C}), \quad k = 1,\dots,\lambda $$ $$ \mathbf{m} \;\leftarrow\; \sum_{k=1}^{\mu} w_k\, \mathbf{x}_{k:\lambda} \qquad (w_1 > w_2 > \dots > w_\mu > 0,\ \ \textstyle\sum_k w_k = 1) $$

$\mathbf{x}_{k:\lambda}$ is the $k$-th best candidate of the generation. Only the ranking matters, never the raw values. That is what makes CMA-ES invariant to any monotone rescaling of $f$.

Learning the landscape's metric

The interesting part is the covariance $\mathbf{C}$. Each generation it is tilted toward the directions that just produced improvement. On a curved, ill-conditioned valley the sampling ellipse stops being a round blob and becomes a long, thin, rotated shape aligned with the valley floor. In effect the algorithm has learned the local second-order metric without ever computing a gradient:

fig. The same Gaussian over three generations on a curved valley. It begins round (gen 1), elongates and rotates as the covariance picks up the productive direction (gen 5), and ends as a long thin ellipse tracking the valley toward the optimum (gen 12).

Concretely, $\mathbf{C}$ is a blend of its old self, a rank-one term that reinforces the path the mean has been travelling, and a rank-$\mu$ term built from this generation's best steps; the overall scale $\sigma$ is adapted separately:

$$ \mathbf{C} \;\leftarrow\; (1-c_1-c_\mu)\,\mathbf{C} \;+\; c_1\,\underbrace{\textcolor{#56d364}{\mathbf{p}_c\,\mathbf{p}_c^{\top}}}_{\text{rank-one}} \;+\; c_\mu \underbrace{\textcolor{#5ec8ff}{\textstyle\sum_k w_k\, \mathbf{y}_k\mathbf{y}_k^{\top}}}_{\text{rank-}\mu} $$ $$ \sigma \;\leftarrow\; \sigma\,\exp\!\Big(\tfrac{c_\sigma}{d_\sigma}\Big(\tfrac{\lVert \mathbf{p}_\sigma\rVert}{\mathbb{E}\lVert\mathcal{N}(0,\mathbf{I})\rVert} - 1\Big)\Big) $$

live. The filled teal ellipse is the $1\sigma$ contour of the sampling Gaussian and the faint ring is $2\sigma$; the dots are one generation's $\lambda$ samples, and the dashed arrow is the step the mean just took. On Rosenbrock the ellipse starts round and slowly turns into a thin diagonal sliver lying along the valley. The panel's axis ratio is the conditioning the covariance has picked up so far.

Drive it Start on Rosenbrock and step through a dozen generations: the ellipse rotates to follow the curved valley while the mean walks down it. A small $\sigma_0$ makes it creep; a large one overshoots for a few generations until the step size adapts back down. More samples $\lambda$ give a steadier covariance estimate, at a higher cost per generation.

Cost and regime The covariance update and its eigendecomposition are $O(d^2)$–$O(d^3)$, so CMA-ES targets low-to-moderate dimension ($d \lesssim 100$) where sample efficiency dominates per-step cost. On smooth, ill-conditioned objectives it is usually the best of the lot: in our benchmarks it drove Rosenbrock to $\sim\!10^{-23}$, far past every population method. On highly multimodal landscapes, though, a single adapting distribution can still get trapped. The full standard $(\mu/\mu_w, \lambda)$ algorithm is about 90 lines in metaheuristics/cma_es.py.

07 Bayesian Optimization

A probabilistic surrogate plus an acquisition function. For when evaluations are the bottleneck.

The surrogate loop

When a single evaluation is genuinely expensive, say a simulation, a training run, or a lab experiment, you can afford only tens of them. Bayesian optimization makes each one count: it fits a cheap probabilistic model of the objective and uses that model to choose where the next evaluation should go.

fig. 1. The loop. Everything except the gold box is cheap arithmetic. You spend that arithmetic on the model and the acquisition so the one expensive call lands somewhere worth it.

The Gaussian-process posterior

Condition a Gaussian process on the observations $(\mathbf{X}, \mathbf{y})$. With kernel $k$ and noise $\sigma_n^2$, the posterior at any new point $x$ is itself Gaussian, with a closed-form mean and variance that share one Gram-matrix inverse:

$$ \mu(x) = \textcolor{#56d364}{\mathbf{k}(x)^{\top}(\mathbf{K}+\sigma_n^2\mathbf{I})^{-1}\mathbf{y}} $$ $$ \sigma^2(x) = k(x,x) - \textcolor{#5ec8ff}{\mathbf{k}(x)^{\top}(\mathbf{K}+\sigma_n^2\mathbf{I})^{-1}\mathbf{k}(x)} $$

$\mathbf{k}(x)_i = k(x, x_i)$. At an observation, $\mu$ is pinned to the data and $\sigma$ collapses to $\sim\!0$. Far from any data $\sigma$ swells, and that swelling is exactly what the acquisition will chase.

The kernel encodes the smoothness prior. The common RBF choice makes nearby inputs strongly correlated, with a length scale $\ell$ that sets how fast that correlation falls off, that is, how wiggly the GP believes $f$ is:

$$ k(x,x') = \exp\!\big(-\lVert x-x'\rVert^2 / 2\ell^2\big) $$

fig. 2. The RBF kernel as a function of distance. A small $\ell$ de-correlates points quickly, so the GP trusts its own interpolation only very locally; a large $\ell$ assumes a smooth function and extrapolates confidently.

Choosing where to look: Expected Improvement

Given the posterior, the next query maximizes an acquisition function. Expected Improvement over the incumbent $f_{\text{best}}$ splits cleanly into an exploit term and an explore term:

$$ \mathrm{EI}(x) = \underbrace{\textcolor{#56d364}{(f_{\text{best}}-\xi-\mu)\,\Phi(z)}}_{\text{exploit}} \;+\; \underbrace{\textcolor{#5ec8ff}{\sigma\,\phi(z)}}_{\text{explore}}, \qquad z=\frac{f_{\text{best}}-\xi-\mu}{\sigma} $$

$\Phi, \phi$ are the standard-normal CDF and PDF; $\xi$ is a margin that tilts the balance. The first term rewards a low predicted mean, the second rewards high uncertainty, so EI keeps probing unexplored regions instead of only the current best.

live. The dashed curve is the hidden objective, blue is the GP mean, the band is $\pm 2\sigma$, and green is EI. Each sample next evaluates $f$ at the EI maximizer (white line) and you watch the band collapse there.

Drive it Hit sample next repeatedly and watch the band tighten where data lands. Raise $\xi$ to bias toward exploration. Drop the length scale $\ell$ and the GP assumes a wigglier function, so it trusts its interpolation less. This is the machinery behind most hyperparameter tuners.

08 Choosing: there is no free lunch

Averaged over all objectives, every optimizer ties. A method wins only by matching structure.

The No Free Lunch theorem (Wolpert & Macready, 1997) is precise about this: no black-box optimizer is better than any other across all possible objectives. So the question is never "which is best" but "what is known about this landscape and this budget".

Match the method to the landscape

Two axes capture most of the decision: how much structure the landscape has (smooth and unimodal versus rugged and multimodal), and how big the evaluation budget is. Each method has a corner where it is the natural pick:

fig. Where each optimizer is the natural choice. Bayesian optimization fits the few-evaluations regime; CMA-ES the smooth, ill-conditioned middle; DE the generic centre; GA and PSO the rugged, multimodal corner; SA and hill-climbing the cheap, single-point baselines. The full decision rule:

Regime	Method	Rationale
expensive evaluations	Bayesian optimization	surrogate model, O(tens) of evals
smooth, ill-conditioned, d ≲ 100	CMA-ES	learns the landscape metric
generic continuous, unknown structure	differential evolution	robust, self-scaling default
rugged, highly multimodal	GA / PSO	population resists local traps
cheap evals, no population budget	simulated annealing	single chain, escapes via Metropolis
baseline to beat	random-restart hill climbing	if nothing beats it, stay simple

Where the real gains come from In practice the largest gains rarely come from the optimizer. They come from a good representation, a warm start from a domain heuristic instead of random initialization, and a cheap surrogate that buys more evaluations. The optimizer is the last ten percent. See docs/04_results.md for the full benchmark table behind this.

09 The twist

We deployed the entire arsenal. The problem turned out to be nearly trivial.

Back to Delta. We threw all of this at the traffic objective: a three-population PSO (exploratory / balanced / exploitative, sharing a global best), GAs over the signal timings, DE and SA on a fast analytical proxy of the simulator, Bayesian optimization with random embeddings to fight the dimensionality, even a Navier-Stokes fluid analogy for green-time allocation.

Then we reverse-engineered the scoring Once the official simulator's exact cost was understood, the verdict was clear: with the given demand the network was so under-saturated that almost any sane plan was within a fraction of a percent of optimal. "All green, never stop anyone" was essentially the answer. The optimizer tournament was contesting the last decimal of a problem that barely was one.

The ceiling was flat

The objective did have an enormous cost range; you can build a genuinely terrible plan. But every reasonable plan landed in a thin band just above the optimum, so all our sophistication moved the score by less than a percent:

fig. The cost range is huge (a bad plan sits far above), but the do-nothing "all-green" baseline already lands in the thin band just over the optimum, and the full optimizer arsenal only shaves a sliver off it. The tournament was fought inside that $<\!1\%$ gap.

Which is the most useful thing the competition could have taught us, and exactly the point of the callout above: characterize the objective before optimizing it. We spent days tuning swarms; an afternoon of demand analysis would have shown the ceiling was flat.

The contest solution was disposable. The toolbox was not. We pulled the algorithms out of the competition code, reimplemented them cleanly against standard benchmark functions where their behavior is legible, tested them, and built this page so the mechanics stick. Everything here runs the real algorithm, live.

See it on the Delta problem There is a separate interactive page that runs the traffic-signal problem itself: a small city, the optimizer working live, and the flat-ceiling finding as a chart you can sweep. Open the interactive Delta demo →

📦 Code, tests, and this site live in one repo: metaheuristics-lab. The optimizers are short, documented modules, and the full Delta application (the traffic-signal problem, the solver, and the page above) lives in delta/. Run python examples/compare_optimizers.py to reproduce the benchmark leaderboard.