Himmelblau's function
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: where x₁² + x₂ = 11 and x₁ + x₂² = 7. At (3, 2), both hold. At the origin, f = 11² + 7² = 170.
What makes it hard
The equations have four solutions, so the function has four global minima, all of value 0: (3, 2), (−2.805118, 3.131313), (−3.779310, −3.283186) and (3.584428, −1.848127). They lie in four separate basins. A search that converges to one point finds one of them, and which one depends on where it starts. Finding all four takes several searches, or a method that keeps several apart.
Representation
A Real genome of 2 genes, each in [−5, 5], the box around the four minima: the point itself. The
fitness is f, to minimize. The function is genoxide's problems::Himmelblau, which also gives the
four minima, computed to full precision.
Algorithm
Twenty independent local searches, each from a random point, with seeds 1 to 20. Each is a hill climber:
- a step makes 10 neighbors of the current point, each with Gaussian noise on both genes, of standard deviation 0.005 (0.0005 of the range of 10);
- the search moves to the best neighbor only if it is strictly better;
- it stops after 3,000 steps, 30,001 evaluations with the starting point.
With steps that small, a search follows its basin down to the minimum at the bottom. Restarting from random points is the simplest way to find several minima: each start lands in some basin. Each search's end point is then assigned to the nearest of the four known minima.
Output
The first line gives the number of searches, the bounds and the length of each search. Then a
table has a row per known minimum: its coordinates, how many searches ended nearest to it, and the
range of values they reached, from the best to the worst. In Python, run evaluates the function
in Rust, so both versions print the same table.
The values are small but not 0. The step size is fixed, so near a minimum few neighbors are better than the current point, and progress slows down to a stop. The smaller the steps, the closer to 0 a search gets, and the longer it takes to reach a minimum from its start: with steps of 0.01 and 1,000 steps, a quarter of the searches stopped between 1e-6 and 6e-6.
The project page plays this run back.
Good results
Each minimum is worth 0. A good result finds all four, each within 1e-6 of 0. The run finds all four, with 3 to 6 searches each, and every search ends between 3.4e-10 and 2.1e-7.
With seeds 1 to 600, in 30 groups of 20 searches, every search ends within 5.3e-7 of 0, and 29 of the 30 groups find all four minima; the other finds three.
Reference: Himmelblau, D. M. (1972). Applied Nonlinear Programming. McGraw-Hill.
Known optimum: 0 (at four points)
Source: examples/himmelblau
Interactive run: tachsin.gr/projects/genoxide/examples/himmelblau
cargo run --release --example himmelblau
//! Himmelblau: find the four global minima of a two-dimensional function by restarting a local
//! search from random points.
//!
//! Each search is a hill climber with Gaussian steps; it ends in the minimum whose basin it
//! started in. 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 himmelblau
//! ```
mod trace;
use genoxide::prelude::*;
use genoxide::problems::{Himmelblau, Problem};
const SEARCHES: u64 = 20;
const STEPS: u64 = 3_000;
fn main() -> Result<()> {
let problem = Himmelblau;
let optimum = problem.optimum().expect("known");
let minima = optimum.solutions();
// per known minimum: the searches that ended nearest it, and the best and worst values they
// reached
let mut found = vec![(0, f64::INFINITY, 0.0f64); minima.len()];
// with GENOXIDE_TRACE=<file>, a trace of the searches for the plot on the example's page
let mut trace = trace::Trace::from_env();
for seed in 1..=SEARCHES {
let search = LocalSearch::builder(problem.representation())
.neighbor(GaussianMutation::per_gene(1.0, 0.0005)?)
.neighbors(10)
.acceptance(Acceptance::Improving)
.minimize()
.seed(seed)
.build()?;
let outcome = Engine::new(search, problem)
.stop_when(Stop::generations(STEPS))
.on_generation(|snapshot| trace.record(snapshot))
.run()?;
let end = outcome.best_genome();
let nearest = (0..minima.len())
.min_by(|&a, &b| distance(end, &minima[a]).total_cmp(&distance(end, &minima[b])))
.expect("four minima");
let value = outcome.best_fitness().score().expect("valid");
let (searches, best, worst) = &mut found[nearest];
*searches += 1;
*best = best.min(value);
*worst = worst.max(value);
}
println!(
"{SEARCHES} local searches from random points in [-5, 5] x [-5, 5], {STEPS} steps each"
);
println!("minimum searches values reached");
for (minimum, (searches, best, worst)) in minima.iter().zip(found) {
println!(
"({:>9.6}, {:>9.6}) {searches:>8} {} to {}",
minimum[0],
minimum[1],
scientific(best),
scientific(worst)
);
}
trace.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/himmelblau/main.py
"""Himmelblau: find the four global minima of a two-dimensional function by restarting a local
search from random points.
Each search is a hill climber with Gaussian steps; it ends in the minimum whose basin it started
in. 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/himmelblau/main.py
"""
import math
import numpy as np
import genoxide as gx
from trace import Trace
SEARCHES = 20
STEPS = 3_000
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
# per known minimum: the searches that ended nearest it, and the best and worst values they
# reached
found = [[0, math.inf, 0.0] for _ in minima]
# with GENOXIDE_TRACE=<file>, a trace of the searches for the plot on the example's page
trace = Trace(minima)
for seed in range(1, SEARCHES + 1):
search = gx.LocalSearch(
problem.genome,
neighbor=gx.GaussianMutation(0.0005, rate=1.0),
neighbors=10,
acceptance=gx.Improving(),
objective="minimize",
seed=seed,
)
result = search.run(problem, generations=STEPS, on_generation=trace.on_generation)
nearest = int(np.argmin(np.linalg.norm(minima - result.best_genome, axis=1)))
found[nearest][0] += 1
found[nearest][1] = min(found[nearest][1], result.best_fitness)
found[nearest][2] = max(found[nearest][2], result.best_fitness)
print(f"{SEARCHES} local searches from random points in [-5, 5] x [-5, 5], {STEPS} steps each")
print("minimum searches values reached")
for minimum, (searches, best, worst) in zip(minima, found):
print(
f"({minimum[0]:>9.6f}, {minimum[1]:>9.6f}) {searches:>8} "
f"{scientific(best)} to {scientific(worst)}"
)
trace.write()
What it prints, from a seeded run:
20 local searches from random points in [-5, 5] x [-5, 5], 3000 steps each
minimum searches values reached
( 3.000000, 2.000000) 3 3.4e-10 to 7.4e-8
(-2.805118, 3.131313) 6 7.4e-10 to 2.1e-7
(-3.779310, -3.283186) 6 9.7e-9 to 2.0e-7
( 3.584428, -1.848127) 5 1.5e-8 to 1.5e-7