Skip to content

Nelder-Mead with restarts on Himmelblau

The problem

Himmelblau's function (1972) is a sum of two squares in two variables:

f(x₁, x₂) = (x₁² + x₂ − 11)² + (x₁ + x₂² − 7)²

It is 0 exactly where both squares are 0, at four points: (3, 2), (−2.805118, 3.131313), (−3.779310, −3.283186) and (3.584428, −1.848127). These are its four global minima.

What makes it hard

The four minima lie in four separate basins. A local method converges to the minimum of the basin it starts in, so one run finds one minimum, and which one depends on the start. Finding all four takes runs from starts in all four basins.

Representation

A Real genome of 2 genes, each in [−5, 5], the box around the four minima. The fitness is f, to minimize: genoxide's problems::Himmelblau, which also gives the four minima, computed to full precision.

Algorithm

NelderMead with local::Restarts::Random { times: 19 }: a Nelder-Mead simplex search (see the Nelder-Mead on Rosenbrock example for its steps) that starts again from a new random point in the box each time its simplex has converged, 19 times. A run has converged when every vertex of its triangle is within 1e-9 of the first step (1e-10 of the range) of the best vertex in each gene. After the last run, the engine stops with StopReason::Converged. A trial point that leaves the box is mirrored back in, so a triangle that meets the edge of the box doesn't flatten against it.

Each run starts from a random point with a triangle of 0.1 of the range on each side, 1 in both genes, and follows its basin down. The best vertex of each run when it converges is assigned to the nearest known minimum.

Output

The first line gives the number of runs. Then a table has a row per known minimum: its coordinates, how many runs ended at it, and the range of values they reached, from the best to the worst. The last line gives the evaluations of all 20 runs. In Python, run evaluates the function in Rust, so both versions print the same table.

The project page plays the runs back: the triangle on the contour of the function, and where each run so far converged.

Good results

Each minimum is worth 0. A good result finds all four, each within 1e-6 of 0. The 20 runs find all four, 4 to 6 runs each, and every run ends between 2.1e-18 and 6.7e-18, after 140 evaluations on average: 2,808 in all.

With seeds 1 to 300, every run of every seed ends at one of the minima, at most 2.4e-17 from 0, and 294 of the 300 seeds find all four. In the other 6, none of the 20 random starts fell in one of the basins. More restarts make that rarer.

The Himmelblau example finds the same four minima with 20 hill climbers of Gaussian steps, 30,001 evaluations each, ending between 3.4e-10 and 2.1e-7: the simplex adapts its size as it closes in, where fixed steps can't.

Reference: Himmelblau, D. M. (1972). Applied Nonlinear Programming. McGraw-Hill.

Known optimum: 0 (at four points)

Source: examples/nelder_mead_himmelblau

Interactive run: tachsin.gr/projects/genoxide/examples/nelder-mead-himmelblau

cargo run --release --example nelder_mead_himmelblau

//! Nelder-Mead with restarts: find the four global minima of Himmelblau's function with one
//! search that starts again from a random point each time its simplex has converged.
//!
//! Each run converges to the minimum of the basin it starts in; twenty runs land in all four. The
//! known minima come from genoxide's `problems::Himmelblau`.
//!
//! 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 nelder_mead_himmelblau
//! ```

mod trace;

use genoxide::prelude::*;
use genoxide::problems::{Himmelblau, Problem};
use std::cell::RefCell;

const RESTARTS: u64 = 19;

