swaits
fa3f2e8fb0
Theme: documentation and project polish. No public-API changes; this is the v0.5 release that elevates heuropt's docs/onboarding/governance to bar-setting status. Adds: - mdbook user guide at docs/book/ with intro, getting-started, defining-problems, choosing-an-algorithm, cookbook (7 recipes), comparison vs other libraries, stability/SemVer, migration guides. Deploys to https://swaits.github.io/heuropt/ via .github/workflows/ docs.yml. - Runnable rustdoc examples on every algorithm (35 of them), all exercised by cargo test --doc. - Three real-world examples: portfolio.rs (multi-obj with budget constraint), hyperparam_tuning.rs (BO + TPE), scheduling.rs (permutation via SA + SwapMutation against Smith's-rule oracle). - Governance: CONTRIBUTING.md, SECURITY.md, CODE_OF_CONDUCT.md (adopting builderscode.org's Builder's Code of Conduct), GitHub issue templates, PR template. Polishes: - README hero with badges + user-guide link. - lib.rs crate-level docs. - CHANGELOG entry for 0.5.0. Bumps Cargo.toml to 0.5.0.
518 lines
19 KiB
Rust
518 lines
19 KiB
Rust
//! CMA-ES — Hansen & Ostermeier 2001 Covariance Matrix Adaptation
|
||
//! Evolution Strategy.
|
||
|
||
use rand_distr::{Distribution, Normal};
|
||
|
||
use crate::algorithms::parallel_eval::evaluate_batch;
|
||
use crate::core::candidate::Candidate;
|
||
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_from_seed;
|
||
use crate::internal::eigen::symmetric_eigen;
|
||
use crate::operators::real::RealBounds;
|
||
use crate::pareto::front::best_candidate;
|
||
use crate::traits::Optimizer;
|
||
|
||
/// Configuration for [`CmaEs`].
|
||
#[derive(Debug, Clone)]
|
||
pub struct CmaEsConfig {
|
||
/// Population size `λ`. Must be at least 4. Hansen recommends
|
||
/// `4 + floor(3 · ln(N))` as a default for `N`-dim problems.
|
||
pub population_size: usize,
|
||
/// Number of generations.
|
||
pub generations: usize,
|
||
/// Initial step size `σ_0`. Often ~ 1/3 of the search range per dim.
|
||
pub initial_sigma: f64,
|
||
/// Recompute the eigendecomposition of `C` every this many generations
|
||
/// to amortize cost. The full algorithm decomposes every generation
|
||
/// (set this to 1); 1–10 is fine for small `N`.
|
||
pub eigen_decomposition_period: usize,
|
||
/// Optional initial mean. If `None`, the mean defaults to the per-axis
|
||
/// midpoint of the bounds. Used by `IpopCmaEs` to inject restart
|
||
/// diversity without shrinking the search box.
|
||
pub initial_mean: Option<Vec<f64>>,
|
||
/// Seed for the deterministic RNG.
|
||
pub seed: u64,
|
||
}
|
||
|
||
impl Default for CmaEsConfig {
|
||
fn default() -> Self {
|
||
Self {
|
||
population_size: 16,
|
||
generations: 200,
|
||
initial_sigma: 0.5,
|
||
eigen_decomposition_period: 1,
|
||
initial_mean: None,
|
||
seed: 42,
|
||
}
|
||
}
|
||
}
|
||
|
||
/// Single-objective real-valued CMA-ES.
|
||
///
|
||
/// Maintains a multivariate Gaussian sampler `mean + σ · N(0, C)`, samples
|
||
/// `λ` offspring from it each generation, selects the `μ` best (weighted),
|
||
/// and updates `mean`, `σ`, and `C` via the standard CMA-ES rules.
|
||
///
|
||
/// `Vec<f64>` decisions only. Bounds come from the embedded `RealBounds`
|
||
/// field; both the initial mean and every offspring are clamped per
|
||
/// dimension.
|
||
///
|
||
/// # Example
|
||
///
|
||
/// ```
|
||
/// use heuropt::prelude::*;
|
||
///
|
||
/// struct Sphere;
|
||
/// impl Problem for Sphere {
|
||
/// type Decision = Vec<f64>;
|
||
/// fn objectives(&self) -> ObjectiveSpace {
|
||
/// ObjectiveSpace::new(vec![Objective::minimize("f")])
|
||
/// }
|
||
/// fn evaluate(&self, x: &Vec<f64>) -> Evaluation {
|
||
/// Evaluation::new(vec![x.iter().map(|v| v * v).sum::<f64>()])
|
||
/// }
|
||
/// }
|
||
///
|
||
/// let mut opt = CmaEs::new(
|
||
/// CmaEsConfig {
|
||
/// population_size: 12,
|
||
/// generations: 100,
|
||
/// initial_sigma: 1.0,
|
||
/// eigen_decomposition_period: 1,
|
||
/// initial_mean: None,
|
||
/// seed: 42,
|
||
/// },
|
||
/// RealBounds::new(vec![(-5.0, 5.0); 5]),
|
||
/// );
|
||
/// let r = opt.run(&Sphere);
|
||
/// // CMA-ES converges aggressively on Sphere.
|
||
/// assert!(r.best.unwrap().evaluation.objectives[0] < 1e-3);
|
||
/// ```
|
||
#[derive(Debug, Clone)]
|
||
pub struct CmaEs {
|
||
/// Algorithm configuration.
|
||
pub config: CmaEsConfig,
|
||
/// Per-variable bounds — used both to seed `mean` (midpoint) and to
|
||
/// clamp every offspring.
|
||
pub bounds: RealBounds,
|
||
}
|
||
|
||
impl CmaEs {
|
||
/// Construct a `CmaEs`.
|
||
pub fn new(config: CmaEsConfig, bounds: RealBounds) -> Self {
|
||
Self { config, bounds }
|
||
}
|
||
}
|
||
|
||
impl<P> Optimizer<P> for CmaEs
|
||
where
|
||
P: Problem<Decision = Vec<f64>> + Sync,
|
||
{
|
||
fn run(&mut self, problem: &P) -> OptimizationResult<P::Decision> {
|
||
assert!(
|
||
self.config.population_size >= 4,
|
||
"CmaEs population_size must be >= 4",
|
||
);
|
||
assert!(
|
||
self.config.initial_sigma > 0.0,
|
||
"CmaEs initial_sigma must be positive",
|
||
);
|
||
assert!(
|
||
self.config.eigen_decomposition_period >= 1,
|
||
"CmaEs eigen_decomposition_period must be >= 1",
|
||
);
|
||
let objectives = problem.objectives();
|
||
assert!(
|
||
objectives.is_single_objective(),
|
||
"CmaEs only supports single-objective problems",
|
||
);
|
||
let direction = objectives.objectives[0].direction;
|
||
|
||
let n = self.bounds.bounds.len();
|
||
let n_f = n as f64;
|
||
let lambda = self.config.population_size;
|
||
let lambda_f = lambda as f64;
|
||
let mu = lambda / 2;
|
||
assert!(mu >= 1, "CmaEs derived mu (= lambda/2) must be >= 1");
|
||
let mut rng = rng_from_seed(self.config.seed);
|
||
|
||
// ---------------------------------------------------------------
|
||
// Selection weights w_i ∝ ln((λ+1)/2) − ln(i) for i = 1..μ,
|
||
// normalized so they sum to 1. Then mu_eff = 1 / Σ w_i².
|
||
// ---------------------------------------------------------------
|
||
let raw_weights: Vec<f64> = (0..mu)
|
||
.map(|i| ((lambda_f + 1.0) / 2.0).ln() - ((i + 1) as f64).ln())
|
||
.collect();
|
||
let sum_w: f64 = raw_weights.iter().sum();
|
||
let weights: Vec<f64> = raw_weights.iter().map(|w| w / sum_w).collect();
|
||
let mu_eff = 1.0 / weights.iter().map(|w| w * w).sum::<f64>();
|
||
|
||
// ---------------------------------------------------------------
|
||
// Standard CMA-ES strategy parameters (Hansen tutorial §7.1).
|
||
// ---------------------------------------------------------------
|
||
let c_sigma = (mu_eff + 2.0) / (n_f + mu_eff + 5.0);
|
||
let d_sigma = 1.0 + 2.0 * ((mu_eff - 1.0) / (n_f + 1.0)).sqrt().max(0.0) + c_sigma;
|
||
let c_c = (4.0 + mu_eff / n_f) / (n_f + 4.0 + 2.0 * mu_eff / n_f);
|
||
let c_1 = 2.0 / ((n_f + 1.3).powi(2) + mu_eff);
|
||
let c_mu = ((1.0 - c_1) * 2.0 * (mu_eff - 2.0 + 1.0 / mu_eff)
|
||
/ ((n_f + 2.0).powi(2) + mu_eff))
|
||
.min(1.0 - c_1);
|
||
// E‖N(0, I)‖ ≈ √n · (1 − 1/(4n) + 1/(21n²))
|
||
let chi_n = n_f.sqrt() * (1.0 - 1.0 / (4.0 * n_f) + 1.0 / (21.0 * n_f * n_f));
|
||
|
||
// ---------------------------------------------------------------
|
||
// Initial state.
|
||
// ---------------------------------------------------------------
|
||
let mut mean: Vec<f64> = if let Some(provided) = self.config.initial_mean.clone() {
|
||
assert_eq!(
|
||
provided.len(),
|
||
self.bounds.bounds.len(),
|
||
"CmaEs initial_mean.len() must equal the bounds dimension",
|
||
);
|
||
// Clamp the user-provided mean into the bounds so the algorithm
|
||
// doesn't start outside the search box.
|
||
provided
|
||
.into_iter()
|
||
.zip(self.bounds.bounds.iter())
|
||
.map(|(v, &(lo, hi))| v.clamp(lo, hi))
|
||
.collect()
|
||
} else {
|
||
self.bounds
|
||
.bounds
|
||
.iter()
|
||
.map(|&(lo, hi)| 0.5 * (lo + hi))
|
||
.collect()
|
||
};
|
||
let mut sigma = self.config.initial_sigma;
|
||
// Covariance C, eigenvectors B, eigenvalues d (square roots of eigenvalues of C).
|
||
let mut c_matrix: Vec<Vec<f64>> = (0..n)
|
||
.map(|i| (0..n).map(|j| if i == j { 1.0 } else { 0.0 }).collect())
|
||
.collect();
|
||
let mut b: Vec<Vec<f64>> = c_matrix.to_vec();
|
||
let mut d: Vec<f64> = vec![1.0; n];
|
||
let mut p_sigma = vec![0.0_f64; n];
|
||
let mut p_c = vec![0.0_f64; n];
|
||
let mut evaluations = 0usize;
|
||
|
||
let normal = Normal::new(0.0, 1.0).expect("Normal::new(0, 1)");
|
||
let mut best_candidate_seen: Option<Candidate<Vec<f64>>> = None;
|
||
|
||
for generation in 0..self.config.generations {
|
||
// Recompute B, d every period generations from C (after symmetrizing).
|
||
if generation % self.config.eigen_decomposition_period == 0 {
|
||
// Force symmetry.
|
||
#[allow(clippy::needless_range_loop)] // body indexes both [i][j] and [j][i].
|
||
for i in 0..n {
|
||
for j in (i + 1)..n {
|
||
let avg = 0.5 * (c_matrix[i][j] + c_matrix[j][i]);
|
||
c_matrix[i][j] = avg;
|
||
c_matrix[j][i] = avg;
|
||
}
|
||
}
|
||
let (eigenvalues, eigenvectors) = symmetric_eigen(&c_matrix, 1e-14, 100);
|
||
// eigenvectors is sorted descending; we don't depend on order
|
||
// for sampling correctness, but we do need positive eigenvalues.
|
||
d = eigenvalues.iter().map(|&v| v.max(1e-20).sqrt()).collect();
|
||
// B is the matrix whose columns are the eigenvectors. The
|
||
// helper returns `eigenvectors[i]` as the i-th *eigenvector*,
|
||
// so b[r][c] should equal eigenvectors[c][r].
|
||
b = (0..n)
|
||
.map(|r| (0..n).map(|c| eigenvectors[c][r]).collect())
|
||
.collect();
|
||
}
|
||
|
||
// ----- Sample λ offspring -----
|
||
let mut z_samples: Vec<Vec<f64>> = Vec::with_capacity(lambda);
|
||
let mut x_samples: Vec<Vec<f64>> = Vec::with_capacity(lambda);
|
||
for _ in 0..lambda {
|
||
let z: Vec<f64> = (0..n).map(|_| normal.sample(&mut rng)).collect();
|
||
// y = B · D · z
|
||
let bd_z: Vec<f64> = (0..n)
|
||
.map(|i| (0..n).map(|j| b[i][j] * d[j] * z[j]).sum::<f64>())
|
||
.collect();
|
||
// x = mean + σ · y, clamped to bounds
|
||
let x: Vec<f64> = (0..n)
|
||
.map(|i| {
|
||
let v = mean[i] + sigma * bd_z[i];
|
||
let (lo, hi) = self.bounds.bounds[i];
|
||
v.clamp(lo, hi)
|
||
})
|
||
.collect();
|
||
z_samples.push(z);
|
||
x_samples.push(x);
|
||
}
|
||
|
||
// Evaluate offspring (parallel-friendly).
|
||
let evaluated = evaluate_batch(problem, x_samples.clone());
|
||
evaluations += evaluated.len();
|
||
|
||
// Track the best candidate ever.
|
||
for c in &evaluated {
|
||
let beats_best = match &best_candidate_seen {
|
||
None => true,
|
||
Some(b) => better_than_so(&c.evaluation, &b.evaluation, direction),
|
||
};
|
||
if beats_best {
|
||
best_candidate_seen = Some(c.clone());
|
||
}
|
||
}
|
||
|
||
// Sort offspring by fitness ascending (best first).
|
||
let mut order: Vec<usize> = (0..lambda).collect();
|
||
order.sort_by(|&a, &b_| {
|
||
compare_so(
|
||
&evaluated[a].evaluation,
|
||
&evaluated[b_].evaluation,
|
||
direction,
|
||
)
|
||
});
|
||
|
||
// ----- Recompute mean from the μ best (weighted average of x) -----
|
||
let old_mean = mean.clone();
|
||
let mut new_mean = vec![0.0_f64; n];
|
||
for k in 0..mu {
|
||
let xk = &x_samples[order[k]];
|
||
let wk = weights[k];
|
||
for i in 0..n {
|
||
new_mean[i] += wk * xk[i];
|
||
}
|
||
}
|
||
mean = new_mean;
|
||
|
||
// ----- Weighted average of z (used for evolution-path updates) -----
|
||
let mut z_weighted = vec![0.0_f64; n];
|
||
for k in 0..mu {
|
||
let zk = &z_samples[order[k]];
|
||
let wk = weights[k];
|
||
for i in 0..n {
|
||
z_weighted[i] += wk * zk[i];
|
||
}
|
||
}
|
||
|
||
// ----- Evolution path for step size: p_σ = (1 - c_σ) p_σ + sqrt(c_σ (2 - c_σ) μ_eff) · B z̄ -----
|
||
let factor_p_sigma = (c_sigma * (2.0 - c_sigma) * mu_eff).sqrt();
|
||
// B · z_weighted (since C^{-1/2} (m_new - m_old) / σ = B · D^{-1} · D · z̄ = B · z̄)
|
||
let bz: Vec<f64> = (0..n)
|
||
.map(|i| (0..n).map(|j| b[i][j] * z_weighted[j]).sum::<f64>())
|
||
.collect();
|
||
for i in 0..n {
|
||
p_sigma[i] = (1.0 - c_sigma) * p_sigma[i] + factor_p_sigma * bz[i];
|
||
}
|
||
|
||
// ----- Step-size update -----
|
||
let p_sigma_norm = p_sigma.iter().map(|x| x * x).sum::<f64>().sqrt();
|
||
sigma *= ((c_sigma / d_sigma) * (p_sigma_norm / chi_n - 1.0)).exp();
|
||
|
||
// Heaviside for h_σ: damp p_c update if the step length is huge.
|
||
let h_sigma = if p_sigma_norm
|
||
/ (1.0 - (1.0 - c_sigma).powi(2 * (generation as i32 + 1))).sqrt()
|
||
< (1.4 + 2.0 / (n_f + 1.0)) * chi_n
|
||
{
|
||
1.0
|
||
} else {
|
||
0.0
|
||
};
|
||
|
||
// ----- Evolution path for C: p_c = (1 - c_c) p_c + h_σ · sqrt(c_c (2 - c_c) μ_eff) · (m_new - m_old)/σ -----
|
||
let factor_p_c = h_sigma * (c_c * (2.0 - c_c) * mu_eff).sqrt();
|
||
for i in 0..n {
|
||
p_c[i] = (1.0 - c_c) * p_c[i] + factor_p_c * (mean[i] - old_mean[i]) / sigma;
|
||
}
|
||
|
||
// ----- Covariance matrix update (rank-1 + rank-μ) -----
|
||
let delta_h = (1.0 - h_sigma) * c_c * (2.0 - c_c);
|
||
#[allow(clippy::needless_range_loop)]
|
||
// body uses both i and j to index c_matrix and offspring.
|
||
for i in 0..n {
|
||
for j in 0..n {
|
||
let mut update = (1.0 - c_1 - c_mu) * c_matrix[i][j]
|
||
+ c_1 * (p_c[i] * p_c[j] + delta_h * c_matrix[i][j]);
|
||
// Rank-μ contribution.
|
||
let mut rank_mu_term = 0.0;
|
||
for k in 0..mu {
|
||
let xk = &x_samples[order[k]];
|
||
let yi = (xk[i] - old_mean[i]) / sigma;
|
||
let yj = (xk[j] - old_mean[j]) / sigma;
|
||
rank_mu_term += weights[k] * yi * yj;
|
||
}
|
||
update += c_mu * rank_mu_term;
|
||
c_matrix[i][j] = update;
|
||
}
|
||
}
|
||
|
||
// Clamp mean to bounds (sigma may push it out otherwise).
|
||
for (i, m) in mean.iter_mut().enumerate() {
|
||
let (lo, hi) = self.bounds.bounds[i];
|
||
*m = m.clamp(lo, hi);
|
||
}
|
||
}
|
||
|
||
// Final population: just the best-seen candidate. Match other
|
||
// single-objective algorithms' convention.
|
||
let best = best_candidate_seen.expect("at least one generation evaluated");
|
||
let final_pop = vec![best.clone()];
|
||
let front = vec![best.clone()];
|
||
let best_opt = best_candidate(&final_pop, &objectives);
|
||
OptimizationResult::new(
|
||
Population::new(final_pop),
|
||
front,
|
||
best_opt,
|
||
evaluations,
|
||
self.config.generations,
|
||
)
|
||
}
|
||
}
|
||
|
||
fn compare_so(
|
||
a: &crate::core::evaluation::Evaluation,
|
||
b: &crate::core::evaluation::Evaluation,
|
||
direction: Direction,
|
||
) -> std::cmp::Ordering {
|
||
match (a.is_feasible(), b.is_feasible()) {
|
||
(true, false) => std::cmp::Ordering::Less,
|
||
(false, true) => std::cmp::Ordering::Greater,
|
||
(false, false) => a
|
||
.constraint_violation
|
||
.partial_cmp(&b.constraint_violation)
|
||
.unwrap_or(std::cmp::Ordering::Equal),
|
||
(true, true) => match direction {
|
||
Direction::Minimize => a.objectives[0]
|
||
.partial_cmp(&b.objectives[0])
|
||
.unwrap_or(std::cmp::Ordering::Equal),
|
||
Direction::Maximize => b.objectives[0]
|
||
.partial_cmp(&a.objectives[0])
|
||
.unwrap_or(std::cmp::Ordering::Equal),
|
||
},
|
||
}
|
||
}
|
||
|
||
fn better_than_so(
|
||
a: &crate::core::evaluation::Evaluation,
|
||
b: &crate::core::evaluation::Evaluation,
|
||
direction: Direction,
|
||
) -> bool {
|
||
compare_so(a, b, direction) == std::cmp::Ordering::Less
|
||
}
|
||
|
||
#[cfg(test)]
|
||
mod tests {
|
||
use super::*;
|
||
use crate::core::evaluation::Evaluation;
|
||
use crate::core::objective::{Objective, ObjectiveSpace};
|
||
use crate::tests_support::{SchafferN1, Sphere1D};
|
||
|
||
/// 5-D Rosenbrock for exercise.
|
||
struct Rosenbrock5D;
|
||
impl Problem for Rosenbrock5D {
|
||
type Decision = Vec<f64>;
|
||
|
||
fn objectives(&self) -> ObjectiveSpace {
|
||
ObjectiveSpace::new(vec![Objective::minimize("f")])
|
||
}
|
||
|
||
fn evaluate(&self, x: &Vec<f64>) -> Evaluation {
|
||
let f: f64 = (0..(x.len() - 1))
|
||
.map(|i| {
|
||
let a = 1.0 - x[i];
|
||
let b = x[i + 1] - x[i] * x[i];
|
||
a * a + 100.0 * b * b
|
||
})
|
||
.sum();
|
||
Evaluation::new(vec![f])
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
fn finds_minimum_of_sphere() {
|
||
let mut opt = CmaEs::new(
|
||
CmaEsConfig {
|
||
population_size: 12,
|
||
generations: 100,
|
||
initial_sigma: 0.5,
|
||
eigen_decomposition_period: 1,
|
||
initial_mean: None,
|
||
seed: 1,
|
||
},
|
||
RealBounds::new(vec![(-5.0, 5.0)]),
|
||
);
|
||
let r = opt.run(&Sphere1D);
|
||
let best = r.best.unwrap();
|
||
assert!(
|
||
best.evaluation.objectives[0] < 1e-8,
|
||
"got f = {}",
|
||
best.evaluation.objectives[0],
|
||
);
|
||
}
|
||
|
||
#[test]
|
||
fn finds_minimum_of_rosenbrock_5d() {
|
||
let mut opt = CmaEs::new(
|
||
CmaEsConfig {
|
||
population_size: 16,
|
||
generations: 400,
|
||
initial_sigma: 0.5,
|
||
eigen_decomposition_period: 1,
|
||
initial_mean: None,
|
||
seed: 1,
|
||
},
|
||
RealBounds::new(vec![(-5.0, 5.0); 5]),
|
||
);
|
||
let r = opt.run(&Rosenbrock5D);
|
||
let best = r.best.unwrap();
|
||
// Rosenbrock is a tough non-convex valley; CMA-ES should still get
|
||
// far closer than random search.
|
||
assert!(
|
||
best.evaluation.objectives[0] < 1.0,
|
||
"got f = {}",
|
||
best.evaluation.objectives[0],
|
||
);
|
||
}
|
||
|
||
#[test]
|
||
fn deterministic_with_same_seed() {
|
||
let cfg = CmaEsConfig {
|
||
population_size: 8,
|
||
generations: 30,
|
||
initial_sigma: 0.5,
|
||
eigen_decomposition_period: 1,
|
||
initial_mean: None,
|
||
seed: 99,
|
||
};
|
||
let mut a = CmaEs::new(cfg.clone(), RealBounds::new(vec![(-5.0, 5.0)]));
|
||
let mut b = CmaEs::new(cfg, RealBounds::new(vec![(-5.0, 5.0)]));
|
||
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 = "single-objective")]
|
||
fn multi_objective_panics() {
|
||
let mut opt = CmaEs::new(CmaEsConfig::default(), RealBounds::new(vec![(-5.0, 5.0)]));
|
||
let _ = opt.run(&SchafferN1);
|
||
}
|
||
|
||
#[test]
|
||
#[should_panic(expected = "population_size must be >= 4")]
|
||
fn small_population_panics() {
|
||
let mut opt = CmaEs::new(
|
||
CmaEsConfig {
|
||
population_size: 3,
|
||
generations: 1,
|
||
initial_sigma: 0.5,
|
||
eigen_decomposition_period: 1,
|
||
initial_mean: None,
|
||
seed: 0,
|
||
},
|
||
RealBounds::new(vec![(-1.0, 1.0)]),
|
||
);
|
||
let _ = opt.run(&Sphere1D);
|
||
}
|
||
}
|