gosh_lbfgs/
line.rs

1// docs
2// Copyright (c) 1990, Jorge Nocedal
3// Copyright (c) 2007-2010 Naoaki Okazaki
4// Copyright (c) 2018-2019 Wenping Guo
5// All rights reserved.
6
7//! # Find a satisfactory step length along predefined search direction
8//!
9//! # Example
10//! ```
11//! use liblbfgs::math::LbfgsMath;
12//! use liblbfgs::Problem;
13//! use liblbfgs::default_evaluate;
14//! use liblbfgs::line::LineSearch;
15//! 
16//! const N: usize = 100;
17//! let mut x = [0.0; N];
18//! for i in (0..N).step_by(2) {
19//!     x[i] = -1.2;
20//!     x[i + 1] = 1.0;
21//! }
22//! 
23//! // construct problem
24//! let mut prb = Problem::new(&mut x, default_evaluate(), None);
25//! prb.evaluate();
26//! // construct initial search direction
27//! prb.update_search_direction();
28//! // Compute the initial step
29//! let mut step = 1.0/prb.search_direction().vec2norm();
30//! 
31//! let ls = LineSearch::default();
32//! let ncall = ls.find(&mut prb, &mut step).expect("line search");
33//! ```
34
35// [[file:~/Workspace/Programming/gosh-rs/lbfgs/lbfgs.note::*base][base:1]]
36use crate::core::*;
37// base:1 ends here
38
39// [[file:~/Workspace/Programming/gosh-rs/lbfgs/lbfgs.note::*algorithm][algorithm:1]]
40/// Line search algorithms.
41#[derive(Debug, Copy, Clone, PartialEq)]
42pub enum LineSearchAlgorithm {
43    /// MoreThuente method proposd by More and Thuente. This is the default for
44    /// regular LBFGS.
45    MoreThuente,
46
47    ///
48    /// Backtracking method with the Armijo condition.
49    ///
50    /// The backtracking method finds the step length such that it satisfies
51    /// the sufficient decrease (Armijo) condition,
52    ///   - f(x + a * d) <= f(x) + ftol * a * g(x)^T d,
53    ///
54    /// where x is the current point, d is the current search direction, and
55    /// a is the step length.
56    ///
57    BacktrackingArmijo,
58
59    /// Backtracking method with strong Wolfe condition.
60    ///
61    /// The backtracking method finds the step length such that it satisfies
62    /// both the Armijo condition (BacktrackingArmijo)
63    /// and the following condition,
64    ///   - |g(x + a * d)^T d| <= gtol * |g(x)^T d|,
65    ///
66    /// where x is the current point, d is the current search direction, and
67    /// a is the step length.
68    ///
69    BacktrackingStrongWolfe,
70
71    ///
72    /// Backtracking method with regular Wolfe condition.
73    ///
74    /// The backtracking method finds the step length such that it satisfies
75    /// both the Armijo condition (BacktrackingArmijo)
76    /// and the curvature condition,
77    ///   - g(x + a * d)^T d >= gtol * g(x)^T d,
78    ///
79    /// where x is the current point, d is the current search direction, and a
80    /// is the step length.
81    ///
82    BacktrackingWolfe,
83}
84
85impl Default for LineSearchAlgorithm {
86    /// The default algorithm (MoreThuente method).
87    fn default() -> Self {
88        LineSearchAlgorithm::MoreThuente
89    }
90}
91// algorithm:1 ends here
92
93// [[file:~/Workspace/Programming/gosh-rs/lbfgs/lbfgs.note::*paramters][paramters:1]]
94#[derive(Debug, Copy, Clone)]
95pub struct LineSearch {
96    /// Various line search algorithms.
97    pub algorithm: LineSearchAlgorithm,
98
99    /// A parameter to control the accuracy of the line search routine.
100    ///
101    /// The default value is 1e-4. This parameter should be greater
102    /// than zero and smaller than 0.5.
103    pub ftol: f64,
104
105    /// A parameter to control the accuracy of the line search routine.
106    ///
107    /// The default value is 0.9. If the function and gradient evaluations are
108    /// inexpensive with respect to the cost of the iteration (which is
109    /// sometimes the case when solving very large problems) it may be
110    /// advantageous to set this parameter to a small value. A typical small
111    /// value is 0.1. This parameter shuold be greater than the ftol parameter
112    /// (1e-4) and smaller than 1.0.
113    pub gtol: f64,
114
115    /// xtol is a nonnegative input variable. termination occurs when the
116    /// relative width of the interval of uncertainty is at most xtol.
117    ///
118    /// The machine precision for floating-point values.
119    ///
120    ///  This parameter must be a positive value set by a client program to
121    ///  estimate the machine precision. The line search routine will terminate
122    ///  with the status code (::LBFGSERR_ROUNDING_ERROR) if the relative width
123    ///  of the interval of uncertainty is less than this parameter.
124    pub xtol: f64,
125
126    /// The minimum step of the line search routine.
127    ///
128    /// The default value is \c 1e-20. This value need not be modified unless
129    /// the exponents are too large for the machine being used, or unless the
130    /// problem is extremely badly scaled (in which case the exponents should be
131    /// increased).
132    pub min_step: f64,
133
134    /// The maximum step of the line search.
135    ///
136    ///  The default value is \c 1e+20. This value need not be modified unless
137    ///  the exponents are too large for the machine being used, or unless the
138    ///  problem is extremely badly scaled (in which case the exponents should
139    ///  be increased).
140    pub max_step: f64,
141
142    /// The maximum number of trials for the line search.
143    ///
144    /// This parameter controls the number of function and gradients evaluations
145    /// per iteration for the line search routine. The default value is 40. Set
146    /// this value to 0, will completely disable line search.
147    ///
148    pub max_linesearch: usize,
149
150    /// Make line search conditions use only gradients.
151    pub gradient_only: bool,
152}
153
154// TODO: better defaults
155impl Default for LineSearch {
156    fn default() -> Self {
157        LineSearch {
158            ftol: 1e-4,
159            gtol: 0.9,
160            xtol: 1e-16,
161            min_step: 1e-20,
162            max_step: 1e20,
163            max_linesearch: 40,
164            gradient_only: false,
165            algorithm: LineSearchAlgorithm::default(),
166        }
167    }
168}
169// paramters:1 ends here
170
171// [[file:~/Workspace/Programming/gosh-rs/lbfgs/lbfgs.note::*adhoc][adhoc:1]]
172impl LineSearch {
173    fn validate_step(&self, step: f64) -> Result<()> {
174        // The step is the minimum value.
175        if step < self.min_step {
176            bail!("The line-search step became smaller than LineSearch::min_step.");
177        }
178        // The step is the maximum value.
179        if step > self.max_step {
180            bail!("The line-search step became larger than LineSearch::max_step.");
181        }
182
183        Ok(())
184    }
185}
186// adhoc:1 ends here
187
188// [[file:~/Workspace/Programming/gosh-rs/lbfgs/lbfgs.note::*entry][entry:1]]
189use crate::lbfgs::Problem;
190
191impl LineSearch {
192    /// Unified entry for line searches
193    ///
194    /// # Arguments
195    ///
196    /// * step: initial step size. On output it will be the optimal step size.
197    /// * direction: proposed searching direction
198    /// * eval_fn: a callback function to evaluate `Problem`
199    ///
200    /// # Return
201    ///
202    /// * On success, return the number of line searching iterations
203    ///
204    pub fn find<E>(&self, prb: &mut Problem<E>, step: &mut f64) -> Result<usize>
205    where
206        E: FnMut(&[f64], &mut [f64]) -> Result<f64>,
207    {
208        // Check the input parameters for errors.
209        ensure!(
210            step.is_sign_positive(),
211            "A logic error (negative line-search step) occurred."
212        );
213
214        // Search for an optimal step.
215        let ls = if self.algorithm == MoreThuente && !prb.orthantwise() {
216            if !self.gradient_only {
217                line_search_morethuente(prb, step, &self)
218            } else {
219                bail!("Gradient only optimization is incompatible with MoreThuente line search.");
220            }
221        } else {
222            line_search_backtracking(prb, step, &self)
223        }
224        .unwrap_or_else(|err| {
225            // Revert to the previous point.
226            error!("line search failed, revert to the previous point!");
227            prb.revert();
228            println!("{:?}", err);
229
230            0
231        });
232
233        Ok(ls)
234    }
235}
236// entry:1 ends here
237
238// [[file:~/Workspace/Programming/gosh-rs/lbfgs/lbfgs.note::*Original documentation by J. Nocera (lbfgs.f)][Original documentation by J. Nocera (lbfgs.f):1]]
239
240// Original documentation by J. Nocera (lbfgs.f):1 ends here
241
242// [[file:~/Workspace/Programming/gosh-rs/lbfgs/lbfgs.note::*core][core:1]]
243use crate::math::*;
244
245pub fn line_search_morethuente<E>(
246    prb: &mut Problem<E>,
247    stp: &mut f64,      // Step size
248    param: &LineSearch, // line search parameters
249) -> Result<usize>
250where
251    E: FnMut(&[f64], &mut [f64]) -> Result<f64>,
252{
253    // Initialize local variables.
254    let dginit = prb.dginit()?;
255    let mut brackt = false;
256    let mut stage1 = 1;
257    let mut uinfo = 0;
258
259    let finit = prb.fx;
260    let dgtest = param.ftol * dginit;
261    let mut width = param.max_step - param.min_step;
262    let mut prev_width = 2.0 * width;
263
264    // The variables stx, fx, dgx contain the values of the step,
265    // function, and directional derivative at the best step.
266    // The variables sty, fy, dgy contain the value of the step,
267    // function, and derivative at the other endpoint of
268    // the interval of uncertainty.
269    // The variables stp, f, dg contain the values of the step,
270    // function, and derivative at the current step.
271    let (mut stx, mut sty) = (0.0, 0.0);
272    let mut fx = finit;
273    let mut fy = finit;
274    let mut dgy = dginit;
275    let mut dgx = dgy;
276
277    for count in 0..param.max_linesearch {
278        // Set the minimum and maximum steps to correspond to the
279        // present interval of uncertainty.
280        let (stmin, stmax) = if brackt {
281            (if stx <= sty { stx } else { sty }, if stx >= sty { stx } else { sty })
282        } else {
283            (stx, *stp + 4.0 * (*stp - stx))
284        };
285
286        // Clip the step in the range of [stpmin, stpmax].
287        if *stp < param.min_step {
288            *stp = param.min_step
289        }
290        if param.max_step < *stp {
291            *stp = param.max_step
292        }
293
294        // If an unusual termination is to occur then let
295        // stp be the lowest point obtained so far.
296        if brackt && (*stp <= stmin || stmax <= *stp || param.max_linesearch <= count + 1 || uinfo != 0)
297            || brackt && stmax - stmin <= param.xtol * stmax
298        {
299            *stp = stx
300        }
301
302        prb.take_line_step(*stp);
303
304        // Evaluate the function and gradient values.
305        prb.evaluate()?;
306        let f = prb.fx;
307        let dg = prb.dg_unchecked();
308        let ftest1 = finit + *stp * dgtest;
309
310        // Test for errors and convergence.
311        if brackt && (*stp <= stmin || stmax <= *stp || uinfo != 0i32) {
312            // Rounding errors prevent further progress.
313            bail!(
314                "A rounding error occurred; alternatively, no line-search step
315satisfies the sufficient decrease and curvature conditions."
316            );
317        }
318
319        if brackt && stmax - stmin <= param.xtol * stmax {
320            bail!("Relative width of the interval of uncertainty is at most xtol.");
321        }
322
323        // FIXME: float == float?
324        if *stp == param.max_step && f <= ftest1 && dg <= dgtest {
325            // The step is the maximum value.
326            bail!("The line-search step became larger than LineSearch::max_step.");
327        }
328        // FIXME: float == float?
329        if *stp == param.min_step && (ftest1 < f || dgtest <= dg) {
330            // The step is the minimum value.
331            bail!("The line-search step became smaller than LineSearch::min_step.");
332        }
333
334        if dg.abs() <= param.gtol * -dginit {
335            // the directional derivative condition hold.
336            return Ok(count);
337        } else if f <= ftest1 && dg.abs() <= param.gtol * -dginit {
338            // The sufficient decrease condition and the directional derivative condition hold.
339            return Ok(count);
340        } else {
341            // In the first stage we seek a step for which the modified
342            // function has a nonpositive value and nonnegative derivative.
343            if 0 != stage1 && f <= ftest1 && param.ftol.min(param.gtol) * dginit <= dg {
344                stage1 = 0;
345            }
346
347            // A modified function is used to predict the step only if
348            // we have not obtained a step for which the modified
349            // function has a nonpositive function value and nonnegative
350            // derivative, and if a lower function value has been
351            // obtained but the decrease is not sufficient.
352            if 0 != stage1 && ftest1 < f && f <= fx {
353                // Define the modified function and derivative values.
354                let fm = f - *stp * dgtest;
355                let mut fxm = fx - stx * dgtest;
356                let mut fym = fy - sty * dgtest;
357                let dgm = dg - dgtest;
358                let mut dgxm = dgx - dgtest;
359                let mut dgym = dgy - dgtest;
360
361                // Call update_trial_interval() to update the interval of
362                // uncertainty and to compute the new step.
363                uinfo = mcstep::update_trial_interval(
364                    &mut stx,
365                    &mut fxm,
366                    &mut dgxm,
367                    &mut sty,
368                    &mut fym,
369                    &mut dgym,
370                    &mut *stp,
371                    fm,
372                    dgm,
373                    stmin,
374                    stmax,
375                    &mut brackt,
376                )?;
377
378                // Reset the function and gradient values for f.
379                fx = fxm + stx * dgtest;
380                fy = fym + sty * dgtest;
381                dgx = dgxm + dgtest;
382                dgy = dgym + dgtest
383            } else {
384                uinfo = mcstep::update_trial_interval(
385                    &mut stx,
386                    &mut fx,
387                    &mut dgx,
388                    &mut sty,
389                    &mut fy,
390                    &mut dgy,
391                    &mut *stp,
392                    f,
393                    dg,
394                    stmin,
395                    stmax,
396                    &mut brackt,
397                )?;
398            }
399
400            // Force a sufficient decrease in the interval of uncertainty.
401            if !(brackt) {
402                continue;
403            }
404
405            if 0.66 * prev_width <= (sty - stx).abs() {
406                *stp = stx + 0.5 * (sty - stx)
407            }
408
409            prev_width = width;
410            width = (sty - stx).abs()
411        }
412    }
413
414    // Maximum number of iteration.
415    info!("The line-search routine reaches the maximum number of evaluations.");
416
417    Ok(param.max_linesearch)
418}
419// core:1 ends here
420
421// [[file:~/Workspace/Programming/gosh-rs/lbfgs/lbfgs.note::*core][core:1]]
422/// Represents the original MCSTEP subroutine by J. Nocera, which is a variant
423/// of More' and Thuente's routine.
424///
425/// The purpose of mcstep is to compute a safeguarded step for a linesearch and
426/// to update an interval of uncertainty for a minimizer of the function.
427///
428/// Documentation is adopted from the original Fortran codes.
429mod mcstep {
430    // dependencies
431    use super::{cubic_minimizer, cubic_minimizer2, quard_minimizer, quard_minimizer2};
432
433    use crate::core::*;
434
435    ///
436    /// Update a safeguarded trial value and interval for line search.
437    ///
438    /// This function assumes that the derivative at the point of x in the
439    /// direction of the step. If the bracket is set to true, the minimizer has
440    /// been bracketed in an interval of uncertainty with endpoints between x
441    /// and y.
442    ///
443    /// # Arguments
444    ///
445    /// * x, fx, and dx: variables which specify the step, the function, and the
446    /// derivative at the best step obtained so far. The derivative must be
447    /// negative in the direction of the step, that is, dx and t-x must have
448    /// opposite signs. On output these parameters are updated appropriately.
449    ///
450    /// * y, fy, and dy: variables which specify the step, the function, and
451    /// the derivative at the other endpoint of the interval of uncertainty. On
452    /// output these parameters are updated appropriately.
453    ///
454    /// * t, ft, and dt: variables which specify the step, the function, and the
455    /// derivative at the current step. If bracket is set true then on input t
456    /// must be between x and y. On output t is set to the new step.
457    ///
458    /// * tmin, tmax: lower and upper bounds for the step.
459    ///
460    /// * `brackt`: Specifies if a minimizer has been bracketed. If the
461    /// minimizer has not been bracketed then on input brackt must be set false.
462    /// If the minimizer is bracketed then on output `brackt` is set true.
463    ///
464    /// # Return
465    /// - Status value. Zero indicates a normal termination.
466    ///
467    pub(crate) fn update_trial_interval(
468        x: &mut f64,
469        fx: &mut f64,
470        dx: &mut f64,
471        y: &mut f64,
472        fy: &mut f64,
473        dy: &mut f64,
474        t: &mut f64,
475        ft: f64,
476        dt: f64,
477        tmin: f64,
478        tmax: f64,
479        brackt: &mut bool,
480    ) -> Result<i32> {
481        // fsigndiff
482        let dsign = dt * (*dx / (*dx).abs()) < 0.0;
483        // minimizer of an interpolated cubic.
484        let mut mc = 0.;
485        // minimizer of an interpolated quadratic.
486        let mut mq = 0.;
487        // new trial value.
488        let mut newt = 0.;
489
490        // Check the input parameters for errors.
491        if *brackt {
492            if *t <= x.min(*y) || x.max(*y) <= *t {
493                // The trival value t is out of the interval.
494                bail!("The line-search step went out of the interval of uncertainty.");
495            } else if 0.0 <= *dx * (*t - *x) {
496                // The function must decrease from x.
497                bail!("The current search direction increases the objective function value.");
498            } else if tmax < tmin {
499                // Incorrect tmin and tmax specified.
500                bail!("A logic error occurred; alternatively, the interval of uncertainty became too small.");
501            }
502        }
503
504        // Trial value selection.
505        let bound = if *fx < ft {
506            // Case 1: a higher function value.
507            // The minimum is brackt. If the cubic minimizer is closer
508            // to x than the quadratic one, the cubic one is taken, else
509            // the average of the minimizers is taken.
510            *brackt = true;
511            cubic_minimizer(&mut mc, *x, *fx, *dx, *t, ft, dt);
512            quard_minimizer(&mut mq, *x, *fx, *dx, *t, ft);
513            if (mc - *x).abs() < (mq - *x).abs() {
514                newt = mc
515            } else {
516                newt = mc + 0.5 * (mq - mc)
517            }
518
519            1
520        } else if dsign {
521            // Case 2: a lower function value and derivatives of
522            // opposite sign. The minimum is brackt. If the cubic
523            // minimizer is closer to x than the quadratic (secant) one,
524            // the cubic one is taken, else the quadratic one is taken.
525            *brackt = true;
526            cubic_minimizer(&mut mc, *x, *fx, *dx, *t, ft, dt);
527            quard_minimizer2(&mut mq, *x, *dx, *t, dt);
528            if (mc - *t).abs() > (mq - *t).abs() {
529                newt = mc
530            } else {
531                newt = mq
532            }
533
534            0
535        } else if dt.abs() < (*dx).abs() {
536            // Case 3: a lower function value, derivatives of the
537            // same sign, and the magnitude of the derivative decreases.
538            // The cubic minimizer is only used if the cubic tends to
539            // infinity in the direction of the minimizer or if the minimum
540            // of the cubic is beyond t. Otherwise the cubic minimizer is
541            // defined to be either tmin or tmax. The quadratic (secant)
542            // minimizer is also computed and if the minimum is brackt
543            // then the the minimizer closest to x is taken, else the one
544            // farthest away is taken.
545            cubic_minimizer2(&mut mc, *x, *fx, *dx, *t, ft, dt, tmin, tmax);
546            quard_minimizer2(&mut mq, *x, *dx, *t, dt);
547            if *brackt {
548                if (*t - mc).abs() < (*t - mq).abs() {
549                    newt = mc
550                } else {
551                    newt = mq
552                }
553            } else if (*t - mc).abs() > (*t - mq).abs() {
554                newt = mc
555            } else {
556                newt = mq
557            }
558
559            1
560        } else {
561            // Case 4: a lower function value, derivatives of the
562            // same sign, and the magnitude of the derivative does
563            // not decrease. If the minimum is not brackt, the step
564            // is either tmin or tmax, else the cubic minimizer is taken.
565            if *brackt {
566                cubic_minimizer(&mut newt, *t, ft, dt, *y, *fy, *dy);
567            } else if *x < *t {
568                newt = tmax
569            } else {
570                newt = tmin
571            }
572
573            0
574        };
575
576        // Update the interval of uncertainty. This update does not
577        // depend on the new step or the case analysis above.
578        // - Case a: if f(x) < f(t),
579        //    x <- x, y <- t.
580        // - Case b: if f(t) <= f(x) && f'(t)*f'(x) > 0,
581        //   x <- t, y <- y.
582        // - Case c: if f(t) <= f(x) && f'(t)*f'(x) < 0,
583        //   x <- t, y <- x.
584        if *fx < ft {
585            /* Case a */
586            *y = *t;
587            *fy = ft;
588            *dy = dt
589        } else {
590            /* Case c */
591            if dsign {
592                *y = *x;
593                *fy = *fx;
594                *dy = *dx
595            }
596            /* Cases b and c */
597            *x = *t;
598            *fx = ft;
599            *dx = dt
600        }
601
602        // Clip the new trial value in [tmin, tmax].
603        if tmax < newt {
604            newt = tmax
605        }
606        if newt < tmin {
607            newt = tmin
608        }
609
610        // Redefine the new trial value if it is close to the upper bound of the
611        // interval.
612        if *brackt && 0 != bound {
613            mq = *x + 0.66 * (*y - *x);
614            if *x < *y {
615                if mq < newt {
616                    newt = mq
617                }
618            } else if newt < mq {
619                newt = mq
620            }
621        }
622
623        // Return the new trial value.
624        *t = newt;
625
626        Ok(0)
627    }
628}
629// core:1 ends here
630
631// [[file:~/Workspace/Programming/gosh-rs/lbfgs/lbfgs.note::*interpolation][interpolation:1]]
632/// Find a minimizer of an interpolated cubic function.
633///
634/// # Arguments
635///  * `cm`: The minimizer of the interpolated cubic.
636///  * `u` : The value of one point, u.
637///  * `fu`: The value of f(u).
638///  * `du`: The value of f'(u).
639///  * `v` : The value of another point, v.
640///  * `fv`: The value of f(v).
641///  * `dv`:  The value of f'(v).
642#[inline]
643fn cubic_minimizer(cm: &mut f64, u: f64, fu: f64, du: f64, v: f64, fv: f64, dv: f64) {
644    let d = v - u;
645    let theta = (fu - fv) * 3.0 / d + du + dv;
646
647    let mut p = theta.abs();
648    let mut q = du.abs();
649    let mut r = dv.abs();
650    let s = (p.max(q)).max(r); // max3(p, q, r)
651    let a = theta / s;
652    let mut gamma = s * (a * a - du / s * (dv / s)).sqrt();
653    if v < u {
654        gamma = -gamma
655    }
656    p = gamma - du + theta;
657    q = gamma - du + gamma + dv;
658    r = p / q;
659    *cm = u + r * d;
660}
661
662/// Find a minimizer of an interpolated cubic function.
663///
664/// # Arguments
665///  * cm  :   The minimizer of the interpolated cubic.
666///  * u   :   The value of one point, u.
667///  * fu  :   The value of f(u).
668///  * du  :   The value of f'(u).
669///  * v   :   The value of another point, v.
670///  * fv  :   The value of f(v).
671///  * dv  :   The value of f'(v).
672///  * xmin:   The minimum value.
673///  * xmax:   The maximum value.
674#[inline]
675fn cubic_minimizer2(cm: &mut f64, u: f64, fu: f64, du: f64, v: f64, fv: f64, dv: f64, xmin: f64, xmax: f64) {
676    // STP - STX
677    let d = v - u;
678    let theta = (fu - fv) * 3.0 / d + du + dv;
679    let mut p = theta.abs();
680    let mut q = du.abs();
681    let mut r = dv.abs();
682    // s = max3(p, q, r);
683    let s = (p.max(q)).max(r); // max3(p, q, r)
684    let a = theta / s;
685
686    let mut gamma = s * (0f64.max(a * a - du / s * (dv / s)).sqrt());
687    // STX < STP
688    if u < v {
689        gamma = -gamma
690    }
691    p = gamma - dv + theta;
692    q = gamma - dv + gamma + du;
693    r = p / q;
694    if r < 0.0 && gamma != 0.0 {
695        *cm = v - r * d;
696    // } else if a < 0 as f64 {
697    //  ELSE IF (STP .GT. STX) THEN
698    } else if v > u {
699        *cm = xmax;
700    } else {
701        *cm = xmin;
702    }
703}
704
705/// Find a minimizer of an interpolated quadratic function.
706///
707/// # Arguments
708/// * qm : The minimizer of the interpolated quadratic.
709/// * u  : The value of one point, u.
710/// * fu : The value of f(u).
711/// * du : The value of f'(u).
712/// * v  : The value of another point, v.
713/// * fv : The value of f(v).
714#[inline]
715fn quard_minimizer(qm: &mut f64, u: f64, fu: f64, du: f64, v: f64, fv: f64) {
716    let a = v - u;
717    *qm = u + du / ((fu - fv) / a + du) / 2.0 * a;
718}
719
720/// Find a minimizer of an interpolated quadratic function.
721///
722/// # Arguments
723/// * `qm` :    The minimizer of the interpolated quadratic.
724/// * `u`  :    The value of one point, u.
725/// * `du` :    The value of f'(u).
726/// * `v`  :    The value of another point, v.
727/// * `dv` :    The value of f'(v).
728#[inline]
729fn quard_minimizer2(qm: &mut f64, u: f64, du: f64, v: f64, dv: f64) {
730    let a = u - v;
731    *qm = v + dv / (dv - du) * a;
732}
733// interpolation:1 ends here
734
735// [[file:~/Workspace/Programming/gosh-rs/lbfgs/lbfgs.note::*core][core:1]]
736use self::LineSearchAlgorithm::*;
737
738/// `prb` holds input variables `x`, gradient `gx` arrays, and function value
739/// `fx`. on input it must contain the base point for the line search. on output
740/// it contains data on x + stp*d.
741pub fn line_search_backtracking<E>(
742    prb: &mut Problem<E>,
743    stp: &mut f64,      // step length
744    param: &LineSearch, // line search parameters
745) -> Result<usize>
746where
747    E: FnMut(&[f64], &mut [f64]) -> Result<f64>,
748{
749    let dginit = prb.dginit()?;
750    let dec: f64 = 0.5;
751    let inc: f64 = 2.1;
752
753    // quick wrapper
754    let orthantwise = prb.orthantwise();
755
756    // The initial value of the objective function.
757    let finit = prb.fx;
758    let dgtest = param.ftol * dginit;
759
760    let mut width: f64;
761    for count in 0..param.max_linesearch {
762        prb.take_line_step(*stp);
763
764        // Evaluate the function and gradient values.
765        prb.evaluate()?;
766
767        if prb.fx > finit + *stp * dgtest {
768            width = dec;
769        } else if param.algorithm == BacktrackingArmijo || orthantwise {
770            // The sufficient decrease condition.
771            // Exit with the Armijo condition.
772            return Ok(count);
773        } else {
774            // Check the Wolfe condition.
775            let dg = prb.dg_unchecked();
776            if dg < param.gtol * dginit {
777                width = inc
778            } else if param.algorithm == BacktrackingWolfe {
779                // Exit with the regular Wolfe condition.
780                return Ok(count);
781            } else if dg > -param.gtol * dginit {
782                width = dec
783            } else {
784                return Ok(count);
785            }
786        }
787
788        // allow energy rises
789        // only check strong wolfe condition
790        if param.gradient_only {
791            info!("allow energy rises");
792            let dg = prb.dg_unchecked();
793            if dg.abs() <= -param.gtol * dginit.abs() {
794                return Ok(count);
795            }
796        }
797
798        param.validate_step(*stp)?;
799        *stp *= width
800    }
801
802    // Maximum number of iteration.
803    info!("The line-search routine reaches the maximum number of evaluations.");
804
805    Ok(param.max_linesearch)
806}
807// core:1 ends here