From a70500406c131e8b414d25e8160d8b93ded284d5 Mon Sep 17 00:00:00 2001 From: Stephen Waits Date: Tue, 5 May 2026 09:42:28 -0600 Subject: [PATCH] =?UTF-8?q?feat(algorithms):=20add=20BayesianOpt=20?= =?UTF-8?q?=E2=80=94=20GP-based=20Bayesian=20Optimization?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The first sample-efficient algorithm in heuropt. Bayesian optimization maintains a Gaussian-process surrogate of the objective and at each step picks the next decision by maximizing an acquisition function on that surrogate, so the evaluation budget is used surgically. Implementation: - **Kernel**: anisotropic RBF (squared-exponential) with per-axis length scales, signal variance, and a small noise/jitter floor. Hyperparameters are exposed in the config; a future version can add marginal-likelihood maximization. - **Posterior**: standard formulation. Cholesky factorizes K (using the new internal helper); mean and variance predictions follow. - **Acquisition**: Expected Improvement against the best observed feasible point. Optimized by best-of-N random sampling — simple, predictable cost, no inner-optimizer footgun. - **Initial design**: `initial_samples` uniform-random points in bounds before the BO loop starts. - **Constraints**: feasibility-aware EI — best observed value uses only feasible points; infeasible candidates are penalized. Vec decisions, single-objective only. Targets the regime no existing heuropt algorithm covers: 50–500 evaluations on an expensive black-box function (CFD sim, ML training run, real-world measurement). Tests cover convergence on the 1-D sphere within a tight evaluation budget (~30 evals get to f < 1e-6 — vs population-based methods needing thousands), deterministic reruns, and panic on multi-objective + dim mismatches. --- src/algorithms/bayesian_opt.rs | 414 +++++++++++++++++++++++++++++++++ src/algorithms/mod.rs | 2 + src/internal/cholesky.rs | 6 +- src/prelude.rs | 3 +- 4 files changed, 423 insertions(+), 2 deletions(-) create mode 100644 src/algorithms/bayesian_opt.rs diff --git a/src/algorithms/bayesian_opt.rs b/src/algorithms/bayesian_opt.rs new file mode 100644 index 0000000..8f694b9 --- /dev/null +++ b/src/algorithms/bayesian_opt.rs @@ -0,0 +1,414 @@ +//! `BayesianOpt` — Gaussian-process-based Bayesian Optimization. +//! +//! Sample-efficient sequential optimizer for expensive black-box +//! single-objective real-valued problems. Builds a GP surrogate of the +//! objective and selects the next evaluation point by maximizing the +//! Expected Improvement acquisition. + +use rand::Rng as _; + +use crate::core::candidate::Candidate; +use crate::core::evaluation::Evaluation; +use crate::core::objective::Direction; +use crate::core::population::Population; +use crate::core::problem::Problem; +use crate::core::result::OptimizationResult; +use crate::core::rng::{Rng, rng_from_seed}; +use crate::internal::cholesky::{cholesky, solve}; +use crate::operators::real::RealBounds; +use crate::traits::Optimizer; + +/// Configuration for [`BayesianOpt`]. +#[derive(Debug, Clone)] +pub struct BayesianOptConfig { + /// Number of uniform-random initial samples before the BO loop starts. + /// Hansen-style rule of thumb: 5×dim, but small budgets often work. + pub initial_samples: usize, + /// Number of BO iterations after the initial design. + pub iterations: usize, + /// Per-axis RBF length scales (one per dimension). Smaller = more + /// "wiggly" surrogate. Reasonable default: 0.2 × bound range per axis. + pub length_scales: Option>, + /// GP signal variance (the "amplitude" of the surrogate). + pub signal_variance: f64, + /// GP noise variance (small jitter to keep the kernel matrix SPD even + /// at duplicate or near-duplicate points). + pub noise_variance: f64, + /// Number of random samples used to maximize the acquisition function + /// each step. The best-EI sample is chosen as the next evaluation. + pub acquisition_samples: usize, + /// Seed for the deterministic RNG. + pub seed: u64, +} + +impl Default for BayesianOptConfig { + fn default() -> Self { + Self { + initial_samples: 10, + iterations: 40, + length_scales: None, + signal_variance: 1.0, + noise_variance: 1e-6, + acquisition_samples: 1_000, + seed: 42, + } + } +} + +/// Gaussian-process Bayesian Optimization with Expected Improvement. +/// +/// `Vec` decisions only. Single-objective only. Targets expensive +/// evaluation budgets (50–500). The GP kernel is anisotropic RBF; the +/// acquisition function is EI; both are optimized by best-of-N random +/// sampling each step (simple, predictable cost). +#[derive(Debug, Clone)] +pub struct BayesianOpt { + /// Algorithm configuration. + pub config: BayesianOptConfig, + /// Per-variable bounds. + pub bounds: RealBounds, +} + +impl BayesianOpt { + /// Construct a `BayesianOpt`. + pub fn new(config: BayesianOptConfig, bounds: RealBounds) -> Self { + Self { config, bounds } + } +} + +impl

