► To cite this version:

Nikolaus Hansen. CMA-ES. Parallel Problem Solving from Nature - PPSN XVIII, 18th International Conference, Sep 2024, Hagenberg, Austria. hal-04709819

HAL Id: hal-04709819

https://inria.hal.science/hal-04709819v1

Submitted on 25 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
CMA-ES

Covariance Matrix Adaptation Evolution Strategy

Nikolaus Hansen

Inria & Ecole Polytechnique, France

September 2024
I am happy to answer questions at any time!

Need the slides? Available under Talks at my homepage:
http://www.cmap.polytechnique.fr/~nikolaus.hansen/talks.html
# Topics

## 0. Problem statement

### 1. What makes an optimization problem difficult to solve?

### 2. How does the CMA-ES work?

* Normal Distribution, Rank-Based Recombination

* Step-Size Adaptation

* Covariance Matrix Adaptation

### 3. What can users do?

* Choice of problem formulation and encoding (not covered)

* Choice of initial solution and initial step-size

* Restarts, Increasing Population Size

* Restricted Covariance Matrix
# Problem Statement

## Continuous Domain Search/Optimization

* Task: **minimize** an **objective function** (fitness function, loss function) in continuous domain

$$f: \mathcal{X} \subseteq \mathbb{R}^n \to \mathbb{R}, \quad \mathbf{x} \mapsto f(\mathbf{x})$$

* **Black Box** scenario (direct search scenario)

<figure><img src="image_2.png" /></figure>

* gradients are not available or not useful

* problem domain specific knowledge is used only within the black box, e.g. within an appropriate encoding

* **Search costs**: number of function evaluations

<figure><img src="image_3.png" /></figure>
# Problem Statement

## Continuous Domain Search/Optimization

* Goal

* fast convergence to the global optimum

... or to a robust solution x
* solution x with **small function value** $f(x)$ with **least search cost**
there are two conflicting objectives

* Typical Examples

* shape optimization (e.g. using CFD)

* model calibration

* parameter calibration

curve fitting, airfoils
biological, physical
controller, plants, images

* Difficulties

* exhaustive search is infeasible

* naive random search takes too long

* deterministic search is not successful / takes too long

**Approach:** stochastic search, Evolutionary Algorithms

<figure><img src="image_4.png" /></figure>
# What Makes a Function Difficult to Solve?

## Why stochastic search?

* non-linear, non-quadratic, non-convex
  on linear and quadratic functions much better
  search policies are available

<figure><img src="image_5.png" /></figure>

* ruggedness
  non-smooth, discontinuous, multimodal, and/or
  noisy function

<figure><img src="image_6.png" /></figure>

* dimensionality (size of search space)
  (considerably) larger than three

<figure><img src="image_7.png" /></figure>

* non-separability
  dependencies between the objective variables

* ill-conditioning

* non-smooth level sets

<figure><img src="image_8.png" /><figcaption>gradient direction Newton direction</figcaption></figure>
# Ruggedness

non-smooth, discontinuous, multimodal, and/or noisy

<figure><img src="image_9.png" /></figure>

cut from a 5-D example, (easily) solvable with evolution strategies

<figure><img src="image_10.png" /></figure>
# Curse of Dimensionality

The term *Curse of dimensionality* (Richard Bellman) refers to effects caused by the **rapid increase in volume** associated with adding extra dimensions to a (mathematical) space.

**Example:** Consider placing 20 points equally spaced onto the interval [0, 1]. Now consider the 10-dimensional space [0, 1]¹⁰. To get **similar coverage** in terms of distance between adjacent points requires 20¹⁰ ≈ 10¹³ points. 20 points appear now as isolated points in a vast empty space.

**Remark:** **distance measures** break down in higher dimensionalities (the central limit theorem kicks in).

**Consequence:** a **search policy** that is valuable in small dimensions **might be useless** in moderate or large dimensional search spaces.

**Example:** exhaustive search.

<figure><img src="image_11.png" /></figure>
# Separable Problems

## Definition (Separable Problem)

A function $f$ is separable if

$$ \arg\min_{(x_1, \dots, x_n)} f(x_1, \dots, x_n) = \left( \arg\min_{x_1} f(x_1, \dots), \dots, \arg\min_{x_n} f(\dots, x_n) \right) $$

⇒ it follows that $f$ can be optimized in a sequence of $n$ independent 1-D optimization processes

## Example: Additively decomposable functions

$$ f(x_1, \dots, x_n) = \sum_{i=1}^{n} f_i(x_i) $$

Rastrigin function

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>
# Non-Separable Problems

## Building a non-separable problem from a separable one (1,2)

### Rotating the coordinate system

* $f : \mathbf{x} \mapsto f(\mathbf{x})$ separable

* $f : \mathbf{x} \mapsto f(\mathbf{Rx})$ non-separable

<figure><img src="image_14.png" /><figcaption>R rotation matrix</figcaption></figure>

<figure><img src="image_15.png" /></figure>

<footnote n="1">¹Hansen, Ostermeier, Gawelczyk (1995). On the adaptation of arbitrary normal mutation distributions in evolution strategies: The generating set adaptation. Sixth ICGA, pp. 57-64, Morgan Kaufmann</footnote>

<footnote n="2">²Salomon (1996). "Reevaluating Genetic Algorithm Performance under Coordinate Rotation of Benchmark Functions; A survey of some theoretical and practical aspects of genetic algorithms." BioSystems, 39(3):263-278</footnote>
# Ill-Conditioned Problems

## Curvature of level sets

Consider the convex-quadratic function

$$f(\mathbf{x}) = \frac{1}{2}(\mathbf{x} - \mathbf{x}^*)^T \mathbf{H} (\mathbf{x} - \mathbf{x}^*) = \frac{1}{2} \sum_{i} h_{i,i} (x_i - x_i^*)^2 + \frac{1}{2} \sum_{i \neq j} h_{i,j} (x_i - x_j^*) (x_j - x_j^*)$$

$\mathbf{H}$ is Hessian matrix of $f$ and symmetric positive definite

<figure><img src="image_16.png" /></figure>

gradient direction $-f'(\mathbf{x})^T$

Newton direction $-\mathbf{H}^{-1}f'(\mathbf{x})^T$

Ill-conditioning means **squeezed level sets** (high curvature). Condition number equals nine here. Condition numbers up to $10^{10}$ are not unusual in real world problems.

If $\mathbf{H} \approx \mathbf{I}$ (small condition number of $\mathbf{H}$) first order information (e.g. the gradient) is sufficient. Otherwise **second order information** (estimation of $\mathbf{H}^{-1}$) **is necessary**.

<figure><img src="image_17.png" /></figure>
# Non-smooth level sets (sharp ridges)

Similar difficulty **but worse** than ill-conditioning

<figure><img src="image_18.png" /><figcaption>1-norm</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>scaled 1-norm</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>1/2-norm</figcaption></figure>

opening angle is the crucial parameter
# What Makes a Function Difficult to Solve?

## ...and what can be done

<table><thead><tr><td>The Problem</td><td>Possible Approaches</td></tr></thead><tbody><tr><td>Dimensionality</td><td>exploiting the problem structure<br>separability, locality/neighborhood, encoding</td></tr><tr><td>Ill-conditioning</td><td>second order approach<br>changes the neighborhood metric</td></tr><tr><td>Ruggedness and<br>non-smooth level<br>sets</td><td>non-local policy, large sampling width (step-size)<br>as large as possible while preserving a<br>reasonable convergence speed<br><br>population-based method, stochastic, non-elitistic<br>recombination operator<br>serves as repair mechanism<br><br>restarts</td></tr></tbody></table>

...metaphors

<figure><img src="image_21.png" /></figure>
# Topics

1. What makes the problem difficult to solve?

2. How does the CMA-ES work?

* Normal Distribution, Rank-Based Recombination

* Step-Size Adaptation

* Covariance Matrix Adaptation

3. What can/should the users do for the CMA-ES to work effectively on their problem?

* Choice of problem formulation and encoding (not covered)

* Choice of initial solution and initial step-size

* Restarts, Increasing Population Size

* Restricted Covariance Matrix
# Stochastic Search

A black box search template to minimize $f: \mathbb{R}^n \to \mathbb{R}$

Initialize distribution parameters $\theta$, set population size $\lambda \in \mathbb{N}$
While not terminate

1. Sample distribution $P(x|\theta) \to x_1, \dots, x_\lambda \in \mathbb{R}^n$

2. Evaluate $x_1, \dots, x_\lambda$ on $f$

3. Update parameters $\theta \leftarrow F_{\theta}(\theta, x_1, \dots, x_\lambda, f(x_1), \dots, f(x_\lambda))$

Everything depends on the definition of $P$ and $F_\theta$

deterministic algorithms are covered as well

In many Evolutionary Algorithms the distribution $P$ is implicitly defined
via operators on a population, in particular, selection, recombination
and mutation

Natural template for (incremental) Estimation of Distribution Algorithms
# Stochastic Search

A black box search template to minimize $f: \mathbb{R}^n \to \mathbb{R}$

Initialize distribution parameters $\theta$, set population size $\lambda \in \mathbb{N}$

While not terminate

1. Sample distribution $P(x|\theta) \to x_1, \dots, x_\lambda \in \mathbb{R}^n$

2. Evaluate $x_1, \dots, x_\lambda$ on $f$

3. Update parameters $\theta \leftarrow F_{\theta}(\theta, x_1, \dots, x_\lambda, f(x_1), \dots, f(x_\lambda))$

