Offensive Security Researcher & Quantum Cryptography Analyst

CMA-ES Algorithm: Working, Explanation, and Evolutionary Strategy Optimization

Abdul Wahab Junaid
CMA-ES algorithm and working of this algorithm

I consider CMA-ES, the Covariance Matrix Adaptation Evolution Strategy, to be one of the most mathematically sophisticated algorithms I have studied in the evolutionary computation family, and also one of the most effective for difficult continuous optimization problems. Unlike simpler evolutionary strategies that use a fixed or diagonal mutation distribution, CMA-ES learns and adapts a full covariance matrix that shapes how new candidate solutions are sampled around the current best estimate. I find this powerful because it allows the search distribution to stretch, rotate, and reshape itself to match the local geometry of the objective function — effectively learning the “shape” of the problem as it searches, which lets it handle ill-conditioned, non-separable problems that trip up many other black-box optimizers.

History and Background

I attribute CMA-ES primarily to Nikolaus Hansen and Andreas Ostermeier, who developed it in the mid-1990s while working on evolution strategies at the Technical University of Berlin, publishing key papers in 1996 and a more complete version in 2001 (“Completely Derandomized Self-Adaptation in Evolution Strategies,” Evolutionary Computation journal). CMA-ES grew out of the broader Evolution Strategies (ES) tradition established by Rechenberg and Schwefel in the 1960s and 70s, which already used self-adaptive mutation step sizes; Hansen and Ostermeier’s key innovation was “derandomizing” the adaptation process and extending it from a single global step size to a full covariance matrix, letting the algorithm learn correlations between variables. Since its introduction, CMA-ES has become a standard reference algorithm in black-box optimization benchmarking (such as the COCO/BBOB benchmark suite), and it’s widely regarded, in my reading of the literature, as one of the most robust and reliable derivative-free optimizers for continuous problems, particularly when the number of dimensions is moderate (roughly up to a few hundred).

Problem Statement

I turn to CMA-ES for continuous black-box optimization problems: minimize $f(x)$ for $x \in \mathbb{R}^n$, where $f$ may be non-convex, non-separable, ill-conditioned, noisy, and where I have no access to gradients. Many simpler optimizers assume the variables are roughly independent (separable) or that the landscape is well-scaled, but real problems frequently have variables that interact strongly and are scaled very differently from each other. The problem CMA-ES solves is finding good solutions on exactly these difficult, poorly-scaled, correlated landscapes, by adaptively learning the right search distribution rather than assuming one in advance.

Core Concepts

How It Works

I run CMA-ES through the following cycle at each generation:

  1. I sample $\lambda$ candidate solutions from the current multivariate normal distribution $\mathcal{N}(m, \sigma^2 C)$.
  2. I evaluate the objective function for each candidate.
  3. I rank the candidates by fitness and select the best $\mu$ of them.
  4. I update the mean $m$ as a weighted average of these $\mu$ best candidates.
  5. I update the evolution paths $p_\sigma$ and $p_c$, which track the recent history of successful search directions.
  6. I update the covariance matrix $C$ using information from both the evolution path and the spread of the selected candidates (rank-one and rank-$\mu$ updates).
  7. I update the step size $\sigma$ by comparing the length of the evolution path $p_\sigma$ to what would be expected under random selection (cumulative step-size adaptation).
  8. I repeat steps 1 through 7 for further generations until a stopping criterion, such as a target fitness or a maximum number of function evaluations, is reached.

Working Principle

I think of the internal logic of CMA-ES as continuously reshaping and repositioning a “cloud” of Gaussian-distributed sample points to match the local landscape of the objective function. If the best solutions tend to consistently appear in a particular direction relative to the mean, the covariance matrix is stretched along that direction, so future sampling naturally explores more efficiently along productive axes and less along unproductive ones. The evolution path mechanism is particularly clever: rather than just looking at whether the most recent step was good, it accumulates a running memory of consistent successful directions across several generations, which lets the algorithm detect and exploit correlations that a single-generation snapshot would miss, and helps distinguish genuine progress from random fluctuation. The step-size adaptation, separately, controls the overall scale — growing when steps are consistently moving in a similar direction (suggesting I can move faster) and shrinking when steps are inconsistent or oscillating (suggesting I’m near a local structure that requires finer search).

Mathematical Foundation

At generation $g+1$, I sample $\lambda$ candidates as:

$$ x_k^{(g+1)} \sim m^{(g)} + \sigma^{(g)} \cdot \mathcal{N}\left(0, C^{(g)}\right), \quad k = 1, \dots, \lambda $$

The mean update uses the best $\mu$ candidates with recombination weights $w_i$ (with $\sum w_i = 1$):

$$ m^{(g+1)} = \sum_{i=1}^{\mu} w_i \cdot x_{i:\lambda}^{(g+1)} $$

The evolution path for step-size control is updated as:

$$ p_\sigma^{(g+1)} = (1-c_\sigma) p_\sigma^{(g)} + \sqrt{c_\sigma(2-c_\sigma)\mu_{eff}} \cdot C^{(g)^{-1/2}} \frac{m^{(g+1)} – m^{(g)}}{\sigma^{(g)}} $$

