# CMA-ES and Advanced Adaptation Mechanisms

Youhei Akimoto, Nikolaus Hansen

► To cite this version:

Youhei Akimoto, Nikolaus Hansen. CMA-ES and Advanced Adaptation Mechanisms. GECCO 2023 - Companion: Companion Conference on Genetic and Evolutionary Computation, Jul 2023, Lisbon, Portugal. ACM, pp.1157-1182, 2023, 10.1145/3583133.3595054. hal-04404028

HAL Id: hal-04404028

https://inria.hal.science/hal-04404028

Submitted on 18 Jan 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>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_1.png" /></figure>

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

# CMA-ES and Advanced Adaptation Mechanisms

## Youhei Akimoto¹ & Nikolaus Hansen²

1. University of Tsukuba, Japan
2. Inria & Ecole Polytechnique, France

akimoto@cs.tsukuba.ac.jp
forename.lastname@inria.fr

Permission to make digital or hard copies of part or all of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the full citation on the first page. Copyrights for third-party components of this work must be honored. For all other uses, contact the owner/author(s).

GECCO '23 Companion, July 15 - 19, 2023, Lisbon, Portugal

© 2023 Copyright held by the owner/author(s).

ACM ISBN 979-8-4007-0120-7...$15.00

https://doi.org/10.1145/3583133.3595054

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

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

# CMA-ES and Advanced Adaptation Mechanisms

Youhei Akimoto¹ & Nikolaus Hansen²

1. University of Tsukuba, Japan

2. Inria & Ecole Polytechnique, France

akimoto@cs.tsukuba.ac.jp
forename.lastname@inria.fr

Permission to make digital or hard copies of part or all of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the full citation on the first page. Copyrights for third-party components of this work must be honored. For all other uses, contact the owner/author(s).

GECCO '23 Companion, July 15 - 19, 2023, Lisbon, Portugal

© 2023 Copyright held by the owner/author(s).

ACM ISBN 979-8-4007-0120-7...$15.00

https://doi.org/10.1145/3583133.3595054

<figure><img src="image_4.png" /></figure>
We are happy to answer questions at any time.
We are happy to answer questions at any time.
# Topics

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/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
# Topics

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/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
# 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_4.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_5.png" /></figure>
# 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_5.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_6.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_6.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) curve fitting, airfoils

* model calibration biological, physical

* parameter calibration 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_7.png" /></figure>
# Topics

## 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/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
# Topics

## 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/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
# 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_7.png" /></figure>

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

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

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

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

* **non-separability**
  dependencies between the objective variables

* **ill-conditioning**

* **non-smooth level sets**

<figure><img src="image_10.png" /><figcaption>gradient direction Newton direction</figcaption></figure>

<figure><img src="image_11.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_8.png" /></figure>

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

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

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

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

* non-separability
  dependencies between the objective variables

* ill-conditioning

* non-smooth level sets

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

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

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

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

<figure><img src="image_13.png" /></figure>
# Ruggedness

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

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

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

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

The term *Curse of dimensionality* (Richard Bellman) refers to problems 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_14.png" /></figure>
# Curse of Dimensionality

The term *Curse of dimensionality* (Richard Bellman) refers to problems 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_14.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_15.png" /></figure>

<figure><img src="image_16.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_15.png" /></figure>

<figure><img src="image_16.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_17.png" /><figcaption>R rotation matrix</figcaption></figure>

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

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

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

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

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

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

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

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

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

<figure><img src="image_27.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>
# 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_17.png" /><figcaption>R rotation matrix</figcaption></figure>

<figure><img src="image_18.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_28.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_29.png" /></figure>
# 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_19.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_20.png" /></figure>
# Non-smooth level sets (sharp ridges)

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

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

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

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

opening angle is the crucial parameter
# Non-smooth level sets (sharp ridges)

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

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

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

<figure><img src="image_23.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>population-based method, stochastic, non-elitistic<br>recombination operator<br>serves as repair mechanism<br>restarts</td></tr></tbody></table>

...metaphors

<figure><img src="image_33.png" /></figure>
# 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>population-based method, stochastic, non-elitistic<br>recombination operator<br>serves as repair mechanism<br>restarts</td></tr></tbody></table>

...metaphors

<figure><img src="image_24.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
# 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
# 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_34.png" /></figure>
# 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_25.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
# 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_35.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 **σ**.

<figure><img src="image_36.png" /></figure>
# 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_26.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 **σ**.

<figure><img src="image_27.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:

