Evolution Strategies

Evolution Strategies

You are an expert in evolution strategies (ES) for continuous, integer, and mixed-integer optimization in operations research. This skill covers the $(\mu/\rho \,\overset{+}{,}\, \lambda)$ framework, the 1/5 success rule, log-normal self-adaptation, cumulative step-size adaptation (CSA), CMA-ES, restart strategies, and integer/mixed-integer handling, plus the main combinatorial entry point: continuous relaxations such as random keys. Use the framework below to choose an ES variant, implement it in clean numpy, and apply it where ES earns its place in OR practice — algorithm-parameter tuning, continuous subproblems, and simulation optimization.

Initial Assessment

Establish these facts before writing any ES code:

• Search-space type. Continuous, integer, mixed-integer, or genuinely combinatorial (permutation, subset, assignment)? ES is native to continuous spaces. For combinatorial structures, decide early between (a) a continuous relaxation with a decoder (random keys) and (b) a different metaheuristic operating on the native encoding — option (b) usually wins (see solution-encodings).
• Dimension n. CMA-ES is the default for $n \lesssim 100$; per-generation cost grows as $O(\lambda n^2)$ with an amortized $O(n^3)$ eigendecomposition. Above a few hundred dimensions, switch to separable/diagonal variants.
• Evaluation cost and budget. Count total affordable evaluations. CMA-ES needs roughly $100n$ to $1000n$ evaluations to show its strength. If the budget is under ~$50n$ (expensive simulations), a model-based tuner (see optuna-hyperparameter-tuning) is usually a better fit.
• Noise. Is the objective deterministic, or stochastic (a simulation, or a randomized algorithm's output)? Noise dictates population sizing, reevaluation policy, and use of common random numbers.
• Gradients. If the objective is differentiable and gradients are cheap, use a gradient method first. ES is for black-box objectives: nonsmooth, noisy, simulation-based, or rugged.
• Bounds and constraints. Box bounds only, or general constraints? Decide per constraint: repair (clip/project), penalty, or resample. Box bounds are routine; general constraints need explicit design.
• Integer or categorical coordinates. Mark every integer coordinate now: rounding inside the objective plus a step-size floor handles them, but only if planned from the start (see the mixed-integer section below). Unordered categorical parameters have no meaningful Gaussian neighborhood — their presence in volume is a signal that irace or Optuna fits better than ES.
• Scaling of variables. Note each variable's natural range and whether it lives on a log scale (rates, temperatures, penalty weights). Plan a normalization map to $[0,1]^n$ or $[0,10]^n$ before optimizing.
• Multimodality expectation. Unimodal-ish (refinement task) suggests a (1+1)-ES or plain CMA-ES; rugged landscapes suggest larger $\lambda$ and IPOP/BIPOP restarts.
• Parallelism. Can $\lambda$ candidates be evaluated concurrently? ES is embarrassingly parallel within a generation; this often decides $\lambda$.
• Quality requirement and baseline. Target precision (e.g., $10^{-8}$ on a benchmark, or "beats default parameters by 2%") and an existing baseline to compare against (random search, default configuration, a local optimizer).
• Reproducibility. Seeds per run, number of repetitions, and the reporting format the results must feed into.

Algorithm Anatomy

The (mu/rho +, lambda) framework

An ES maintains $\mu$ parents. Each generation it creates $\lambda$ offspring; each offspring recombines $\rho$ parents and mutates the result by adding Gaussian noise:

$$ x' = \text{recombine}(x{i1}, \dots, x{i\rho}) + \sigma \odot z, \qquad z \sim \mathcal{N}(0, C). $$

Selection is the defining choice (Beyer & Schwefel 2002, "Evolution strategies — a comprehensive introduction"):

| Scheme | Survivor pool | Character | |---|---|---| | $(\mu, \lambda)$ — comma | best $\mu$ of the $\lambda$ offspring only | Non-elitist; forgets parents. Required for reliable self-adaptation — mis-adapted step sizes die out. Needs $\lambda \gtrsim 5\mu$. | | $(\mu + \lambda)$ — plus | best $\mu$ of parents $\cup$ offspring | Elitist; never loses the incumbent. Safer with tiny budgets, but step sizes can lock up at a local optimum. | | $(1+1)$ | better of parent and child | Minimal ES; pairs with the 1/5 success rule. Strong cheap local refiner. |

What the theory buys you: progress rates

On the sphere model $f(x) = \|x\|^2$ at distance $R$ from the optimum, the (1+1)-ES achieves maximal expected progress at $\sigma^* \approx 1.224\, R/n$, and the resulting convergence is linear with rate $\Theta(1/n)$: each factor-of-ten improvement in $f$ costs $O(n)$ evaluations (Beyer 2001, "The Theory of Evolution Strategies"). Recombination helps: a $(\mu/\mu, \lambda)$-ES with intermediate recombination gains a speed-up of order $\mu$ from genetic-repair averaging of mutation noise — the reason CMA-ES recombines all $\mu$ selected parents with weights instead of pairs. Two practical consequences:

1. Budgets scale linearly with dimension on smooth problems — a 50-D problem needs roughly five times the evaluations of a 10-D one for the same precision. Use this to sanity-check whether a requested precision is affordable before running anything.
2. **Progress is sharply peaked around $\sigma^$.* A step size off by a factor of 10 cuts progress by roughly two orders of magnitude. This is why every serious ES adapts $\sigma$ online, and why a flat convergence curve almost always means broken step-size control rather than a hard problem.

Step-size control — the heart of ES

A fixed $\sigma$ fails: progress on a sphere-like region requires $\sigma$ proportional to the distance from the optimum. Three control mechanisms, in increasing sophistication:

1. 1/5 success rule (Rechenberg 1973, "Evolutionsstrategie"). On the sphere, the optimal success probability of a (1+1)-ES is ≈ 0.2. Count successes over a window; if the rate exceeds 1/5 increase $\sigma$, if below decrease it (factor $c \approx 0.85$ per window of $n$ trials; Schwefel 1981).
2. Log-normal self-adaptation (Schwefel 1995, "Evolution and Optimum Seeking"). Each individual carries its own step-size vector, mutated before the object variables:

$$ \sigmai' = \sigmai \cdot \exp(\tau' z0 + \tau zi), \qquad \tau' = \frac{1}{\sqrt{2n}}, \quad \tau = \frac{1}{\sqrt{2\sqrt{n}}}, $$