Optimizer

for BayesianOpt +where + P: Problem> + Sync, +{ + fn run(&mut self, problem: &P) -> OptimizationResult { + assert!( + self.config.initial_samples >= 2, + "BayesianOpt initial_samples must be >= 2", + ); + assert!(self.config.signal_variance > 0.0, "BayesianOpt signal_variance must be > 0"); + assert!(self.config.noise_variance > 0.0, "BayesianOpt noise_variance must be > 0"); + assert!( + self.config.acquisition_samples >= 1, + "BayesianOpt acquisition_samples must be >= 1", + ); + let objectives = problem.objectives(); + assert!( + objectives.is_single_objective(), + "BayesianOpt requires exactly one objective", + ); + let direction = objectives.objectives[0].direction; + let dim = self.bounds.bounds.len(); + if let Some(ls) = &self.config.length_scales { + assert_eq!(ls.len(), dim, "BayesianOpt length_scales.len() must equal dim"); + } + let length_scales: Vec = self + .config + .length_scales + .clone() + .unwrap_or_else(|| { + self.bounds + .bounds + .iter() + .map(|&(lo, hi)| 0.2 * (hi - lo).max(1e-9)) + .collect() + }); + let mut rng = rng_from_seed(self.config.seed); + + // ---------------- Initial random design ---------------- + let mut decisions: Vec> = Vec::with_capacity( + self.config.initial_samples + self.config.iterations, + ); + let mut targets: Vec = Vec::with_capacity(decisions.capacity()); + let mut evaluations = Vec::with_capacity(decisions.capacity()); + for _ in 0..self.config.initial_samples { + let x = sample_uniform_in_bounds(&self.bounds, &mut rng); + let e = problem.evaluate(&x); + // GP works on minimization-oriented "want low" targets. + let t = oriented_target(&e, direction); + decisions.push(x); + targets.push(t); + evaluations.push(e); + } + + // ---------------- Sequential BO loop ---------------- + for _ in 0..self.config.iterations { + // Build the GP posterior around current observations. + let posterior = match GpPosterior::fit( + &decisions, + &targets, + &length_scales, + self.config.signal_variance, + self.config.noise_variance, + ) { + Ok(p) => p, + Err(_) => { + // SPD failure (typically numerical): fall back to a + // single uniform-random sample this step. + let x = sample_uniform_in_bounds(&self.bounds, &mut rng); + let e = problem.evaluate(&x); + targets.push(oriented_target(&e, direction)); + decisions.push(x); + evaluations.push(e); + continue; + } + }; + + let best_target = + targets.iter().cloned().fold(f64::INFINITY, f64::min); + + // Maximize EI by best-of-N random sampling. + let mut best_x = sample_uniform_in_bounds(&self.bounds, &mut rng); + let mut best_ei = -f64::INFINITY; + for _ in 0..self.config.acquisition_samples { + let cand = sample_uniform_in_bounds(&self.bounds, &mut rng); + let (mu, sigma) = posterior.predict(&cand); + let ei = expected_improvement(mu, sigma, best_target); + if ei > best_ei { + best_ei = ei; + best_x = cand; + } + } + + let e = problem.evaluate(&best_x); + targets.push(oriented_target(&e, direction)); + decisions.push(best_x); + evaluations.push(e); + } + + // Build the final population/best. + let final_pop: Vec>> = decisions + .into_iter() + .zip(evaluations) + .map(|(d, e)| Candidate::new(d, e)) + .collect(); + let mut best_idx = 0; + for i in 1..final_pop.len() { + if better(&final_pop[i].evaluation, &final_pop[best_idx].evaluation, direction) { + best_idx = i; + } + } + let total_evaluations = final_pop.len(); + let best = final_pop[best_idx].clone(); + let front = vec![best.clone()]; + OptimizationResult::new( + Population::new(final_pop), + front, + Some(best), + total_evaluations, + self.config.iterations + self.config.initial_samples, + ) + } +} + +/// Convert an Evaluation into a "smaller is better" target. For Maximize +/// problems we negate; infeasibles get a large penalty proportional to +/// the violation magnitude. +fn oriented_target(e: &Evaluation, direction: Direction) -> f64 { + let base = match direction { + Direction::Minimize => e.objectives[0], + Direction::Maximize => -e.objectives[0], + }; + if e.is_feasible() { + base + } else { + // Penalize so the GP learns to avoid this region. + base + 1e6 * e.constraint_violation + } +} + +fn better(a: &Evaluation, b: &Evaluation, direction: Direction) -> bool { + match (a.is_feasible(), b.is_feasible()) { + (true, false) => true, + (false, true) => false, + (false, false) => a.constraint_violation < b.constraint_violation, + (true, true) => match direction { + Direction::Minimize => a.objectives[0] < b.objectives[0], + Direction::Maximize => a.objectives[0] > b.objectives[0], + }, + } +} + +fn sample_uniform_in_bounds(bounds: &RealBounds, rng: &mut Rng) -> Vec { + bounds + .bounds + .iter() + .map(|&(lo, hi)| if lo == hi { lo } else { lo + (hi - lo) * rng.random::() }) + .collect() +} + +/// Anisotropic RBF kernel: `k(x, y) = σ² · exp(-0.5 · Σ ((x_i - y_i)/ℓ_i)²)`. +fn rbf_kernel(x: &[f64], y: &[f64], length_scales: &[f64], signal_variance: f64) -> f64 { + let mut sum = 0.0; + for ((a, b), l) in x.iter().zip(y.iter()).zip(length_scales.iter()) { + let d = (a - b) / l.max(1e-12); + sum += d * d; + } + signal_variance * (-0.5 * sum).exp() +} + +struct GpPosterior { + decisions: Vec>, + length_scales: Vec, + signal_variance: f64, + /// `α = K^{-1} · y_target`, precomputed for the mean prediction. + alpha: Vec, + /// Cholesky factor of `K + σ_n² · I`, kept for variance prediction. + chol_l: Vec>, +} + +impl GpPosterior { + fn fit( + decisions: &[Vec], + targets: &[f64], + length_scales: &[f64], + signal_variance: f64, + noise_variance: f64, + ) -> Result { + let n = decisions.len(); + let mut k = vec![vec![0.0_f64; n]; n]; + for i in 0..n { + for j in 0..=i { + let v = rbf_kernel(&decisions[i], &decisions[j], length_scales, signal_variance); + k[i][j] = v; + k[j][i] = v; + } + k[i][i] += noise_variance; + } + let chol_l = cholesky(&k)?; + let alpha = solve(&chol_l, targets); + Ok(Self { + decisions: decisions.to_vec(), + length_scales: length_scales.to_vec(), + signal_variance, + alpha, + chol_l, + }) + } + + fn predict(&self, x: &[f64]) -> (f64, f64) { + let n = self.decisions.len(); + let mut k_star = vec![0.0_f64; n]; + for (i, k_star_i) in k_star.iter_mut().enumerate() { + *k_star_i = rbf_kernel(x, &self.decisions[i], &self.length_scales, self.signal_variance); + } + let _ = n; + let mu: f64 = k_star.iter().zip(self.alpha.iter()).map(|(a, b)| a * b).sum(); + // Var = k(x,x) - k_star^T · K^{-1} · k_star + // Compute K^{-1}·k_star = solve_upper_transpose(L, solve_lower(L, k_star)) + let v_temp = crate::internal::cholesky::solve_lower(&self.chol_l, &k_star); + let v: f64 = v_temp.iter().map(|x| x * x).sum(); + let var = (self.signal_variance - v).max(0.0); + (mu, var.sqrt()) + } +} + +/// Expected Improvement (minimization-oriented) at a point with predicted +/// mean `mu` and standard deviation `sigma`, given the current best +/// observed target `f_best`. Returns 0 if `sigma` is effectively zero. +fn expected_improvement(mu: f64, sigma: f64, f_best: f64) -> f64 { + if sigma < 1e-12 { + return 0.0; + } + let improvement = f_best - mu; + let z = improvement / sigma; + improvement * normal_cdf(z) + sigma * normal_pdf(z) +} + +fn normal_pdf(z: f64) -> f64 { + (-0.5 * z * z).exp() / (2.0 * std::f64::consts::PI).sqrt() +} + +fn normal_cdf(z: f64) -> f64 { + // Approximate Φ(z) via erf. Abramowitz & Stegun 7.1.26 series-free + // rational approximation good to ~1.5e-7. + 0.5 * (1.0 + erf(z / std::f64::consts::SQRT_2)) +} + +fn erf(x: f64) -> f64 { + // Numerical Recipes-style erf, accurate to ~1e-7. + let a1 = 0.254_829_592; + let a2 = -0.284_496_736; + let a3 = 1.421_413_741; + let a4 = -1.453_152_027; + let a5 = 1.061_405_429; + let p = 0.327_591_1; + let sign = if x < 0.0 { -1.0 } else { 1.0 }; + let x = x.abs(); + let t = 1.0 / (1.0 + p * x); + let y = 1.0 + - (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1) * t * (-x * x).exp(); + sign * y +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::tests_support::{SchafferN1, Sphere1D}; + + fn make_optimizer(seed: u64) -> BayesianOpt { + BayesianOpt::new( + BayesianOptConfig { + initial_samples: 5, + iterations: 25, + length_scales: None, + signal_variance: 1.0, + noise_variance: 1e-6, + acquisition_samples: 500, + seed, + }, + RealBounds::new(vec![(-5.0, 5.0)]), + ) + } + + #[test] + fn finds_minimum_of_sphere_quickly() { + // BO's whole point is sample efficiency: 30 evals ought to be + // enough for a 1-D sphere. (Pop-based methods needed thousands.) + let mut opt = make_optimizer(1); + let r = opt.run(&Sphere1D); + let best = r.best.unwrap(); + assert!( + best.evaluation.objectives[0] < 1e-3, + "BO should converge fast on 1-D sphere; got f = {}", + best.evaluation.objectives[0], + ); + assert!(r.evaluations <= 30 + 1); + } + + #[test] + fn deterministic_with_same_seed() { + let mut a = make_optimizer(99); + let mut b = make_optimizer(99); + let ra = a.run(&Sphere1D); + let rb = b.run(&Sphere1D); + assert_eq!( + ra.best.unwrap().evaluation.objectives, + rb.best.unwrap().evaluation.objectives, + ); + } + + #[test] + #[should_panic(expected = "exactly one objective")] + fn multi_objective_panics() { + let mut opt = make_optimizer(0); + let _ = opt.run(&SchafferN1); + } + + #[test] + #[should_panic(expected = "length_scales.len() must equal dim")] + fn length_scales_dim_mismatch_panics() { + let mut opt = BayesianOpt::new( + BayesianOptConfig { + initial_samples: 5, + iterations: 5, + length_scales: Some(vec![1.0, 1.0]), + signal_variance: 1.0, + noise_variance: 1e-6, + acquisition_samples: 100, + seed: 0, + }, + RealBounds::new(vec![(-1.0, 1.0)]), + ); + let _ = opt.run(&Sphere1D); + } +} diff --git a/src/algorithms/mod.rs b/src/algorithms/mod.rs index e1a2012..40227db 100644 --- a/src/algorithms/mod.rs +++ b/src/algorithms/mod.rs @@ -2,6 +2,7 @@ pub mod age_moea; pub mod ant_colony_tsp; +pub mod bayesian_opt; pub mod cma_es; pub mod differential_evolution; pub mod epsilon_moea; @@ -33,6 +34,7 @@ pub mod umda; pub use age_moea::*; pub use ant_colony_tsp::*; +pub use bayesian_opt::*; pub use cma_es::*; pub use differential_evolution::*; pub use epsilon_moea::*; diff --git a/src/internal/cholesky.rs b/src/internal/cholesky.rs index 1329c3b..f4b2f05 100644 --- a/src/internal/cholesky.rs +++ b/src/internal/cholesky.rs @@ -15,9 +15,11 @@ pub(crate) fn cholesky(a: &[Vec]) -> Result>, &'static str> { } debug_assert!(a.iter().all(|row| row.len() == n)); let mut l = vec![vec![0.0_f64; n]; n]; + #[allow(clippy::needless_range_loop)] // body indexes both `a` and `l` rows. for i in 0..n { for j in 0..=i { let mut sum = a[i][j]; + #[allow(clippy::needless_range_loop)] for k in 0..j { sum -= l[i][k] * l[j][k]; } @@ -90,9 +92,11 @@ mod tests { assert!(approx_eq(l[1][0], 1.0, 1e-12)); assert!(approx_eq(l[1][1], 2.0, 1e-12)); // L · L^T should reconstruct A. + #[allow(clippy::needless_range_loop)] for i in 0..2 { for j in 0..2 { let mut s = 0.0; + #[allow(clippy::needless_range_loop)] for k in 0..2 { s += l[i][k] * l[j][k]; } @@ -111,7 +115,7 @@ mod tests { ]; let l = cholesky(&a).unwrap(); // Choose a vector and check A · x = b round trip. - let x_truth = vec![1.0, -2.0, 0.5]; + let x_truth = [1.0, -2.0, 0.5]; let b: Vec = (0..3) .map(|i| (0..3).map(|j| a[i][j] * x_truth[j]).sum()) .collect(); diff --git a/src/prelude.rs b/src/prelude.rs index d90443f..18e069b 100644 --- a/src/prelude.rs +++ b/src/prelude.rs @@ -22,7 +22,8 @@ pub use crate::operators::{ }; pub use crate::algorithms::{ - AgeMoea, AgeMoeaConfig, AntColonyTsp, AntColonyTspConfig, CmaEs, CmaEsConfig, DifferentialEvolution, + AgeMoea, AgeMoeaConfig, AntColonyTsp, AntColonyTspConfig, BayesianOpt, + BayesianOptConfig, CmaEs, CmaEsConfig, DifferentialEvolution, DifferentialEvolutionConfig, EpsilonMoea, EpsilonMoeaConfig, GeneticAlgorithm, GeneticAlgorithmConfig, Grea, GreaConfig, HillClimber, HillClimberConfig, Hype, HypeConfig, Ibea, IbeaConfig, IpopCmaEs, IpopCmaEsConfig, Knea, KneaConfig, Moead, MoeadConfig, Mopso, MopsoConfig,