$$ \mathcal{N}(\mathbf{x}, \mathbf{A}) + \mathcal{N}(\mathbf{y}, \mathbf{B}) \sim \mathcal{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_37.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_28.png" /></figure>
# Normal Distribution

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

probability density of the 1-D standard normal distribution

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

probability density of a 2-D normal distribution

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

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

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

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

probability density of the 1-D standard normal distribution

probability density of
a 2-D normal
distribution

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

<figure><img src="image_32.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_42.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_43.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_33.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_34.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_44.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_45.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_35.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 $\mathbf{A} \times \mathcal{N}(\mathbf{0}, \mathbf{I}) \sim \mathcal{N}(\mathbf{0}, \mathbf{A}\mathbf{A}^{\mathrm{T}})$ holds for all $\mathbf{A}$.

<figure><img src="image_36.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$$

*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)

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

<figure><img src="image_46.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$$

*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)

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

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

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

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

<figure><img src="image_49.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})||$$

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

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

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

<figure><img src="image_40.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})|| $$

<figure><img src="image_41.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
# 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} + \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.
# 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} + \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_51.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_52.png" /></figure>
# 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_42.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_43.png" /></figure>
# Basic Invariance in Search Space

* translation invariance

is true for most optimization algorithms

<figure><img src="image_53.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_54.png" /></figure>
# Basic Invariance in Search Space

* translation invariance

is true for most optimization algorithms

<figure><img src="image_44.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_45.png" /></figure>
# Summary

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

On 20D Sphere Function: $f(x) = \sum_{i=1}^{N} [x]_{i}^{2}$

* ES without adaptation can't approach the optimum → **adaptation required**
# Summary

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

On 20D Sphere Function: $f(x) = \sum_{i=1}^{N} [x]_{i}^{2}$

* ES without adaptation can't approach the optimum → **adaptation required**
# 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_56.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_57.png" /></figure>
# 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>
# 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

* **two-point adaptation**$^f$ (TPA)
  akin to a rudimentary line search along the mean shift

<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>
# 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

* **two-point adaptation**$^f$ (TPA)
  akin to a rudimentary line search along the mean shift

<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>
# Path Length Control (CSA)

## The Concept of Cumulative Step-Size Adaptation

$$x_i = m + \sigma y_i$$

$$m \leftarrow m + \sigma y_w$$

## Measure the length of the evolution path

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

<figure><img src="image_58.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_59.png" /></figure>
# Path Length Control (CSA)

## The Concept of Cumulative Step-Size Adaptation

$$x_i = m + \sigma y_i$$

$$m \leftarrow m + \sigma 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 $m \in \mathbb{R}^n$, $\sigma \in \mathbb{R}_+$, evolution path $p_\sigma = \mathbf{0}$, set $c_\sigma \approx 4/n$, $d_\sigma \approx 1$.

$m \leftarrow m + \sigma y_w$ where $y_w = \sum_{i=1}^{\mu} w_i 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} y_w$

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

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

## The Equations

Initialize **m** ∈ ℝⁿ, **σ** ∈ ℝ⁺, evolution path **p**_σ = **0**,
set **c**_σ ≈ 4/n, **d**_σ ≈ 1.

**m** ← **m** + **σ**y<sub>w</sub> where y<sub>w</sub> = Σ<sub>i=1</sub><sup>μ</sup> w<sub>i</sub> y<sub>i:</lambda</sub> update mean

**p**_σ ← (1 - c_σ) **p**_σ + √<u>1 - (1 - c_σ)<sup>2</sup></u> √<u>μ_w</u> y_w
accounts for 1-c_σ accounts for w_i

**σ** ← **σ** × exp(<u>(c_σ / d_σ) (||p_σ|| / E||N(0, I)||) - 1</u>) update step-size
>1 ⇐⇒ ||p_σ|| is greater than its expectation

<figure><img src="image_51.png" /></figure>
<figure><img src="image_61.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_62.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)$$

<figure><img src="image_63.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.
# 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}(\mathbf{0}, \mathbf{I})\|}{\|m^{(g)} - m^{(g-1)}\|_{\mathbf{C}^{(g)}}} (m^{(g)} - m^{(g-1)}) \quad \|x\|_{\mathbf{C}} := x^{\mathrm{T}} \mathbf{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)$$

<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.
# 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_64.png" /></figure>

<figure><img src="image_65.png" /></figure>
# 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_55.png" /></figure>

<figure><img src="image_56.png" /></figure>
# Alternatives: Success-Based Step-Size Control

comparing the fitness distributions of current and previous iterations