with one global draw $z0 \sim \mathcal{N}(0,1)$ shared by all coordinates and independent $zi$ per coordinate. Selection then implicitly favors individuals whose step sizes generated good moves. Works only under comma selection with adequate $\lambda/\mu$.

1. Cumulative step-size adaptation (CSA) (Ostermeier, Gawelczyk & Hansen 1994; Hansen & Ostermeier 2001, "Completely derandomized self-adaptation in evolution strategies"). Derandomized: instead of letting selection judge step sizes indirectly, accumulate the path of mean shifts $p_\sigma$ and compare its length against the expectation under random selection, $E\|\mathcal{N}(0, I)\| \approx \sqrt{n}\,(1 - \tfrac{1}{4n} + \tfrac{1}{21n^2})$:

$$ \sigma \leftarrow \sigma \cdot \exp\!\left(\frac{c\sigma}{d\sigma}\left(\frac{\|p_\sigma\|}{E\|\mathcal{N}(0,I)\|} - 1\right)\right). $$

A longer-than-random path means consecutive steps point the same way — increase $\sigma$; a shorter path means steps cancel — decrease it.

CMA-ES essentials

CMA-ES (Hansen & Ostermeier 2001; Hansen 2016, "The CMA Evolution Strategy: A Tutorial") additionally adapts a full covariance matrix $C$, learning the local metric of the problem (variable scaling and correlations). Per generation: sample $xk = m + \sigma\, B D zk$ with $C = B D^2 B^\top$, rank by objective, move the mean by weighted recombination of the $\mu$ best, and update $C$ by a rank-one term (evolution path $p_c$) plus a rank-$\mu$ term:

$$ C \leftarrow (1 - c1 - c\mu)\, C + c1\, pc pc^\top + c\mu \sum{i=1}^{\mu} wi\, y{i:\lambda} y{i:\lambda}^\top . $$

CMA-ES is invariant to order-preserving transformations of the objective (it only uses ranks) and to affine transformations of the search space (given matching initialization) — the reason it solves ill-conditioned, non-separable problems that defeat isotropic or per-coordinate step-size ES.

Choosing a variant

