echidna_optim/
convergence.rs

1use num_traits::Float;
2
3/// Parameters controlling convergence checks.
4#[derive(Debug, Clone)]
5pub struct ConvergenceParams<F> {
6    /// Maximum number of iterations (default: 100).
7    pub max_iter: usize,
8    /// Gradient norm tolerance: stop when `||g|| < grad_tol` (default: 1e-8).
9    pub grad_tol: F,
10    /// Step size tolerance: stop when `||x_{k+1} - x_k|| < step_tol` (default: 1e-12).
11    pub step_tol: F,
12    /// Function change tolerance: stop when `|f_{k+1} - f_k| < func_tol` (default: 0, disabled).
13    pub func_tol: F,
14}
15
16impl Default for ConvergenceParams<f64> {
17    fn default() -> Self {
18        ConvergenceParams {
19            max_iter: 100,
20            grad_tol: 1e-8,
21            step_tol: 1e-12,
22            func_tol: 0.0,
23        }
24    }
25}
26
27impl Default for ConvergenceParams<f32> {
28    fn default() -> Self {
29        ConvergenceParams {
30            max_iter: 100,
31            grad_tol: 1e-5,
32            step_tol: 1e-7,
33            func_tol: 0.0,
34        }
35    }
36}
37
38/// Compute the L2 norm of a vector.
39///
40/// For `len() >= KAHAN_THRESHOLD`, uses Neumaier/Kahan compensated
41/// summation. Naive recursive summation of `len` terms accumulates
42/// `O(len·eps)` relative error in the worst case; at `len = 10^4`
43/// (f32) or `10^{13}` (f64) the ULP-level noise can leak into the
44/// convergence test and make the optimizer oscillate. Kahan drops
45/// the error to `O(eps)` independent of `len`.
46pub fn norm<F: Float>(v: &[F]) -> F {
47    kahan_sum(v.iter().map(|&x| x * x)).sqrt()
48}
49
50/// Compute the dot product of two vectors.
51pub fn dot<F: Float>(a: &[F], b: &[F]) -> F {
52    debug_assert_eq!(a.len(), b.len());
53    kahan_sum(a.iter().zip(b.iter()).map(|(&x, &y)| x * y))
54}
55
56/// The shared post-step convergence gate: gradient norm, then step size,
57/// then relative function change — the order every solver used.
58///
59/// Relative func_tol: absolute `|f_prev - f_val| < tol` is scale-blind — a
60/// tolerance of 1e-8 means ULP-precision on large-magnitude objectives
61/// (|f| ≈ 1e8) and impossibly tight on tiny ones. Scaling by `(1 + |f|)`
62/// makes the criterion track the problem. NOT for the solvers' pre-loop
63/// gradient-only checks: with no step taken, `f_prev == f_val` would fire
64/// the function-change predicate spuriously.
65pub(crate) fn check_convergence<F: Float>(
66    grad_norm: F,
67    step_norm: F,
68    f_prev: F,
69    f_val: F,
70    p: &ConvergenceParams<F>,
71) -> Option<crate::result::TerminationReason> {
72    use crate::result::TerminationReason;
73    if grad_norm < p.grad_tol {
74        return Some(TerminationReason::GradientNorm);
75    }
76    if step_norm < p.step_tol {
77        return Some(TerminationReason::StepSize);
78    }
79    if p.func_tol > F::zero() && (f_prev - f_val).abs() < p.func_tol * (F::one() + f_val.abs()) {
80        return Some(TerminationReason::FunctionChange);
81    }
82    None
83}
84
85/// Threshold above which compensated summation beats naive summation.
86/// Below this, naive summation's runtime advantage dominates and the
87/// precision gap is negligible.
88const KAHAN_THRESHOLD: usize = 64;
89
90/// Neumaier's improved Kahan summation. Handles arbitrary input
91/// magnitudes (unlike plain Kahan, which struggles when a term is
92/// larger than the running sum). Falls back to naive summation for
93/// short sequences where the overhead isn't worth it.
94#[inline]
95fn kahan_sum<F: Float, I: Iterator<Item = F>>(iter: I) -> F {
96    let mut it = iter;
97    let mut s = F::zero();
98    let mut c = F::zero();
99    let mut n = 0usize;
100    // Gather a small prefix to decide whether to use compensated summation.
101    let mut prefix: [F; KAHAN_THRESHOLD] = [F::zero(); KAHAN_THRESHOLD];
102    for slot in prefix.iter_mut() {
103        if let Some(x) = it.next() {
104            *slot = x;
105            n += 1;
106        } else {
107            break;
108        }
109    }
110    if n < KAHAN_THRESHOLD {
111        // Short case: naive accumulation. Cheaper and precision loss is
112        // bounded by `n·eps` — negligible for n < 64.
113        for &x in prefix.iter().take(n) {
114            s = s + x;
115        }
116        return s;
117    }
118    // Long case: Neumaier compensated summation.
119    for x in prefix.iter().copied().chain(it) {
120        let t = s + x;
121        if s.abs() >= x.abs() {
122            c = c + ((s - t) + x);
123        } else {
124            c = c + ((x - t) + s);
125        }
126        s = t;
127    }
128    s + c
129}
130
131#[cfg(test)]
132mod tests {
133    use super::*;
134
135    // M36: Kahan summation — very long vectors should yield tight norms.
136    #[test]
137    fn m36_norm_kahan_tight_for_long_vector() {
138        // 10_000 copies of 1.0 — exact norm is sqrt(10_000) = 100.0.
139        let v: Vec<f64> = (0..10_000).map(|_| 1.0).collect();
140        let n = norm(&v);
141        assert!(
142            (n - 100.0).abs() < 1e-12,
143            "norm of 10k ones not near 100: got {}",
144            n
145        );
146    }
147
148    #[test]
149    fn m36_norm_short_vector_still_works() {
150        let v = vec![3.0_f64, 4.0];
151        let n = norm(&v);
152        assert!((n - 5.0).abs() < 1e-15);
153    }
154}
echidna_optim/convergence.rs

echidna_optim/
convergence.rs