Everything depends on the definition of $P$ and $F_\theta$

deterministic algorithms are covered as well

In many Evolutionary Algorithms the distribution $P$ is implicitly defined via *operators on a population*, in particular, selection, recombination and mutation

Natural template for (incremental) *Estimation of Distribution Algorithms*
# Stochastic Search

A black box search template to minimize $f: \mathbb{R}^n \to \mathbb{R}$

Initialize distribution parameters $\theta$, set population size $\lambda \in \mathbb{N}$

**While not terminate**

1. Sample distribution $P(x|\theta) \to x_1, \dots, x_\lambda \in \mathbb{R}^n$

2. Evaluate $x_1, \dots, x_\lambda$ on $f$

3. Update parameters $\theta \leftarrow F_{\theta}(\theta, x_1, \dots, x_\lambda, f(x_1), \dots, f(x_\lambda))$

Everything depends on the definition of $P$ and $F_\theta$

deterministic algorithms are covered as well

In many Evolutionary Algorithms the distribution $P$ is implicitly defined
via **operators on a population**, in particular, selection, recombination
and mutation

Natural template for (incremental) Estimation of Distribution Algorithms
# Stochastic Search

A black box search template to minimize $f: \mathbb{R}^n \to \mathbb{R}$

Initialize distribution parameters $\theta$, set population size $\lambda \in \mathbb{N}$
While not terminate

1. Sample distribution $P(x|\theta) \to x_1, \dots, x_\lambda \in \mathbb{R}^n$

2. Evaluate $x_1, \dots, x_\lambda$ on $f$

3. Update parameters $\theta \leftarrow F_{\theta}(\theta, x_1, \dots, x_\lambda, f(x_1), \dots, f(x_\lambda))$

Everything depends on the definition of $P$ and $F_\theta$

deterministic algorithms are covered as well

In many Evolutionary Algorithms the distribution $P$ is implicitly defined
via operators on a population, in particular, selection, recombination
and mutation

Natural template for (incremental) Estimation of Distribution Algorithms
# Stochastic Search

A black box search template to minimize $f: \mathbb{R}^n \to \mathbb{R}$

Initialize distribution parameters $\theta$, set population size $\lambda \in \mathbb{N}$
While not terminate

1. Sample distribution $P(x|\theta) \to x_1, \dots, x_\lambda \in \mathbb{R}^n$

2. Evaluate $x_1, \dots, x_\lambda$ on $f$

3. Update parameters $\theta \leftarrow F_{\theta}(\theta, x_1, \dots, x_\lambda, f(x_1), \dots, f(x_\lambda))$

Everything depends on the definition of $P$ and $F_\theta$

deterministic algorithms are covered as well

In many Evolutionary Algorithms the distribution $P$ is implicitly defined
via operators on a population, in particular, selection, recombination
and mutation

Natural template for (incremental) Estimation of Distribution Algorithms
# Stochastic Search

A black box search template to minimize $f: \mathbb{R}^n \to \mathbb{R}$

Initialize distribution parameters $\theta$, set population size $\lambda \in \mathbb{N}$
While not terminate

1. Sample distribution $P(x|\theta) \to x_1, \dots, x_\lambda \in \mathbb{R}^n$

2. Evaluate $x_1, \dots, x_\lambda$ on $f$

3. Update parameters $\theta \leftarrow F_{\theta}(\theta, x_1, \dots, x_\lambda, f(x_1), \dots, f(x_\lambda))$

Everything depends on the definition of $P$ and $F_\theta$

deterministic algorithms are covered as well

In many Evolutionary Algorithms the distribution $P$ is implicitly defined
via operators on a population, in particular, selection, recombination
and mutation

Natural template for (incremental) Estimation of Distribution Algorithms
# Stochastic Search

A black box search template to minimize $f: \mathbb{R}^n \to \mathbb{R}$

Initialize distribution parameters $\theta$, set population size $\lambda \in \mathbb{N}$
While not terminate

1. Sample distribution $P(x|\theta) \to x_1, \dots, x_\lambda \in \mathbb{R}^n$

2. Evaluate $x_1, \dots, x_\lambda$ on $f$

3. Update parameters $\theta \leftarrow F_{\theta}(\theta, x_1, \dots, x_\lambda, f(x_1), \dots, f(x_\lambda))$

Everything depends on the definition of $P$ and $F_\theta$

deterministic algorithms are covered as well

In many Evolutionary Algorithms the distribution $P$ is implicitly defined
via operators on a population, in particular, selection, recombination
and mutation

Natural template for (incremental) Estimation of Distribution Algorithms
# Stochastic Search

A black box search template to minimize $f: \mathbb{R}^n \to \mathbb{R}$

Initialize distribution parameters $\theta$, set population size $\lambda \in \mathbb{N}$
While not terminate

1. Sample distribution $P(x|\theta) \to x_1, \dots, x_\lambda \in \mathbb{R}^n$

2. Evaluate $x_1, \dots, x_\lambda$ on $f$

3. Update parameters $\theta \leftarrow F_{\theta}(\theta, x_1, \dots, x_\lambda, f(x_1), \dots, f(x_\lambda))$

Everything depends on the definition of $P$ and $F_\theta$

deterministic algorithms are covered as well

In many Evolutionary Algorithms the distribution $P$ is implicitly defined
via **operators on a population**, in particular, selection, recombination
and mutation

Natural template for (incremental) *Estimation of Distribution Algorithms*
# The CMA-ES

**Input:** $m \in \mathbb{R}^n$; $\sigma \in \mathbb{R}_+$; $\lambda \in \mathbb{N}_{\ge 2}$, usually $\lambda \ge 5$, default $4 + \lfloor 3 \log n \rfloor$

**Set** $c_m = 1$; $c_1 \approx 2/n^2$; $c_\mu \approx \mu_w/n^2$; $c_c \approx 4/n$; $c_\sigma \approx 1/\sqrt{n}$; $d_\sigma \approx 1$; $w_{i=1...\lambda}$ decreasing in $i$ and $\sum_i^\mu w_i = 1$, $w_\mu > 0 \ge w_{\mu+1}$, $\mu_w^{-1} := \sum_{i=1}^\mu w_i^2 \approx 3/\lambda$

**Initialize** $C = I$, and $p_c = 0$, $p_\sigma = 0$

**While** not terminate

$x_i = m + \sigma y_i$, where $y_i \sim N_i(0, C)$ for $i = 1, ..., \lambda$ sampling

$m \leftarrow m + c_m \sigma y_w$, where $y_w = \sum_{i=1}^\mu w_i y_{rk^{-1}(i)}$ update mean

$p_\sigma \leftarrow (1 - c_\sigma) p_\sigma + \sqrt{1 - (1 - c_\sigma)^2} \sqrt{\mu_w} C^{-\frac{1}{2}} y_w$ path for $\sigma$

$p_c \leftarrow (1 - c_c) p_c + \mathbf{1}_{[0,2n]}\{\|p_\sigma\|^2\} \sqrt{1 - (1 - c_c)^2} \sqrt{\mu_w} y_w$ path for $C$

$\sigma \leftarrow \sigma \times \exp\left(\frac{c_\sigma}{d_\sigma} \left(\frac{\|p_\sigma\|}{E\|N(0,I)\|} - 1\right)\right)$ update of $\sigma$

$C \leftarrow C + c_\mu \sum_{i=1}^\lambda w_{rk(i)} (y_i y_i^\top - C) + c_1 (p_c p_c^\top - C)$ update $C$

**Not covered:** termination, restarts, useful output, search boundaries and encoding, corrections for: positive definiteness guaranty, $p_c$ variance loss, $c_\sigma$ and $d_\sigma$ for large $\lambda$

<figure><img src="image_22.png" /></figure>
# Stochastic Search

A black box search template to minimize $f: \mathbb{R}^n \to \mathbb{R}$

Initialize distribution parameters $\theta$, set population size $\lambda \in \mathbb{N}$
While not terminate

1. Sample distribution $P(x|\theta) \to x_1, \dots, x_\lambda \in \mathbb{R}^n$

2. Evaluate $x_1, \dots, x_\lambda$ on $f$

3. Update parameters $\theta \leftarrow F_{\theta}(\theta, x_1, \dots, x_\lambda, f(x_1), \dots, f(x_\lambda))$

Everything depends on the definition of $P$ and $F_\theta$

deterministic algorithms are covered as well

In many Evolutionary Algorithms the distribution $P$ is implicitly defined
via operators on a population, in particular, selection, recombination
and mutation

Natural template for (incremental) Estimation of Distribution Algorithms
# Evolution Strategies

New search points are sampled normally distributed

$$x_i \sim m + \sigma N_i(0, C) \quad \text{for } i = 1, \dots, \lambda$$

as perturbations of **m**, where $x_i, m \in \mathbb{R}^n, \sigma \in \mathbb{R}_+, C \in \mathbb{R}^{n \times n}$ where

<figure><img src="image_23.png" /></figure>

* the **mean** vector $m \in \mathbb{R}^n$ represents the favorite solution

* the so-called **step-size** $\sigma \in \mathbb{R}_+$ controls the **step length**

* the **covariance matrix** $C \in \mathbb{R}^{n \times n}$ determines the **shape** of the distribution ellipsoid

here, all new points are sampled with the same parameters

The question remains how to update **m**, **C**, and **\sigma**.

<figure><img src="image_24.png" /></figure>
# Why Normal Distributions?

1. widely observed in nature, for example as phenotypic traits