## Generalizations of 1/5th-success-rule for non-elitist and multi-recombinant ES

* Median Success Rule [Ait Elhara et al., 2013]

* Population Success Rule [Loshchilov, 2014]

controls a *success probability*

## An advantage over CSA and TPA: Cheap Computation

* It depends only on $\lambda$.

* cf. CSA and TPA require a computation of $C^{-1/2}x$ and $C^{-1}x$, respectively.

[Ait Elhara et al., 2013] 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.

[Loshchilov, 2014] Loshchilov, I. (2014). A computationally efficient limited memory cma-es for large scale optimization. In Proc. of the GECCO, pages 397–404.
# Alternatives: Success-Based Step-Size Control

comparing the fitness distributions of current and previous iterations

## Generalizations of 1/5th-success-rule for non-elitist and multi-recombinant ES

* Median Success Rule [Ait Elhara et al., 2013]

* Population Success Rule [Loshchilov, 2014]

controls a *success probability*

## An advantage over CSA and TPA: Cheap Computation

* It depends only on $\lambda$.

* cf. CSA and TPA require a computation of $C^{-1/2}x$ and $C^{-1}x$, respectively.

[Ait Elhara et al., 2013] 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.

[Loshchilov, 2014] Loshchilov, I. (2014). A computationally efficient limited memory cma-es for large scale optimization. In Proc. of the GECCO, pages 397–404.
Step-Size Control Is Not Enough

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

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

On 20D TwoAxes Function: $f(x) = \sum_{i=1}^{N/2} [Rx]_i^2 + a^2 \sum_{i=N/2+1}^{N} [Rx]_i^2$, R: orthogonal

* 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_57.png" /></figure>

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

On 20D TwoAxes Function: $f(x) = \sum_{i=1}^{N/2} [Rx]_i^2 + a^2 \sum_{i=N/2+1}^{N} [Rx]_i^2$, R: orthogonal

* 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_68.png" /></figure>

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

On 20D TwoAxes Function: $f(x) = \sum_{i=1}^{N/2} [Rx]_i^2 + a^2 \sum_{i=N/2+1}^{N} [Rx]_i^2$, R: orthogonal

* 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_59.png" /></figure>

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

On 20D TwoAxes Function: $f(x) = \sum_{i=1}^{N/2} [Rx]_i^2 + a^2 \sum_{i=N/2+1}^{N} [Rx]_i^2$, R: orthogonal

* convergence speed of CSA-ES becomes lower as the function becomes ill conditioned ($a^2$ becomes greater) ⇒ **covariance matrix adaptation required**
# 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_70.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_71.png" /></figure>
# 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_61.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_62.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_72.png" /></figure>

initial distribution, C = I

... equations

<figure><img src="image_73.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_63.png" /></figure>

initial distribution, **C** = **I**

... equations

<figure><img src="image_64.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_74.png" /></figure>

initial distribution, **C** = **I**

... equations

<figure><img src="image_75.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_65.png" /></figure>

initial distribution, C = I

... equations

<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_76.png" /></figure>

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

<figure><img src="image_77.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>

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

<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_78.png" /></figure>

mixture of distribution **C** and step $y_w$,

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

<figure><img src="image_79.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>

mixture of distribution **C** and step $y_w$,

$$C \leftarrow 0.8 \times C + 0.2 \times y_w y_w^T$$

<figure><img src="image_70.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_80.png" /></figure>

new distribution (disregarding $\sigma$)

<figure><img src="image_81.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>

new distribution (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_82.png" /></figure>

new distribution (disregarding $\sigma$)

<figure><img src="image_83.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>

new distribution (disregarding $\sigma$)

<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_84.png" /></figure>

movement of the population mean $m$

<figure><img src="image_85.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>

movement of the population mean $m$

<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_86.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_87.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>

$$ \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_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_88.png" /></figure>

new distribution,

$$ C \leftarrow 0.8 \times C + 0.2 \times y_w y_w^T $$

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

**successful steps**, $y_w$, to appear again

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

approximation of the expected fitness

...equations

<figure><img src="image_89.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>

new distribution,

$$ C \leftarrow 0.8 \times C + 0.2 \times y_w y_w^T $$

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

**successful steps**, $y_w$, to appear again

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

approximation of the expected fitness

...equations

<figure><img src="image_80.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),$$

$$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_90.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>
# 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),$$

$$m \leftarrow m + \sigma y_w \quad \text{where} \quad 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} \quad \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_81.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 principle components / the principle axes of the mutation ellipsoid

* learns a new **rotated problem representation**

