rustyqlib/core/optimization/
conjugate_gradient.rs

1//! Nonlinear conjugate gradient (Polak-Ribiere+): steepest descent's
2//! cost per iteration with far better search directions on ill-
3//! conditioned valleys, no matrix storage.
4
5use super::line_search::backtracking;
6use super::numerics::{dot, norm_inf, numeric_gradient};
7use super::{OptimConfig, OptimResult};
8
9/// Minimize `f` from `x0` by Polak-Ribiere conjugate gradient with the
10/// `beta >= 0` safeguard (automatic restart to steepest descent).
11/// `grad` falls back to central finite differences when absent.
12pub fn conjugate_gradient(
13    cfg: &OptimConfig,
14    f: &dyn Fn(&[f64]) -> f64,
15    grad: Option<&dyn Fn(&[f64]) -> Vec<f64>>,
16    x0: &[f64],
17) -> OptimResult {
18    let g_of = |x: &[f64]| match grad {
19        Some(g) => g(x),
20        None => numeric_gradient(f, x),
21    };
22    let mut x = x0.to_vec();
23    let mut fx = f(&x);
24    let mut g = g_of(&x);
25    let mut dir: Vec<f64> = g.iter().map(|gi| -gi).collect();
26    for it in 0..cfg.max_iter {
27        if norm_inf(&g) <= cfg.tol {
28            return OptimResult { x, value: fx, iterations: it, converged: true };
29        }
30        let mut slope = dot(&g, &dir);
31        if slope >= 0.0 {
32            // not a descent direction: restart with steepest descent
33            dir = g.iter().map(|gi| -gi).collect();
34            slope = dot(&g, &dir);
35        }
36        let alpha0 = 1.0_f64.min(1.0 / norm_inf(&g).max(1e-12));
37        let (x_new, f_new) = match backtracking(f, &x, fx, &dir, slope, alpha0) {
38            Some(step) => step,
39            None => return OptimResult { x, value: fx, iterations: it, converged: false },
40        };
41        let g_new = g_of(&x_new);
42        // Polak-Ribiere+ conjugacy factor
43        let beta = (dot(&g_new, &g_new) - dot(&g_new, &g)) / dot(&g, &g).max(1e-300);
44        let beta = beta.max(0.0);
45        dir = g_new.iter().zip(&dir).map(|(gi, di)| -gi + beta * di).collect();
46        x = x_new;
47        fx = f_new;
48        g = g_new;
49    }
50    OptimResult { x, value: fx, iterations: cfg.max_iter, converged: norm_inf(&g) <= cfg.tol }
51}
52
53#[cfg(test)]
54mod tests {
55    use super::*;
56
57    fn rosenbrock(x: &[f64]) -> f64 {
58        (1.0 - x[0]).powi(2) + 100.0 * (x[1] - x[0] * x[0]).powi(2)
59    }
60
61    #[test]
62    fn minimizes_the_rosenbrock_valley() {
63        let f = |x: &[f64]| rosenbrock(x);
64        let r = conjugate_gradient(&OptimConfig::new(1e-6, 20_000), &f, None, &[-1.2, 1.0]);
65        assert!((r.x[0] - 1.0).abs() < 1e-3 && (r.x[1] - 1.0).abs() < 1e-3, "{r:?}");
66        assert!(r.value < 1e-6, "{r:?}");
67    }
68
69    #[test]
70    fn converges_on_an_ill_conditioned_quadratic() {
71        // condition number 10^6 across four dimensions
72        let f = |x: &[f64]| {
73            x.iter()
74                .enumerate()
75                .map(|(i, xi)| 10.0_f64.powi(2 * i as i32) * xi * xi)
76                .sum::<f64>()
77        };
78        let cfg = OptimConfig::new(1e-8, 50_000);
79        let cg = conjugate_gradient(&cfg, &f, None, &[1.0, 1.0, 1.0, 1.0]);
80        assert!(cg.converged, "{cg:?}");
81        assert!(cg.value < 1e-10, "{cg:?}");
82    }
83}
rustyqlib/core/optimization/conjugate_gradient.rs

rustyqlib/core/optimization/
conjugate_gradient.rs
