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
- Search distribution: A multivariate normal (Gaussian) distribution $\mathcal{N}(m, \sigma^2 C)$ that candidate solutions are sampled from at each generation.
- Mean vector (m): The current estimate of the optimum, updated each generation as a weighted average of the best-performing sampled points.
- Step size (σ): A global scalar controlling the overall scale of the search distribution.
- Covariance matrix (C): A matrix capturing the shape, scale, and orientation of the search distribution, learned over generations to reflect the local structure of the objective function.
- Evolution paths (p_σ, p_c): Accumulated history of successive steps, used to detect correlations between generations and to speed up adaptation.
- Recombination weights: Weights assigned to the best $\mu$ individuals in a generation, used to compute the weighted mean update, giving more influence to better-performing points.
How It Works
I run CMA-ES through the following cycle at each generation:
- I sample $\lambda$ candidate solutions from the current multivariate normal distribution $\mathcal{N}(m, \sigma^2 C)$.
- I evaluate the objective function for each candidate.
- I rank the candidates by fitness and select the best $\mu$ of them.
- I update the mean $m$ as a weighted average of these $\mu$ best candidates.
- I update the evolution paths $p_\sigma$ and $p_c$, which track the recent history of successful search directions.
- 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).
- 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).
- 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$.
- I sample $\lambda=6$ candidates around $(5,5)$ using the identity covariance, so initially the spread is circular.
- Evaluating fitness, I notice candidates that moved more in the $x_1$ direction tend to improve fitness more per unit distance than those moving in $x_2$, because the $10x_2^2$ term penalizes $x_2$ movement much more heavily.
- I select the best $\mu=3$ candidates and compute the new mean, which shifts noticeably in the direction of reduced fitness — likely more toward smaller $x_1$ initially, since that axis is “cheaper” to explore.
- As generations proceed, the covariance matrix $C$ gradually becomes anisotropic: it shrinks along the $x_2$ direction (since large steps there are usually penalized) and can remain relatively larger along $x_1$, reshaping the sampling ellipse to match the elongated bowl shape of the objective function.
- This adaptation lets CMA-ES converge efficiently toward $(0,0)$ even though the problem is poorly scaled between dimensions — something a fixed, isotropic-step algorithm would struggle with, since it would need very small steps to avoid overshooting in $x_2$ while making frustratingly slow progress in $x_1$.
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
- It adapts to correlations and different scalings between variables automatically, handling ill-conditioned, non-separable problems well.
- It has strong theoretical convergence properties, including affine invariance, and quasi-parameter-free operation (very few settings need manual tuning).
- It performs excellently on the standard black-box optimization benchmarks (like COCO/BBOB), often outperforming other derivative-free methods on medium-dimensional problems.
- The step-size adaptation mechanism intelligently speeds up or slows down the search based on recent progress patterns.
- It requires no gradient information and works well on noisy objective functions.
Disadvantages
- Its per-generation computational cost grows cubically with dimensionality due to the covariance matrix operations, limiting practicality for very high-dimensional problems (thousands of variables) without special large-scale variants.
- It requires storing and manipulating an $n \times n$ covariance matrix, which becomes memory-intensive at high dimensions.
- It’s more complex to implement correctly than simpler evolutionary algorithms, given the intricate update equations for paths, step size, and covariance.
- Like all stochastic optimizers, it can converge to a local rather than global optimum on highly multimodal problems without restart strategies.
- It’s designed for continuous domains and needs adaptation or a different algorithm entirely for discrete or combinatorial problems.
Applications
- Robotics, particularly for optimizing control policies and reinforcement learning agent parameters in continuous action spaces.
- Hyperparameter optimization for machine learning models, especially when the parameter space has strong interactions.
- Aerodynamic and structural engineering design optimization.
- Direct policy search in reinforcement learning, where CMA-ES is used as a strong baseline or component in black-box policy optimization.
- Computer graphics and animation, tuning physical simulation or character control parameters.
- Any black-box simulation-based optimization problem with moderate dimensionality and no gradient access.
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
- Restart strategies (IPOP-CMA-ES, BIPOP-CMA-ES): Restarting with increasing (or alternating small/large) population sizes when the search stagnates helps escape local optima on multimodal problems.
- Active covariance matrix adaptation: Incorporating information from the worst-performing individuals (not just the best) to actively shrink the covariance in unpromising directions, speeding convergence further.
- Efficient covariance updates: Using rank-one and rank-$\mu$ updates together, and only performing full eigen-decomposition periodically rather than every generation, balances accuracy with computational cost.
- Boundary handling: For bounded problems, techniques like resampling out-of-bound candidates or penalizing them preserve the integrity of the covariance adaptation.
- Large-scale variants: For high-dimensional problems, using limited-memory or separable approximations of CMA-ES (like sep-CMA-ES or VD-CMA) avoids the full $O(n^2)$ to $O(n^3)$ costs.
Common Mistakes
- Using an inappropriately small initial step size $\sigma$, which can cause the algorithm to converge prematurely to a poor local region before it has explored enough.
- Neglecting to periodically re-decompose or numerically stabilize the covariance matrix, leading to numerical issues (loss of positive-definiteness) over many generations.
- Choosing a population size $\lambda$ that’s too small for the problem’s dimensionality, which starves the covariance adaptation of enough information to learn useful correlations.
- Ignoring recommended default parameter formulas (for $c_1$, $c_\mu$, $c_\sigma$, etc.) and guessing arbitrary values, which can destabilize the adaptation dynamics.
- Applying vanilla CMA-ES directly to very high-dimensional or highly multimodal problems without considering restarts or large-scale variants, leading to poor scalability or premature convergence.
Further Reading
- Hansen, N., Ostermeier, A. “Completely Derandomized Self-Adaptation in Evolution Strategies.” Evolutionary Computation, 2001: https://direct.mit.edu/evco/article/9/2/159/892/Completely-Derandomized-Self-Adaptation-in
- Hansen, N. “The CMA Evolution Strategy: A Tutorial.” 2016 (updated): https://arxiv.org/abs/1604.00772
- Hansen, N., Auger, A. “CMA-ES: Evolution Strategies and Covariance Matrix Adaptation.” GECCO Companion, 2011: https://dl.acm.org/doi/10.1145/2001858.2002123
- COCO/BBOB benchmarking platform: https://numbbo.github.io/coco/
- Wikipedia overview: https://en.wikipedia.org/wiki/CMA-ES