components are independent (only) in the new representation

* 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_91.png" /></figure>
$$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 principle components / the principle axes of the mutation ellipsoid

* learns a new **rotated problem representation**

components are independent (only) in the new representation

<figure><img src="image_82.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_83.png" /></figure>
# Invariance Under Rigid Search Space Transformation

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

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

<figure><img src="image_84.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_93.png" /></figure>

<figure><img src="image_94.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_85.png" /></figure>

for example, invariance under search space rotation (separable ⇔ non-separable)
# 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_95.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}$, $c_c \ll 1$. **History information** is accumulated in the evolution path.
# 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_86.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.

<figure><img src="image_87.png" /></figure>
# 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_96.png" /></figure>

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

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

The recursive construction of the evolution path (cumulation):

$$ \mathbf{p}_c \leftarrow \underbrace{(1 - c_e)}_{\text{decay factor}} \mathbf{p}_c + \underbrace{\sqrt{1 - (1 - c_e)^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}$, $c_e \ll 1$. **History information** is accumulated in the evolution path.

<figure><img src="image_97.png" /></figure>
# 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 - c_e)^{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 - c_e)}_{\text{decay factor}} \mathbf{p}_c + \underbrace{\sqrt{1 - (1 - c_e)^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}$, $c_e \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_98.png" /></figure>
“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_99.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_100.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_101.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_102.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

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

## 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.

<figure><img src="image_103.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_104.png" /></figure>
# Cumulation

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

## 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.

<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(\mathbf{x}) = x_1^2 + 10^6 \sum_{i=2}^n x_i^2$

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

$c_c = 1$ (no cumulation)

$c_c = 1/\sqrt{n}$

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

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>
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(\mathbf{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)$

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>
# Rank-μ Update

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

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

$$ 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.

with μ = λ weights can be negative <fnref n="10" />

The rank-μ update then reads

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

where $c_{\mathrm{cov}} \approx \mu_w/n^2$ and $c_{\mathrm{cov}} \le 1$.

<footnote n="10">¹⁰ Jastrebski and Arnold (2006). Improving evolution strategies through active covariance matrix adaptation. CEC.</footnote>
# Rank-μ Update

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

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

$$ 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.

with μ = λ weights can be negative 10

The rank-μ update then reads

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

where $c_{\mathrm{cov}} \approx \mu_w/n^2$ and $c_{\mathrm{cov}} \le 1$.

<footnote n="10">10 Jastrebski and Arnold (2006). Improving evolution strategies through active covariance matrix adaptation. CEC.</footnote>
<figure><img src="image_106.png" /></figure>

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

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

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

<figure><img src="image_108.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_109.png" /></figure>
<figure><img src="image_97.png" /></figure>

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

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

$$
\begin{aligned}
C_{\mu} &= \frac{1}{\mu} \sum 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_110.png" /></figure>

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

<figure><img src="image_111.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_112.png" /></figure>

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

rank-μ CMA
conducts a
PCA of
steps

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

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

<figure><img src="image_114.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_115.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>
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, 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, 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>
## 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>
## Rank-one update

best f-value

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

step-size

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

mean coordinates

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

principle axis length

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

## Rank-μ update

best f-value

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

step-size

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

mean coordinates

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

principle axis length

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

## Hybrid (combined) update

best f-value

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

step-size

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

mean coordinates

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

principle axis length

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

$$f_{\text{TwoAxes}}(x) = \sum_{i=1}^{5} x_i^2 + 10^6 \sum_{i=6}^{10} x_i^2$$

$\lambda = 10$ (default for $N = 10$)
## Rank-one update

best f-value

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

step-size

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

mean coordinates

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

principle axis length

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

## Rank-μ update

best f-value

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

step-size

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

mean coordinates

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

principle axis length

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

## Hybrid (combined) update

best f-value

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

step-size

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

mean coordinates

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

principle axis length

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

$$f_{\text{TwoAxes}}(x) = \sum_{i=1}^{5} x_i^2 + 10^6 \sum_{i=6}^{10} x_i^2$$

$\lambda = 10$ (default for $N = 10$)
Rank-one update

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

Rank-μ update

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

Hybrid (combined) update

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

$$f_{\text{TwoAxes}}(x) = \sum_{i=1}^{5} x_i^2 + 10^6 \sum_{i=6}^{10} x_i^2$$

$$\lambda = 50$$
Rank-one update

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

Rank-μ update

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

Hybrid (combined) update

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

$$f_{\text{TwoAxes}}(x) = \sum_{i=1}^{5} x_i^2 + 10^6 \sum_{i=6}^{10} x_i^2$$

$$\lambda = 50$$
# Different Types of Ill-Conditioning

(α: Axes Ratio = 10)

Cigar Type:

1 long axis = n-1 short axes

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

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

Discus Type:

1 short axis = n-1 long axes

$$f(x) = \alpha \cdot x_1^2 + \sum_{i=1}^{n} x_i^2$$

<figure><img src="image_132.png" /></figure>
# Different Types of Ill-Conditioning

(α: Axes Ratio = 10)

Cigar Type:

1 long axis = n-1 short axes

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

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

Discus Type:

1 short axis = n-1 long axes

$$f(x) = \alpha \cdot x_1^2 + \sum_{i=1}^{n} x_i^2$$

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

utilize negative weights [Jastrebski and Arnold, 2006]

## Active Update (rewriting)

decreasing the variances in unpromising directions

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

* 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.
# Active Update

utilize negative weights [Jastrebski and Arnold, 2006]

## Active Update (rewriting)

decreasing the variances in unpromising directions

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

* 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.
# On 10D Discus Function

10D Discus Function (axis ratio: $\alpha = 10^3$)

$$f(x) = \alpha^2 \cdot x_1^2 + \sum_{i=1}^{n} x_i^2$$

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

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

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

* Positive: wait for the smallest eig(C) decreasing

* Active: decrease the smallest eig(C) actively
# On 10D Discus Function

10D Discus Function (axis ratio: $\alpha = 10^3$)

$$f(x) = \alpha^2 \cdot x_1^2 + \sum_{i=1}^{n} x_i^2$$

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

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

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

* Positive: wait for the smallest eig(C) decreasing

* Active: decrease the smallest eig(C) actively
# Summary

## Active Covariance Matrix Adaptation + Cumulation

$$ 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-\lfloor \lambda/2 \rfloor+1}^{\lambda} |w_i| y_{i:\lambda} y_{i:\lambda}^T $$

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

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

These components compensate each other

* cumulation: excels to learn a long axis, but inefficient for a large $\lambda$

* rank-$\mu$ update: efficient for a large $\lambda$

* active update: effective to learn short axes

An important yet solvable issue of active update

* The positive definiteness of C will be violated if $c_{\mu}^{-}$ is not small enough

* The positive definiteness can be guaranteed w.p.1 by controlling $c_{\mu}^{-} w_i$
# Summary

## Active Covariance Matrix Adaptation + Cumulation

$$ 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-\lfloor \lambda/2 \rfloor+1}^{\lambda} |w_i| y_{i:\lambda} y_{i:\lambda}^T $$

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

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

