Weierstrass
The problem
The Weierstrass function sums cosines of growing frequency and shrinking amplitude:
f(x) = Σᵢ Σₖ aᵏ cos(2π bᵏ (xᵢ + 0.5)) − n Σₖ aᵏ cos(π bᵏ), a = 0.5, b = 3, k from 0 to 20
each xᵢ in [−0.5, 0.5]
Its minimum is 0, at the origin: each gene's sum is at least −Σ aᵏ, reached where every cosine is −1, at the integers, and the second term is −n Σ aᵏ, since every bᵏ is odd. Here n = 10. It's the function F11 of the CEC 2005 report (Suganthan et al. 2005), shifted and rotated there, with its constants and bounds; Liang et al. (2006) and the CEC 2014 report have the same form, and BBOB's f16 another one. It's named after Weierstrass's (1872) continuous, nowhere-differentiable function.
What makes it hard
Each gene's sum is a fractal: ripples on ripples, each three times faster and half as high as the last, down to 3²⁰ ≈ 3.5·10⁹ periods per unit. The sum stops there, so the function is smooth, but it's steep and has local minima at every scale. A search that has found the right valley at one scale still has to find it at every finer one.
Representation
A Real genome of 10 genes, each in [−0.5, 0.5]: the point x itself. The fitness is f(x), to
minimize. The function is genoxide's problems::Weierstrass, which brings its bounds and its
minimum.
Algorithm
Five algorithms, each from seeds 1 to 10, with a budget of 10,000 evaluations per dimension, 100,000 per run, and a target of 1e-8:
- CMA-ES (Hansen and Ostermeier, 2001, Evolutionary Computation 9(2): 159-195), which samples a
population of 10 from a normal distribution and adapts its mean, its step size and its covariance
matrix, from a step size of 0.3 of each gene's range and a random start, and stops once it has
converged (
cmaes::Restarts::Stop): sampling on around its point to the end of the budget wouldn't change its best; - the same with IPOP restarts (Auger and Hansen, 2005, IEEE CEC 2005: 1769-1776): a run that has converged starts again from a random point with twice the population;
- differential evolution with genoxide's defaults, SHADE (Tanabe and Fukunaga, CEC 2013), with a population of 100 and its restarts on stagnation;
- particle swarm optimization (Kennedy and Eberhart, 1995), 40 particles with Clerc and Kennedy's constriction coefficients and a global topology;
- a real-coded genetic algorithm: a population of 100, tournaments of 3, simulated binary crossover (Deb and Agrawal, 1995) with η = 15 and polynomial mutation with η = 20 at a rate of 1/10 per gene.
An evaluation is slow, 21 cosines per gene, some of arguments up to 2·10¹⁰, so each generation's points are evaluated in parallel, which gives the same results on any number of threads.
Output
The first line gives the dimension, the seeds and the budget. Then a row per algorithm: how many of
its 10 runs reached the minimum, to within 1e-8, the median of their evaluations (a dash if none
did), and the median of every run's best error, to two significant digits. The function is evaluated
with genoxide's portable math, so the runs are the same on every platform, and in Python, run
evaluates it in Rust, so both versions print the same.
The project page plays back another run: CMA-ES with IPOP restarts on the function in 2 dimensions, so that the population can be drawn on its contour. It meets the target after 816 evaluations.
Good results
The minimum is 0. CMA-ES with IPOP restarts reaches it in all 10 runs, after a median of 10,195 evaluations, and SHADE in all 10, after 42,300. PSO reaches it 8 times, after 17,720. CMA-ES without restarts reaches it once, and its median run ends at 2.6e-3; the genetic algorithm never does, its median run ending at 1.2e-2.
Known optimum: 0 (at the origin)
Source: examples/weierstrass
Interactive run: tachsin.gr/projects/genoxide/examples/weierstrass
cargo run --release --example weierstrass
//! Weierstrass: minimize the Weierstrass function, rippled at every scale, in 10 dimensions.
//!
//! Runs CMA-ES without and with IPOP restarts (a population that doubles at each restart),
//! differential evolution (SHADE), particle swarm optimization and a real-coded genetic algorithm
//! from 10 seeds each, and counts the runs that reach the minimum, 0 at the origin, to within
//! 1e-8. The function is genoxide's `problems::Weierstrass`. Its 21 cosines per gene, some of
//! arguments up to 2·10¹⁰, make each evaluation slow: the runs evaluate in parallel, with the same
//! results on any number of threads, and CMA-ES without restarts stops once it has converged.
//!
//! With `GENOXIDE_TRACE=<file>`, it also writes a trace of a run for the plot on the example's
//! page, with `trace.rs`.
//!
//! ```text
//! cargo run --release --example weierstrass
//! ```
mod trace;
use genoxide::prelude::*;
use genoxide::problems::{Problem, Weierstrass};
const DIMENSIONS: usize = 10;
const SEEDS: u64 = 10;
const BUDGET: u64 = 10_000 * DIMENSIONS as u64;
// a run stops once its error to the minimum is at most this
const ERROR: f64 = 1e-8;
const ALGORITHMS: [&str; 5] = ["CMA-ES", "CMA-ES with IPOP", "DE", "PSO", "GA"];
fn main() -> Result<()> {
let problem = Weierstrass::new(DIMENSIONS);
let minimum = problem.optimum().expect("known").value();
println!(
"Weierstrass in {DIMENSIONS} dimensions, {SEEDS} seeds, {BUDGET} evaluations at most per run"
);
println!("algorithm at min evaluations median error");
for algorithm in ALGORITHMS {
// the evaluations of the runs that reach the minimum, and every run's best error
let mut evaluations = Vec::new();
let mut errors = Vec::new();
for seed in 1..=SEEDS {
let outcome = run(algorithm, problem, seed, minimum + ERROR)?;
if outcome.stop_reason() == StopReason::Target {
evaluations.push(outcome.evaluations() as f64);
}
// rounding can put a solution a few ulps below the minimum
let best = outcome.best_fitness().score().expect("valid");
errors.push((best - minimum).max(0.0));
}
let reached = format!("{}/{SEEDS}", evaluations.len());
let evaluations = median(evaluations).map_or("-".to_string(), |e| format!("{e:.0}"));
let error = median(errors).expect("a run");
println!(
"{algorithm:<16} {reached:>6} {evaluations:>11} {:>12}",
format!("{error:.1e}")
);
}
println!("evaluations: the median of the runs that reach the minimum");
// with GENOXIDE_TRACE=<file>, a trace for the plot on the example's page, of a separate
// run in 2 dimensions: the plot is the function's contour
trace::record_small()?;
Ok(())
}
// a run of `algorithm` from `seed`, until its best is at most `target` or it has used BUDGET
// evaluations
fn run(algorithm: &str, problem: Weierstrass, seed: u64, target: f64) -> Result<Outcome<Reals>> {
let real = problem.representation();
let stop = Stop::target(target).or(Stop::evaluations(BUDGET));
match algorithm {
"CMA-ES" | "CMA-ES with IPOP" => {
// without restarts, the run ends once it has converged: sampling on around its point
// wouldn't change its best
let restarts = if algorithm == "CMA-ES" {
cmaes::Restarts::Stop
} else {
cmaes::Restarts::Ipop
};
let cmaes = Cmaes::builder(real)
.restarts(restarts)
.minimize()
.seed(seed)
.build()?;
Engine::new(cmaes, problem)
.stop_when(stop)
.parallel(true)
.run()
}
"DE" => {
let de = De::builder(real).minimize().seed(seed).build()?;
Engine::new(de, problem)
.stop_when(stop)
.parallel(true)
.run()
}
"PSO" => {
let pso = Pso::builder(real)
.population_size(40)
.minimize()
.seed(seed)
.build()?;
Engine::new(pso, problem)
.stop_when(stop)
.parallel(true)
.run()
}
_ => {
let ga = Ga::builder(real)
.population_size(100)
.select(Tournament::new(3)?)
.crossover(SimulatedBinaryCrossover::new(15.0)?)
.mutate(PolynomialMutation::per_gene(1.0 / DIMENSIONS as f64, 20.0)?)
.minimize()
.seed(seed)
.build()?;
Engine::new(ga, problem)
.stop_when(stop)
.parallel(true)
.run()
}
}
}
// the median of `values`, None without any
fn median(mut values: Vec<f64>) -> Option<f64> {
values.sort_by(f64::total_cmp);
let middle = values.len() / 2;
match values.len() {
0 => None,
n if n % 2 == 1 => Some(values[middle]),
_ => Some((values[middle - 1] + values[middle]) / 2.0),
}
}
python examples/weierstrass/main.py
"""Weierstrass: minimize the Weierstrass function, rippled at every scale, in 10 dimensions.
Runs CMA-ES without and with IPOP restarts (a population that doubles at each restart), differential
evolution (SHADE), particle swarm optimization and a real-coded genetic algorithm from 10 seeds
each, and counts the runs that reach the minimum, 0 at the origin, to within 1e-8. The function is
genoxide's `problems::Weierstrass`, which run evaluates in Rust. Its 21 cosines per gene, some of
arguments up to 2·10¹⁰, make each evaluation slow: the runs evaluate in parallel, with the same
results on any number of threads, and CMA-ES without restarts stops once it has converged.
With ``GENOXIDE_TRACE=<file>``, it also writes a trace of a run for the plot on the example's page,
with trace.py.
python examples/weierstrass/main.py
"""
import genoxide as gx
from trace import record_small
DIMENSIONS = 10
SEEDS = 10
BUDGET = 10_000 * DIMENSIONS
# a run stops once its error to the minimum is at most this
ERROR = 1e-8
ALGORITHMS = ["CMA-ES", "CMA-ES with IPOP", "DE", "PSO", "GA"]
def build(name, genome, seed):
"""The algorithm called ``name``, on ``genome``, from ``seed``."""
if name == "CMA-ES":
# without restarts, the run ends once it has converged: sampling on around its point
# wouldn't change its best
return gx.Cmaes(genome, restarts="stop", objective="minimize", seed=seed)
if name == "CMA-ES with IPOP":
return gx.Cmaes(genome, restarts="ipop", objective="minimize", seed=seed)
if name == "DE":
return gx.De(genome, objective="minimize", seed=seed)
if name == "PSO":
return gx.Pso(genome, population_size=40, objective="minimize", seed=seed)
return gx.Ga(
genome,
population_size=100,
select=gx.Tournament(3),
crossover=gx.SimulatedBinaryCrossover(15.0),
mutation=gx.PolynomialMutation(20.0, rate=1 / DIMENSIONS),
objective="minimize",
seed=seed,
)
def median(values):
"""The median of ``values``, None without any."""
values = sorted(values)
middle = len(values) // 2
if not values:
return None
return values[middle] if len(values) % 2 else (values[middle - 1] + values[middle]) / 2
def error_text(error):
"""An error to two significant digits, as Rust writes it: 9.9e-9."""
mantissa, exponent = f"{error:.1e}".split("e")
return f"{mantissa}e{int(exponent)}"
problem = gx.problems.Weierstrass(DIMENSIONS)
minimum = problem.optimum.value
print(
f"Weierstrass in {DIMENSIONS} dimensions, {SEEDS} seeds, "
f"{BUDGET} evaluations at most per run"
)
print("algorithm at min evaluations median error")
for name in ALGORITHMS:
# the evaluations of the runs that reach the minimum, and every run's best error
evaluations, errors = [], []
for seed in range(1, SEEDS + 1):
result = build(name, problem.genome, seed).run(
problem, target=minimum + ERROR, evaluations=BUDGET, parallel=True
)
if result.stop_reason == "target":
evaluations.append(float(result.evaluations))
# rounding can put a solution a few ulps below the minimum
errors.append(max(result.best_fitness - minimum, 0.0))
reached = f"{len(evaluations)}/{SEEDS}"
middle = median(evaluations)
evaluations_text = "-" if middle is None else f"{middle:.0f}"
error = error_text(median(errors))
print(f"{name:<16} {reached:>6} {evaluations_text:>11} {error:>12}")
print("evaluations: the median of the runs that reach the minimum")
# with GENOXIDE_TRACE=<file>, a trace for the plot on the example's page, of a separate run in
# 2 dimensions: the plot is the function's contour
record_small()
What it prints, from a seeded run:
Weierstrass in 10 dimensions, 10 seeds, 100000 evaluations at most per run
algorithm at min evaluations median error
CMA-ES 1/10 4070 2.6e-3
CMA-ES with IPOP 10/10 10195 8.3e-9
DE 10/10 42300 8.8e-9
PSO 8/10 17720 9.1e-9
GA 0/10 - 1.2e-2
evaluations: the median of the runs that reach the minimum