Introduction

Imagine you are trying to find the highest point on a foggy hillside. You cannot see the slope — you can only ask "is this spot higher than the last one?" Standard optimization algorithms demand that you know the gradient, the direction of steepest ascent, at every step. But what if the function has no gradient, or is too noisy, or is a black box you cannot peek inside?

CMA-ES — the Covariance Matrix Adaptation Evolution Strategy — was designed for exactly this situation. Introduced by Nikolaus Hansen and Andreas Ostermeier in 1996 and refined into its modern form by 2001, it treats optimization as evolution: a population of candidate solutions is sampled from a multivariate Gaussian distribution, the better solutions are selected, and the distribution itself is updated to concentrate probability mass near the good region.

The key insight is what gets updated. Naïve evolution strategies only move the mean. CMA-ES also updates the full covariance matrix — the parameter that controls the shape, size and orientation of the search cloud. If the landscape stretches along a diagonal, the ellipse tilts to match it. If the optimum is near, the ellipse shrinks. The algorithm learns its own geometry from the history of successful steps, without ever computing a derivative.

That makes CMA-ES a member of a family of methods called non-convex optimization strategies, but one that can handle landscapes far wilder than gradient descent can manage.

Try It

The canvas below shows a 2-D fitness landscape (brighter = higher fitness). A Gaussian search cloud — drawn as an ellipse — samples candidate points, keeps the best, and updates its mean and covariance matrix each generation.

<!-- CMA-ES demo — evolving search ellipse -->
<div class="controls">
  <button id="btn-step" type="button">Step</button>
  <button id="btn-run" type="button">Run</button>
  <button id="btn-reset" type="button" class="ghost">Reset</button>
</div>
<canvas id="canvas" width="360" height="270"></canvas>
<div class="info-row">
  <span id="gen-stat">Generation: 0</span>
  <span class="sep">|</span>
  <span id="sigma-stat">Step size σ: —</span>
</div>
<div id="status" class="status"></div>

