Skip to content

Bayesian optimization in batches on Hartmann 6-D

The problem

Hartmann's function in 6 dimensions, a sum of four Gaussian wells on [0, 1]⁶, whose global minimum is −3.32237; the Hartmann 6-D example gives its formula and constants. Here it stands for an expensive function that can be evaluated several times at once: a simulation that runs for an hour, on 4 machines. What counts is how many rounds of evaluations the search takes, each as long as one evaluation when its points run side by side, and how many evaluations in all.

What makes it hard

Its second deepest minimum, −3.2032, is nearly as deep as the global one and far from it, with a basin as wide: a search that finds it first can stay there. The Bayesian optimization of one point at a time already has to choose between exploring and refining; a batch has to choose 4 points before it knows the value of any.

Representation

A Real genome of 6 genes in [0, 1]: genoxide's problems::Hartmann6, evaluated in Rust in both languages.

Algorithm

Bo with batches of 4 (.batch(4)), its other settings the defaults: a Latin hypercube of 2(n + 1) = 14 points, then a Gaussian process with Matérn's 5/2 kernel fitted to every evaluation, and the log expected improvement maximized for each point.

The 4 points of a round are chosen one after the other, as Ginsbourger, Le Riche and Carraro (2010) propose: after each, the model is told the point with a fantasized value, without fitting its hyperparameters again, and the next point maximizes the acquisition of that model. The fantasy here is the Kriging believer (bo::Fantasy::KrigingBeliever, the default): the model's own mean at the point. It leaves the model's mean where it was and removes its uncertainty at the point, so the expected improvement there vanishes and the next point goes elsewhere. The constant liar (bo::Fantasy::ConstantLiar) tells the model a fixed value instead, the lowest, mean or highest value so far: the higher the lie, the farther the next points go.

Engine::parallel(true) evaluates the 4 points of a round at once; a seed gives the same points on any number of threads. The run stops within 1e-4 of the minimum, or after 200 evaluations.

As a contrast, the same search one point a round, with the same seed.

Output

The first lines give the problem and the method. Then a row per round: its number (0 is the initial design), the evaluations so far, the best value and its distance above the global minimum. The last lines give the evaluations and rounds each search took to come within 1e-4 of the minimum, and the distance of the batches' best point from the minimum's. Both versions print the same rows: the problem is evaluated in Rust.

The project page plots the best value's distance above the minimum after each round, for both searches.

Good results

A good result is within 1e-4 of −3.32237. The batches get there in round 15, after 74 evaluations, 3.9e-5 above it; one point a round needs 50 evaluations, but 36 rounds: with an evaluation of an hour and 4 machines, 15 hours against 36. A batch spends more evaluations, since each of its points is chosen knowing less than a point chosen after the others' results.

Not every seed finds the global minimum. Over seeds 1 to 20, within 200 evaluations: 13 batch searches came within 1e-4, after 62 to 90 evaluations (a median of 70, 14 rounds); the other 7 stayed at the second minimum, −3.2032. One point a round: 12 of 20 (a median of 50 evaluations and 36 rounds), the others at −3.2032 too. The seed of this example, 3, is the first to reach the global minimum both ways. With the constant liar's lowest value and a larger initial design (30 or 60 points), or the log transform of the values, the share stayed at 10 to 13 of 20 (to 1e-3); another run from another design is the way out, as the Hartmann 6-D example shows for CMA-ES with restarts.

The constant liars do about as well on this problem, to 1e-3 within 200 evaluations: 13 of 20 with the lowest value, 12 with the mean, 11 with the highest, which explores most and needed a median of 86 evaluations, against the believer's 13. On Branin and Hartmann 3, every seed of 20 reached 1e-3 within 80 evaluations with each fantasy: the believer and the lowest lie in a median of 34 evaluations on Branin and 30 on Hartmann 3, the higher lies in 42 to 54 and 40.

Reference: Ginsbourger, D., Le Riche, R. and Carraro, L. (2010). Kriging is well-suited to parallelize optimization. In Computational Intelligence in Expensive Optimization Problems, Springer: 131-162.

Known optimum: −3.32237 at (0.20169, 0.15001, 0.47687, 0.27533, 0.31165, 0.65730) (best known)

Source: examples/bo_hartmann6

Interactive run: tachsin.gr/projects/genoxide/examples/bo-hartmann6

cargo run --release --example bo_hartmann6
//! Bayesian optimization of Hartmann's 6-D function in batches: 4 points a round, chosen one after
//! the other with the Kriging believer and evaluated in parallel, to within 1e-4 of the global
//! minimum. Then, as a contrast, the same search one point a round: fewer evaluations, more rounds.
//!
//! With `GENOXIDE_TRACE=<file>`, it also writes a trace of both runs for the plot on the
//! example's page, with `trace.rs`.
//!
//! ```text
//! cargo run --release --example bo_hartmann6
//! ```

mod trace;

use genoxide::prelude::*;
use genoxide::problems::{Hartmann6, Problem};

// how close to the global minimum, and the evaluations each run may take at most
const TOLERANCE: f64 = 1e-4;
const BUDGET: u64 = 200;
const SEED: u64 = 3;