### These components compensate each other

* cumulation: excels to learn a long axis, but inefficient for a large $\lambda$

* rank-$\mu$ update: efficient for a large $\lambda$

* active update: effective to learn short axes

### An important yet solvable issue of active update

* The positive definiteness of $C$ will be violated if $c_{\mu}^{-}$ is not small enough

* The positive definiteness can be guaranteed w.p.1 by controlling $c_{\mu}^{-} w_i$
**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_136.png" /></figure>
**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_127.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
# 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
# Default Parameter Values

CMA-ES + (B)IPOP Restart Strategy = Quasi-Parameter Free Optimizer

The following parameters were identified in carefully chosen experimental set ups.

* related to selection and recombination

* λ: offspring number, new solutions sampled, population size

* μ: parent number, solutions involved in mean update

* wᵢ: recombination weights

* related to C-update

* 1 - cₑ: decay rate for the evolution path, cumulation factor

* c₁: learning rate for rank-one update of C

* cμ: learning rate for rank-μ update of C

* related to σ-update

* 1 - cσ: decay rate of the evolution path

* dσ: damping for σ-change

The default values depends only on the **dimension**. They do in the first place **not depend on the objective function**.
# Default Parameter Values

CMA-ES + (B)IPOP Restart Strategy = Quasi-Parameter Free Optimizer

The following parameters were identified in carefully chosen experimental set ups.

* related to selection and recombination

* λ: offspring number, new solutions sampled, population size

* μ: parent number, solutions involved in mean update

* wᵢ: recombination weights

* related to C-update

* 1 - cₑ: decay rate for the evolution path, cumulation factor

* c₁: learning rate for rank-one update of C

* cμ: learning rate for rank-μ update of C

* related to σ-update

* 1 - cσ: decay rate of the evolution path

* dσ: damping for σ-change

The default values depends only on the **dimension**. They do in the first place **not depend on the objective function**.
# Parameters to be set depending on the problem

