Coverage for /usr/lib/python3/dist-packages/scipy/optimize/_linesearch.py: 6%
305 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1"""
2Functions
3---------
4.. autosummary::
5 :toctree: generated/
7 line_search_armijo
8 line_search_wolfe1
9 line_search_wolfe2
10 scalar_search_wolfe1
11 scalar_search_wolfe2
13"""
14from warnings import warn
16from scipy.optimize import _minpack2 as minpack2
17import numpy as np
19__all__ = ['LineSearchWarning', 'line_search_wolfe1', 'line_search_wolfe2',
20 'scalar_search_wolfe1', 'scalar_search_wolfe2',
21 'line_search_armijo']
23class LineSearchWarning(RuntimeWarning):
24 pass
27#------------------------------------------------------------------------------
28# Minpack's Wolfe line and scalar searches
29#------------------------------------------------------------------------------
31def line_search_wolfe1(f, fprime, xk, pk, gfk=None,
32 old_fval=None, old_old_fval=None,
33 args=(), c1=1e-4, c2=0.9, amax=50, amin=1e-8,
34 xtol=1e-14):
35 """
36 As `scalar_search_wolfe1` but do a line search to direction `pk`
38 Parameters
39 ----------
40 f : callable
41 Function `f(x)`
42 fprime : callable
43 Gradient of `f`
44 xk : array_like
45 Current point
46 pk : array_like
47 Search direction
49 gfk : array_like, optional
50 Gradient of `f` at point `xk`
51 old_fval : float, optional
52 Value of `f` at point `xk`
53 old_old_fval : float, optional
54 Value of `f` at point preceding `xk`
56 The rest of the parameters are the same as for `scalar_search_wolfe1`.
58 Returns
59 -------
60 stp, f_count, g_count, fval, old_fval
61 As in `line_search_wolfe1`
62 gval : array
63 Gradient of `f` at the final point
65 """
66 if gfk is None:
67 gfk = fprime(xk, *args)
69 gval = [gfk]
70 gc = [0]
71 fc = [0]
73 def phi(s):
74 fc[0] += 1
75 return f(xk + s*pk, *args)
77 def derphi(s):
78 gval[0] = fprime(xk + s*pk, *args)
79 gc[0] += 1
80 return np.dot(gval[0], pk)
82 derphi0 = np.dot(gfk, pk)
84 stp, fval, old_fval = scalar_search_wolfe1(
85 phi, derphi, old_fval, old_old_fval, derphi0,
86 c1=c1, c2=c2, amax=amax, amin=amin, xtol=xtol)
88 return stp, fc[0], gc[0], fval, old_fval, gval[0]
91def scalar_search_wolfe1(phi, derphi, phi0=None, old_phi0=None, derphi0=None,
92 c1=1e-4, c2=0.9,
93 amax=50, amin=1e-8, xtol=1e-14):
94 """
95 Scalar function search for alpha that satisfies strong Wolfe conditions
97 alpha > 0 is assumed to be a descent direction.
99 Parameters
100 ----------
101 phi : callable phi(alpha)
102 Function at point `alpha`
103 derphi : callable phi'(alpha)
104 Objective function derivative. Returns a scalar.
105 phi0 : float, optional
106 Value of phi at 0
107 old_phi0 : float, optional
108 Value of phi at previous point
109 derphi0 : float, optional
110 Value derphi at 0
111 c1 : float, optional
112 Parameter for Armijo condition rule.
113 c2 : float, optional
114 Parameter for curvature condition rule.
115 amax, amin : float, optional
116 Maximum and minimum step size
117 xtol : float, optional
118 Relative tolerance for an acceptable step.
120 Returns
121 -------
122 alpha : float
123 Step size, or None if no suitable step was found
124 phi : float
125 Value of `phi` at the new point `alpha`
126 phi0 : float
127 Value of `phi` at `alpha=0`
129 Notes
130 -----
131 Uses routine DCSRCH from MINPACK.
133 """
135 if phi0 is None:
136 phi0 = phi(0.)
137 if derphi0 is None:
138 derphi0 = derphi(0.)
140 if old_phi0 is not None and derphi0 != 0:
141 alpha1 = min(1.0, 1.01*2*(phi0 - old_phi0)/derphi0)
142 if alpha1 < 0:
143 alpha1 = 1.0
144 else:
145 alpha1 = 1.0
147 phi1 = phi0
148 derphi1 = derphi0
149 isave = np.zeros((2,), np.intc)
150 dsave = np.zeros((13,), float)
151 task = b'START'
153 maxiter = 100
154 for i in range(maxiter):
155 stp, phi1, derphi1, task = minpack2.dcsrch(alpha1, phi1, derphi1,
156 c1, c2, xtol, task,
157 amin, amax, isave, dsave)
158 if task[:2] == b'FG':
159 alpha1 = stp
160 phi1 = phi(stp)
161 derphi1 = derphi(stp)
162 else:
163 break
164 else:
165 # maxiter reached, the line search did not converge
166 stp = None
168 if task[:5] == b'ERROR' or task[:4] == b'WARN':
169 stp = None # failed
171 return stp, phi1, phi0
174line_search = line_search_wolfe1
177#------------------------------------------------------------------------------
178# Pure-Python Wolfe line and scalar searches
179#------------------------------------------------------------------------------
181# Note: `line_search_wolfe2` is the public `scipy.optimize.line_search`
183def line_search_wolfe2(f, myfprime, xk, pk, gfk=None, old_fval=None,
184 old_old_fval=None, args=(), c1=1e-4, c2=0.9, amax=None,
185 extra_condition=None, maxiter=10):
186 """Find alpha that satisfies strong Wolfe conditions.
188 Parameters
189 ----------
190 f : callable f(x,*args)
191 Objective function.
192 myfprime : callable f'(x,*args)
193 Objective function gradient.
194 xk : ndarray
195 Starting point.
196 pk : ndarray
197 Search direction. The search direction must be a descent direction
198 for the algorithm to converge.
199 gfk : ndarray, optional
200 Gradient value for x=xk (xk being the current parameter
201 estimate). Will be recomputed if omitted.
202 old_fval : float, optional
203 Function value for x=xk. Will be recomputed if omitted.
204 old_old_fval : float, optional
205 Function value for the point preceding x=xk.
206 args : tuple, optional
207 Additional arguments passed to objective function.
208 c1 : float, optional
209 Parameter for Armijo condition rule.
210 c2 : float, optional
211 Parameter for curvature condition rule.
212 amax : float, optional
213 Maximum step size
214 extra_condition : callable, optional
215 A callable of the form ``extra_condition(alpha, x, f, g)``
216 returning a boolean. Arguments are the proposed step ``alpha``
217 and the corresponding ``x``, ``f`` and ``g`` values. The line search
218 accepts the value of ``alpha`` only if this
219 callable returns ``True``. If the callable returns ``False``
220 for the step length, the algorithm will continue with
221 new iterates. The callable is only called for iterates
222 satisfying the strong Wolfe conditions.
223 maxiter : int, optional
224 Maximum number of iterations to perform.
226 Returns
227 -------
228 alpha : float or None
229 Alpha for which ``x_new = x0 + alpha * pk``,
230 or None if the line search algorithm did not converge.
231 fc : int
232 Number of function evaluations made.
233 gc : int
234 Number of gradient evaluations made.
235 new_fval : float or None
236 New function value ``f(x_new)=f(x0+alpha*pk)``,
237 or None if the line search algorithm did not converge.
238 old_fval : float
239 Old function value ``f(x0)``.
240 new_slope : float or None
241 The local slope along the search direction at the
242 new value ``<myfprime(x_new), pk>``,
243 or None if the line search algorithm did not converge.
246 Notes
247 -----
248 Uses the line search algorithm to enforce strong Wolfe
249 conditions. See Wright and Nocedal, 'Numerical Optimization',
250 1999, pp. 59-61.
252 The search direction `pk` must be a descent direction (e.g.
253 ``-myfprime(xk)``) to find a step length that satisfies the strong Wolfe
254 conditions. If the search direction is not a descent direction (e.g.
255 ``myfprime(xk)``), then `alpha`, `new_fval`, and `new_slope` will be None.
257 Examples
258 --------
259 >>> import numpy as np
260 >>> from scipy.optimize import line_search
262 A objective function and its gradient are defined.
264 >>> def obj_func(x):
265 ... return (x[0])**2+(x[1])**2
266 >>> def obj_grad(x):
267 ... return [2*x[0], 2*x[1]]
269 We can find alpha that satisfies strong Wolfe conditions.
271 >>> start_point = np.array([1.8, 1.7])
272 >>> search_gradient = np.array([-1.0, -1.0])
273 >>> line_search(obj_func, obj_grad, start_point, search_gradient)
274 (1.0, 2, 1, 1.1300000000000001, 6.13, [1.6, 1.4])
276 """
277 fc = [0]
278 gc = [0]
279 gval = [None]
280 gval_alpha = [None]
282 def phi(alpha):
283 fc[0] += 1
284 return f(xk + alpha * pk, *args)
286 fprime = myfprime
288 def derphi(alpha):
289 gc[0] += 1
290 gval[0] = fprime(xk + alpha * pk, *args) # store for later use
291 gval_alpha[0] = alpha
292 return np.dot(gval[0], pk)
294 if gfk is None:
295 gfk = fprime(xk, *args)
296 derphi0 = np.dot(gfk, pk)
298 if extra_condition is not None:
299 # Add the current gradient as argument, to avoid needless
300 # re-evaluation
301 def extra_condition2(alpha, phi):
302 if gval_alpha[0] != alpha:
303 derphi(alpha)
304 x = xk + alpha * pk
305 return extra_condition(alpha, x, phi, gval[0])
306 else:
307 extra_condition2 = None
309 alpha_star, phi_star, old_fval, derphi_star = scalar_search_wolfe2(
310 phi, derphi, old_fval, old_old_fval, derphi0, c1, c2, amax,
311 extra_condition2, maxiter=maxiter)
313 if derphi_star is None:
314 warn('The line search algorithm did not converge', LineSearchWarning)
315 else:
316 # derphi_star is a number (derphi) -- so use the most recently
317 # calculated gradient used in computing it derphi = gfk*pk
318 # this is the gradient at the next step no need to compute it
319 # again in the outer loop.
320 derphi_star = gval[0]
322 return alpha_star, fc[0], gc[0], phi_star, old_fval, derphi_star
325def scalar_search_wolfe2(phi, derphi, phi0=None,
326 old_phi0=None, derphi0=None,
327 c1=1e-4, c2=0.9, amax=None,
328 extra_condition=None, maxiter=10):
329 """Find alpha that satisfies strong Wolfe conditions.
331 alpha > 0 is assumed to be a descent direction.
333 Parameters
334 ----------
335 phi : callable phi(alpha)
336 Objective scalar function.
337 derphi : callable phi'(alpha)
338 Objective function derivative. Returns a scalar.
339 phi0 : float, optional
340 Value of phi at 0.
341 old_phi0 : float, optional
342 Value of phi at previous point.
343 derphi0 : float, optional
344 Value of derphi at 0
345 c1 : float, optional
346 Parameter for Armijo condition rule.
347 c2 : float, optional
348 Parameter for curvature condition rule.
349 amax : float, optional
350 Maximum step size.
351 extra_condition : callable, optional
352 A callable of the form ``extra_condition(alpha, phi_value)``
353 returning a boolean. The line search accepts the value
354 of ``alpha`` only if this callable returns ``True``.
355 If the callable returns ``False`` for the step length,
356 the algorithm will continue with new iterates.
357 The callable is only called for iterates satisfying
358 the strong Wolfe conditions.
359 maxiter : int, optional
360 Maximum number of iterations to perform.
362 Returns
363 -------
364 alpha_star : float or None
365 Best alpha, or None if the line search algorithm did not converge.
366 phi_star : float
367 phi at alpha_star.
368 phi0 : float
369 phi at 0.
370 derphi_star : float or None
371 derphi at alpha_star, or None if the line search algorithm
372 did not converge.
374 Notes
375 -----
376 Uses the line search algorithm to enforce strong Wolfe
377 conditions. See Wright and Nocedal, 'Numerical Optimization',
378 1999, pp. 59-61.
380 """
382 if phi0 is None:
383 phi0 = phi(0.)
385 if derphi0 is None:
386 derphi0 = derphi(0.)
388 alpha0 = 0
389 if old_phi0 is not None and derphi0 != 0:
390 alpha1 = min(1.0, 1.01*2*(phi0 - old_phi0)/derphi0)
391 else:
392 alpha1 = 1.0
394 if alpha1 < 0:
395 alpha1 = 1.0
397 if amax is not None:
398 alpha1 = min(alpha1, amax)
400 phi_a1 = phi(alpha1)
401 #derphi_a1 = derphi(alpha1) evaluated below
403 phi_a0 = phi0
404 derphi_a0 = derphi0
406 if extra_condition is None:
407 def extra_condition(alpha, phi):
408 return True
410 for i in range(maxiter):
411 if alpha1 == 0 or (amax is not None and alpha0 == amax):
412 # alpha1 == 0: This shouldn't happen. Perhaps the increment has
413 # slipped below machine precision?
414 alpha_star = None
415 phi_star = phi0
416 phi0 = old_phi0
417 derphi_star = None
419 if alpha1 == 0:
420 msg = 'Rounding errors prevent the line search from converging'
421 else:
422 msg = "The line search algorithm could not find a solution " + \
423 "less than or equal to amax: %s" % amax
425 warn(msg, LineSearchWarning)
426 break
428 not_first_iteration = i > 0
429 if (phi_a1 > phi0 + c1 * alpha1 * derphi0) or \
430 ((phi_a1 >= phi_a0) and not_first_iteration):
431 alpha_star, phi_star, derphi_star = \
432 _zoom(alpha0, alpha1, phi_a0,
433 phi_a1, derphi_a0, phi, derphi,
434 phi0, derphi0, c1, c2, extra_condition)
435 break
437 derphi_a1 = derphi(alpha1)
438 if (abs(derphi_a1) <= -c2*derphi0):
439 if extra_condition(alpha1, phi_a1):
440 alpha_star = alpha1
441 phi_star = phi_a1
442 derphi_star = derphi_a1
443 break
445 if (derphi_a1 >= 0):
446 alpha_star, phi_star, derphi_star = \
447 _zoom(alpha1, alpha0, phi_a1,
448 phi_a0, derphi_a1, phi, derphi,
449 phi0, derphi0, c1, c2, extra_condition)
450 break
452 alpha2 = 2 * alpha1 # increase by factor of two on each iteration
453 if amax is not None:
454 alpha2 = min(alpha2, amax)
455 alpha0 = alpha1
456 alpha1 = alpha2
457 phi_a0 = phi_a1
458 phi_a1 = phi(alpha1)
459 derphi_a0 = derphi_a1
461 else:
462 # stopping test maxiter reached
463 alpha_star = alpha1
464 phi_star = phi_a1
465 derphi_star = None
466 warn('The line search algorithm did not converge', LineSearchWarning)
468 return alpha_star, phi_star, phi0, derphi_star
471def _cubicmin(a, fa, fpa, b, fb, c, fc):
472 """
473 Finds the minimizer for a cubic polynomial that goes through the
474 points (a,fa), (b,fb), and (c,fc) with derivative at a of fpa.
476 If no minimizer can be found, return None.
478 """
479 # f(x) = A *(x-a)^3 + B*(x-a)^2 + C*(x-a) + D
481 with np.errstate(divide='raise', over='raise', invalid='raise'):
482 try:
483 C = fpa
484 db = b - a
485 dc = c - a
486 denom = (db * dc) ** 2 * (db - dc)
487 d1 = np.empty((2, 2))
488 d1[0, 0] = dc ** 2
489 d1[0, 1] = -db ** 2
490 d1[1, 0] = -dc ** 3
491 d1[1, 1] = db ** 3
492 [A, B] = np.dot(d1, np.asarray([fb - fa - C * db,
493 fc - fa - C * dc]).flatten())
494 A /= denom
495 B /= denom
496 radical = B * B - 3 * A * C
497 xmin = a + (-B + np.sqrt(radical)) / (3 * A)
498 except ArithmeticError:
499 return None
500 if not np.isfinite(xmin):
501 return None
502 return xmin
505def _quadmin(a, fa, fpa, b, fb):
506 """
507 Finds the minimizer for a quadratic polynomial that goes through
508 the points (a,fa), (b,fb) with derivative at a of fpa.
510 """
511 # f(x) = B*(x-a)^2 + C*(x-a) + D
512 with np.errstate(divide='raise', over='raise', invalid='raise'):
513 try:
514 D = fa
515 C = fpa
516 db = b - a * 1.0
517 B = (fb - D - C * db) / (db * db)
518 xmin = a - C / (2.0 * B)
519 except ArithmeticError:
520 return None
521 if not np.isfinite(xmin):
522 return None
523 return xmin
526def _zoom(a_lo, a_hi, phi_lo, phi_hi, derphi_lo,
527 phi, derphi, phi0, derphi0, c1, c2, extra_condition):
528 """Zoom stage of approximate linesearch satisfying strong Wolfe conditions.
530 Part of the optimization algorithm in `scalar_search_wolfe2`.
532 Notes
533 -----
534 Implements Algorithm 3.6 (zoom) in Wright and Nocedal,
535 'Numerical Optimization', 1999, pp. 61.
537 """
539 maxiter = 10
540 i = 0
541 delta1 = 0.2 # cubic interpolant check
542 delta2 = 0.1 # quadratic interpolant check
543 phi_rec = phi0
544 a_rec = 0
545 while True:
546 # interpolate to find a trial step length between a_lo and
547 # a_hi Need to choose interpolation here. Use cubic
548 # interpolation and then if the result is within delta *
549 # dalpha or outside of the interval bounded by a_lo or a_hi
550 # then use quadratic interpolation, if the result is still too
551 # close, then use bisection
553 dalpha = a_hi - a_lo
554 if dalpha < 0:
555 a, b = a_hi, a_lo
556 else:
557 a, b = a_lo, a_hi
559 # minimizer of cubic interpolant
560 # (uses phi_lo, derphi_lo, phi_hi, and the most recent value of phi)
561 #
562 # if the result is too close to the end points (or out of the
563 # interval), then use quadratic interpolation with phi_lo,
564 # derphi_lo and phi_hi if the result is still too close to the
565 # end points (or out of the interval) then use bisection
567 if (i > 0):
568 cchk = delta1 * dalpha
569 a_j = _cubicmin(a_lo, phi_lo, derphi_lo, a_hi, phi_hi,
570 a_rec, phi_rec)
571 if (i == 0) or (a_j is None) or (a_j > b - cchk) or (a_j < a + cchk):
572 qchk = delta2 * dalpha
573 a_j = _quadmin(a_lo, phi_lo, derphi_lo, a_hi, phi_hi)
574 if (a_j is None) or (a_j > b-qchk) or (a_j < a+qchk):
575 a_j = a_lo + 0.5*dalpha
577 # Check new value of a_j
579 phi_aj = phi(a_j)
580 if (phi_aj > phi0 + c1*a_j*derphi0) or (phi_aj >= phi_lo):
581 phi_rec = phi_hi
582 a_rec = a_hi
583 a_hi = a_j
584 phi_hi = phi_aj
585 else:
586 derphi_aj = derphi(a_j)
587 if abs(derphi_aj) <= -c2*derphi0 and extra_condition(a_j, phi_aj):
588 a_star = a_j
589 val_star = phi_aj
590 valprime_star = derphi_aj
591 break
592 if derphi_aj*(a_hi - a_lo) >= 0:
593 phi_rec = phi_hi
594 a_rec = a_hi
595 a_hi = a_lo
596 phi_hi = phi_lo
597 else:
598 phi_rec = phi_lo
599 a_rec = a_lo
600 a_lo = a_j
601 phi_lo = phi_aj
602 derphi_lo = derphi_aj
603 i += 1
604 if (i > maxiter):
605 # Failed to find a conforming step size
606 a_star = None
607 val_star = None
608 valprime_star = None
609 break
610 return a_star, val_star, valprime_star
613#------------------------------------------------------------------------------
614# Armijo line and scalar searches
615#------------------------------------------------------------------------------
617def line_search_armijo(f, xk, pk, gfk, old_fval, args=(), c1=1e-4, alpha0=1):
618 """Minimize over alpha, the function ``f(xk+alpha pk)``.
620 Parameters
621 ----------
622 f : callable
623 Function to be minimized.
624 xk : array_like
625 Current point.
626 pk : array_like
627 Search direction.
628 gfk : array_like
629 Gradient of `f` at point `xk`.
630 old_fval : float
631 Value of `f` at point `xk`.
632 args : tuple, optional
633 Optional arguments.
634 c1 : float, optional
635 Value to control stopping criterion.
636 alpha0 : scalar, optional
637 Value of `alpha` at start of the optimization.
639 Returns
640 -------
641 alpha
642 f_count
643 f_val_at_alpha
645 Notes
646 -----
647 Uses the interpolation algorithm (Armijo backtracking) as suggested by
648 Wright and Nocedal in 'Numerical Optimization', 1999, pp. 56-57
650 """
651 xk = np.atleast_1d(xk)
652 fc = [0]
654 def phi(alpha1):
655 fc[0] += 1
656 return f(xk + alpha1*pk, *args)
658 if old_fval is None:
659 phi0 = phi(0.)
660 else:
661 phi0 = old_fval # compute f(xk) -- done in past loop
663 derphi0 = np.dot(gfk, pk)
664 alpha, phi1 = scalar_search_armijo(phi, phi0, derphi0, c1=c1,
665 alpha0=alpha0)
666 return alpha, fc[0], phi1
669def line_search_BFGS(f, xk, pk, gfk, old_fval, args=(), c1=1e-4, alpha0=1):
670 """
671 Compatibility wrapper for `line_search_armijo`
672 """
673 r = line_search_armijo(f, xk, pk, gfk, old_fval, args=args, c1=c1,
674 alpha0=alpha0)
675 return r[0], r[1], 0, r[2]
678def scalar_search_armijo(phi, phi0, derphi0, c1=1e-4, alpha0=1, amin=0):
679 """Minimize over alpha, the function ``phi(alpha)``.
681 Uses the interpolation algorithm (Armijo backtracking) as suggested by
682 Wright and Nocedal in 'Numerical Optimization', 1999, pp. 56-57
684 alpha > 0 is assumed to be a descent direction.
686 Returns
687 -------
688 alpha
689 phi1
691 """
692 phi_a0 = phi(alpha0)
693 if phi_a0 <= phi0 + c1*alpha0*derphi0:
694 return alpha0, phi_a0
696 # Otherwise, compute the minimizer of a quadratic interpolant:
698 alpha1 = -(derphi0) * alpha0**2 / 2.0 / (phi_a0 - phi0 - derphi0 * alpha0)
699 phi_a1 = phi(alpha1)
701 if (phi_a1 <= phi0 + c1*alpha1*derphi0):
702 return alpha1, phi_a1
704 # Otherwise, loop with cubic interpolation until we find an alpha which
705 # satisfies the first Wolfe condition (since we are backtracking, we will
706 # assume that the value of alpha is not too small and satisfies the second
707 # condition.
709 while alpha1 > amin: # we are assuming alpha>0 is a descent direction
710 factor = alpha0**2 * alpha1**2 * (alpha1-alpha0)
711 a = alpha0**2 * (phi_a1 - phi0 - derphi0*alpha1) - \
712 alpha1**2 * (phi_a0 - phi0 - derphi0*alpha0)
713 a = a / factor
714 b = -alpha0**3 * (phi_a1 - phi0 - derphi0*alpha1) + \
715 alpha1**3 * (phi_a0 - phi0 - derphi0*alpha0)
716 b = b / factor
718 alpha2 = (-b + np.sqrt(abs(b**2 - 3 * a * derphi0))) / (3.0*a)
719 phi_a2 = phi(alpha2)
721 if (phi_a2 <= phi0 + c1*alpha2*derphi0):
722 return alpha2, phi_a2
724 if (alpha1 - alpha2) > alpha1 / 2.0 or (1 - alpha2/alpha1) < 0.96:
725 alpha2 = alpha1 / 2.0
727 alpha0 = alpha1
728 alpha1 = alpha2
729 phi_a0 = phi_a1
730 phi_a1 = phi_a2
732 # Failed to find a suitable step length
733 return None, phi_a1
736#------------------------------------------------------------------------------
737# Non-monotone line search for DF-SANE
738#------------------------------------------------------------------------------
740def _nonmonotone_line_search_cruz(f, x_k, d, prev_fs, eta,
741 gamma=1e-4, tau_min=0.1, tau_max=0.5):
742 """
743 Nonmonotone backtracking line search as described in [1]_
745 Parameters
746 ----------
747 f : callable
748 Function returning a tuple ``(f, F)`` where ``f`` is the value
749 of a merit function and ``F`` the residual.
750 x_k : ndarray
751 Initial position.
752 d : ndarray
753 Search direction.
754 prev_fs : float
755 List of previous merit function values. Should have ``len(prev_fs) <= M``
756 where ``M`` is the nonmonotonicity window parameter.
757 eta : float
758 Allowed merit function increase, see [1]_
759 gamma, tau_min, tau_max : float, optional
760 Search parameters, see [1]_
762 Returns
763 -------
764 alpha : float
765 Step length
766 xp : ndarray
767 Next position
768 fp : float
769 Merit function value at next position
770 Fp : ndarray
771 Residual at next position
773 References
774 ----------
775 [1] "Spectral residual method without gradient information for solving
776 large-scale nonlinear systems of equations." W. La Cruz,
777 J.M. Martinez, M. Raydan. Math. Comp. **75**, 1429 (2006).
779 """
780 f_k = prev_fs[-1]
781 f_bar = max(prev_fs)
783 alpha_p = 1
784 alpha_m = 1
785 alpha = 1
787 while True:
788 xp = x_k + alpha_p * d
789 fp, Fp = f(xp)
791 if fp <= f_bar + eta - gamma * alpha_p**2 * f_k:
792 alpha = alpha_p
793 break
795 alpha_tp = alpha_p**2 * f_k / (fp + (2*alpha_p - 1)*f_k)
797 xp = x_k - alpha_m * d
798 fp, Fp = f(xp)
800 if fp <= f_bar + eta - gamma * alpha_m**2 * f_k:
801 alpha = -alpha_m
802 break
804 alpha_tm = alpha_m**2 * f_k / (fp + (2*alpha_m - 1)*f_k)
806 alpha_p = np.clip(alpha_tp, tau_min * alpha_p, tau_max * alpha_p)
807 alpha_m = np.clip(alpha_tm, tau_min * alpha_m, tau_max * alpha_m)
809 return alpha, xp, fp, Fp
812def _nonmonotone_line_search_cheng(f, x_k, d, f_k, C, Q, eta,
813 gamma=1e-4, tau_min=0.1, tau_max=0.5,
814 nu=0.85):
815 """
816 Nonmonotone line search from [1]
818 Parameters
819 ----------
820 f : callable
821 Function returning a tuple ``(f, F)`` where ``f`` is the value
822 of a merit function and ``F`` the residual.
823 x_k : ndarray
824 Initial position.
825 d : ndarray
826 Search direction.
827 f_k : float
828 Initial merit function value.
829 C, Q : float
830 Control parameters. On the first iteration, give values
831 Q=1.0, C=f_k
832 eta : float
833 Allowed merit function increase, see [1]_
834 nu, gamma, tau_min, tau_max : float, optional
835 Search parameters, see [1]_
837 Returns
838 -------
839 alpha : float
840 Step length
841 xp : ndarray
842 Next position
843 fp : float
844 Merit function value at next position
845 Fp : ndarray
846 Residual at next position
847 C : float
848 New value for the control parameter C
849 Q : float
850 New value for the control parameter Q
852 References
853 ----------
854 .. [1] W. Cheng & D.-H. Li, ''A derivative-free nonmonotone line
855 search and its application to the spectral residual
856 method'', IMA J. Numer. Anal. 29, 814 (2009).
858 """
859 alpha_p = 1
860 alpha_m = 1
861 alpha = 1
863 while True:
864 xp = x_k + alpha_p * d
865 fp, Fp = f(xp)
867 if fp <= C + eta - gamma * alpha_p**2 * f_k:
868 alpha = alpha_p
869 break
871 alpha_tp = alpha_p**2 * f_k / (fp + (2*alpha_p - 1)*f_k)
873 xp = x_k - alpha_m * d
874 fp, Fp = f(xp)
876 if fp <= C + eta - gamma * alpha_m**2 * f_k:
877 alpha = -alpha_m
878 break
880 alpha_tm = alpha_m**2 * f_k / (fp + (2*alpha_m - 1)*f_k)
882 alpha_p = np.clip(alpha_tp, tau_min * alpha_p, tau_max * alpha_p)
883 alpha_m = np.clip(alpha_tm, tau_min * alpha_m, tau_max * alpha_m)
885 # Update C and Q
886 Q_next = nu * Q + 1
887 C = (nu * Q * (C + eta) + fp) / Q_next
888 Q = Q_next
890 return alpha, xp, fp, Fp, C, Q