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/// Threshold above which compensated summation beats naive summation.
57/// Below this, naive summation's runtime advantage dominates and the
58/// precision gap is negligible.
59const KAHAN_THRESHOLD: usize = 64;
60
61/// Neumaier's improved Kahan summation. Handles arbitrary input
62/// magnitudes (unlike plain Kahan, which struggles when a term is
63/// larger than the running sum). Falls back to naive summation for
64/// short sequences where the overhead isn't worth it.
65#[inline]
66fn kahan_sum<F: Float, I: Iterator<Item = F>>(iter: I) -> F {
67    let mut it = iter;
68    let mut s = F::zero();
69    let mut c = F::zero();
70    let mut n = 0usize;
71    // Gather a small prefix to decide whether to use compensated summation.
72    let mut prefix: [F; KAHAN_THRESHOLD] = [F::zero(); KAHAN_THRESHOLD];
73    for slot in prefix.iter_mut() {
74        if let Some(x) = it.next() {
75            *slot = x;
76            n += 1;
77        } else {
78            break;
79        }
80    }
81    if n < KAHAN_THRESHOLD {
82        // Short case: naive accumulation. Cheaper and precision loss is
83        // bounded by `n·eps` — negligible for n < 64.
84        for &x in prefix.iter().take(n) {
85            s = s + x;
86        }
87        return s;
88    }
89    // Long case: Neumaier compensated summation.
90    for &x in prefix.iter() {
91        let t = s + x;
92        if s.abs() >= x.abs() {
93            c = c + ((s - t) + x);
94        } else {
95            c = c + ((x - t) + s);
96        }
97        s = t;
98    }
99    for x in it {
100        let t = s + x;
101        if s.abs() >= x.abs() {
102            c = c + ((s - t) + x);
103        } else {
104            c = c + ((x - t) + s);
105        }
106        s = t;
107    }
108    s + c
109}
echidna_optim/convergence.rs

echidna_optim/
convergence.rs