## Initialization and termination conditions

The following should be set or implemented depending on the problem.

* related to the initial search distribution

* $m^{(0)}$: initial mean vector

* $\sigma^{(0)}$ (or $\sqrt{C_{i,i}^{(0)}}$): initial (coordinate-wise) standard deviation

* related to stopping conditions

* max. func. evals.

* max. iterations

* function value tolerance

* min. axis length

* stagnation

## *Practical Hints:*

* start with an initial guess $m^{(0)}$ with a relatively small step-size $\sigma^{(0)}$ to *locally* improve the current guess;

* then increase the step-size, e.g., by factor of 10, to *globally* search for a better solution.
# Parameters to be set depending on the problem

## Initialization and termination conditions

The following should be set or implemented depending on the problem.

* related to the initial search distribution

* $m^{(0)}$: initial mean vector

* $\sigma^{(0)}$ (or $\sqrt{C_{i,i}^{(0)}}$): initial (coordinate-wise) standard deviation

* related to stopping conditions

* max. func. evals.

* max. iterations

* function value tolerance

* min. axis length

* stagnation

## *Practical Hints:*

* start with an initial guess $m^{(0)}$ with a relatively small step-size $\sigma^{(0)}$ to *locally* improve the current guess;

* then increase the step-size, e.g., by factor of 10, to *globally* search for a better solution.
# Python CMA-ES Implementation

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

## pycma

A Python implementation of CMA-ES and a few related numerical optimization tools.

The Covariance Matrix Adaptation Evolution Strategy (CMA-ES) is a stochastic derivative-free numerical optimization algorithm for difficult (non-convex, ill-conditioned, multi-modal, rugged, noisy) optimization problems in continuous search spaces.

**Useful links:**

* A quick start guide with a few usage examples

* The API Documentation

* Hints for how to use this (kind of) optimization module in practice

## Installation of the latest release

### Type

```
python -m pip install cma
```

in a system shell to install the latest *release* from the Python Package Index (PyPI). The release link also provides more installation hints and a quick start guide.
# Python CMA-ES Implementation

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

## pycma

A Python implementation of CMA-ES and a few related numerical optimization tools.

The Covariance Matrix Adaptation Evolution Strategy (CMA-ES) is a stochastic derivative-free numerical optimization algorithm for difficult (non-convex, ill-conditioned, multi-modal, rugged, noisy) optimization problems in continuous search spaces.

**Useful links:**

* A quick start guide with a few usage examples

* The API Documentation

* Hints for how to use this (kind of) optimization module in practice

## Installation of the latest release

### Type

```
python -m pip install cma
```

in a system shell to install the latest *release* from the Python Package Index (PyPI). The release link also provides more installation hints and a quick start guide.
# Python CMA-ES Demo

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

## Optimizing 10D Rosenbrock Function

In [1]:
```python
import cma
opts = cma.CMAOptions()
opts['ftarget'] = 1e-4
opts['maxfevals'] = 1e6

cma.fmin(cma ff.rosen,
         x0=[0.0] * 10,
         sigma0=0.1,
         options=opts)
```

(5_w,10)-aCMA-ES (mu_w=3.2,w_1=45%) in dimension 10 (seed=909490, Mon Apr 16 13:39:57 2018)

Iterat #Fevals function value axis ratio sigma min&max std t[m:s]
1 10 1.169928472214858e+01 1.0e+00 9.12e-02 9e-02 9e-02 0:00.0
2 20 1.363303277917634e+01 1.1e+00 8.33e-02 8e-02 8e-02 0:00.0
3 30 1.232089008099892e+01 1.2e+00 7.55e-02 7e-02 8e-02 0:00.0
100 1000 5.724977739870999e+00 9.1e+00 1.65e-02 7e-03 2e-02 0:00.1
200 2000 2.550841127554589e+00 1.5e+01 3.97e-02 1e-02 4e-02 0:00.2
300 3000 3.674986141687857e-01 1.5e+01 2.76e-02 3e-03 2e-02 0:00.4
400 4000 1.266345464781239e-03 5.0e+01 1.18e-02 8e-04 2e-02 0:00.5
420 4200 7.039461687999381e-05 5.5e+01 4.04e-03 2e-04 5e-03 0:00.5

termination on ftarget=0.0001 (Mon Apr 16 13:39:58 2018)

final/bestever f-value = 2.804423e-05 2.804423e-05

incumbent solution: [ 0.9998542 0.99996219 0.9999681 1.00000445 0.99998977 0.99968537 0.99954974 0.99918266 ...]