/* Layout and control styles */
* { box-sizing: border-box; margin: 0; padding: 0; }
body { font-family: system-ui, sans-serif; color: #222; background: #fff; }
.controls { display: flex; gap: .5rem; flex-wrap: wrap; margin-bottom: .5rem; }
button { font: 600 13px system-ui; padding: .4rem .85rem; border: 1px solid #1d3557;
         background: #1d3557; color: #fff; border-radius: 7px; cursor: pointer; }
button.ghost { background: #fff; color: #1d3557; }
button:disabled { opacity: .4; cursor: default; }
canvas { display: block; border: 1px solid #d0d7df; border-radius: 8px; width: 100%; }
.info-row { font-size: .85rem; color: #444; margin-top: .35rem; display: flex;
            gap: .35rem; align-items: center; }
.sep { color: #aaa; }
.status { font-size: .9rem; font-weight: 600; min-height: 1.3em; margin-top: .3rem; color: #0a7d33; }

const i18nFmt = (s, v) => s.replace(/%(\w+)%/g, (_, k) => String(v[k] ?? ''));
/* CMA-ES core — adapts mean, covariance matrix and step size each generation */

// Canvas setup
const canvas = document.getElementById('canvas');
const ctx = canvas.getContext('2d');
const W = canvas.width, H = canvas.height;

// Diagonal-ridge fitness landscape (brighter = higher fitness)
const PEAK = { x: 0.65, y: 0.4 };
function fitness(x, y) {
  // Rotated Gaussian ridge — peak at (0.65, 0.4), stretched along 30-degree axis
  const dx = x - PEAK.x, dy = y - PEAK.y;
  const rot30 = 0.866 * dx + 0.5 * dy;
  const perp30 = -0.5 * dx + 0.866 * dy;
  return Math.exp(-(rot30 * rot30) / 0.03 - (perp30 * perp30) / 0.12);
}

// CMA-ES state variables
let m = [0.2, 0.75];
let C = [[1, 0], [0, 1]];
let sigma = 0.25;
let gen = 0;
let running = false;
let rafId = null;
const LAMBDA = 12;
const MU = 6;

// Log-linear recombination weights for mu elite samples
const rawW = Array.from({ length: MU }, (_, i) => Math.log(MU + 0.5) - Math.log(i + 1));
const sumW = rawW.reduce((a, b) => a + b, 0);
const weights = rawW.map(w => w / sumW);
const muEff = 1 / weights.reduce((a, w) => a + w * w, 0);

// Evolution paths for sigma and C adaptation
let ps = [0, 0];
let pc = [0, 0];

// CMA-ES learning-rate constants (standard formulas, n=2)
const n = 2;
const csig = (muEff + 2) / (n + muEff + 5);
const dsig = 1 + 2 * Math.max(0, Math.sqrt((muEff - 1) / (n + 1)) - 1) + csig;
const cc = (4 + muEff / n) / (n + 4 + 2 * muEff / n);
const c1 = 2 / ((n + 1.3) * (n + 1.3) + muEff);
const cmu = Math.min(1 - c1, 2 * (muEff - 2 + 1 / muEff) / ((n + 2) * (n + 2) + muEff));
const chiN = Math.sqrt(n) * (1 - 1 / (4 * n) + 1 / (21 * n * n));

function matAdd(A, B) { return [[A[0][0]+B[0][0], A[0][1]+B[0][1]], [A[1][0]+B[1][0], A[1][1]+B[1][1]]]; }
function matScale(A, s) { return [[A[0][0]*s, A[0][1]*s], [A[1][0]*s, A[1][1]*s]]; }
function outerProduct(v) { return [[v[0]*v[0], v[0]*v[1]], [v[1]*v[0], v[1]*v[1]]]; }

// Cholesky decomposition of 2x2 positive-definite matrix
function cholesky(C) {
  const a = C[0][0], b = C[0][1], d = C[1][1];
  const L00 = Math.sqrt(Math.max(a, 1e-12));
  const L10 = b / L00;
  const L11 = Math.sqrt(Math.max(d - L10 * L10, 1e-12));
  return [[L00, 0], [L10, L11]];
}

function sampleNormal() { // Box-Muller transform to sample from N(0,1)
  const u = 1 - Math.random(), v = Math.random();
  return Math.sqrt(-2 * Math.log(u)) * Math.cos(2 * Math.PI * v);
}

function cmaStep() {
  // Sample lambda points from the current Gaussian
  const L = cholesky(C);
  const samples = [];
  for (let k = 0; k < LAMBDA; k++) {
    const z0 = sampleNormal(), z1 = sampleNormal();
    const d0 = sigma * (L[0][0] * z0 + L[0][1] * z1);
    const d1 = sigma * (L[1][0] * z0 + L[1][1] * z1);
    const x = m[0] + d0, y = m[1] + d1;
    samples.push({ x, y, f: fitness(x, y), d: [d0, d1] });
  }
  samples.sort((a, b) => b.f - a.f);
  const elite = samples.slice(0, MU);

  // Weighted mean of elite samples (new distribution center)
  const newM = [0, 0];
  elite.forEach((s, i) => { newM[0] += weights[i] * s.x; newM[1] += weights[i] * s.y; });
  const step = [(newM[0] - m[0]) / sigma, (newM[1] - m[1]) / sigma];

  // Approximate C^(-1/2) step for sigma path (use identity approximation here)
  const det = C[0][0] * C[1][1] - C[0][1] * C[1][0];
  const Cinv = matScale([[C[1][1], -C[0][1]], [-C[1][0], C[0][0]]], 1 / Math.max(det, 1e-12));
  const tmp = [Cinv[0][0]*step[0]+Cinv[0][1]*step[1], Cinv[1][0]*step[0]+Cinv[1][1]*step[1]];
  const Cinv_half_step = [tmp[0], tmp[1]]; // simplified; full CMA uses C^(-1/2)

  // Update sigma path and covariance path
  const sqrtMuEff = Math.sqrt(muEff);
  ps = [
    (1 - csig) * ps[0] + Math.sqrt(csig * (2 - csig)) * sqrtMuEff * step[0],
    (1 - csig) * ps[1] + Math.sqrt(csig * (2 - csig)) * sqrtMuEff * step[1]
  ];
  const hSig = Math.sqrt(ps[0]*ps[0]+ps[1]*ps[1]) / (chiN * Math.sqrt(1-(1-csig)**(2*(gen+1)))) < 1.4+2/(n+1) ? 1 : 0;
  pc = [
    (1 - cc) * pc[0] + hSig * Math.sqrt(cc * (2 - cc)) * sqrtMuEff * step[0],
    (1 - cc) * pc[1] + hSig * Math.sqrt(cc * (2 - cc)) * sqrtMuEff * step[1]
  ];

  // Rank-1 + rank-mu covariance update
  let newC = matScale(C, 1 - c1 - cmu);
  newC = matAdd(newC, matScale(outerProduct(pc), c1));
  let rankMu = [[0,0],[0,0]];
  elite.forEach((s, i) => {
    const dy = [(s.x - m[0])/sigma, (s.y - m[1])/sigma];
    rankMu = matAdd(rankMu, matScale(outerProduct(dy), weights[i]));
  });
  newC = matAdd(newC, matScale(rankMu, cmu));

  // Cumulative step-size adaptation
  const psNorm = Math.sqrt(ps[0]*ps[0]+ps[1]*ps[1]);
  sigma = sigma * Math.exp((csig / dsig) * (psNorm / chiN - 1));
  sigma = Math.max(0.005, Math.min(sigma, 1.5));

  m = newM;
  C = newC;
  gen++;
}

// Draw the fitness landscape as a color image
function drawLandscape() {
  const id = ctx.createImageData(W, H);
  for (let py = 0; py < H; py++) {
    for (let px = 0; px < W; px++) {
      const fx = px / W, fy = py / H;
      const f = fitness(fx, fy);
      const v = Math.round(255 * (0.2 + 0.8 * f));
      const idx = (py * W + px) * 4;
      id.data[idx] = Math.round(v * 0.55);
      id.data[idx+1] = Math.round(v * 0.72);
      id.data[idx+2] = v;
      id.data[idx+3] = 255;
    }
  }
  ctx.putImageData(id, 0, 0);
}

function drawEllipse(cx, cy, C, sigma, color, alpha) {
  // Draw the search ellipse via eigenvalue decomposition
  const px = cx * W, py = cy * H;
  const sc = sigma * Math.min(W, H);
  const a = C[0][0], b = C[0][1], d = C[1][1];
  const tr = a + d, det = a * d - b * b;
  const disc = Math.sqrt(Math.max(0, tr * tr / 4 - det));
  const l1 = Math.sqrt(Math.max(0, tr / 2 + disc));
  const l2 = Math.sqrt(Math.max(0, tr / 2 - disc));
  const angle = Math.atan2(l1 * l1 - a, b) || 0;
  ctx.save();
  ctx.translate(px, py);
  ctx.rotate(angle);
  ctx.beginPath();
  ctx.ellipse(0, 0, l1 * sc * 1.5, l2 * sc * 1.5, 0, 0, 2 * Math.PI);
  ctx.globalAlpha = alpha;
  ctx.strokeStyle = color;
  ctx.lineWidth = 2;
  ctx.stroke();
  ctx.globalAlpha = 1;
  ctx.restore();
}

function drawSamples(samples) {
  samples.forEach((s, i) => {
    const px = s.x * W, py = s.y * H;
    ctx.beginPath();
    ctx.arc(px, py, i < MU ? 4 : 2.5, 0, 2 * Math.PI);
    ctx.fillStyle = i < MU ? '#e63946' : 'rgba(200,200,200,0.7)';
    ctx.fill();
  });
}

function drawMean(mx, my) {
  const px = mx * W, py = my * H;
  ctx.beginPath();
  ctx.arc(px, py, 5, 0, 2 * Math.PI);
  ctx.fillStyle = '#1d3557';
  ctx.fill();
}

let lastSamples = [];

function render() {
  drawLandscape();
  drawEllipse(m[0], m[1], C, sigma, '#ff9f1c', 0.9);
  if (lastSamples.length) drawSamples(lastSamples);
  drawMean(m[0], m[1]);
  document.getElementById('gen-stat').textContent = i18nFmt('Generation: %gen%', { gen });
  document.getElementById('sigma-stat').textContent = i18nFmt('Step size σ: %sigma%', { sigma: sigma.toFixed(4) });
}

function stepOnce() {
  // Run one CMA-ES generation and redraw
  const L = cholesky(C);
  const samples = [];
  for (let k = 0; k < LAMBDA; k++) {
    const z0 = sampleNormal(), z1 = sampleNormal();
    const d0 = sigma * (L[0][0] * z0 + L[0][1] * z1);
    const d1 = sigma * (L[1][0] * z0 + L[1][1] * z1);
    const x = m[0] + d0, y = m[1] + d1;
    samples.push({ x, y, f: fitness(x, y), d: [d0, d1] });
  }
  samples.sort((a, b) => b.f - a.f);
  lastSamples = samples;
  cmaStep();
  render();
  const peakFit = fitness(m[0], m[1]);
  if (peakFit > 0.98) {
    document.getElementById('status').textContent = i18nFmt('Converged at generation %gen% — ellipse locked onto the peak!', { gen });
  }
}

function reset() {
  if (rafId) { cancelAnimationFrame(rafId); rafId = null; }
  running = false;
  m = [0.2, 0.75];
  C = [[1, 0], [0, 1]];
  sigma = 0.25;
  gen = 0;
  ps = [0, 0];
  pc = [0, 0];
  lastSamples = [];
  document.getElementById('status').textContent = '';
  document.getElementById('btn-run').textContent = 'Run';
  document.getElementById('gen-stat').textContent = 'Generation: 0';
  document.getElementById('sigma-stat').textContent = 'Step size σ: —';
  render();
}

function runLoop() {
  stepOnce();
  if (running && gen < 200 && fitness(m[0], m[1]) < 0.999) {
    rafId = requestAnimationFrame(runLoop);
  } else {
    running = false;
    document.getElementById('btn-run').textContent = 'Run';
  }
}

document.getElementById('btn-step').addEventListener('click', () => {
  if (!running) stepOnce();
});
document.getElementById('btn-run').addEventListener('click', () => {
  running = !running;
  document.getElementById('btn-run').textContent = running ? 'Pause' : 'Run';
  if (running) runLoop();
  else if (rafId) { cancelAnimationFrame(rafId); rafId = null; }
});
document.getElementById('btn-reset').addEventListener('click', reset);

render();
Result
Console

Press Step to advance one generation at a time and observe how the ellipse tilts, stretches and shrinks toward the peak. Press Run to let it converge automatically. Notice that on a diagonal ridge the ellipse rotates to align with it — that rotation is the covariance matrix at work.

The Real Complexity

At each generation, CMA-ES maintains three core state variables: the mean m\mathbf{m}, the covariance matrix C\mathbf{C}, and the step size σ\sigma. Updating them is where the magic — and the cost — live.

Sampling. λ\lambda new candidates are drawn as xk=m+σ C1/2 zk\mathbf{x}_k = \mathbf{m} + \sigma\,\mathbf{C}^{1/2}\,\mathbf{z}_k where each zk∼N(0,I)\mathbf{z}_k \sim \mathcal{N}(\mathbf{0}, \mathbf{I}). Computing C1/2\mathbf{C}^{1/2} via eigendecomposition costs O(n3)O(n^3) but only needs refreshing every O(n)O(n) generations.

Mean update. The new mean is the weighted average of the μ\mu best candidates (the elite half):

m′=∑i=1μwi xi:λ\mathbf{m}' = \sum_{i=1}^{\mu} w_i\, \mathbf{x}_{i:\lambda}

Covariance update. This is the heart of CMA-ES. Two complementary mechanisms adapt C\mathbf{C}:

  • Rank-μ\mu update — accumulates the outer products (xi:λ−m)(xi:λ−m)⊤(\mathbf{x}_{i:\lambda} - \mathbf{m})(\mathbf{x}_{i:\lambda} - \mathbf{m})^{\top} from all elite samples, giving a low-rank estimate of the local curvature.
  • Rank-1 update — maintains an evolution path pc\mathbf{p}_c, a running sum of consecutive mean shifts that captures correlations across generations, and adds pc pc⊤\mathbf{p}_c\,\mathbf{p}_c^{\top} to C\mathbf{C}.

Together they give the update C′=(1−c1−cμ)C+c1 pcpc⊤+cμ∑iwi(xi:λ−m)(xi:λ−m)⊤\mathbf{C}' = (1-c_1-c_\mu)\mathbf{C} + c_1\,\mathbf{p}_c\mathbf{p}_c^{\top} + c_\mu\sum_i w_i(\mathbf{x}_{i:\lambda}-\mathbf{m})(\mathbf{x}_{i:\lambda}-\mathbf{m})^{\top}.

Step-size control. A separate evolution path pσ\mathbf{p}_\sigma tracks whether consecutive steps are correlated; if they are, σ\sigma increases (steps are too small); if they zigzag, σ\sigma shrinks. This cumulative step-size adaptation keeps the algorithm in a productive regime.

Computational cost. The full algorithm is O(n2)O(n^2) per function evaluation in nn dimensions once the eigendecomposition is amortized — expensive compared to gradient descent's O(n)O(n), but gradient-free methods that only scale the axes (diagonal covariance) fail on rotated or ill-conditioned problems. CMA-ES trades cost for robustness.

The algorithm is not NP-hard to run; the covariance update is polynomial. What is hard is the underlying black-box optimization problem itself: no algorithm can find the global optimum of an arbitrary function in polynomial time with only function-value queries.

Where It Matters

CMA-ES shines wherever gradients are absent, noisy or misleading:

  • Neuroevolution and reinforcement learning: OpenAI's 2017 paper showed CMA-ES and related evolution strategies matching deep-RL performance on Atari and MuJoCo without backpropagation. The policy network weights are the search space.
  • Robot locomotion: physical simulations are noisy and sometimes non-differentiable (contact forces, joint limits). CMA-ES tunes gait parameters directly from episode returns.
  • Aerodynamic shape optimization: computational fluid dynamics solvers are expensive black boxes. CMA-ES finds wing profiles that reduce drag without needing adjoint solvers.
  • Hyperparameter search: neural network training is sensitive to learning rate, batch size, architecture choices — CMA-ES explores this space more efficiently than random search when the number of hyperparameters is moderate (n≲100n \lesssim 100).
  • Protein and drug design: scoring functions for molecular fitness are expensive, non-differentiable and full of local optima. CMA-ES navigates them with far fewer evaluations than grid search.

CMA-ES is considered one of the best general-purpose optimizers for problems with nn up to a few hundred dimensions and expensive-to-evaluate objective functions. For very high-dimensional problems (n≫1000n \gg 1000), variants like sep-CMA-ES (diagonal covariance) or VD-CMA trade some power for scalability. Compare this landscape with non-convex optimization more broadly, or with discrete methods like simulated annealing which tackle similar rugged landscapes.

Conclusion

CMA-ES does something quietly remarkable: it treats the search distribution itself as the thing to be learned. By updating a full covariance matrix from the history of successful steps, it builds an internal model of the landscape's geometry — its ridges, valleys and correlations — without ever computing a gradient.

The price is quadratic cost in dimension, but the payoff is an optimizer that degrades gracefully on the hardest real-world problems: noisy, non-differentiable, multi-modal black boxes that defeat gradient-based methods entirely. Three decades after its introduction, CMA-ES remains one of the most reliable tools in the optimization toolkit.

So the next time you see a robot learning to walk or a wing shape being tuned by simulation, there is a good chance a search ellipse is somewhere in the loop — rotating, shrinking, chasing the peak.

Share this article

Pick a channel — or use your device's native share sheet.

Comments

Loading comments...

https://www.kipuhub.com/en/article/cma-es/Content licensed under CC BY-NC 4.0.