2. only stable distribution with finite variance

stable means that the sum of normal variates is again normal:

$$N(\mathbf{x}, \mathbf{A}) + N(\mathbf{y}, \mathbf{B}) \sim N(\mathbf{x} + \mathbf{y}, \mathbf{A} + \mathbf{B})$$

helpful in **design and analysis** of algorithms related to the **central limit theorem**

3. most convenient way to generate **isotropic** search points

the isotropic distribution does **not favor any direction**, rotational invariant

4. maximum entropy distribution with finite variance

the least possible assumptions on $f$ in the distribution shape

<figure><img src="image_25.png" /></figure>
# Normal Distribution

<figure><img src="image_26.png" /></figure>

probability density of the 1-D standard normal distribution

<figure><img src="image_27.png" /></figure>

probability density of a 2-D normal distribution

<figure><img src="image_28.png" /></figure>

<figure><img src="image_29.png" /></figure>
# The Multi-Variate (n-Dimensional) Normal Distribution

Any multi-variate normal distribution $N(m, C)$ is uniquely determined by its mean value $m \in \mathbb{R}^n$ and its symmetric positive definite $n \times n$ covariance matrix $C$.

## The mean value $m$

* determines the displacement (translation)

* value with the largest density (modal value)

* the distribution is symmetric about the distribution mean

<figure><img src="image_30.png" /></figure>

## The covariance matrix $C$

* determines the shape

* geometrical interpretation: any covariance matrix can be uniquely identified with the iso-density ellipsoid $\{x \in \mathbb{R}^n | (x - m)^T C^{-1} (x - m) = n\}$

<figure><img src="image_31.png" /></figure>
Effect of Dimensionality

<figure><img src="image_32.png" /></figure>

<figure><img src="image_33.png" /></figure>

<figure><img src="image_34.png" /></figure>

$$ ||\mathcal{N}(\mathbf{0}, \mathbf{I})|| \rightarrow \mathcal{N}(\sqrt{n-1/2}, \frac{1}{2}) \text{ with modal value } \sqrt{n-1} $$

yet: maximum entropy distribution

also consider a difference between two vectors:

$$ ||\mathcal{N}(\mathbf{0}, \mathbf{I}) - \mathcal{N}(\mathbf{0}, \mathbf{I})|| \sim ||\mathcal{N}(\mathbf{0}, \mathbf{I}) + \mathcal{N}(\mathbf{0}, \mathbf{I})|| \sim \sqrt{2} ||\mathcal{N}(\mathbf{0}, \mathbf{I})|| $$
Effect of Dimensionality

<figure><img src="image_35.png" /></figure>

<figure><img src="image_36.png" /></figure>

<figure><img src="image_37.png" /></figure>

$$ ||\mathcal{N}(\mathbf{0}, \mathbf{I})|| \rightarrow \mathcal{N}(\sqrt{n-1/2}, \frac{1}{2}) \text{ with modal value } \sqrt{n-1} $$

yet: maximum entropy distribution

also consider a difference between two vectors:

$$ ||\mathcal{N}(\mathbf{0}, \mathbf{I}) - \mathcal{N}(\mathbf{0}, \mathbf{I})|| \sim ||\mathcal{N}(\mathbf{0}, \mathbf{I}) + \mathcal{N}(\mathbf{0}, \mathbf{I})|| \sim \sqrt{2} ||\mathcal{N}(\mathbf{0}, \mathbf{I})|| $$
# Stochastic Search

A black box search template to minimize $f: \mathbb{R}^n \to \mathbb{R}$

Initialize distribution parameters $\theta$, set population size $\lambda \in \mathbb{N}$
While not terminate

1. Sample distribution $P(x|\theta) \to x_1, \dots, x_\lambda \in \mathbb{R}^n$

2. Evaluate $x_1, \dots, x_\lambda$ on $f$

3. Update parameters $\theta \leftarrow F_{\theta}(\theta, x_1, \dots, x_\lambda, f(x_1), \dots, f(x_\lambda))$

Everything depends on the definition of $P$ and $F_\theta$

deterministic algorithms are covered as well

In many Evolutionary Algorithms the distribution $P$ is implicitly defined
via operators on a population, in particular, selection, recombination
and mutation

Natural template for (incremental) Estimation of Distribution Algorithms
# The $(\mu/\mu, \lambda)$-ES, Update of the Distribution Mean

## Non-elitist selection and intermediate (weighted) recombination

Given the i-th solution point $\mathbf{x}_i = \mathbf{m} + \sigma \underbrace{\mathcal{N}_i(\mathbf{0}, \mathbf{C})}_{=:\mathbf{y}_i} = \mathbf{m} + \sigma \mathbf{y}_i$

Let $\mathbf{x}_{i:\lambda}$ the i-th ranked solution point, such that $f(\mathbf{x}_{1:\lambda}) \le \dots \le f(\mathbf{x}_{\lambda:\lambda})$. The new mean reads

$$ \mathbf{m} \leftarrow \sum_{i=1}^{\mu} w_i \mathbf{x}_{i:\lambda} = \mathbf{m} + c_m \sigma \underbrace{\sum_{i=1}^{\mu} w_i \mathbf{y}_{i:\lambda}}_{=:\mathbf{y}_w} $$

where

$$ w_1 \ge \dots \ge w_\mu > 0, \quad \sum_{i=1}^{\mu} w_i = 1, \quad \frac{1}{\sum_{i=1}^{\mu} w_i^2} =: \mu_w \approx \frac{\lambda}{4} $$

The best $\mu$ points are selected from the new solutions (non-elitistic) and weighted intermediate recombination is applied.
# Invariance Under Monotonically Increasing Functions

## Rank-based algorithms

Update of all parameters uses only the ranks

$$f(x_{1:\lambda}) \le f(x_{2:\lambda}) \le \dots \le f(x_{\lambda:\lambda})$$

<figure><img src="image_38.png" /></figure>

$$g(f(x_{1:\lambda})) \le g(f(x_{2:\lambda})) \le \dots \le g(f(x_{\lambda:\lambda})) \quad \forall g$$

g is strictly monotonically increasing

g preserves ranks

3

Whitley 1989. The GENITOR algorithm and selection pressure: Why rank-based allocation of reproductive trials is best,

<figure><img src="image_39.png" /></figure>
# Basic Invariance in Search Space

* translation invariance

is true for most optimization algorithms

<figure><img src="image_40.png" /></figure>

Identical behavior on $f$ and $f_a$

$$f: \mathbf{x} \mapsto f(\mathbf{x}), \quad \mathbf{x}^{(t=0)} = \mathbf{x}_0$$

$$f_a: \mathbf{x} \mapsto f(\mathbf{x} - \mathbf{a}), \quad \mathbf{x}^{(t=0)} = \mathbf{x}_0 + \mathbf{a}$$

No difference can be observed w.r.t. the argument of $f$

<figure><img src="image_41.png" /></figure>
# Invariance Under Rigid Search Space Transformation

<figure><img src="image_42.png" /></figure>

for example, invariance under search space rotation (separable ⇔ non-separable)
# Invariance Under Rigid Search Space Transformation

$$f = h_{\text{Rast}} \circ R$$

f-level sets in dimension 2

$$f = h \circ R$$

<figure><img src="image_43.png" /></figure>

<figure><img src="image_44.png" /></figure>

for example, invariance under search space rotation (separable ⇔ non-separable)
# Invariance

> The grand aim of all science is to cover the greatest number of empirical facts by logical deduction from the smallest number of hypotheses or axioms.

> — Albert Einstein

* Empirical performance results

* from benchmark functions

* from solved real world problems

are only useful if they do **generalize** to other problems

* **Invariance** is a strong **non-empirical** statement about generalization

* generalizing (identical) performance from a single function to a whole class of functions

consequently, invariance is important for the evaluation of search algorithms

<figure><img src="image_45.png" /></figure>
# Summary

<figure><img src="image_46.png" /></figure>

On 20D Sphere Function: $f(\mathbf{x}) = \sum_{i=1}^{n} x_i^2 = \|\mathbf{x}\|^2$

* ES without adaptation can't approach the optimum $\Rightarrow$ **adaptation required**
# Topics

## 0. Problem statement

## 1. What makes an optimization problem difficult to solve?

## 2. How does the CMA-ES work?

* Normal Distribution, Rank-Based Recombination

* Step-Size Adaptation

* Covariance Matrix Adaptation

## 3. What can users do?

* Choice of problem formulation and encoding (not covered)

* Choice of initial solution and initial step-size

* Restarts, Increasing Population Size

* Restricted Covariance Matrix
# Methods for Step-Size Control

* **1/5-th success rule**$^{ab}$, often applied with “+”-selection

increase step-size if more than 20% of the new solutions are successful, decrease otherwise

* **$\sigma$-self-adaptation**$^c$, applied with “,”-selection

mutation is applied to the step-size and the better, according to the objective function value, is selected

simplified “global” self-adaptation

* **path length control**$^d$ (Cumulative Step-size Adaptation, CSA)$^e$

self-adaptation derandomized and non-localized

f

<footnote>$^a$Rechenberg 1973, *Evolutionsstrategie, Optimierung technischer Systeme nach Prinzipien der biologischen Evolution*, Frommann-Holzboog</footnote>

<footnote>$^b$Schumer and Steiglitz 1968. Adaptive step size random search. *IEEE TAC*</footnote>

<footnote>$^c$Schwefel 1981, *Numerical Optimization of Computer Models*, Wiley</footnote>