fn main() -> Result<()> {
    let optimum = Hartmann6.optimum().expect("known");
    let minimum = optimum.value();
    println!("Hartmann's 6-D function in [0, 1]^6: global minimum {minimum:.6}");
    println!(
        "14 points of a Latin hypercube, then 4 points a round by log-EI and the Kriging believer"
    );
    println!("round  evaluations          best     f - f*");
    let (batched, rounds) = search(4, minimum, true)?;
    let best = batched.best_fitness().score().expect("valid");
    let x = batched.best_genome();
    let distance = x
        .iter()
        .zip(&optimum.solutions()[0][..])
        .map(|(a, b)| (a - b) * (a - b))
        .sum::<f64>()
        .sqrt();
    println!(
        "4 points a round: within {TOLERANCE:.0e} of the minimum after {} evaluations in {} \
         rounds, {distance:.1e} from its point",
        batched.evaluations(),
        batched.generations()
    );
    assert!(best - minimum <= TOLERANCE);

    // the contrast: one point a round, the same seed
    let (single, single_rounds) = search(1, minimum, false)?;
    let reached = single.best_fitness().score().expect("valid") - minimum <= TOLERANCE;
    println!(
        "1 point a round:  {} after {} evaluations in {} rounds",
        if reached {
            format!("within {TOLERANCE:.0e} of the minimum")
        } else {
            "not within the tolerance".to_string()
        },
        single.evaluations(),
        single.generations()
    );
    trace::write(&rounds, &single_rounds);
    Ok(())
}

// a search in batches of `batch` points evaluated in parallel, to within the tolerance of
// `minimum` or the budget, printing a row per round if `print`; its outcome, and its best value's
// distance above the minimum after each round
fn search(batch: usize, minimum: f64, print: bool) -> Result<(Outcome<Reals>, Vec<f64>)> {
    let bo = Bo::builder(Hartmann6.representation())
        .batch(batch)
        .minimize()
        .seed(SEED)
        .build()?;
    let mut rounds = Vec::new();
    let outcome = Engine::new(bo, Hartmann6)
        .parallel(true)
        .stop_when(Stop::target(minimum + TOLERANCE).or(Stop::evaluations(BUDGET)))
        .on_generation(|snapshot| {
            let progress = snapshot.progress();
            let best = progress.best().and_then(Fitness::score).expect("valid");
            rounds.push(best - minimum);
            if print {
                println!(
                    "{:>5} {:>12} {:>13.6} {:>10}",
                    progress.generation(),
                    progress.evaluations(),
                    best,
                    format!("{:.1e}", best - minimum)
                );
            }
        })
        .run()?;
    Ok((outcome, rounds))
}
python examples/bo_hartmann6/main.py
"""Bayesian optimization of Hartmann's 6-D function in batches: 4 points a round, chosen one after
the other with the Kriging believer and evaluated in parallel, to within 1e-4 of the global
minimum. Then, as a contrast, the same search one point a round: fewer evaluations, more rounds.

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

    python examples/bo_hartmann6/main.py
"""

import numpy as np

import genoxide as gx

import trace

# how close to the global minimum, and the evaluations each run may take at most
TOLERANCE = 1e-4
BUDGET = 200
SEED = 3


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.Hartmann6()
minimum = problem.optimum.value


def search(batch, show):
    """A search in batches of ``batch`` points evaluated in parallel, to within the tolerance of
    the minimum or the budget, printing a row per round if ``show``; its result, and its best
    value's distance above the minimum after each round."""
    rounds = []

    def on_generation(progress):
        rounds.append(progress.best_fitness - minimum)
        if show:
            print(
                f"{progress.generation:>5} {progress.evaluations:>12} "
                f"{progress.best_fitness:>13.6f} {scientific(progress.best_fitness - minimum):>10}"
            )

    bo = gx.Bo(problem.genome, batch=batch, objective="minimize", seed=SEED)
    result = bo.run(
        problem,
        target=minimum + TOLERANCE,
        evaluations=BUDGET,
        parallel=True,
        on_generation=on_generation,
    )
    return result, rounds


print(f"Hartmann's 6-D function in [0, 1]^6: global minimum {minimum:.6f}")
print("14 points of a Latin hypercube, then 4 points a round by log-EI and the Kriging believer")
print("round  evaluations          best     f - f*")
batched, rounds = search(4, True)
distance = float(np.linalg.norm(batched.best_genome - problem.optimum.solutions[0]))
print(
    f"4 points a round: within 1e-4 of the minimum after {batched.evaluations} evaluations in "
    f"{batched.generations} rounds, {scientific(distance)} from its point"
)
assert batched.best_fitness - minimum <= TOLERANCE

# the contrast: one point a round, the same seed
single, single_rounds = search(1, False)
reached = single.best_fitness - minimum <= TOLERANCE
print(
    f"1 point a round:  {'within 1e-4 of the minimum' if reached else 'not within the tolerance'}"
    f" after {single.evaluations} evaluations in {single.generations} rounds"
)
trace.write(rounds, single_rounds)

What it prints, from a seeded run:

Hartmann's 6-D function in [0, 1]^6: global minimum -3.322368
14 points of a Latin hypercube, then 4 points a round by log-EI and the Kriging believer
round  evaluations          best     f - f*
    0           14     -0.885455      2.4e0
    1           18     -0.942238      2.4e0
    2           22     -0.942238      2.4e0
    3           26     -1.154749      2.2e0
    4           30     -1.170522      2.2e0
    5           34     -1.381891      1.9e0
    6           38     -1.381891      1.9e0
    7           42     -1.875810      1.4e0
    8           46     -2.431634     8.9e-1
    9           50     -2.756340     5.7e-1
   10           54     -2.756340     5.7e-1
   11           58     -2.969910     3.5e-1
   12           62     -3.241253     8.1e-2
   13           66     -3.317829     4.5e-3
   14           70     -3.321314     1.1e-3
   15           74     -3.322329     3.9e-5
4 points a round: within 1e-4 of the minimum after 74 evaluations in 15 rounds, 1.4e-3 from its point
1 point a round:  within 1e-4 of the minimum after 50 evaluations in 36 rounds