fn main() -> Result<()> {
    let problem = Himmelblau;
    let optimum = problem.optimum().expect("known");
    let minima = optimum.solutions();
    let nelder_mead = NelderMead::builder(problem.representation())
        .restarts(local::Restarts::Random { times: RESTARTS })
        .minimize()
        .seed(1)
        .build()?;
    // the best vertex of each run when it converged
    let mut ends: Vec<Individual<Reals>> = Vec::new();
    // with GENOXIDE_TRACE=<file>, a trace of the runs for the plot on the example's page
    let trace = RefCell::new(trace::Trace::from_env());
    let outcome = Engine::new(nelder_mead, problem)
        .stop_when(Stop::evaluations(100_000))
        .on_generation(|snapshot| trace.borrow_mut().record(snapshot))
        .control(|nelder_mead, _| {
            if nelder_mead.converged() {
                let end = nelder_mead.population()[0].clone();
                trace.borrow_mut().end(end.genome());
                ends.push(end);
            }
            Ok(())
        })
        .run()?;
    assert_eq!(outcome.stop_reason(), StopReason::Converged);

    // per known minimum: the runs that ended nearest it, and the best and worst values they
    // reached
    let mut found = vec![(0, f64::INFINITY, 0.0f64); minima.len()];
    for end in &ends {
        let nearest = (0..minima.len())
            .min_by(|&a, &b| {
                let (a, b) = (
                    distance(end.genome(), &minima[a]),
                    distance(end.genome(), &minima[b]),
                );
                a.total_cmp(&b)
            })
            .expect("four minima");
        let value = end.fitness().and_then(Fitness::score).expect("valid");
        let (runs, best, worst) = &mut found[nearest];
        *runs += 1;
        *best = best.min(value);
        *worst = worst.max(value);
    }

    println!(
        "{} runs of Nelder-Mead in [-5, 5] x [-5, 5]: one, then {RESTARTS} restarts from random points",
        ends.len()
    );
    println!("minimum                  runs  values reached");
    for (minimum, (runs, best, worst)) in minima.iter().zip(found) {
        println!(
            "({:>9.6}, {:>9.6})  {runs:>4}  {} to {}",
            minimum[0],
            minimum[1],
            scientific(best),
            scientific(worst)
        );
    }
    println!(
        "{} evaluations in all, {:.1} per run",
        outcome.evaluations(),
        outcome.evaluations() as f64 / ends.len() as f64
    );
    trace.into_inner().write();
    Ok(())
}

fn distance(a: &Reals, b: &Reals) -> f64 {
    a.iter()
        .zip(b.iter())
        .map(|(x, y)| (x - y) * (x - y))
        .sum::<f64>()
        .sqrt()
}

// two significant digits, e.g. 1.2e-7
fn scientific(value: f64) -> String {
    format!("{value:.1e}")
}

python examples/nelder_mead_himmelblau/main.py

"""Nelder-Mead with restarts: find the four global minima of Himmelblau's function with one search
that starts again from a random point each time its simplex has converged.

Each run converges to the minimum of the basin it starts in; twenty runs land in all four. The
known minima come from genoxide's problems.Himmelblau.

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/nelder_mead_himmelblau/main.py
"""

import math

import numpy as np

import genoxide as gx

from trace import Trace

RESTARTS = 19


def scientific(value):
    """Two significant digits, e.g. 1.2e-7."""
    mantissa, exponent = f"{value:.1e}".split("e")
    return f"{mantissa}e{int(exponent)}"


problem = gx.problems.Himmelblau()
minima = problem.optimum.solutions
nelder_mead = gx.NelderMead(problem.genome, restarts=RESTARTS, objective="minimize", seed=1)
# the best vertex of each run when it converged, and its value
ends = []
# with GENOXIDE_TRACE=<file>, a trace of the runs for the plot on the example's page
trace = Trace(minima)


def control(algorithm, progress):
    if algorithm.converged:
        end = progress.population[0]
        trace.end(end)
        ends.append((end, float(progress.scores[0])))


result = nelder_mead.run(
    problem, evaluations=100_000, on_generation=trace.on_generation, control=control
)
assert result.stop_reason == "converged"

# per known minimum: the runs that ended nearest it, and the best and worst values they reached
found = [[0, math.inf, 0.0] for _ in minima]
for end, value in ends:
    nearest = int(np.argmin(np.linalg.norm(minima - end, axis=1)))
    found[nearest][0] += 1
    found[nearest][1] = min(found[nearest][1], value)
    found[nearest][2] = max(found[nearest][2], value)

print(
    f"{len(ends)} runs of Nelder-Mead in [-5, 5] x [-5, 5]: one, then {RESTARTS} restarts from "
    "random points"
)
print("minimum                  runs  values reached")
for minimum, (runs, best, worst) in zip(minima, found):
    print(
        f"({minimum[0]:>9.6f}, {minimum[1]:>9.6f})  {runs:>4}  "
        f"{scientific(best)} to {scientific(worst)}"
    )
print(f"{result.evaluations} evaluations in all, {result.evaluations / len(ends):.1f} per run")
trace.write()

What it prints, from a seeded run:

20 runs of Nelder-Mead in [-5, 5] x [-5, 5]: one, then 19 restarts from random points
minimum                  runs  values reached
( 3.000000,  2.000000)     5  2.1e-18 to 3.5e-18
(-2.805118,  3.131313)     5  2.5e-18 to 6.7e-18
(-3.779310, -3.283186)     6  3.0e-18 to 6.6e-18
( 3.584428, -1.848127)     4  2.8e-18 to 4.8e-18
2808 evaluations in all, 140.4 per run