<footnote>$^d$Hansen & Ostermeier 2001, *Completely Derandomized Self-Adaptation in Evolution Strategies*, *Evol. Comput.* 9 (2)</footnote>

<footnote>$^e$Ostermeier *et al* 1994, Step-size adaptation based on non-local use of selection information, *PPSN IV*</footnote>

<footnote>$^f$Hansen *et al* 2014, How to assess step-size adaptation mechanisms in randomised search, *PPSN IX*</footnote>
# Evolution Strategies

## Recalling

New search points are sampled normally distributed

$$x_i \sim m + \sigma N_i(0, C) \quad \text{for } i = 1, \dots, \lambda$$

as perturbations of $m$, where $x_i, m \in \mathbb{R}^n, \sigma \in \mathbb{R}_+, C \in \mathbb{R}^{n \times n}$ where

<figure><img src="image_47.png" /></figure>

* the **mean** vector $m \in \mathbb{R}^n$ represents the favorite solution and $m \leftarrow \sum_{i=1}^{\mu} w_i x_{i:\lambda}$

* the so-called **step-size** $\sigma \in \mathbb{R}_+$ controls the **step length**

* the **covariance matrix** $C \in \mathbb{R}^{n \times n}$ determines the **shape** of the distribution ellipsoid

The remaining question is how to update $\sigma$ and $C$.

<figure><img src="image_48.png" /></figure>
# Path Length Control (CSA)

## The Concept of Cumulative Step-Size Adaptation

$$ \boldsymbol{x}_i = \boldsymbol{m} + \sigma \boldsymbol{y}_i $$

$$ \boldsymbol{m} \leftarrow \boldsymbol{m} + \sigma \boldsymbol{y}_w $$

## Measure the length of the evolution path

the pathway of the mean vector **m** in the generation sequence

<figure><img src="image_49.png" /></figure>

loosely speaking steps are

* perpendicular under random selection (in expectation)

* perpendicular in the desired situation (to be most efficient)

<figure><img src="image_50.png" /></figure>
# Path Length Control (CSA)

## The Equations

Initialize $\textcolor{green}{m} \in \mathbb{R}^n, \sigma \in \mathbb{R}_+, \text{evolution path } p_\sigma = \textcolor{blue}{0},$ set $c_\sigma \approx \mu_w/(n + \mu_w), d_\sigma \approx 1$.

$\textcolor{green}{m} \leftarrow \textcolor{green}{m} + \sigma \textcolor{blue}{y}_w \quad \text{where } \textcolor{blue}{y}_w = \sum_{i=1}^{\mu} w_i \textcolor{blue}{y}_{i:\lambda}$ update mean

$p_\sigma \leftarrow (1 - c_\sigma) p_\sigma + \underbrace{\sqrt{1 - (1 - c_\sigma)^2}}_{\text{accounts for } 1-c_\sigma} \sqrt{\mu_w} \textcolor{blue}{y}_w$

$\sigma \leftarrow \sigma \times \underbrace{\exp\left(\frac{c_\sigma}{d_\sigma}\left(\frac{\|p_\sigma\|}{\mathbb{E}\|\mathcal{N}(\textcolor{blue}{0}, \textcolor{blue}{I})\|} - 1\right)\right)}_{>1 \iff \|p_\sigma\| \text{ is greater than its expectation}}$ update step-size

<figure><img src="image_51.png" /></figure>
<figure><img src="image_52.png" /><figcaption>(5/5, 10)-CSA-ES, default parameters</figcaption></figure>

$$f(\boldsymbol{x}) = \sum_{i=1}^{n} x_i^2$$

in [-0.2, 0.8]^n

for n = 30

<figure><img src="image_53.png" /></figure>
# Two-Point Step-Size Adaptation (TPA)

* Sample a pair of symmetric points along the previous mean shift

$$x_{1/2} = m^{(g)} \pm \sigma^{(g)} \frac{\|\mathcal{N}(0, I)\|}{\|m^{(g)} - m^{(g-1)}\|_{C^{(g)}}} (m^{(g)} - m^{(g-1)}) \quad \|x\|_{C} := x^{\mathrm{T}} C^{-1} x$$

* Compare the ranking of $x_1$ and $x_2$ among $\lambda$ current population

$$s^{(g+1)} = (1 - c_s)s^{(g)} + c_s \underbrace{\frac{\mathrm{rank}(x_2) - \mathrm{rank}(x_1)}{\lambda - 1}}_{> 0 \text{ if forward direction ranks better}}$$

* Update the step-size

$$\sigma^{(g+1)} = \sigma^{(g)} \exp \left( \frac{s^{(g+1)}}{d_{\sigma}} \right)$$

where $c_s \approx 0.3$ and $d_\sigma \approx 1 + 2\ln(n)$

<figure><img src="image_54.png" /></figure>

[Salomon, 1998] Salomon, R. (1998). Evolutionary algorithms and gradient search: similarities and differences. IEEE Transactions on Evolutionary Computation, 2(2), pages 45–55.

[Hansen, 2008] Hansen, N. (2008). CMA-ES with two-point step-size adaptation. [research report] rr-6527, 2008. Inria-00276854v5.

[Hansen et al, 2014] Hansen, N., Atamna, A., and Auger, A. (2014). How to assess step-size adaptation mechanisms in randomised search. In *Parallel Problem Solving from Nature–PPSN XIII*, pages 60–69. Springer.

[Akimoto and Hansen, 2020] Akimoto, Y. and Hansen, N. (2020). Diagonal acceleration for covariance matrix adaptation evolution strategies. *Evolutionary computation*, 28(3), pages 405–435.
# Step-Size Control: Summary

<figure><img src="image_55.png" /></figure>

## Why Step-Size Control?

* to achieve linear convergence at near-optimal rate

## Cumulative Step-Size Adaptation

* efficient and robust for $\lambda \le N$

* inefficient on functions with (many) ineffective axes

## Alternative Step-Size Adaptation Mechanisms

* Two-Point Step-Size Adaptation

* Median Success Rule,[b] Population Success Rule[c]

the effective adaptation of the overall population diversity seems yet to pose open questions, in particular with recombination or without entire control over the realized distribution. a

<footnote><sup>a</sup>Hansen et al. How to Assess Step-Size Adaptation Mechanisms in Randomised Search. PPSN 2014</footnote>

<footnote>[b] Ait Elhara, O., Auger, A., and Hansen, N. (2013). A median success rule for non-elitist evolution strategies: Study of feasibility. In Proc. of the GECCO, pages 415–422.</footnote>

<footnote>[c] Loshchilov, I. (2014). A computationally efficient limited memory CMA-ES for large scale optimization. In Proc. of the GECCO, pages 397–404.</footnote>
Step-Size Control Is Not Enough

<figure><img src="image_56.png" /></figure>

<figure><img src="image_57.png" /></figure>

On 20D TwoAxes Function: $f(\mathbf{x}) = \sum_{i=1}^{\lfloor n/2 \rfloor} x_i^2 + a \sum_{i=\lfloor n/2 \rfloor+1}^{n} x_i^2$

* convergence speed of CSA-ES becomes lower as the function becomes ill conditioned ($a^2$ becomes greater) ⇒ covariance matrix adaptation required
Step-Size Control Is Not Enough

<figure><img src="image_58.png" /></figure>

<figure><img src="image_59.png" /></figure>

On 20D TwoAxes Function: $f(\mathbf{x}) = \sum_{i=1}^{\lfloor n/2 \rfloor} x_i^2 + a \sum_{i=\lfloor n/2 \rfloor+1}^{n} x_i^2$

* convergence speed of CSA-ES becomes lower as the function becomes ill conditioned ($a^2$ becomes greater) ⇒ covariance matrix adaptation required
# Topics

## 0. Problem statement

## 1. What makes an optimization problem difficult to solve?

## 2. How does the CMA-ES work?

* Normal Distribution, Rank-Based Recombination

* Step-Size Adaptation

* Covariance Matrix Adaptation

## 3. What can users do?

* Choice of problem formulation and encoding (not covered)

* Choice of initial solution and initial step-size

* Restarts, Increasing Population Size

* Restricted Covariance Matrix
# Topics

## 0. Problem statement

### 1. What makes an optimization problem difficult to solve?

### 2. How does the CMA-ES work?

* Normal Distribution, Rank-Based Recombination

* Step-Size Adaptation

* Covariance Matrix Adaptation

* Rank-One Update and Cumulation

* Rank-$\mu$ Update

* Active Covariance Update

* Choice of problem formulation and encoding (not covered)

* Choice of initial solution and initial step-size

* Restarts, Increasing Population Size

* Restricted Covariance Matrix
# The Multi-Variate (n-Dimensional) Normal Distribution

Any multi-variate normal distribution $N(m, C)$ is uniquely determined by its mean value $m \in \mathbb{R}^n$ and its symmetric positive definite $n \times n$ covariance matrix $C$.

## The mean value $m$

* determines the displacement (translation)

* value with the largest density (modal value)

* the distribution is symmetric about the distribution mean

<figure><img src="image_60.png" /></figure>

## The covariance matrix $C$

* determines the shape

* geometrical interpretation: any covariance matrix can be uniquely identified with the iso-density ellipsoid $\{x \in \mathbb{R}^n | (x - m)^T C^{-1} (x - m) = n\}$

<figure><img src="image_61.png" /></figure>
... any **covariance matrix** can be uniquely identified with the iso-density ellipsoid
$$ \{ \mathbf{x} \in \mathbb{R}^n \mid (\mathbf{x} - \mathbf{m})^{\mathrm{T}} \mathbf{C}^{-1} (\mathbf{x} - \mathbf{m}) = \mathbf{n} \} $$