std deviations: [ 0.00023937 0.00022203 0.00024836 0.00024782 0.0003 1258 0.00043481 0.00078261 0.0014964 ...]
# Python CMA-ES Demo

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

## Optimizing 10D Rosenbrock Function

In [1]:
```python
import cma
opts = cma.CMAOptions()
opts['ftarget'] = 1e-4
opts['maxfevals'] = 1e6

cma.fmin(cma ff.rosen,
         x0=[0.0] * 10,
         sigma0=0.1,
         options=opts)
```

(5_w,10)-aCMA-ES (mu_w=3.2,w_1=45%) in dimension 10 (seed=909490, Mon Apr 16 13:39:57 2018)

Iterat #Fevals function value axis ratio sigma min&max std t[m:s]
1 10 1.169928472214858e+01 1.0e+00 9.12e-02 9e-02 9e-02 0:00.0
2 20 1.363303277917634e+01 1.1e+00 8.33e-02 8e-02 8e-02 0:00.0
3 30 1.232089008099892e+01 1.2e+00 7.55e-02 7e-02 8e-02 0:00.0
100 1000 5.724977739870999e+00 9.1e+00 1.65e-02 7e-03 2e-02 0:00.1
200 2000 2.550841127554589e+00 1.5e+01 3.97e-02 1e-02 4e-02 0:00.2
300 3000 3.674986141687857e-01 1.5e+01 2.76e-02 3e-03 2e-02 0:00.4
400 4000 1.266345464781239e-03 5.0e+01 1.18e-02 8e-04 2e-02 0:00.5
420 4200 7.039461687999381e-05 5.5e+01 4.04e-03 2e-04 5e-03 0:00.5
termination on ftarget=0.0001 (Mon Apr 16 13:39:58 2018)
final/bestever f-value = 2.804423e-05 2.804423e-05
incumbent solution: [ 0.9998542 0.99996219 0.9999681 1.00000445 0.
99998977 0.99968537
0.99954974 0.99918266 ...]
std deviations: [ 0.00023937 0.00022203 0.00024836 0.00024782 0.0003
1258 0.00043481
0.00078261 0.0014964 ...]
# Python CMA-ES Demo

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

## Optimizing 10D Rosenbrock Function

cma.plot()

<figure><img src="image_137.png" /><figcaption>Figure 325</figcaption></figure>

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

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

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

## Optimizing 10D Rosenbrock Function

cma.plot()

<figure><img src="image_128.png" /><figcaption>Figure 325</figcaption></figure>

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

Two approaches for multimodal functions: Try again with

- a larger population size

- a smaller initial step-size (and random initial mean vector)
# Multimodality

Two approaches for multimodal functions: Try again with

- a larger population size

- a smaller initial step-size (and random initial mean vector)
# 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_140.png" /></figure>

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

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

<figure><img src="image_143.png" /></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 **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_130.png" /></figure>

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

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

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

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

<figure><img src="image_144.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

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

<figure><img src="image_134.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_145.png" /></figure>

a large population size can have a negative effect
# 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_135.png" /></figure>

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

a large population size can have a negative effect
# 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
# 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
# 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
# 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
# Motivation of the Restricted Covariance Matrix

## Bottlenecks of the CMA-ES on high dimensional problems

1. $O(N^2)$ Time and Space Complexities

* to store and update $C \in \mathbb{R}^{N \times N}$

* to compute the eigen decomposition of $C$

2. $O(1/N^2)$ Learning Rates for C-Update

* $c_{\mu} \approx \mu_{w}/N^2$

* $c_1 \approx 2/N^2$

## Exploit prior knowledge on the problem structure such as separability

⇒ decrease the degrees of freedom of the covariance matrix for

* less time and space complexities

* a higher learning rates that potentially accelerate the adaptation
# Motivation of the Restricted Covariance Matrix

## Bottlenecks of the CMA-ES on high dimensional problems

1. $O(N^2)$ Time and Space Complexities

* to store and update $C \in \mathbb{R}^{N \times N}$

* to compute the eigen decomposition of $C$

2. $O(1/N^2)$ Learning Rates for C-Update

* $c_{\mu} \approx \mu_{w}/N^2$

* $c_1 \approx 2/N^2$

## Exploit prior knowledge on the problem structure such as separability

⇒ decrease the degrees of freedom of the covariance matrix for

* less time and space complexities

* a higher learning rates that potentially accelerate the adaptation
# Variants with Restricted Covariance Matrix

