Skip to main content

rustyqlib/core/optimization/
steepest_descent.rs

1//! Steepest descent: follow `-gradient` with a backtracking line
2//! search. Linear convergence — the robust baseline, and the reference
3//! the fancier methods are tested against.
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 steepest descent. `grad` falls back to
10/// central finite differences when absent.
11pub fn steepest_descent(
12    cfg: &OptimConfig,
13    f: &dyn Fn(&[f64]) -> f64,
14    grad: Option<&dyn Fn(&[f64]) -> Vec<f64>>,
15    x0: &[f64],
16) -> OptimResult {
17    let g_of = |x: &[f64]| match grad {
18        Some(g) => g(x),
19        None => numeric_gradient(f, x),
20    };
21    let mut x = x0.to_vec();
22    let mut fx = f(&x);
23    for it in 0..cfg.max_iter {
24        let g = g_of(&x);
25        if norm_inf(&g) <= cfg.tol {
26            return OptimResult { x, value: fx, iterations: it, converged: true };
27        }
28        let dir: Vec<f64> = g.iter().map(|gi| -gi).collect();
29        let slope = dot(&g, &dir);
30        // scale the first trial step to a unit-size move
31        let alpha0 = 1.0_f64.min(1.0 / norm_inf(&g).max(1e-12));
32        match backtracking(f, &x, fx, &dir, slope, alpha0) {
33            Some((x_new, f_new)) => {
34                x = x_new;
35                fx = f_new;
36            }
37            None => return OptimResult { x, value: fx, iterations: it, converged: false },
38        }
39    }
40    let g = g_of(&x);
41    OptimResult { x, value: fx, iterations: cfg.max_iter, converged: norm_inf(&g) <= cfg.tol }
42}
43
44#[cfg(test)]
45mod tests {
46    use super::*;
47
48    #[test]
49    fn minimizes_a_convex_quadratic() {
50        let f = |x: &[f64]| 2.0 * (x[0] - 3.0).powi(2) + (x[1] + 1.0).powi(2);
51        let r = steepest_descent(&OptimConfig::new(1e-8, 5000), &f, None, &[10.0, -10.0]);
52        assert!(r.converged, "{r:?}");
53        assert!((r.x[0] - 3.0).abs() < 1e-6 && (r.x[1] + 1.0).abs() < 1e-6, "{:?}", r.x);
54    }
55
56    #[test]
57    fn analytic_gradient_is_used_when_given() {
58        let f = |x: &[f64]| (x[0] - 1.0).powi(2);
59        let g = |x: &[f64]| vec![2.0 * (x[0] - 1.0)];
60        let r = steepest_descent(&OptimConfig::default(), &f, Some(&g), &[7.0]);
61        assert!(r.converged && (r.x[0] - 1.0).abs() < 1e-7, "{r:?}");
62    }
63}