Rotated hyper-ellipsoid
The problem
The rotated hyper-ellipsoid, as Molga and Smutnicki (2005, section 2.3) define it, sums the sums of squares of the first genes:
f(x) = Σᵢ Σⱼ≤ᵢ xⱼ², each xᵢ in [−65.536, 65.536]
Its minimum is 0, at the origin. Here n = 30. Its origin is unknown; genoxide takes it, and its bounds, from Molga and Smutnicki, and it's still to be checked against an original (issue #168).
Despite its name, it isn't rotated. Gene j appears in the n − j + 1 sums from i = j on, so
f(x) = Σⱼ (n − j + 1) xⱼ²
an ellipsoid along the axes, with weights from 30 down to 1: genoxide's axis-parallel ellipsoid with its genes reversed. Molga and Smutnicki describe Schwefel's problem 1.2, Σᵢ (Σⱼ≤ᵢ xⱼ)², whose ellipsoids are rotated, but write this formula, which other collections repeat.
What makes it hard
As written, little: it's a separable quadratic with a condition number of 30, which a search can solve a gene at a time.
The example then rotates it, with genoxide's problems::Rotated and seed 1: an orthogonal matrix,
drawn from normal numbers made orthonormal by Gram-Schmidt as BBOB draws its rotations, turns the
function about its minimum. The ellipsoid's axes are no longer the genes' axes, and the genes
interact: a scale per gene no longer fits it, nor do steps along the axes. The plan of genoxide's
test problems suggested this page for comparing CMA-ES with a full and with a diagonal covariance
matrix.
Representation
A Real genome of 30 genes, each in [−65.536, 65.536]: the point x itself. The fitness is f(x), to
minimize. The function is genoxide's problems::RotatedHyperEllipsoid, and its rotation
problems::Rotated, which keeps the bounds and the minimum.
Algorithm
Five algorithms, each with a budget of 10,000 evaluations per dimension, 300,000 in all, and a target of 1e-8, from seed 1:
- CMA-ES (Hansen and Ostermeier, 2001, Evolutionary Computation 9(2): 159-195), which samples a population of 14 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;
- sep-CMA-ES (Ros and Hansen, 2008, PPSN X: 296-305), the same with a diagonal covariance matrix: a scale per gene but no correlations;
- differential evolution with genoxide's defaults, SHADE (Tanabe and Fukunaga, CEC 2013), with a population of 100;
- 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/30 per gene.
Output
The first line gives the dimension and the budget. Then two tables, the function as written and
rotated by an orthogonal matrix: a row per algorithm, the evaluations it had used when its best
error first reached each value of the heading, and the best error it found, to two significant
digits. A dash is an error not reached. 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 on the function in 2 dimensions, 2x₁² + x₂², rotated with seed 1, so that the population can be drawn on its contour. It meets the target after 342 evaluations.
Good results
The minimum is 0. As written, sep-CMA-ES reaches 1e-8 first, after 4,886 evaluations, then CMA-ES (7,140), PSO (27,520) and SHADE (31,900); the genetic algorithm ends at 1.1e-3.
Rotated, CMA-ES takes about as long, 7,294 evaluations: its full covariance matrix learns the ellipsoid's axes, whichever they are. sep-CMA-ES takes 2.5 times as long, 12,180, since its diagonal matrix can only scale the genes; with a condition number of 30, it still gets there. SHADE and PSO slow down by two and three and a half times (59,000 and 96,600), and the genetic algorithm, which recombines genes position by position, ends at 4.2.
Reference: Molga, M. and Smutnicki, C. (2005). Test functions for optimization needs.
Known optimum: 0 (at the origin)
Source: examples/rotated_hyper_ellipsoid
Interactive run: tachsin.gr/projects/genoxide/examples/rotated-hyper-ellipsoid
cargo run --release --example rotated_hyper_ellipsoid
//! Rotated hyper-ellipsoid: minimize the "rotated" hyper-ellipsoid in 30 dimensions, which is in fact
//! axis-parallel, and then the same function truly rotated.
//!
//! Compares how fast CMA-ES, with a full and with a diagonal covariance matrix (sep-CMA-ES),
//! differential evolution, particle swarm optimization and a real-coded genetic algorithm close in
//! on the minimum, 0 at the origin: the evaluations each takes until its error is at most 1, 1e-2, 1e-4, 1e-6 and 1e-8.
//! The function is genoxide's `problems::RotatedHyperEllipsoid`. Then the same on the function rotated by an orthogonal matrix, with genoxide's
//! `problems::Rotated`.
//!
//! 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 rotated_hyper_ellipsoid
//! ```
mod trace;
use genoxide::observer::Snapshot;
use genoxide::prelude::*;
use genoxide::problems::{Problem, Rotated, RotatedHyperEllipsoid};
const DIMENSIONS: usize = 30;
const BUDGET: u64 = 10_000 * DIMENSIONS as u64;
// the errors at which the table gives each run's evaluations
const ERRORS: [f64; 5] = [1e0, 1e-2, 1e-4, 1e-6, 1e-8];
const COLUMNS: [&str; 5] = ["1", "1e-2", "1e-4", "1e-6", "1e-8"];
fn main() -> Result<()> {
println!("Rotated hyper-ellipsoid in {DIMENSIONS} dimensions, {BUDGET} evaluations at most");
compare(
"Rotated hyper-ellipsoid",
&RotatedHyperEllipsoid::new(DIMENSIONS),
)?;
// the same function rotated by an orthogonal matrix: now the genes interact
let rotated = Rotated::new(RotatedHyperEllipsoid::new(DIMENSIONS), 1);
compare("Rotated by an orthogonal matrix (seed 1)", &rotated)?;
// 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(())
}
// the table of the five algorithms on `problem`, after a line that names it
fn compare<P>(name: &str, problem: &P) -> Result<()>
where
P: Problem<Representation = Real> + FitnessFunction<Reals, Output = f64> + Clone,
{
let minimum = problem.optimum().expect("known").value();
let stop = || Stop::target(minimum + 1e-8).or(Stop::evaluations(BUDGET));
println!("{name}: evaluations until the error is at most");
print!("{:<10}", "algorithm");
COLUMNS.iter().for_each(|column| print!("{column:>9}"));
println!("{:>9}", "best");
for (name, covariance) in [
("CMA-ES", cmaes::Covariance::Full),
("sep-CMA-ES", cmaes::Covariance::Diagonal),
] {
let cmaes = Cmaes::builder(problem.representation())
.covariance(covariance)
.minimize()
.seed(1)
.build()?;
let mut reached = Reached::new(minimum);
let outcome = Engine::new(cmaes, problem.clone())
.stop_when(stop())
.on_generation(|snapshot| reached.record(snapshot))
.run()?;
reached.print(name, &outcome);
}
let de = De::builder(problem.representation())
.minimize()
.seed(1)
.build()?;
let mut reached = Reached::new(minimum);
let outcome = Engine::new(de, problem.clone())
.stop_when(stop())
.on_generation(|snapshot| reached.record(snapshot))
.run()?;
reached.print("DE", &outcome);
let pso = Pso::builder(problem.representation())
.population_size(40)
.minimize()
.seed(1)
.build()?;
let mut reached = Reached::new(minimum);
let outcome = Engine::new(pso, problem.clone())
.stop_when(stop())
.on_generation(|snapshot| reached.record(snapshot))
.run()?;
reached.print("PSO", &outcome);
let ga = Ga::builder(problem.representation())
.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(1)
.build()?;
let mut reached = Reached::new(minimum);
let outcome = Engine::new(ga, problem.clone())
.stop_when(stop())
.on_generation(|snapshot| reached.record(snapshot))
.run()?;
reached.print("GA", &outcome);
Ok(())
}
// the evaluations after the first generation whose best error was at most each of ERRORS, for a
// function whose minimum is `minimum`
struct Reached {
minimum: f64,
evaluations: [Option<u64>; 5],
}
impl Reached {
fn new(minimum: f64) -> Self {
let evaluations = [None; 5];
Self {
minimum,
evaluations,
}
}
fn record(&mut self, snapshot: &Snapshot<'_, Reals>) {
let progress = snapshot.progress();
let Some(best) = progress.best().and_then(Fitness::score) else {
return;
};
let error = best - self.minimum;
for (reached, bound) in self.evaluations.iter_mut().zip(ERRORS) {
if reached.is_none() && error <= bound {
*reached = Some(progress.evaluations());
}
}
}
// a row of the table: the evaluations, "-" for an error not reached, and the best error
fn print(&self, name: &str, outcome: &Outcome<Reals>) {
print!("{name:<10}");
for reached in self.evaluations {
let reached = reached.map_or("-".to_string(), |evaluations| evaluations.to_string());
print!("{reached:>9}");
}
// rounding can put a solution a few ulps below the minimum
let best = outcome.best_fitness().score().expect("valid");
let error = (best - self.minimum).max(0.0);
println!("{:>9}", format!("{error:.1e}"));
}
}
python examples/rotated_hyper_ellipsoid/main.py
"""Rotated hyper-ellipsoid: minimize the "rotated" hyper-ellipsoid in 30 dimensions, which is in
fact axis-parallel, and then the same function truly rotated.
Compares how fast CMA-ES, with a full and with a diagonal covariance matrix (sep-CMA-ES),
differential evolution, particle swarm optimization and a real-coded genetic algorithm close in on
the minimum, 0 at the origin: the evaluations each takes until its error is at most 1, 1e-2, 1e-4,
1e-6 and 1e-8. The function is genoxide's `problems::RotatedHyperEllipsoid`. Then the same on the
function rotated by an orthogonal matrix, with genoxide's `problems::Rotated`.
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/rotated_hyper_ellipsoid/main.py
"""
import genoxide as gx
from trace import record_small
DIMENSIONS = 30
BUDGET = 10_000 * DIMENSIONS
# the errors at which the table gives each run's evaluations
ERRORS = [1e0, 1e-2, 1e-4, 1e-6, 1e-8]
COLUMNS = ["1", "1e-2", "1e-4", "1e-6", "1e-8"]
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)}"
class Reached:
"""The evaluations after the first generation whose best error was at most each of ERRORS,
for a function whose minimum is ``minimum``."""
def __init__(self, minimum):
self.minimum = minimum
self.evaluations = [None] * len(ERRORS)
def record(self, progress):
if progress.best_fitness is None:
return
error = progress.best_fitness - self.minimum
for i, bound in enumerate(ERRORS):
if self.evaluations[i] is None and error <= bound:
self.evaluations[i] = progress.evaluations
def print(self, name, result):
"""A row of the table: the evaluations, "-" for an error not reached, and the best
error."""
cells = ["-" if reached is None else str(reached) for reached in self.evaluations]
# rounding can put a solution a few ulps below the minimum
cells.append(error_text(max(result.best_fitness - self.minimum, 0.0)))
print(f"{name:<10}" + "".join(f"{cell:>9}" for cell in cells))
def compare(name, problem):
"""The table of the five algorithms on ``problem``, after a line that names it."""
minimum = problem.optimum.value
print(f"{name}: evaluations until the error is at most")
print(f"{'algorithm':<10}" + "".join(f"{column:>9}" for column in COLUMNS) + f"{'best':>9}")
genome = problem.genome
for label, algorithm in (
("CMA-ES", gx.Cmaes(genome, objective="minimize", seed=1)),
("sep-CMA-ES", gx.Cmaes(genome, covariance="diagonal", objective="minimize", seed=1)),
("DE", gx.De(genome, objective="minimize", seed=1)),
("PSO", gx.Pso(genome, population_size=40, objective="minimize", seed=1)),
(
"GA",
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=1,
),
),
):
reached = Reached(minimum)
result = algorithm.run(
problem, target=minimum + 1e-8, evaluations=BUDGET, on_generation=reached.record
)
reached.print(label, result)
print(f"Rotated hyper-ellipsoid in {DIMENSIONS} dimensions, {BUDGET} evaluations at most")
compare("Rotated hyper-ellipsoid", gx.problems.RotatedHyperEllipsoid(DIMENSIONS))
# the same function rotated by an orthogonal matrix: now the genes interact
rotated = gx.problems.Rotated(gx.problems.RotatedHyperEllipsoid(DIMENSIONS), seed=1)
compare("Rotated by an orthogonal matrix (seed 1)", rotated)
# 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:
Rotated hyper-ellipsoid in 30 dimensions, 300000 evaluations at most
Rotated hyper-ellipsoid: evaluations until the error is at most
algorithm 1 1e-2 1e-4 1e-6 1e-8 best
CMA-ES 3262 4410 5362 6160 7140 8.7e-9
sep-CMA-ES 2142 2772 3556 4354 4886 9.6e-9
DE 13500 18100 23100 27100 31900 7.9e-9
PSO 10960 16320 20280 24160 27520 9.4e-9
GA 44250 138965 - - - 1.1e-3
Rotated by an orthogonal matrix (seed 1): evaluations until the error is at most
algorithm 1 1e-2 1e-4 1e-6 1e-8 best
CMA-ES 3430 4536 5348 6314 7294 8.8e-9
sep-CMA-ES 3920 6342 8540 10668 12180 9.5e-9
DE 18300 30500 41100 50500 59000 9.9e-9
PSO 25120 42640 57360 79760 96600 1.0e-8
GA - - - - - 4.2e0