## CMA-ES Variants with Restricted Covariance Matrices

* Sep-CMA [Ros and Hansen, 2008]
  ► C = D. D: Diagonal

* VD-CMA [Akimoto et al., 2014]
  ► C = D(I + v v^T)D. D: Diagonal, v ∈ R^N.

* LM-CMA [Loshchilov, 2014]
  ► C = I + Σ_k = 1 v_i v_i^T. v_i ∈ R^N.

* VkD-CMA [Akimoto and Hansen, 2016]
  ► C = D(I + Σ_k = 1 v_i v_i^T)D. v_i ∈ R^N.

[Ros and Hansen, 2008] Ros, R. and Hansen, N. (2008). A simple modification in CMA-ES achieving linear time and space complexity. In Parallel Problem Solving from Nature - PPSN X, pages 296–305. Springer.

[Akimoto et al., 2014] Akimoto, Y., Auger, A., and Hansen, N. (2014). Comparison-based natural gradient optimization in high dimension. In Proceedings of Genetic and Evolutionary Computation Conference, pages 373–380, Vancouver, BC, Canada.

[Loshchilov, 2014] Loshchilov, I. (2014). A computationally efficient limited memory cma-es for large scale optimization. In Proceedings of Genetic and Evolutionary Computation Conference, pages 397–404.

[Akimoto and Hansen, 2016] Akimoto, Y. and Hansen, N. (2016). Projection-based restricted covariance matrix adaptation for high dimension. In Genetic and Evolutionary Computation Conference, GECCO 2016, Denver, Colorado, USA, July 20-24, 2016, page (accepted). ACM.
# Variants with Restricted Covariance Matrix

## CMA-ES Variants with Restricted Covariance Matrices

* Sep-CMA [Ros and Hansen, 2008]
  $C = D. D$: Diagonal

* C = D. D: Diagonal

* VD-CMA [Akimoto et al., 2014]
  $C = D(I + v\mathbf{v}^T)D. D$: Diagonal, $v \in \mathbb{R}^N$.

* C = D(I + v\mathbf{v}^T)D. D: Diagonal, $v \in \mathbb{R}^N$.

* LM-CMA [Loshchilov, 2014]
  $C = I + \sum_{i=1}^{k} v_i v_i^T. v_i \in \mathbb{R}^N$.

* C = I + \sum_{i=1}^{k} v_i v_i^T. v_i \in \mathbb{R}^N$.

* VkD-CMA [Akimoto and Hansen, 2016]
  $C = D(I + \sum_{i=1}^{k} v_i v_i^T)D. v_i \in \mathbb{R}^N$.

* C = D(I + \sum_{i=1}^{k} v_i v_i^T)D. v_i \in \mathbb{R}^N$.

[Ros and Hansen, 2008] Ros, R. and Hansen, N. (2008). A simple modification in CMA-ES achieving linear time and space complexity. In Parallel Problem Solving from Nature - PPSN X, pages 296-305. Springer.

[Akimoto et al., 2014] Akimoto, Y., Auger, A., and Hansen, N. (2014). Comparison-based natural gradient optimization in high dimension. In Proceedings of Genetic and Evolutionary Computation Conference, pages 373-380, Vancouver, BC, Canada.

[Loshchilov, 2014] Loshchilov, I. (2014). A computationally efficient limited memory cma-es for large scale optimization. In Proceedings of Genetic and Evolutionary Computation Conference, pages 397-404.

[Akimoto and Hansen, 2016] Akimoto, Y. and Hansen, N. (2016). Projection-based restricted covariance matrix adaptation for high dimension. In Genetic and Evolutionary Computation Conference, GECCO 2016, Denver, Colorado, USA, July 20-24, 2016, page (accepted). ACM.
# Demo: On 100D Separable Ellipsoid Function

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

Separable-CMA

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

CMA

* CMA needed 10 times more FEs + more CPU time

* However, Sep-CMA won't be able to solve rotated ellipsoid function as efficiently as it solves separable ellipsoid
# Demo: On 100D Separable Ellipsoid Function

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

Separable-CMA

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

CMA

* CMA needed 10 times more FEs + more CPU time

* However, Sep-CMA won't be able to solve rotated ellipsoid function as efficiently as it solves separable ellipsoid
Summary and Final Remarks
# 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 cave: 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_148.png" /></figure>
# The Continuous Search Problem

Difficulties of a non-linear optimization problem are

* dimensionality and non-separability

demands to exploit problem structure, e.g. neighborhood cave: 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_139.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_149.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_140.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_150.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_141.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/
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/