Rosenbrock

The problem

Rosenbrock's function sums the same term over each pair of neighboring genes:

f(x) = Σᵢ₌₁ⁿ⁻¹ [100 (xᵢ₊₁ − xᵢ²)² + (xᵢ − 1)²],   each xᵢ in [−30, 30]

Rosenbrock (1960) defined it in two dimensions, with no bounds, to test his method of minimization. This is the chained form in n dimensions, with the bounds of Yao, Liu and Lin (1999, IEEE Transactions on Evolutionary Computation 3(2): 82-102, function f5). It isn't the "extended" form, which sums over the separate pairs x₁x₂, x₃x₄, and so on. The minimum is 0, at (1, …, 1). Here n = 30, Yao, Liu and Lin's dimension.

In two dimensions, the first part of the term, 100 (x₂ − x₁²)², is 0 on the parabola x₂ = x₁². That parabola is the floor of a valley. Along the floor, f is (x₁ − 1)², which falls by only 4 from x₁ = −1 to x₁ = 1. Across the floor, f rises fast: a step of 0.1 away from (1, 1) in x₂ costs 1. Rosenbrock started his searches from (−1.2, 1), where f is 24.2.

What makes it hard

Finding the valley is easy: the steep walls lead to it. Following it is hard. The floor is narrow and bends, so a direction that goes along the floor at one point leaves it a little further on. An algorithm that steps along the axes, or in all directions alike, must take steps small enough to stay inside the valley, and needs many of them to travel along it.

In n dimensions the terms form a chain: each gene wants to be the square of the one before it, xᵢ₊₁ ≈ xᵢ², and x₁ wants to be 1. The genes can't be solved one at a time. Yao, Liu and Lin count the function as unimodal, but for n from 4 to 30 it also has a local minimum (Shang and Qiu, 2006, Evolutionary Computation 14(1): 119-126), with x₁ near −1 from n = 6 on.

Representation

A Real genome of 30 genes, each in [−30, 30]: the point x itself. The fitness is f(x), to minimize. The function is genoxide's problems::Rosenbrock, which brings its bounds and its minimum.

Algorithm

Three algorithms, each with a budget of 10,000 evaluations per dimension, 300,000 in all, and a target of 1e-8 above the minimum.

CMA-ES (Hansen and Ostermeier, 2001, Evolutionary Computation 9(2): 159-195) samples a population from a normal distribution, and adapts its mean, step size and covariance matrix. The covariance matrix is what suits this function: the distribution stretches along the directions of past successful steps, which here are the directions of the valley, and it keeps turning as the valley bends. It uses genoxide's defaults: a population of 4 + ⌊3 ln 30⌋ = 14, a step size of 0.3 of each gene's range and a random start, with IPOP restarts (Auger and Hansen, 2005, IEEE CEC 2005: 1769-1776): a run that converges to the local minimum starts again from a random point with twice the population.

L-SHADE (Tanabe and Fukunaga, 2014, IEEE CEC 2014: 1658-1665) is a differential evolution that adapts its scale factor and crossover rate from successful trials. Its steps are differences between members of the population, so they line up with the valley once the population does. Its population starts at 18 times the number of genes, 540, and shrinks linearly to 4 over the budget.