<figure><img src="image_62.png" /><figcaption>Lines of Equal Density</figcaption></figure>

where **I** is the identity matrix (isotropic case) and **D** is a diagonal matrix (reasonable for separable problems) and **A** × $\mathcal{N}(\mathbf{0}, \mathbf{I}) \sim \mathcal{N}(\mathbf{0}, \mathbf{A}\mathbf{A}^{\mathrm{T}})$ holds for all **A**.

<figure><img src="image_63.png" /></figure>
# Multivariate Normal Distribution and Eigenvalues

For any positive definite symmetric C,

$$C = d_1^2 \mathbf{b}_1 \mathbf{b}_1^T + \cdots + d_N^2 \mathbf{b}_N \mathbf{b}_N^T \quad (1)$$

*d*<sub>*i*</sub>: square root of the eigenvalue of C

*b*<sub>*i*</sub>: eigenvector of C, corresponding to d<sub>*i*</sub>

The multivariate normal distribution N(m, C)

$$\mathcal{N}(m, C) \sim m + \mathcal{N}(0, d_1^2)\mathbf{b}_1 + \cdots + \mathcal{N}(0, d_N^2)\mathbf{b}_N \quad (2)$$

<figure><img src="image_64.png" /></figure>

(1) ⇔ (2) holds for any **b**<sub>*i*</sub> ∈ ℝ<sup>*n*</sup> and d<sub>*i*</sub> ∈ ℝ
# Evolution Strategies

## Recalling

New search points are sampled normally distributed

$$x_i \sim m + \sigma N_i(0, C) \quad \text{for } i = 1, \dots, \lambda$$

as perturbations of $m$, where $x_i, m \in \mathbb{R}^n, \sigma \in \mathbb{R}_+, C \in \mathbb{R}^{n \times n}$ where

<figure><img src="image_65.png" /></figure>

* the **mean** vector $m \in \mathbb{R}^n$ represents the favorite solution

* the so-called **step-size** $\sigma \in \mathbb{R}_+$ controls the **step length**

* the **covariance matrix** $C \in \mathbb{R}^{n \times n}$ determines the **shape** of the distribution ellipsoid

The remaining question is how to update $C$.

<figure><img src="image_66.png" /></figure>
# Covariance Matrix Adaptation

## Rank-One Update

$$ m \leftarrow m + \sigma y_w, \quad y_w = \sum_{i=1}^{\mu} w_i y_{i:\lambda}, \quad y_i \sim N_i(0, C) $$

<figure><img src="image_67.png" /></figure>

initial distribution, C = I

... equations

<figure><img src="image_68.png" /></figure>
# Covariance Matrix Adaptation

## Rank-One Update

$$m \leftarrow m + \sigma y_w, \quad y_w = \sum_{i=1}^{\mu} w_i y_{i:\lambda}, \quad y_i \sim N_i(0, C)$$

<figure><img src="image_69.png" /></figure>

initial distribution, C = I

... equations

<figure><img src="image_70.png" /></figure>
# Covariance Matrix Adaptation

## Rank-One Update

$$m \leftarrow m + \sigma y_w, \quad y_w = \sum_{i=1}^{\mu} w_i y_{i:\lambda}, \quad y_i \sim N_i(0, C)$$

<figure><img src="image_71.png" /></figure>

$y_w$, movement of the population mean $m$ (disregarding $\sigma$)

<figure><img src="image_72.png" /></figure>
# Covariance Matrix Adaptation

## Rank-One Update

$$ m \leftarrow m + \sigma y_w, \quad y_w = \sum_{i=1}^{\mu} w_i y_{i:\lambda}, \quad y_i \sim N_i(0, C) $$

<figure><img src="image_73.png" /></figure>

$$ \text{mixture of distribution C and step } y_w, \quad C \leftarrow 0.8 \times C + 0.2 \times y_w y_w^T $$

<figure><img src="image_74.png" /></figure>
# Covariance Matrix Adaptation

## Rank-One Update

$$ m \leftarrow m + \sigma y_w, \quad y_w = \sum_{i=1}^{\mu} w_i y_{i:\lambda}, \quad y_i \sim N_i(0, C) $$

<figure><img src="image_75.png" /></figure>

new distribution (disregarding $\sigma$)

<figure><img src="image_76.png" /></figure>
# Covariance Matrix Adaptation

## Rank-One Update

$$m \leftarrow m + \sigma y_w, \quad y_w = \sum_{i=1}^{\mu} w_i y_{i:\lambda}, \quad y_i \sim N_i(0, C)$$

<figure><img src="image_77.png" /></figure>

new distribution (disregarding $\sigma$)

<figure><img src="image_78.png" /></figure>
# Covariance Matrix Adaptation

## Rank-One Update

$$ m \leftarrow m + \sigma y_w, \quad y_w = \sum_{i=1}^{\mu} w_i y_{i:\lambda}, \quad y_i \sim N_i(0, C) $$

<figure><img src="image_79.png" /></figure>

movement of the population mean $m$

<figure><img src="image_80.png" /></figure>
# Covariance Matrix Adaptation

## Rank-One Update

$$ m \leftarrow m + \sigma y_w, \quad y_w = \sum_{i=1}^{\mu} w_i y_{i:\lambda}, \quad y_i \sim N_i(0, C) $$

<figure><img src="image_81.png" /></figure>

$$ \text{mixture of distribution C and step } y_w, \quad C \leftarrow 0.8 \times C + 0.2 \times y_w y_w^T $$

<figure><img src="image_82.png" /></figure>
# Covariance Matrix Adaptation

## Rank-One Update

$$ \mathbf{m} \leftarrow \mathbf{m} + \sigma \mathbf{y}_w, \quad \mathbf{y}_w = \sum_{i=1}^{\mu} w_i \mathbf{y}_{i:\lambda}, \quad \mathbf{y}_i \sim \mathcal{N}_i(\mathbf{0}, \mathbf{C}) $$

<figure><img src="image_83.png" /></figure>

new distribution,

$$ \mathbf{C} \leftarrow 0.8 \times \mathbf{C} + 0.2 \times \mathbf{y}_w \mathbf{y}_w^T $$

the ruling principle: the adaptation **increases the likelihood of**

**successful steps**, $\mathbf{y}_w$, to appear again

another viewpoint: the adaptation **follows a natural gradient**

approximation of the expected fitness

...equations

<figure><img src="image_84.png" /></figure>
# Covariance Matrix Adaptation

## Rank-One Update

Initialize $m \in \mathbb{R}^n$, and $C = I$, set $\sigma = 1$, learning rate $c_{cov} \approx 2/n^2$

While not terminate

$$x_i = m + \sigma y_i, \quad y_i \sim N_i(0, C), \quad i = 1, \dots, \lambda$$

$$m \leftarrow m + \sigma y_w \quad \text{where } y_w = \sum_{i=1}^{\mu} w_i y_{i:\lambda}$$

$$C \leftarrow (1 - c_{cov}) C + c_{cov} \mu_w \underbrace{y_w y_w^T}_{\text{rank-one}} \quad \text{where } \mu_w = \frac{1}{\sum_{i=1}^{\mu} w_i^2} \ge 1$$

The rank-one update has been found independently in several domains<fnref n="6" /> <fnref n="7" /> <fnref n="8" /> <fnref n="9" />

<figure><img src="image_85.png" /></figure>

<footnote n="6">⁶ Kjellström&Taxén 1981. Stochastic Optimization in System Design, IEEE TCS</footnote>

<footnote n="7">⁷ Hansen&Ostermeier 1996. Adapting arbitrary normal mutation distributions in evolution strategies: The covariance matrix adaptation, ICEC</footnote>

<footnote n="8">⁸ Ljung 1999. System Identification: Theory for the User</footnote>

<footnote n="9">⁹ Haario et al 2001. An adaptive Metropolis algorithm, JSTOR</footnote>
$$ C \leftarrow (1 - c_{\text{cov}}) C + c_{\text{cov}} \mu_w y_w y_w^T $$

covariance matrix adaptation

* learns all **pairwise dependencies** between variables

off-diagonal entries in the covariance matrix reflect the dependencies

* conducts a **principle component analysis** (PCA) of steps $y_w$, sequentially in time and space

eigenvectors of the covariance matrix C are the principal components / the principal axes of the mutation ellipsoid

* learns a new **rotated problem representation** [a]

components are independent (only)
in the new representation

<figure><img src="image_86.png" /></figure>

* learns a **new** (Mahalanobis) **metric**

variable metric method

* approximates the **inverse Hessian** on quadratic functions

transformation into the sphere function

* for $\mu = 1$: conducts a **natural gradient ascent** on the distribution $\mathcal{N}$ entirely independent of the given coordinate system

<figure><img src="image_87.png" /></figure>

<footnote>[a] invariant under (independently of) search space rotations</footnote>
# Topics

## 0. Problem statement

### 1. What makes an optimization problem difficult to solve?

### 2. How does the CMA-ES work?

* Normal Distribution, Rank-Based Recombination

* Step-Size Adaptation

* Covariance Matrix Adaptation

* Rank-One Update and Cumulation

* Rank-$\mu$ Update

* Active Covariance Update

* Choice of problem formulation and encoding (not covered)

* Choice of initial solution and initial step-size

* Restarts, Increasing Population Size

* Restricted Covariance Matrix
# Cumulation

