Constrained Bayesian optimization
The problem
The toy problem of Gramacy et al. (2016, section 1): minimize a linear function of two variables in the unit square, subject to two nonlinear constraints,
minimize f(x) = x₁ + x₂, x in [0, 1]²
subject to c₁(x) = 3/2 − x₁ − 2x₂ − ½ sin(2π(x₁² − 2x₂)) ≤ 0
c₂(x) = x₁² + x₂² − 3/2 ≤ 0
The paper gives its global minimizer as about (0.1954, 0.4044), where f is about 0.5998, c₁ is
active and c₂ isn't, and two local minimizers, about (0.7197, 0.1411), f ≈ 0.8609, and (0, 0.75)
on the bound. tests/reference/gramacy_toy.py solves c₁ = 0 with the Lagrange condition, which
says the two partial derivatives of c₁ are equal there, with mpmath from the paper's points:
the minimum is 0.5997880520 at (0.1951226835, 0.4046653685), c₂ = −1.298 there, and a scan of the
box on a grid of 1/2000 confirms it's the global one. The paper's coordinates differ from these in
the fourth digit, and its f agrees: along the boundary, f changes slowly.
The objective is known and cheap here, but the method treats it as a black box, as for a simulation whose output and constraints all come from one expensive run.
What makes it hard
The minimum lies on the boundary of c₁, whose sine makes the feasible region's edge a wave: the objective falls toward the infeasible corner (0, 0), so the best points are exactly where the constraint stops them, and a search that ignores the constraint goes the wrong way. Bayesian optimization has only the evaluations to learn where that boundary is.
Representation
A Real genome of 2 genes in [0, 1]. The fitness function is a constraint::Constrained with 2
constraints: it writes c₁ and c₂ into a slice and returns x₁ + x₂, and genoxide makes the fitness
(score, violation) of it, the violation the sum of the positive values, while Bo reads the
values one by one. In Python, the function returns (value, g) and run takes constraints=2.
Both use genoxide's portable sine, so both versions give the same run.
Algorithm
Bo with its defaults. With the constraints' values, it fits a Gaussian process to the objective
and one to each constraint, and maximizes the log expected improvement over the best feasible
point plus the logarithm of the probability that a point is feasible, P(c₁ ≤ 0) P(c₂ ≤ 0) under
the constraints' models: the expected constrained improvement of Gardner et al. (2014), the
product of the expected improvement and the probability of feasibility, through its logarithm.
Until a feasible point is evaluated, the search maximizes the probability of feasibility alone.
The run stops within 1e-5 of the minimum, or after 60 evaluations.
As a contrast, the same search with a fitness function that returns only (score, violation):
without the values, Bo models the score alone, and only the best point found is chosen by Deb's
rules.
Output
The first lines give the problem and the method. Then a row per evaluation: its number, the point, x₁ + x₂, both constraints' values (feasible at 0 or below) and the best feasible value's distance above the minimum so far. Evaluations 1 to 6 are the initial design; the rest are the points the models chose. The last lines give the best feasible point found and how far it is from the minimum, and the contrast's feasible points and best.
The project page plays the run back: at each step, the probability of feasibility over the square with the points so far, and the acquisition that chose the next point.
Good results
A good result is feasible and within 1e-5 of 0.599788. Three points of the design are feasible, the best 0.96; the models then probe the line x₁ = 0, where the objective is lowest, and from evaluation 11 on, they walk the boundary c₁ = 0 to the minimum: 1.3e-2 above it at evaluation 11, 3.9e-4 at 15, 8.4e-6 at 21, (0.194733, 0.405064), with c₁ just below 0. That point is 5.6e-4 from the minimizer: along the boundary, f changes slowly.
Over seeds 1 to 20, every search came within 1e-5, after 18 to 39 evaluations, half of them within 26; within 1e-3, after 13 to 35, half within 19.
The contrast doesn't come close: with only the violation, the search models x₁ + x₂ alone and heads for the infeasible corner, where the function is lowest. In its 21 evaluations, 6 points are feasible, and the best of them is 0.944530, 3.4e-1 above the minimum.
Known optimum: 0.599788 at (0.195123, 0.404665)
Source: examples/bo_constrained
Interactive run: tachsin.gr/projects/genoxide/examples/bo-constrained
cargo run --release --example bo_constrained
//! Constrained Bayesian optimization of the toy problem of Gramacy et al. (2016): a linear
//! objective on [0, 1]² with two constraints whose values the fitness function gives one by one,
//! each modeled by a Gaussian process. The search maximizes the log expected improvement over the
//! best feasible point plus the logarithm of the probability of feasibility, and reaches the
//! global minimum, on the boundary of a wavy constraint, to within 1e-5. Then, as a contrast, the
//! same search told only the total violation, which it ignores.
//!
//! 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 bo_constrained
//! ```
mod trace;
use genoxide::constraint::Constrained;
use genoxide::math::sin;
use genoxide::prelude::*;
use std::cell::RefCell;
use std::f64::consts::PI;
// the global minimum and its point (tests/reference/gramacy_toy.py, with mpmath)
const MINIMUM: f64 = 0.599_788_052_010_067_6;
const MINIMIZER: [f64; 2] = [0.195_122_683_472_071_76, 0.404_665_368_537_995_8];
// how close to the minimum, and the evaluations the search may take at most
const TOLERANCE: f64 = 1e-5;
const BUDGET: u64 = 60;
// x₁ + x₂, and the values of the two constraints g(x) <= 0 into `g`
fn toy(x: &Reals, g: &mut [f64]) -> f64 {
g[0] = 1.5 - x[0] - 2.0 * x[1] - 0.5 * sin(2.0 * PI * (x[0] * x[0] - 2.0 * x[1]));
g[1] = x[0] * x[0] + x[1] * x[1] - 1.5;
x[0] + x[1]
}
fn main() -> Result<()> {
println!("Gramacy et al.'s toy problem: minimize x1 + x2 on [0, 1]^2 subject to");
println!(
" c1 = 1.5 - x1 - 2 x2 - sin(2 pi (x1^2 - 2 x2)) / 2 <= 0, c2 = x1^2 + x2^2 - 1.5 <= 0"
);
println!(
"global minimum {MINIMUM:.6} at ({:.6}, {:.6}), on the boundary c1 = 0",
MINIMIZER[0], MINIMIZER[1]
);
println!("6 points of a Latin hypercube, then a point per step by log-EI x P(feasible)");
println!(
"evaluation x1 x2 f c1 c2 best feasible f - f*"
);
let bo = Bo::builder(Real::uniform(2, 0.0..=1.0)?)
.minimize()
.seed(1)
.build()?;
// with GENOXIDE_TRACE=<file>, a trace of the run for the plot on the example's page
let trace = RefCell::new(trace::Trace::from_env());
let mut printed = 0;
let mut best = f64::INFINITY;
let mut engine = Engine::new(bo, Constrained::new(2, toy))
.stop_when(Stop::target(MINIMUM + TOLERANCE).or(Stop::evaluations(BUDGET)))
.control(|bo, progress| {
// the points evaluated in this generation
for index in printed..bo.population().len() {
let individual = &bo.population().as_slice()[index];
let (x, fitness) = (
individual.genome(),
individual.fitness().expect("evaluated"),
);
let g = bo.constraint_values(index);
if fitness.is_feasible() {
best = best.min(fitness.score().expect("valid"));
}
let gap = if best.is_finite() {
format!("{:.1e}", best - MINIMUM)
} else {
"none yet".to_string()
};
println!(
"{:>10} {:>9.6} {:>9.6} {:>9.6} {:>10.6} {:>10.6} {gap:>20}",
index + 1,
x[0],
x[1],
fitness.score().expect("valid"),
g[0],
g[1]
);
}
printed = bo.population().len();
trace.borrow_mut().record(bo, progress);
Ok(())
});
let outcome = engine.run()?;
drop(engine);
let x = outcome.best_genome();
let value = outcome.best_fitness().score().expect("valid");
assert!(outcome.best_fitness().is_feasible());
let distance = (x[0] - MINIMIZER[0]).hypot(x[1] - MINIMIZER[1]);
println!(
"{} evaluations: the best feasible point ({:.6}, {:.6}), {value:.6}, {:.1e} above the \
minimum, {distance:.1e} from its point",
outcome.evaluations(),
x[0],
x[1],
value - MINIMUM
);
assert!(value - MINIMUM <= TOLERANCE);
trace.into_inner().write();
// the contrast: the same problem as (score, violation), without the constraints' values
let violation = |x: &Reals| {
let mut g = [0.0; 2];
let score = toy(x, &mut g);
(score, g[0].max(0.0) + g[1].max(0.0))
};
let blind = Bo::builder(Real::uniform(2, 0.0..=1.0)?)
.minimize()
.seed(1)
.build()?;
let mut engine =
Engine::new(blind, violation).stop_when(Stop::evaluations(outcome.evaluations()));
let contrast = engine.run()?;
let feasible = engine
.algorithm()
.population()
.iter()
.filter(|individual| individual.fitness().is_some_and(Fitness::is_feasible))
.count();
let fitness = contrast.best_fitness();
println!(
"without the constraints' values: {feasible} feasible points of {}, the best {}",
contrast.evaluations(),
if fitness.is_feasible() {
format!(
"{:.6}, {:.1e} above the minimum",
fitness.score().expect("valid"),
fitness.score().expect("valid") - MINIMUM
)
} else {
format!("infeasible, by {:.1e}", fitness.violation())
}
);
Ok(())
}
python examples/bo_constrained/main.py
"""Constrained Bayesian optimization of the toy problem of Gramacy et al. (2016): a linear objective
on [0, 1]^2 with two constraints whose values the fitness function gives one by one, each modeled
by a Gaussian process. The search maximizes the log expected improvement over the best feasible
point plus the logarithm of the probability of feasibility, and reaches the global minimum, on the
boundary of a wavy constraint, to within 1e-5. Then, as a contrast, the same search told only the
total violation, which it ignores.
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/bo_constrained/main.py
"""
import math
import numpy as np
import genoxide as gx
from trace import Trace
# the global minimum and its point (tests/reference/gramacy_toy.py, with mpmath)
MINIMUM = 0.5997880520100676
MINIMIZER = [0.19512268347207176, 0.4046653685379958]
# how close to the minimum, and the evaluations the search may take at most
TOLERANCE = 1e-5
BUDGET = 60
def scientific(value):
"""Two significant digits, e.g. 1.2e-7."""
mantissa, exponent = f"{value:.1e}".split("e")
return f"{mantissa}e{int(exponent)}"
def constraints(x):
"""The values of the two constraints g(x) <= 0, with genoxide's portable sine."""
wave = gx.math.sin(2.0 * math.pi * (x[0] * x[0] - 2.0 * x[1]))
return np.array([1.5 - x[0] - 2.0 * x[1] - 0.5 * wave, x[0] * x[0] + x[1] * x[1] - 1.5])
def toy(x):
"""x1 + x2, and the constraints' values."""
return x[0] + x[1], constraints(x)
def feasible(g):
return bool(np.all(g <= 0.0))
print("Gramacy et al.'s toy problem: minimize x1 + x2 on [0, 1]^2 subject to")
print(" c1 = 1.5 - x1 - 2 x2 - sin(2 pi (x1^2 - 2 x2)) / 2 <= 0, c2 = x1^2 + x2^2 - 1.5 <= 0")
print(
f"global minimum {MINIMUM:.6f} at ({MINIMIZER[0]:.6f}, {MINIMIZER[1]:.6f}), on the boundary "
"c1 = 0"
)
print("6 points of a Latin hypercube, then a point per step by log-EI x P(feasible)")
print("evaluation x1 x2 f c1 c2 best feasible f - f*")
# with GENOXIDE_TRACE=<file>, a trace of the run for the plot on the example's page
trace = Trace(MINIMIZER, MINIMUM, constraints)
shown = {"printed": 0, "best": math.inf}
def on_generation(progress):
# the points evaluated in this generation
for index in range(shown["printed"], len(progress.population)):
x = progress.population[index]
value, g = toy(x)
if feasible(g):
shown["best"] = min(shown["best"], float(value))
gap = scientific(shown["best"] - MINIMUM) if math.isfinite(shown["best"]) else "none yet"
print(
f"{index + 1:>10} {x[0]:>9.6f} {x[1]:>9.6f} {value:>9.6f} {g[0]:>10.6f} "
f"{g[1]:>10.6f} {gap:>20}"
)
shown["printed"] = len(progress.population)
space = gx.Real((0.0, 1.0), length=2)
bo = gx.Bo(space, objective="minimize", seed=1)
result = bo.run(
toy,
constraints=2,
target=MINIMUM + TOLERANCE,
evaluations=BUDGET,
on_generation=on_generation,
control=trace.record,
)
x = result.best_genome
value = float(x[0] + x[1])
assert feasible(constraints(x))
distance = math.hypot(x[0] - MINIMIZER[0], x[1] - MINIMIZER[1])
print(
f"{result.evaluations} evaluations: the best feasible point ({x[0]:.6f}, {x[1]:.6f}), "
f"{value:.6f}, {scientific(value - MINIMUM)} above the minimum, {scientific(distance)} from "
"its point"
)
assert value - MINIMUM <= TOLERANCE
trace.write()
# the contrast: the same problem as (score, violation), without the constraints' values
def violation(x):
value, g = toy(x)
return value, max(g[0], 0.0) + max(g[1], 0.0)
evaluated = {}
def keep(progress):
evaluated["population"] = progress.population
blind = gx.Bo(space, objective="minimize", seed=1)
contrast = blind.run(violation, evaluations=result.evaluations, on_generation=keep)
count = sum(feasible(constraints(x)) for x in evaluated["population"])
best = contrast.best_genome
g = constraints(best)
if feasible(g):
score = float(best[0] + best[1])
summary = f"{score:.6f}, {scientific(score - MINIMUM)} above the minimum"
else:
summary = f"infeasible, by {scientific(max(g[0], 0.0) + max(g[1], 0.0))}"
print(
f"without the constraints' values: {count} feasible points of {contrast.evaluations}, the "
f"best {summary}"
)
What it prints, from a seeded run:
Gramacy et al.'s toy problem: minimize x1 + x2 on [0, 1]^2 subject to
c1 = 1.5 - x1 - 2 x2 - sin(2 pi (x1^2 - 2 x2)) / 2 <= 0, c2 = x1^2 + x2^2 - 1.5 <= 0
global minimum 0.599788 at (0.195123, 0.404665), on the boundary c1 = 0
6 points of a Latin hypercube, then a point per step by log-EI x P(feasible)
evaluation x1 x2 f c1 c2 best feasible f - f*
1 0.193732 0.765545 0.959277 -0.204590 -0.876409 3.6e-1
2 0.971102 0.545241 1.516343 -0.161850 -0.259672 3.6e-1
3 0.451976 0.061188 0.513164 0.679542 -1.291974 3.6e-1
4 0.541108 0.210487 0.751595 0.898457 -1.162897 3.6e-1
5 0.039182 0.427164 0.466345 0.207200 -1.315996 3.6e-1
6 0.674825 0.892300 1.567125 -0.520084 -0.248411 3.6e-1
7 0.000000 0.508147 0.508147 0.534808 -1.241786 3.6e-1
8 0.000000 0.000000 0.000000 1.500000 -1.500000 3.6e-1
9 0.000000 0.341614 0.341614 0.360133 -1.383300 3.6e-1
10 0.000000 0.767732 0.767732 -0.145956 -0.910588 1.7e-1
11 0.212447 0.400488 0.612935 -0.013086 -1.294476 1.3e-2
12 0.190897 0.387819 0.578716 0.034616 -1.313155 1.3e-2
13 0.205325 0.398267 0.603592 -0.001670 -1.299225 3.8e-3
14 0.191724 0.406987 0.598711 0.001595 -1.297603 3.8e-3
15 0.193398 0.406780 0.600178 -0.000220 -1.297127 3.9e-4
16 0.560359 0.000000 0.560359 0.479528 -1.185998 3.9e-4
17 0.000000 0.201622 0.201622 1.382342 -1.459348 3.9e-4
18 0.172931 0.425271 0.598202 0.024969 -1.289239 3.9e-4
19 0.194758 0.405116 0.599874 -0.000089 -1.297951 8.6e-5
20 0.194735 0.405068 0.599803 -0.000009 -1.297998 1.5e-5
21 0.194733 0.405064 0.599796 -0.000001 -1.298003 8.4e-6
21 evaluations: the best feasible point (0.194733, 0.405064), 0.599796, 8.4e-6 above the minimum, 5.6e-4 from its point
without the constraints' values: 6 feasible points of 21, the best 0.944530, 3.4e-1 above the minimum