Particle swarm optimization (Kennedy and Eberhart, 1995, Proceedings of ICNN'95: 1942-1948) moves 40 particles, each pulled towards its own best point and the swarm's best, with Clerc and Kennedy's constriction coefficients (2002, IEEE Transactions on Evolutionary Computation 6(1): 58-73). It has no model of the valley's direction.

Output

The first line gives the dimension, the minimum and the budget. Then one line per algorithm: the error of its best point, which is its value minus the minimum, to 4 decimals; the evaluations it took; and how many of its 30 genes are more than 0.01 from 1. A run stops as soon as it is within 1e-8 of the minimum, so an error of 0.0000 means that it met the target. In Python, run evaluates the function in Rust, so both versions print the same.

The project page plays back another run: CMA-ES on Rosenbrock in 2 dimensions, so that the population can be drawn on the function's contour.

Good results

The minimum is 0. CMA-ES reaches the target after about 56,000 evaluations, in its first run, without restarting, and L-SHADE after about 240,000: the adapted covariance matrix makes CMA-ES four times faster.

Both reach it from every seed tried. With seeds 1 to 100, CMA-ES with IPOP restarts reaches the target every time, after at most 124,000 evaluations. Without restarts, 13 of the 100 runs end in the local minimum, at 3.99 with x₁ ≈ −1, and stay there. With seeds 1 to 20, L-SHADE reaches the target every time, after at most 254,000 evaluations.

PSO ends at an error of 6.79 with the budget spent. Its best point is on the valley floor, but far along it from the minimum: x₁ to x₁₂ are within 0.01 of 1, and from there each gene is about the square of the one before, down to x₂₁ ≈ 0.75, x₂₄ ≈ 0.12 and x₃₀ ≈ 0. The swarm found the valley and moved along it too slowly to reach the end.

Reference: Rosenbrock, H. H. (1960). An automatic method for finding the greatest or least value of a function. The Computer Journal 3(3): 175-184.

Known optimum: 0 (at (1, …, 1))

Source: examples/rosenbrock

Interactive run: tachsin.gr/projects/genoxide/examples/rosenbrock

RustPython

cargo run --release --example rosenbrock

//! Rosenbrock: minimize a function whose minimum lies at the end of a narrow curved valley, in 30
//! dimensions.
//!
//! Compares CMA-ES with IPOP restarts (which learns the valley's direction), L-SHADE (differential
//! evolution with a population that shrinks over the budget) and particle swarm optimization. The
//! global minimum is 0, at (1, …, 1). The function is genoxide's `problems::Rosenbrock`.
//!
//! With `GENOXIDE_TRACE=<file>`, it also writes a trace of its run for the plot on the example's
//! page, with `trace.rs`.
//!
//! ```text
//! cargo run --release --example rosenbrock
//! ```

mod trace;

use genoxide::prelude::*;
use genoxide::problems::{Optimum, Problem, Rosenbrock};

const DIMENSIONS: usize = 30;
const BUDGET: u64 = 10_000 * DIMENSIONS as u64;

fn main() -> Result<()> {
    let problem = Rosenbrock::new(DIMENSIONS);
    let optimum = problem.optimum().expect("known");
    let target = optimum.value() + 1e-8;
    let stop = || Stop::target(target).or(Stop::evaluations(BUDGET));
    println!(
        "Rosenbrock in {DIMENSIONS} dimensions: minimum {:.4}, {BUDGET} evaluations at most",
        optimum.value()
    );

    // IPOP restarts: a run that ends in the local minimum near x₁ = −1 starts again
    let cmaes = Cmaes::builder(problem.representation())
        .restarts(cmaes::Restarts::Ipop)
        .minimize()
        .seed(1)
        .build()?;
    let outcome = Engine::new(cmaes, problem).stop_when(stop()).run()?;
    report("CMA-ES with IPOP", &outcome, &optimum);

    let l_shade = De::l_shade(problem.representation(), BUDGET)
        .minimize()
        .seed(1)
        .build()?;
    let outcome = Engine::new(l_shade, problem).stop_when(stop()).run()?;
    report("L-SHADE", &outcome, &optimum);

    let pso = Pso::builder(problem.representation())
        .population_size(40)
        .minimize()
        .seed(1)
        .build()?;
    let outcome = Engine::new(pso, problem).stop_when(stop()).run()?;
    report("PSO", &outcome, &optimum);

    // 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_2d()?;
    Ok(())
}

// the error to the minimum, the evaluations, and how many genes are more than 0.01 from the
// minimum's
fn report(name: &str, outcome: &Outcome<Reals>, optimum: &Optimum<Reals>) {
    let best = outcome.best_fitness().score().expect("valid");
    // rounding can put a solution a few ulps below the minimum
    let error = (best - optimum.value()).max(0.0);
    let solution = &optimum.solutions()[0];
    let genes = outcome.best().genome().iter().zip(solution.iter());
    let off = genes
        .filter(|(x, minimum)| (*x - *minimum).abs() > 0.01)
        .count();
    let evaluations = outcome.evaluations();
    print!("{name}: error {error:.4} after {evaluations} evaluations, ");
    println!("{off} of {DIMENSIONS} genes off by more than 0.01");
}

python examples/rosenbrock/main.py

"""Rosenbrock: minimize a function whose minimum lies at the end of a narrow curved valley, in 30
dimensions.

Compares CMA-ES with IPOP restarts (which learns the valley's direction), L-SHADE (differential
evolution with a population that shrinks over the budget) and particle swarm optimization. The
global minimum is 0, at (1, …, 1). The function is genoxide's problems.Rosenbrock, which run
evaluates in Rust.

With ``GENOXIDE_TRACE=<file>``, it also writes a trace of its run for the plot on the example's
page, with trace.py.

    python examples/rosenbrock/main.py
"""

import genoxide as gx
import numpy as np

from trace import record_2d

DIMENSIONS = 30
BUDGET = 10_000 * DIMENSIONS

problem = gx.problems.Rosenbrock(DIMENSIONS)
optimum = problem.optimum
target = optimum.value + 1e-8
print(
    f"Rosenbrock in {DIMENSIONS} dimensions: minimum {optimum.value:.4f}, "
    f"{BUDGET} evaluations at most"
)
for name, algorithm in (
    # IPOP restarts: a run that ends in the local minimum near x1 = -1 starts again
    (
        "CMA-ES with IPOP",
        gx.Cmaes(problem.genome, restarts="ipop", objective="minimize", seed=1),
    ),
    ("L-SHADE", gx.De(problem.genome, l_shade=BUDGET, objective="minimize", seed=1)),
    ("PSO", gx.Pso(problem.genome, population_size=40, objective="minimize", seed=1)),
):
    result = algorithm.run(problem, target=target, evaluations=BUDGET)
    # the error to the minimum (rounding can put a solution a few ulps below it), and how many
    # genes are more than 0.01 from the minimum's
    error = max(result.best_fitness - optimum.value, 0.0)
    off = int(np.sum(np.abs(result.best_genome - optimum.solutions[0]) > 0.01))
    print(
        f"{name}: error {error:.4f} after {result.evaluations} evaluations, "
        f"{off} of {DIMENSIONS} genes off by more than 0.01"
    )

# 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_2d()

What it prints, from a seeded run:

Rosenbrock in 30 dimensions: minimum 0.0000, 300000 evaluations at most
CMA-ES with IPOP: error 0.0000 after 55524 evaluations, 0 of 30 genes off by more than 0.01
L-SHADE: error 0.0000 after 240480 evaluations, 0 of 30 genes off by more than 0.01
PSO: error 6.7871 after 300000 evaluations, 18 of 30 genes off by more than 0.01