## The Evolution Path

## Evolution Path

Conceptually, the evolution path is the **search path** the strategy takes **over a number of generation steps**. It can be expressed as a sum of consecutive **steps** of the mean **m**.

<figure><img src="image_88.png" /></figure>

An exponentially weighted sum of steps $\mathbf{y}_w$ is used

$$ \mathbf{p}_c \propto \sum_{i=0}^{g} \underbrace{(1 - \mathbf{c}_c)^{g-i}}_{\substack{\text{exponentially} \\ \text{fading weights}}} \mathbf{y}_w^{(i)} $$

The recursive construction of the evolution path (cumulation):

$$ \mathbf{p}_c \leftarrow \underbrace{(1 - \mathbf{c}_c)}_{\text{decay factor}} \mathbf{p}_c + \underbrace{\sqrt{1 - (1 - \mathbf{c}_c)^2} \sqrt{\mu_w}}_{\text{normalization factor}} \underbrace{\mathbf{y}_w}_{\text{input}} = \frac{m - m_{old}}{\sigma} $$

where $\mu_w = \frac{1}{\sum w_i^2}$, $\mathbf{c}_c \ll 1$. **History information** is accumulated in the evolution path.
“Cumulation” is a widely used technique and also know as

* exponential smoothing in time series, forecasting

* exponentially weighted *moving average*

* iterate averaging in stochastic approximation

* momentum in the back-propagation algorithm for ANNs

* ...

“Cumulation” conducts a *low-pass filtering*, but there is more to it...

...why?

<figure><img src="image_89.png" /></figure>
# Cumulation

## Utilizing the Evolution Path

We used $y_w y_w^T$ for updating C. Because $y_w y_w^T = -y_w (-y_w)^T$ the sign of $y_w$ is lost.

$$C \leftarrow (1 - c_{\text{cov}})C + c_{\text{cov}} \mu_w y_w y_w^T$$

<figure><img src="image_90.png" /></figure>

The **sign information** (signifying correlation *between* steps) is (re-)introduced by using the **evolution path**.

$$p_c \leftarrow \underbrace{(1 - c_c)}_{\text{decay factor}} p_c + \underbrace{\sqrt{1 - (1 - c_c)^2} \sqrt{\mu_w}}_{\text{normalization factor}} y_w$$

$$C \leftarrow (1 - c_{\text{cov}})C + c_{\text{cov}} \underbrace{p_c p_c^T}_{\text{rank-one}}$$

where $\mu_w = \frac{1}{\sum w_i^2}$, $c_{\text{cov}} \ll c_c \ll 1$ such that $1/c_c$ is the “backward time horizon”.

<figure><img src="image_91.png" /></figure>
# Cumulation

## Utilizing the Evolution Path

We used $y_w y_w^T$ for updating C. Because $y_w y_w^T = -y_w (-y_w)^T$ the sign of $y_w$ is lost.

$$C \leftarrow (1 - c_{\text{cov}})C + c_{\text{cov}} \mu_w y_w y_w^T$$

<figure><img src="image_92.png" /></figure>

The **sign information** (signifying correlation between steps) is (re-)introduced by using the **evolution path**.

$$p_c \leftarrow \underbrace{(1 - c_c)}_{\text{decay factor}} p_c + \underbrace{\sqrt{1 - (1 - c_c)^2} \sqrt{\mu_w}}_{\text{normalization factor}} y_w$$

$$C \leftarrow (1 - c_{\text{cov}})C + c_{\text{cov}} \underbrace{p_c p_c^T}_{\text{rank-one}}$$

where $\mu_w = \frac{1}{\sum w_i^2}$, $c_{\text{cov}} \ll c_c \ll 1$ such that $1/c_c$ is the “backward time horizon”.

<figure><img src="image_93.png" /></figure>
# Cumulation

## Utilizing the Evolution Path

We used $y_w y_w^T$ for updating C. Because $y_w y_w^T = -y_w (-y_w)^T$ the sign of $y_w$ is lost.

$$C \leftarrow (1 - c_{\text{cov}})C + c_{\text{cov}} \mu_w y_w y_w^T$$

<figure><img src="image_94.png" /></figure>

The **sign information** (signifying correlation *between* steps) is (re-)introduced by using the **evolution path**.

$$p_c \leftarrow \underbrace{(1 - c_c)}_{\text{decay factor}} p_c + \underbrace{\sqrt{1 - (1 - c_c)^2} \sqrt{\mu_w}}_{\text{normalization factor}} y_w$$

$$C \leftarrow (1 - c_{\text{cov}})C + c_{\text{cov}} \underbrace{p_c p_c^T}_{\text{rank-one}}$$

where $\mu_w = \frac{1}{\sum w_i^2}$, $c_{\text{cov}} \ll c_c \ll 1$ such that $1/c_c$ is the “backward time horizon”.

<figure><img src="image_95.png" /></figure>
Using an **evolution path** for the **rank-one update** of the covariance matrix reduces the number of function evaluations to adapt to a straight ridge **from about $O(n^2)$ to $O(n)$**. (a)

Number of $f$-evaluations divided by dimension on the cigar function $f(x) = x_1^2 + 10^6 \sum_{i=2}^n x_i^2$

<figure><img src="image_96.png" /></figure>

$c_c = 1$ (no cumulation)

$c_c = 1/\sqrt{n}$

$c_c = 1/n, 3/(n+3)$

with CSA-ES (without CMA): $\approx 3 \times 10^7 n$ evaluations

The overall model complexity is $n^2$ but important parts of the model can be learned in time of order $n$

<footnote>a Hansen & Auger 2013. Principled design of continuous stochastic search: From theory to practice.</footnote>
# Topics

## 0. Problem statement

### 1. What makes an optimization problem difficult to solve?

### 2. How does the CMA-ES work?

* Normal Distribution, Rank-Based Recombination

* Step-Size Adaptation

* Covariance Matrix Adaptation

* Rank-One Update and Cumulation

* Rank-$\mu$ Update

* Active Covariance Update

* Choice of problem formulation and encoding (not covered)

* Choice of initial solution and initial step-size

* Restarts, Increasing Population Size

* Restricted Covariance Matrix
Rank-μ Update

$$ \mathbf{x}_i = \mathbf{m} + \sigma \mathbf{y}_i, \quad \mathbf{y}_i \sim \mathcal{N}_i(\mathbf{0}, \mathbf{C}), $$

$$ \mathbf{m} \leftarrow \mathbf{m} + \sigma \mathbf{y}_w \quad \mathbf{y}_w = \sum_{i=1}^{\mu} w_i \mathbf{y}_{i:\lambda} $$

The rank-μ update extends the update rule for large population sizes λ using μ > 1 vectors to update C at each generation step.

The weighted empirical covariance matrix

$$ \mathbf{C}_{\mu} = \sum_{i=1}^{\mu} w_i \mathbf{y}_{i:\lambda} \mathbf{y}_{i:\lambda}^{\mathrm{T}} $$

computes a weighted mean of the outer products of the best μ steps and has rank min(μ, n) with probability one.

The rank-μ update then reads

$$ \mathbf{C} \leftarrow (1 - c_{\mathrm{cov}}) \mathbf{C} + c_{\mathrm{cov}} \mathbf{C}_{\mu} $$

where $c_{\mathrm{cov}} \approx \mu_w/n^2$ and $c_{\mathrm{cov}} \le 1$.
<figure><img src="image_97.png" /></figure>

$x_i = m + \sigma y_i, y_i \sim \mathcal{N}(0, C)$

<figure><img src="image_98.png" /></figure>

$$
\begin{aligned}
C_{\mu} &= \frac{1}{\mu} \sum_{i:\lambda} y_{i:\lambda} y_{i:\lambda}^{\mathrm{T}} \\
C &\leftarrow (1-1) \times C + 1 \times C_{\mu}
\end{aligned}
$$

<figure><img src="image_99.png" /></figure>

$m_{\text{new}} \leftarrow m + \frac{1}{\mu} \sum y_i : \lambda$

sampling of $\lambda = 150$
solutions where
$\mathbf{C} = \mathbf{I}$ and $\sigma = 1$

calculating **C** where
$\mu = 50$,
$w_1 = \cdots = w_\mu = \frac{1}{\mu}$,
and $c_{cov} = 1$

new distribution

<figure><img src="image_100.png" /></figure>
Rank-μ CMA versus Estimation of Multivariate Normal Algorithm EMNA<sub>global</sub><sup>11</sup>

<figure><img src="image_101.png" /></figure>

$x_i = m_{old} + y_i, \quad y_i \sim \mathcal{N}(\mathbf{0}, \mathbf{C})$

<figure><img src="image_102.png" /></figure>

$\mathbf{C} \leftarrow \frac{1}{\mu} \sum (x_{i:\lambda} - m_{old})(x_{i:\lambda} - m_{old})^\mathrm{T}$

<figure><img src="image_103.png" /></figure>

$m_{new} = m_{old} + \frac{1}{\mu} \sum y_{i:\lambda}$

rank-μ CMA
conducts a
PCA of
steps

<figure><img src="image_104.png" /></figure>

$x_i = m_{old} + y_i, \quad y_i \sim \mathcal{N}(\mathbf{0}, \mathbf{C})$

<figure><img src="image_105.png" /></figure>

$\mathbf{C} \leftarrow \frac{1}{\mu} \sum (x_{i:\lambda} - m_{new})(x_{i:\lambda} - m_{new})^\mathrm{T}$

<figure><img src="image_106.png" /></figure>