| Situation | Use | |---|---| | $n \le 100$, non-separable or ill-conditioned, budget ≥ $100n$ evaluations | CMA-ES (default choice) | | Cheap local refinement of a good solution; near-unimodal region | (1+1)-ES with 1/5 rule | | Very cheap evaluations, simple code wanted, massive parallel evaluation | $(\mu, \lambda)$ self-adaptive ES | | Rugged/multimodal landscape | CMA-ES + IPOP/BIPOP restarts | | $n$ in the hundreds-thousands | sep-CMA-ES, diagonal variants, OpenAI-ES-style NES | | Budget under ~$50n$ evaluations | Bayesian/TPE tuning instead (see optuna-hyperparameter-tuning) | | Native permutation/subset structure with good local moves | A discrete metaheuristic, not ES; keep ES via random keys only as a baseline |

Parameter guidance

| Parameter | Typical setting | Increasing it buys | At the cost of | |---|---|---|---| | $\lambda$ | $4 + \lfloor 3\ln n\rfloor$ (CMA-ES); $\ge 5\mu$ (self-adaptive comma-ES) | Global search, noise robustness | Evaluations per generation; slower convergence per evaluation on unimodal $f$ | | $\mu$ | $\lambda/2$ with log-decreasing weights (CMA); $\lambda/7$–$\lambda/4$ (SA-ES) | Smoother mean updates, robustness | Lower selection pressure | | $\rho$ | 2, or $\mu$ (global intermediate) | Averaging cancels mutation noise | Loss of population diversity | | $\sigma_0$ | 0.2–0.5 × variable range; start point in the domain core | Early exploration, escape from bad init | Overshooting; wasted early evaluations | | $\tau, \tau'$ | $1/\sqrt{2\sqrt{n}}$, $1/\sqrt{2n}$ | Faster step-size learning | Step-size noise, premature collapse | | Selection (+ vs ,) | comma for self-adaptation; plus for tiny budgets | (+): monotone incumbent | (+): stagnating step sizes on multimodal $f$ | | Restart multiplier | ×2 population per restart (IPOP) | Systematic global search | Budget split across restarts |

Mutation-operator internals (Gaussian vs Cauchy, correlated mutations, discrete perturbations) are covered in mutation-and-perturbation-operators; representation choices and decoder design in solution-encodings. Use those skills rather than re-deriving operators here.

Reusable ES Engine

The engine is problem-independent: the objective maps an $(m, n)$ array of candidate rows to an $(m,)$ array of values (minimization), so fitness evaluation is one vectorized call per generation.

