rust_poly/poly/roots/single_root/
newton.rs

1use super::{line_search_accelerate, line_search_decelerate, LazyDerivatives};
2use crate::{
3    num::{Complex, Zero},
4    poly::roots,
5    roots::initial_guess::initial_guess_smallest,
6    scalar::SafeConstants,
7    util::{
8        self,
9        complex::c_from_f64,
10        doc_macros::{errors_no_converge, panic_t_from_f64},
11    },
12    Poly, RealScalar,
13};
14use num::One;
15
16/// Find a single root using Newton's method
17///
18/// # Errors
19#[doc = errors_no_converge!()]
20///
21/// # Panics
22#[doc = panic_t_from_f64!()]
23#[allow(clippy::similar_names)]
24#[allow(clippy::items_after_statements)]
25#[allow(clippy::type_complexity)]
26pub fn newton<T: RealScalar>(
27    poly: &Poly<T>,
28    epsilon: T,
29    max_iter: Option<usize>,
30    //min_iter: Option<usize>,
31    initial_guess: Option<Complex<T>>,
32) -> std::result::Result<(Vec<Complex<T>>, u128), roots::Error<Vec<Complex<T>>>> {
33    log::trace!("starting with arguments: {{poly: \"{poly}\", epsilon: {epsilon}, max_iter: \"{max_iter:?}\", initial_guess: \"{initial_guess:?}\"}}");
34
35    let mut eval_counter = 0;
36    let mut guess = initial_guess.unwrap_or_else(|| initial_guess_smallest(poly));
37    let mut guess_old = guess.clone();
38    let mut guess_old_old = guess.clone();
39    let mut guess_delta = Complex::one();
40
41    // number of iterations without improvements after which we assume we're in
42    // a cycle set it to a prime number to avoid meta-cycles forming
43    const CYCLE_COUNT_THRESHOLD: usize = 17;
44    let mut cycle_counter = 0;
45    let mut best_guess = guess.clone();
46    let mut best_px_norm = guess.norm_sqr();
47
48    let mut diffs = LazyDerivatives::new(poly);
49
50    // until convergence
51    for i in util::iterator::saturating_counter() {
52        let px = poly.eval(guess.clone());
53        eval_counter += 1;
54
55        log::trace!("best_guess: \"{best_guess:?}\"");
56
57        // stopping criterion 1: converged
58        if px.norm_sqr() <= epsilon {
59            log::trace!("stopping because target precision reached");
60            return Ok((vec![best_guess], eval_counter));
61        }
62
63        // stopping criterion 2: no improvement predicted due to numeric precision
64        if i > 3
65            && super::stopping_criterion_garwick(
66                guess.clone(),
67                guess_old.clone(),
68                guess_old_old.clone(),
69            )
70        {
71            log::trace!("stopping because garwick heuristic says no improvement is possible");
72            return Ok((vec![best_guess], eval_counter));
73        }
74
75        // max iter exceeded
76        if max_iter.is_some_and(|max| i >= max) {
77            log::trace!("did not converge {{best_guess: {best_guess}, poly: {poly}}}");
78            return Err(roots::Error::NoConverge(vec![best_guess]));
79        }
80
81        // check for cycles
82        if px.norm_sqr() >= best_px_norm {
83            cycle_counter += 1;
84            if cycle_counter > CYCLE_COUNT_THRESHOLD {
85                // arbitrary constants
86                const ROTATION_RADIANS: f64 = 0.925_024_5;
87                const SCALE: f64 = 5.0;
88
89                cycle_counter = 0;
90                log::trace!("cycle detected, backing off {{current_guess: {guess}, best_guess: {best_guess}}}");
91
92                // TODO: when const trait methods are supported, this should be
93                //       made fully const.
94                let backoff = c_from_f64(Complex::from_polar(SCALE, ROTATION_RADIANS));
95                // reverting to older base guess, but offset
96                guess = best_guess.clone() - guess_delta.clone() * backoff;
97            }
98        } else {
99            cycle_counter = 0;
100            best_guess = guess.clone();
101            best_px_norm = px.norm_sqr();
102        }
103
104        let pdx = diffs.get_nth_derivative(1).eval(guess.clone());
105        eval_counter += 1;
106
107        guess_delta = compute_delta(px.clone(), pdx.clone(), guess_delta);
108
109        // naive newton step, we're gonna improve this.
110        let mut guess_new = guess.clone() - guess_delta.clone();
111        let px_new = poly.eval(guess_new.clone());
112        let pdx_new = diffs.get_nth_derivative(1).eval(guess_new.clone());
113
114        if !check_will_converge(
115            guess_new.clone(),
116            guess.clone(),
117            px_new.clone(),
118            pdx_new,
119            pdx,
120        ) {
121            // if the current guess isn't captured, we adjust our guess more
122            // aggressively until it is captured, i.e. it is close enough that
123            // it "falls" towards the correct guess instead of jumping somewhere
124            // else.
125            if px_new.norm_sqr() > px.norm_sqr() {
126                let res = line_search_decelerate(poly, guess.clone(), guess_delta.clone());
127                guess_new = res.0;
128                eval_counter += res.1;
129            } else {
130                let res = line_search_accelerate(poly, guess.clone(), guess_delta.clone());
131                guess_new = res.0;
132                eval_counter += res.1;
133            }
134        }
135
136        guess_old_old = guess_old;
137        guess_old = guess;
138        guess = guess_new;
139    }
140    unreachable!()
141}
142
143fn compute_delta<T: RealScalar>(
144    px: Complex<T>,
145    pdx: Complex<T>,
146    delta_old: Complex<T>,
147) -> Complex<T> {
148    const EXPLODE_THRESHOLD: f64 = 5.0;
149
150    if pdx.is_zero() {
151        // these are arbitrary, originally chosen by Madsen.
152        const ROTATION_RADIANS: f64 = 0.925_024_5;
153        const SCALE: f64 = 5.0;
154        // TODO: when const trait methods are supported, this should be
155        //       made fully const.
156        let backoff = c_from_f64(Complex::from_polar(SCALE, ROTATION_RADIANS));
157        return delta_old * backoff;
158    }
159
160    let delta = px / pdx;
161
162    if delta.norm_sqr() > delta_old.norm_sqr() * T::from_f64(EXPLODE_THRESHOLD).expect("overflow") {
163        // these are arbitrary, originally chosen by Madsen.
164        const ROTATION_RADIANS: f64 = 0.925_024_5;
165        const SCALE: f64 = 5.0;
166        // TODO: when const trait methods are supported, this should be
167        //       made fully const.
168
169        let backoff = c_from_f64(Complex::from_polar(SCALE, ROTATION_RADIANS))
170            .scale(delta_old.norm_sqr() / delta.norm_sqr());
171
172        return delta * backoff;
173    }
174
175    delta
176}
177
178/// Heuristic that tries to determine if the current guess is close enough to the
179/// real root to be "captured", i.e. if the current guess is guaranteed to
180/// converge to this root, rather than jumping off to a different root.
181///
182/// This condition is based on [Ostrowski 1966](https://doi.org/10.2307/2005025),
183/// but uses an approximation by [Henrik Vestermark 2020](http://dx.doi.org/10.13140/RG.2.2.30423.34728).
184fn check_will_converge<T: RealScalar>(
185    guess: Complex<T>,
186    guess_old: Complex<T>,
187    px: Complex<T>,
188    pdx: Complex<T>,
189    pdx_old: Complex<T>,
190) -> bool {
191    !(px.clone() * pdx.clone()).is_small()
192        && ((px / pdx.clone()).norm_sqr()
193            * T::from_u8(2).expect("overflow")
194            * ((pdx_old - pdx.clone()) / (guess_old - guess)).norm_sqr()
195            <= pdx.norm_sqr())
196}
197
198#[cfg(test)]
199mod test {
200    use crate::{num::One, roots::newton_deflate, util::__testing::check_roots, Poly64};
201
202    #[test]
203    pub fn degree_0() {
204        let mut p = Poly64::one();
205        let roots = newton_deflate(&mut p, Some(1E-14), Some(100), &[]).unwrap();
206        assert!(roots.is_empty());
207        assert!(p.is_one());
208    }
209
210    #[test]
211    fn degree_1() {
212        let roots_expected = vec![complex!(1.0)];
213        let mut p = crate::Poly::from_roots(&roots_expected);
214        let roots = newton_deflate(&mut p, Some(1E-14), Some(100), &[]).unwrap();
215        assert!(check_roots(roots, roots_expected, 1E-12));
216    }
217
218    #[test]
219    fn degree_2() {
220        let roots_expected = vec![complex!(1.0), complex!(2.0)];
221        let mut p = crate::Poly::from_roots(&roots_expected);
222        let roots = newton_deflate(&mut p, Some(1E-14), Some(100), &[]).unwrap();
223        assert!(check_roots(roots, roots_expected, 1E-12));
224    }
225
226    #[test]
227    fn degree_3() {
228        let roots_expected = vec![complex!(1.0), complex!(2.0), complex!(3.0)];
229        let mut p = crate::Poly::from_roots(&roots_expected);
230        let roots = newton_deflate(&mut p, Some(1E-14), Some(100), &[]).unwrap();
231        assert!(check_roots(roots, roots_expected, 1E-5));
232    }
233
234    #[test]
235    fn degree_3_complex() {
236        let roots_expected = vec![complex!(1.0), complex!(0.0, 1.0), complex!(0.0, -1.0)];
237        let mut p = crate::Poly::from_roots(&roots_expected);
238        let roots = newton_deflate(&mut p, Some(1E-14), Some(100), &[]).unwrap();
239        assert!(check_roots(roots, roots_expected, 1E-4));
240    }
241
242    #[test]
243    fn degree_5_multiplicity_3() {
244        let roots_expected = vec![
245            complex!(1.0),
246            complex!(2.0),
247            complex!(2.0),
248            complex!(2.0),
249            complex!(3.0),
250        ];
251        let mut p = crate::Poly::from_roots(&roots_expected);
252        let roots = newton_deflate(&mut p, Some(1E-14), Some(100), &[]).unwrap();
253        assert!(check_roots(roots.clone(), roots_expected, 0.2), "{roots:?}");
254    }
255
256    #[test]
257    fn degree_5_2_zeros() {
258        let roots_expected = vec![
259            complex!(0.0),
260            complex!(0.0),
261            complex!(1.0),
262            complex!(2.0),
263            complex!(3.0),
264        ];
265        let mut p = crate::Poly::from_roots(&roots_expected);
266        let roots = newton_deflate(&mut p, Some(1E-14), Some(100), &[]).unwrap();
267        assert!(check_roots(roots, roots_expected, 1E-5));
268    }
269}
rust_poly/poly/roots/single_root/newton.rs

rust_poly/poly/roots/single_root/
newton.rs