$m_{new} = m_{old} + \frac{1}{\mu} \sum y_{i:\lambda}$

EMNA<sub>global</sub>
conducts a
PCA of
points

sampling of $\lambda = 150$
solutions (dots)

calculating **C** from $\mu = 50$
solutions

new distribution

$m_{new}$ is the minimizer for the variances when calculating **C**

<footnote n="11">11 Hansen, N. (2006). The CMA Evolution Strategy: A Comparing Review. In J.A. Lozano, P. Larranga, I. Inza and E. Bengoetxea (Eds.). Towards a new evolutionary computation. Advances in estimation of distribution algorithms. pp. 75-102</footnote>
## The rank-μ update

* increases the possible learning rate in large populations roughly from $2/n^2$ to $\mu_w/n^2$

* can reduce the number of necessary generations roughly from $O(n^2)$ to $O(n)^{(12)}$ given $\mu_w \propto \lambda \propto n$

Therefore the rank-μ update is the primary mechanism whenever a large population size is used

say $\lambda \ge 3n + 10$

## The rank-one update

* uses the evolution path and reduces the number of necessary function evaluations to learn straight ridges from $O(n^2)$ to $O(n)$.

Rank-one update and rank-μ update can be combined

...all equations

<footnote n="12">¹² Hansen, Müller, and Koumoutsakos 2003. Reducing the Time Complexity of the Derandomized Evolution Strategy with Covariance Matrix Adaptation (CMA-ES). *Evolutionary Computation*, 11(1), pp. 1-18</footnote>
## The rank-μ update

* increases the possible learning rate in large populations roughly from $2/n^2$ to $\mu_w/n^2$

* can reduce the number of necessary generations roughly from $O(n^2)$ to $O(n)^{(12)}$ given $\mu_w \propto \lambda \propto n$

Therefore the rank-μ update is the primary mechanism whenever a large population size is used

say $\lambda \ge 3n + 10$

## The rank-one update

* uses the evolution path and reduces the number of necessary function evaluations to learn straight ridges from $O(n^2)$ to $O(n)$.

Rank-one update and rank-μ update can be combined

...all equations

<footnote n="12">¹² Hansen, Müller, and Koumoutsakos 2003. Reducing the Time Complexity of the Derandomized Evolution Strategy with Covariance Matrix Adaptation (CMA-ES). *Evolutionary Computation*, 11(1), pp. 1-18</footnote>
## The rank-μ update

* increases the possible learning rate in large populations roughly from $2/n^2$ to $\mu_w/n^2$

* can reduce the number of necessary generations roughly from $O(n^2)$ to $O(n)^{(12)}$ given $\mu_w \propto \lambda \propto n$

Therefore the rank-μ update is the primary mechanism whenever a large population size is used

say $\lambda \ge 3n + 10$

## The rank-one update

* uses the evolution path and reduces the number of necessary function evaluations to learn straight ridges from $O(n^2)$ to $O(n)$.

Rank-one update and rank-μ update can be combined

...all equations

<footnote n="12">¹² Hansen, Müller, and Koumoutsakos 2003. Reducing the Time Complexity of the Derandomized Evolution Strategy with Covariance Matrix Adaptation (CMA-ES). *Evolutionary Computation*, 11(1), pp. 1-18</footnote>
# Active Update

utilize negative weights [Jastrebski and Arnold, 2006]

$$ C \leftarrow (1 - c_1 - c_\mu + c_\mu^-) C + c_1 p_c p_c^T + c_\mu \sum_{i=1}^{\lfloor \lambda/2 \rfloor} w_i y_{i:\lambda} y_{i:\lambda}^T - c_\mu^- \sum_{i=\lambda-\mu}^{\lambda} |w_i| y_{i:\lambda} y_{i:\lambda}^T $$
$$ \text{decreasing the variances in unpromising directions} $$

* $-|w_i| < 0$ (for $i \ge \lambda - \mu$): negative weight assigned to $y_{i:\lambda}$, $\sum_{i=\lambda-\mu}^{\lambda} |w_i| \approx 1$

* $c_\mu^- > 0$: learning rate for the active update

* increases the variance in the directions of $p_c$ and promising steps $y_{i:\lambda}$ ($i \le \lfloor \lambda/2 \rfloor$)

* decrease the variance in the directions of unpromising steps $y_{i:\lambda}$ ($i \ge \lambda - \lfloor \lambda/2 \rfloor + 1$)

* keep the variance in the subspace orthogonal to the above

[Jastrebski and Arnold, 2006] Jastrebski, G. and Arnold, D. V. (2006). Improving Evolution Strategies through Active Covariance Matrix Adaptation. In 2006 IEEE Congress on Evolutionary Computation, pages 9719–9726.
**Input:** $m \in \mathbb{R}^n$; $\sigma \in \mathbb{R}_+$; $\lambda \in \mathbb{N}_{\ge 2}$, usually $\lambda \ge 5$, default $4 + \lfloor 3 \log n \rfloor$

**Set** $c_m = 1$; $c_1 \approx 2/n^2$; $c_\mu \approx \mu_w/n^2$; $c_c \approx 4/n$; $c_\sigma \approx 1/\sqrt{n}$; $d_\sigma \approx 1$; $w_{i=1...\lambda}$ decreasing in $i$ and $\sum_i^\mu w_i = 1$, $w_\mu > 0 \ge w_{\mu+1}$, $\mu_w^{-1} := \sum_{i=1}^\mu w_i^2 \approx 3/\lambda$

**Initialize** $C = I$, and $p_c = 0$, $p_\sigma = 0$

**While** not terminate

$x_i = m + \sigma y_i$, where $y_i \sim N_i(0, C)$ for $i = 1, ..., \lambda$ sampling

$m \leftarrow m + c_m \sigma y_w$, where $y_w = \sum_{i=1}^\mu w_i y_{rk^{-1}(i)}$ update mean

$p_\sigma \leftarrow (1 - c_\sigma) p_\sigma + \sqrt{1 - (1 - c_\sigma)^2} \sqrt{\mu_w} C^{-\frac{1}{2}} y_w$ path for $\sigma$

$p_c \leftarrow (1 - c_c) p_c + \mathbf{1}_{[0,2n]} \{\|p_\sigma\|^2\} \sqrt{1 - (1 - c_c)^2} \sqrt{\mu_w} y_w$ path for $C$

$\sigma \leftarrow \sigma \times \exp \left( \frac{c_\sigma}{d_\sigma} \left( \frac{\|p_\sigma\|}{E\|N(0,I)\|} - 1 \right) \right)$ update of $\sigma$

$C \leftarrow C + c_\mu \sum_{i=1}^\lambda w_{rk(i)} (y_i y_i^\top - C) + c_1 (p_c p_c^\top - C)$ update $C$

**Not covered:** termination, restarts, useful output, search boundaries and encoding, corrections for: positive definiteness guaranty, $p_c$ variance loss, $c_\sigma$ and $d_\sigma$ for large $\lambda$

<figure><img src="image_107.png" /></figure>
# Python CMA-ES Implementation

https://github.com/

<figure><img src="image_108.png" /></figure>
# Python CMA-ES Demo

https://github.com/CMA-ES/pycma

## Optimizing the 11D Rosenbrock Function

```python
1 import cma
2
3 x, es = cma.fmin2(cma ff.rosen, 11 * [0.1], 0.1, {'verb disp overwrite': 200})
4 es .plot(); # cma .plot() # works too
```

<figure><img src="image_109.png" /></figure>

1.2s

Python

(5_w,11)-aCMA-ES (mu_w=3.4,w_1=42%) in dimension 11 (seed=329070, Fri Sep 13 13:55:03 2024)

<table><thead><tr><th>Iterat</th><th>#Fevals</th><th>function value</th><th>axis ratio</th><th>sigma</th><th>min&amp;max</th><th>std</th><th>t[m:s]</th></tr></thead><tbody><tr><td>1</td><td>11</td><td>1.807461200564364e+01</td><td>1.0e+00</td><td>9.01e-02</td><td>9e-02</td><td>9e-02</td><td>0:00.0</td></tr><tr><td>2</td><td>22</td><td>1.609456060870233e+01</td><td>1.1e+00</td><td>8.18e-02</td><td>8e-02</td><td>8e-02</td><td>0:00.0</td></tr><tr><td>3</td><td>33</td><td>1.501834332056772e+01</td><td>1.2e+00</td><td>7.39e-02</td><td>7e-02</td><td>7e-02</td><td>0:00.0</td></tr><tr><td>100</td><td>1100</td><td>7.225343889168434e+00</td><td>5.6e+00</td><td>1.84e-02</td><td>8e-03</td><td>2e-02</td><td>0:00.1</td></tr><tr><td>200</td><td>2200</td><td>3.475523910468565e+00</td><td>1.1e+01</td><td>2.94e-02</td><td>8e-03</td><td>3e-02</td><td>0:00.1</td></tr><tr><td>619</td><td>6809</td><td>2.441368464931154e-14</td><td>6.4e+01</td><td>2.14e-07</td><td>4e-09</td><td>1e-07</td><td>0:00.3</td></tr></tbody></table>

termination on tolfun=1e-11 (Fri Sep 13 13:55:04 2024)

final/bestever f-value = 1.235402e-14 1.235402e-14 after 6810/6810 evaluations

incumbent solution: [1. 1. 1. 1. 1. 1. 1. 1. ...]