SELF-ADAPTIVE (mu/rho +, lambda)-ES  — minimization
  initialize mu parents: x ~ U(bounds), per-coordinate sigma = 0.3 * range
  evaluate parents
  repeat until evaluation budget exhausted:
      for the lambda offspring (vectorized):
          pick rho distinct parents; intermediate recombination of x and sigma
          sigma  ESResult:
    """(mu/rho +, lambda)-ES with log-normal self-adaptation of per-coordinate step sizes.

    objective maps an (m, n) array of candidates to an (m,) array of values
    (minimization). bounds is an (n, 2) array of [lower, upper] per coordinate.
    """
    rng = np.random.default_rng(seed)
    n = bounds.shape[0]
    lower, upper = bounds[:, 0], bounds[:, 1]
    span = upper - lower
    tau_global = 1.0 / np.sqrt(2.0 * n)
    tau_coord = 1.0 / np.sqrt(2.0 * np.sqrt(n))
    sigma_floor = 1e-12

    X = lower + rng.random((mu, n)) * span          # parent object variables
    S = np.tile(0.3 * span, (mu, 1))                # parent step-size vectors
    F = objective(X)
    evals = mu

    best = int(np.argmin(F))
    x_best, f_best = X[best].copy(), float(F[best])
    history = [f_best]

    while evals  tuple[np.ndarray, float]:
    """(1+1)-ES with Rechenberg's 1/5 success rule, window of n mutations.

    objective uses the batch signature (m, n) -> (m,) for consistency with
    the population engine; here m == 1.
    """
    rng = np.random.default_rng(seed)
    x = np.asarray(x0, dtype=float).copy()
    n = x.size
    f = float(objective(x[None, :])[0])
    sigma = float(sigma0)
    successes = 0
    for k in range(1, max_evaluations):
        y = x + sigma * rng.standard_normal(n)
        fy = float(objective(y[None, :])[0])
        if fy  0.2 else (sigma * c if rate  np.ndarray:
    """f(x) = sum x_i^2; minimum 0 at the origin. Separable, unimodal."""
    return np.sum(X**2, axis=1)

def rosenbrock(X: np.ndarray) -> np.ndarray:
    """Banana valley; minimum 0 at (1, ..., 1). Non-separable, ill-conditioned."""
    return np.sum(
        100.0 * (X[:, 1:] - X[:, :-1] ** 2) ** 2 + (1.0 - X[:, :-1]) ** 2, axis=1
    )

def rastrigin(X: np.ndarray) -> np.ndarray:
    """10n + sum(x_i^2 - 10 cos(2 pi x_i)); ~10^n local optima. Multimodal."""
    return 10.0 * X.shape[1] + np.sum(X**2 - 10.0 * np.cos(2 * np.pi * X), axis=1)

Run the engine and contrast plus vs comma selection — the experiment every ES user should do once:

import numpy as np

# Uses self_adaptive_es, sphere, rosenbrock, rastrigin from the blocks above.

bounds_10 = np.tile(np.array([-5.12, 5.12]), (10, 1))

for name, f in [("sphere", sphere), ("rosenbrock", rosenbrock), ("rastrigin", rastrigin)]:
    for plus in (False, True):
        vals = [
            self_adaptive_es(
                f, bounds_10, mu=15, lam=100, plus_selection=plus,
                max_evaluations=50_000, seed=s,
            ).f_best
            for s in range(5)
        ]
        tag = "(15+100)" if plus else "(15,100)"
        print(f"{name:10s} {tag}: median {np.median(vals):.3e}  best {min(vals):.3e}")
# Expected: sphere solved to ~1e-8 by both; on rosenbrock both crawl along the
# valley (final f roughly 1e-1 to 1e1) because per-coordinate step sizes cannot
# represent the rotating curved metric — the motivation for CMA-ES; on rastrigin
# the comma version reaches lower medians than plus, which stagnates earlier
# because elitist step sizes shrink at the first decent local optimum.

Interpretation guidance: if comma and plus perform identically, $\lambda/\mu$ is too small for self-adaptation to matter. If sphere convergence is not log-linear (a straight line on a semi-log convergence plot), the step-size mechanism is broken — debug that before trusting any multimodal result.

Worked Example 2: Tuning a Simulated-Annealing Solver with CMA-ES

Parameter tuning is the highest-value ES application in combinatorial optimization: the tuning space is small ($n$ = 2–10), continuous or mixed, noisy, and each evaluation is expensive (full solver runs). CMA-ES handles all four properties. Alternatives: irace (López-Ibáñez et al. 2016, "The irace package") for categorical-heavy spaces, and Optuna's TPE for pruning-friendly setups (see optuna-hyperparameter-tuning).

First, the CMA-ES itself — minimal but complete, following Hansen (2016):

CMA-ES — one generation
  sample    z_k ~ N(0, I);  y_k = B D z_k;  x_k = m + sigma * y_k   (k = 1..lambda)
  rank by f;  y_w = sum_{i=1..mu} w_i * y_{i:lambda}
  mean      m  CMAResult:
    """CMA-ES (Hansen 2016 tutorial parameterization, positive weights only).

    objective maps an (m, n) array of candidates to an (m,) array of values.
    """
    rng = np.random.default_rng(seed)
    m = np.asarray(x0, dtype=float).copy()
    n = m.size
    if lam is None:
        lam = 4 + int(3 * np.log(n))
    mu = lam // 2
    w = np.log(mu + 0.5) - np.log(np.arange(1, mu + 1))
    w /= w.sum()
    mu_eff = 1.0 / np.sum(w**2)

    c_sigma = (mu_eff + 2) / (n + mu_eff + 5)
    d_sigma = 1 + 2 * max(0.0, np.sqrt((mu_eff - 1) / (n + 1)) - 1) + c_sigma
    c_c = (4 + mu_eff / n) / (n + 4 + 2 * mu_eff / n)
    c_1 = 2 / ((n + 1.3) ** 2 + mu_eff)
    c_mu = min(1 - c_1, 2 * (mu_eff

…

## Source & license

This open-source skill is cataloged on AgentStack and links to its original source — we do not rehost the code.

- **Author:** [hajibabaie](https://github.com/hajibabaie)
- **Source:** [hajibabaie/combinatorial-optimization-skills](https://github.com/hajibabaie/combinatorial-optimization-skills)
- **License:** MIT

Install and usage instructions live in the source repository linked above.