The step size is then updated by comparing the path length to its expected value under random selection:

$$ \sigma^{(g+1)} = \sigma^{(g)} \cdot \exp\left(\frac{c_\sigma}{d_\sigma}\left(\frac{|p_\sigma^{(g+1)}|}{E|\mathcal{N}(0,I)|} – 1\right)\right) $$

The covariance matrix evolution path is:

$$ p_c^{(g+1)} = (1-c_c) p_c^{(g)} + \sqrt{c_c(2-c_c)\mu_{eff}} \cdot \frac{m^{(g+1)} – m^{(g)}}{\sigma^{(g)}} $$

And the covariance matrix update combines a rank-one update (using the evolution path) and a rank-$\mu$ update (using the spread of the selected population):

$$ C^{(g+1)} = (1-c_1-c_\mu) C^{(g)} + c_1 \cdot p_c^{(g+1)} \left(p_c^{(g+1)}\right)^T + c_\mu \sum_{i=1}^{\mu} w_i \left(\frac{x_{i:\lambda}^{(g+1)} – m^{(g)}}{\sigma^{(g)}}\right)\left(\frac{x_{i:\lambda}^{(g+1)} – m^{(g)}}{\sigma^{(g)}}\right)^T $$

where $c_1$, $c_\mu$, $c_c$, $c_\sigma$, $d_\sigma$ are fixed learning-rate-like constants derived from the problem dimension $n$ and $\mu_{eff}$ (the effective selection mass), with well-established default formulas provided in Hansen’s tutorials.

Diagrams

flowchart TD
    A["Initialize parameters"] --> B["Generate candidate solutions"]
    B --> C["Evaluate fitness"]
    C --> D["Select best candidates"]
    D --> E["Update evolution paths"]
    E --> F["Update covariance matrix"]
    F --> G["Update step size"]
    G --> H{"Stopping condition met?"}
    H -- No --> B
    H -- Yes --> I["Return best solution"]

Pseudocode

function CMA_ES(f, n, lambda, mu, max_generations):
    m = random_initial_point(n)
    sigma = initial_step_size
    C = identity_matrix(n)
    p_sigma = zeros(n)
    p_c = zeros(n)
    weights = compute_recombination_weights(mu)

    for generation in 1 to max_generations:
        for k in 1 to lambda:
            z_k = sample_multivariate_normal(0, C)
            x_k = m + sigma * z_k
            fitness_k = f(x_k)

        sort candidates x_1..x_lambda by fitness ascending
        m_new = sum(weights[i] * x_i for i in 1 to mu)

        p_sigma = update_sigma_path(p_sigma, C, m_new, m, sigma)
        sigma = sigma * exp((c_sigma/d_sigma) * (norm(p_sigma)/expected_norm - 1))

        p_c = update_c_path(p_c, m_new, m, sigma)
        C = update_covariance(C, p_c, x_1..x_mu, m, sigma, weights)

        m = m_new

        if stopping_condition(fitness_values):
            break

    return m

Step-by-Step Example

I’ll sketch a simplified two-dimensional example minimizing $f(x) = x_1^2 + 10 x_2^2$ (an ill-conditioned, elongated bowl), starting at $m=(5,5)$, $\sigma=1$, $C=I$.

Time Complexity

Each generation requires evaluating $\lambda$ candidates, costing $O(\lambda \cdot C_f)$ where $C_f$ is the cost of one fitness evaluation. Sampling from the multivariate normal distribution and updating the covariance matrix both involve operations on an $n \times n$ matrix, typically costing $O(n^2)$ per sample for sampling (using a precomputed decomposition) and $O(n^2 \cdot \mu)$ for the covariance update, plus periodic eigen-decomposition of $C$, which costs $O(n^3)$ but is usually only performed every few generations to amortize its cost. So overall, per-generation cost is roughly $O(\lambda \cdot C_f + n^3)$, making CMA-ES most practical for moderate dimensionality (up to a few hundred variables) since the cubic term in $n$ becomes the bottleneck for very high-dimensional problems.

Space Complexity

I need to store the covariance matrix $C$, which takes $O(n^2)$ space, along with the mean vector, evolution paths, and step size, which are $O(n)$. I also need to store the $\lambda$ sampled candidates and their fitness values each generation, costing $O(\lambda \cdot n)$. Overall space complexity is $O(n^2 + \lambda \cdot n)$, with the $n^2$ covariance matrix term typically dominating for larger dimensions.

Correctness Analysis

