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