std deviations: [3.73610667e-09 3.79778919e-09 3.80249505e-09 3.97976791e-09 4.15284285e-09 5.01249494e-09 8.22295518e-09 1.37738208e-08 ...]

<figure><img src="image_110.png" /><figcaption>Figure 374</figcaption></figure>
100 1100 7.225343889168434e+00 5.6e+00 1.84e-02 8e-03 2e-02 0:00.1
200 2200 3.475523910468565e+00 1.1e+01 2.94e-02 8e-03 3e-02 0:00.1
619 6809 2.441368464931154e-14 6.4e+01 2.14e-07 4e-09 1e-07 0:00.3

termination on tolfun=1e-11 (Fri Sep 13 13:55:04 2024)

final/bestever f-value = 1.235402e-14 1.235402e-14 after 6810/6810 evaluations

incumbent solution: [1. 1. 1. 1. 1. 1. 1. 1. ...]

std deviations: [3.73610667e-09 3.79778919e-09 3.80249505e-09 3.97976791e-09
4.15284285e-09 5.01249494e-09 8.22295518e-09 1.37738208e-08 ...]

<figure><img src="image_111.png" /><figcaption>Figure 374</figcaption></figure>

<figure><img src="image_112.png" /></figure>
<figure><img src="image_113.png" /></figure>

1 cma.CMAOptions('tol')

<figure><img src="image_114.png" /></figure>

0.0s

Python

{'mindx': '0' #v minimal std in any arbitrary direction, cave interference with tol*',
'minstd': '0' #v minimal std (scalar or vector) in any coordinate direction, cave interference with tol*',
'tolconditioncov': '1e14' #v stop if the condition of the covariance matrix is above `tolconditioncov`',
'tolfacupx': '1e3' #v termination when step-size increases by tolfacupx (diverges). That is, the initial step-size
'tolflatfitness': '1' #v iterations tolerated with flat fitness before termination',
'tolfun': '1e-11' #v termination criterion: tolerance in function value, quite useful',
'tolfunhist': '1e-12' #v termination criterion: tolerance in function value history',
'tolfunrel': '0' #v termination criterion: relative tolerance in function value: Delta f current < tolfunrel * (mean of last 10 generations)
'tolstagnation': 'int(100 + 100 * N**1.5) // popsize' #v termination if no improvement over tolstagnation iteration
'tolupsigma': '1e20' #v sigma/sigma0 > tolupsigma * max(eivenvals(C)**0.5) indicates "creeping behavior" with usual tolerance
'tolx': '1e-11' #v termination criterion: tolerance in x-changes',
'tolxstagnation': '[1e-9, 20, 0.1]' #v termination thresholds for Delta of [mean, iterations, iterations fraction]
# Topics

## 0. Problem statement

1. What makes an optimization problem difficult to solve?

2. How does the CMA-ES work?

* Normal Distribution, Rank-Based Recombination

* Step-Size Adaptation

* Covariance Matrix Adaptation

3. What can users do?

* Choice of problem formulation and encoding (not covered)

* Choice of initial solution and initial step-size

* Restarts, Increasing Population Size

* Restricted Covariance Matrix
# Multimodality

## Approaches for multimodal functions: Try again with

* the final solution as initial solution (non-elitist) and small step-size

* a larger population size

* a different initial mean vector (and a smaller initial step-size)

## A restart with a **large population size** helps if the objective function has a **well global structure**

* functions such as Schaffer, Rastrigin, BBOB function 15~19

* loosely, unimodal global structure + deterministic noise

<figure><img src="image_115.png" /></figure>

<figure><img src="image_116.png" /></figure>

<figure><img src="image_117.png" /></figure>

<figure><img src="image_118.png" /></figure>
# Multimodality

Hansen and Kern. Evaluating the CMA Evolution Strategy on Multimodal Test Functions, PPSN 2004.

<figure><img src="image_119.png" /><figcaption>**Fig. 1.** Success rate to reach $f_{stop} = 10^{-10}$ versus population size for (a) Rastrigin function (b) Griewank function for dimensions $n = 2$ ('--○--'), $n = 5$ ('-.-×-.-'), $n = 10$ ('—□—'), $n = 20$ ('- - + - -'), $n = 40$ ('- · - ◊ - · -'), and $n = 80$ ('—∇—').</figcaption></figure>
# Multimodality

## Approaches for multimodal functions: Try again with

* the final solution as initial solution (non-elitist) and small step-size

* a larger population size

* a different initial mean vector (and a smaller initial step-size)

## A restart with a **small initial step-size** helps if the objective function has a **weak global structure**

* functions such as Schwefel, Bi-Sphere, BBOB function 20~24

<figure><img src="image_120.png" /></figure>

<figure><img src="image_121.png" /></figure>

occasionally, a large population size is even
counterproductive
# Restart Strategy

It makes the CMA-ES parameter free

IPOP: Restart with increasing the population size

* start with the default population size

* double the population size after each trial (parameter sweep)

* may be considered as gold standard for automated restarts

BIPOP: IPOP regime + Local search regime

* IPOP regime: restart with increasing population size

* Local search regime: restart with a smaller step-size and a smaller population size than the IPOP regime
# Summary and Final Remarks
# The Continuous Search Problem

Difficulties of a non-linear optimization problem are

* dimensionality and non-separability

demands to exploit problem structure, e.g. neighborhood
caveat: design of benchmark functions

* ill-conditioning

demands to acquire a second order model

* ruggedness

demands a non-local (stochastic? population based?) approach

<figure><img src="image_122.png" /></figure>
# Main Characteristics of (CMA) Evolution Strategies

1. Multivariate normal distribution to generate new search points
   follows the maximum entropy principle

2. Rank-based selection
   implies invariance, same performance on $g(f(x))$ for any increasing $g$
   more invariance properties are featured

3. Step-size control facilitates fast (log-linear) convergence and
   possibly linear scaling with the dimension
   in CMA-ES based on an **evolution path** (a non-local trajectory)

4. Covariance matrix adaptation (CMA) **increases the likelihood of**
   **previously successful steps** and can improve performance by
   orders of magnitude

the update follows the natural gradient

C $\propto H^{-1} \Leftrightarrow$ adapts a variable metric

$\Leftrightarrow$ new (rotated) problem representation

$\Rightarrow$ $f: x \mapsto g(x^T H x)$ reduces to $x \mapsto x^T x$

<figure><img src="image_123.png" /></figure>
# Limitations

## of CMA Evolution Strategies

* **internal CPU-time**: $10^{-8}n^2$ seconds per function evaluation on a 2GHz PC, tweaks are available

1 000 000 *f*-evaluations in 100-D take 100 seconds *internal* CPU-time

variants with restricted covariance matrix such as Sep-CMA

* better methods are presumably available in case of

* partly separable problems

* specific problems, for example with cheap gradients
    * specific methods

* small dimension ($n \ll 10$)

* for example Nelder-Mead

* small running times (number of *f*-evaluations < $100n$)
    * model-based methods

<figure><img src="image_124.png" /></figure>
Thank you

**Source code** for CMA-ES in C, C++, Java, Matlab, Octave, Python, R, Scilab and

**Practical hints** for problem formulation, variable encoding, parameter setting
are available (or linked to) at
https://cma-es.github.io/
# On Sphere with Low Effective Dimension

On a function with low effective dimension

* $f(\mathbf{x}) = \sum_{i=1}^{M} [\mathbf{x}]_i^2, \mathbf{x} \in \mathbb{R}^N, M \le N$.

* $N - M$ variables do not affect the function value

<figure><img src="image_125.png" /></figure>

<figure><img src="image_126.png" /></figure>
# Maximum Likelihood Update

The new distribution mean **m** maximizes the log-likelihood

$$m_{\text{new}} = \arg \max_{m} \sum_{i=1}^{\mu} w_i \log p_{\mathcal{N}}(\mathbf{x}_{i:\lambda} | m)$$

independently of the given covariance matrix

The rank-$\mu$ update matrix $C_\mu$ maximizes the log-likelihood

$$C_{\mu} = \arg \max_{C} \sum_{i=1}^{\mu} w_i \log p_{\mathcal{N}} \left( \frac{\mathbf{x}_{i:\lambda} - m_{\text{old}}}{\sigma} \middle| m_{\text{old}}, C \right)$$

$\log p_{\mathcal{N}}(\mathbf{x}|m, C) = -\frac{1}{2} \log \det(2\pi C) - \frac{1}{2}(\mathbf{x} - m)^{\mathrm{T}} C^{-1} (\mathbf{x} - m)$
$p_{\mathcal{N}}$ is the density of the multi-variate normal distribution

<figure><img src="image_127.png" /></figure>
# Strategy Internal Parameters

* related to selection and recombination

* $\lambda$, offspring number, new solutions sampled, population size

* $\mu$, parent number, solutions involved in updates of *m*, *C*, and $\sigma$

* $w_{i=1,...,\mu}$, recombination weights

* related to C-update

* $c_c$, decay rate for the evolution path

* $c_1$, learning rate for rank-one update of C

* $c_\mu$, learning rate for rank-$\mu$ update of C

* related to $\sigma$-update

* $c_\sigma$, decay rate of the evolution path

* $d_\sigma$, damping for $\sigma$-change

Parameters were identified in carefully chosen experimental set ups. **Parameters do not in the first place depend on the objective function** and are not meant to be in the users choice.

Only(?) the population size $\lambda$ (and the initial $\sigma$) might be reasonably varied in a wide range, *depending on the objective function*

Useful: restarts with increasing population size (IPOP)

<figure><img src="image_128.png" /></figure>