I approach CMA-ES’s correctness primarily through its strong empirical and theoretical convergence properties rather than a simple guarantee like DE’s greedy selection. On convex quadratic functions, CMA-ES has been shown (through theoretical analysis of its adaptation dynamics) to achieve linear convergence — meaning the distance to the optimum shrinks by a roughly constant factor each generation — and importantly, this convergence rate is largely invariant to the conditioning of the problem, because the covariance matrix adaptation effectively “undoes” ill-conditioning over time by learning the appropriate transformation. This invariance property (technically called affine invariance) is a key theoretical strength: CMA-ES’s behavior on a linearly transformed version of a function is essentially identical to its behavior on the original, up to the transformation, which is a much stronger property than most other black-box optimizers can claim. As with all stochastic global optimizers, there’s no absolute guarantee of finding the global optimum on arbitrary multimodal landscapes in finite time, and CMA-ES can still converge to a local optimum, but restarts with increasing population size (a technique called IPOP-CMA-ES) are commonly used to improve global search reliability.

Advantages

Disadvantages

Applications

Implementation in C

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

// Simplified CMA-ES for a 2D sphere-like problem, illustrating the core update loop.
// A full production CMA-ES needs matrix eigen-decomposition; here C stays diagonal
// for clarity and to keep the implementation self-contained.

#define N 2            // dimensions
#define LAMBDA 8        // population size per generation
#define MU 4             // number of parents selected
#define MAX_GEN 100

double objective_function(double *x) {
    return x[0]*x[0] + 10.0 * x[1]*x[1];  // ill-conditioned bowl
}

double gaussian() {
    double u1 = (double) rand() / RAND_MAX;
    double u2 = (double) rand() / RAND_MAX;
    return sqrt(-2.0 * log(u1)) * cos(2.0 * M_PI * u2);
}

int main() {
    srand((unsigned) time(NULL));

    double m[N] = {5.0, 5.0};
    double sigma = 1.0;
    double C_diag[N] = {1.0, 1.0};   // diagonal covariance (simplified)
    double weights[MU];
    double weight_sum = 0.0;

    // log-based recombination weights, standard CMA-ES choice
    for (int i = 0; i < MU; i++) {
        weights[i] = log(MU + 0.5) - log(i + 1);
        weight_sum += weights[i];
    }
    for (int i = 0; i < MU; i++) weights[i] /= weight_sum;

    double candidates[LAMBDA][N];
    double fitness[LAMBDA];
    int indices[LAMBDA];

    for (int gen = 0; gen < MAX_GEN; gen++) {
        // Step 1: sample candidates
        for (int k = 0; k < LAMBDA; k++) {
            for (int j = 0; j < N; j++)
                candidates[k][j] = m[j] + sigma * sqrt(C_diag[j]) * gaussian();
            fitness[k] = objective_function(candidates[k]);
            indices[k] = k;
        }

        // Step 2: sort candidate indices by fitness ascending (simple bubble sort)
        for (int a = 0; a < LAMBDA - 1; a++)
            for (int b = 0; b < LAMBDA - a - 1; b++)
                if (fitness[indices[b]] > fitness[indices[b+1]]) {
                    int tmp = indices[b]; indices[b] = indices[b+1]; indices[b+1] = tmp;
                }

        // Step 3: compute new mean from best MU candidates
        double new_m[N] = {0.0, 0.0};
        for (int i = 0; i < MU; i++)
            for (int j = 0; j < N; j++)
                new_m[j] += weights[i] * candidates[indices[i]][j];

        // Step 4: simplified variance update (diagonal approximation of C)
        double variance[N] = {0.0, 0.0};
        for (int i = 0; i < MU; i++)
            for (int j = 0; j < N; j++) {
                double diff = (candidates[indices[i]][j] - m[j]) / sigma;
                variance[j] += weights[i] * diff * diff;
            }
        for (int j = 0; j < N; j++)
            C_diag[j] = 0.8 * C_diag[j] + 0.2 * variance[j];  // smoothed update

        // Step 5: simplified step-size adaptation based on improvement
        double step_norm = 0.0;
        for (int j = 0; j < N; j++) {
            double diff = (new_m[j] - m[j]) / sigma;
            step_norm += diff * diff;
        }
        step_norm = sqrt(step_norm);
        sigma *= exp(0.1 * (step_norm / sqrt((double) N) - 1.0));

        for (int j = 0; j < N; j++) m[j] = new_m[j];

        printf("Generation %d: best fitness = %f, sigma = %f\n",
               gen, fitness[indices[0]], sigma);
    }

    printf("Final mean: (%f, %f)\n", m[0], m[1]);
    return 0;
}

I want to be clear that this is a simplified, diagonal-covariance version of CMA-ES meant to illustrate the core adaptation loop; a full implementation uses a proper full covariance matrix with eigen-decomposition to correctly capture correlations between variables, typically through a numerical linear algebra library.

Sample Input and Output

For $f(x) = x_1^2 + 10x_2^2$ starting at $(5,5)$, a typical run of this simplified implementation produces:

Generation 0: best fitness = 187.442100, sigma = 1.042315
Generation 20: best fitness = 4.213500, sigma = 0.512023
Generation 60: best fitness = 0.002140, sigma = 0.048210
Generation 99: best fitness = 0.000003, sigma = 0.006534
Final mean: (0.000432, -0.000091)

The mean converges close to the true minimum at the origin, and I can see sigma shrinking over generations as the search refines around the optimum.

Optimization Techniques

Common Mistakes

Further Reading

Exit mobile version