Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/optimize/_trustregion.py: 16%
154 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
1"""Trust-region optimization."""
2import math
3import warnings
5import numpy as np
6import scipy.linalg
7from ._optimize import (_check_unknown_options, _status_message,
8 OptimizeResult, _prepare_scalar_function)
9from scipy.optimize._hessian_update_strategy import HessianUpdateStrategy
10from scipy.optimize._differentiable_functions import FD_METHODS
11__all__ = []
14def _wrap_function(function, args):
15 # wraps a minimizer function to count number of evaluations
16 # and to easily provide an args kwd.
17 ncalls = [0]
18 if function is None:
19 return ncalls, None
21 def function_wrapper(x, *wrapper_args):
22 ncalls[0] += 1
23 # A copy of x is sent to the user function (gh13740)
24 return function(np.copy(x), *(wrapper_args + args))
26 return ncalls, function_wrapper
29class BaseQuadraticSubproblem:
30 """
31 Base/abstract class defining the quadratic model for trust-region
32 minimization. Child classes must implement the ``solve`` method.
34 Values of the objective function, Jacobian and Hessian (if provided) at
35 the current iterate ``x`` are evaluated on demand and then stored as
36 attributes ``fun``, ``jac``, ``hess``.
37 """
39 def __init__(self, x, fun, jac, hess=None, hessp=None):
40 self._x = x
41 self._f = None
42 self._g = None
43 self._h = None
44 self._g_mag = None
45 self._cauchy_point = None
46 self._newton_point = None
47 self._fun = fun
48 self._jac = jac
49 self._hess = hess
50 self._hessp = hessp
52 def __call__(self, p):
53 return self.fun + np.dot(self.jac, p) + 0.5 * np.dot(p, self.hessp(p))
55 @property
56 def fun(self):
57 """Value of objective function at current iteration."""
58 if self._f is None:
59 self._f = self._fun(self._x)
60 return self._f
62 @property
63 def jac(self):
64 """Value of Jacobian of objective function at current iteration."""
65 if self._g is None:
66 self._g = self._jac(self._x)
67 return self._g
69 @property
70 def hess(self):
71 """Value of Hessian of objective function at current iteration."""
72 if self._h is None:
73 self._h = self._hess(self._x)
74 return self._h
76 def hessp(self, p):
77 if self._hessp is not None:
78 return self._hessp(self._x, p)
79 else:
80 return np.dot(self.hess, p)
82 @property
83 def jac_mag(self):
84 """Magnitude of jacobian of objective function at current iteration."""
85 if self._g_mag is None:
86 self._g_mag = scipy.linalg.norm(self.jac)
87 return self._g_mag
89 def get_boundaries_intersections(self, z, d, trust_radius):
90 """
91 Solve the scalar quadratic equation ||z + t d|| == trust_radius.
92 This is like a line-sphere intersection.
93 Return the two values of t, sorted from low to high.
94 """
95 a = np.dot(d, d)
96 b = 2 * np.dot(z, d)
97 c = np.dot(z, z) - trust_radius**2
98 sqrt_discriminant = math.sqrt(b*b - 4*a*c)
100 # The following calculation is mathematically
101 # equivalent to:
102 # ta = (-b - sqrt_discriminant) / (2*a)
103 # tb = (-b + sqrt_discriminant) / (2*a)
104 # but produce smaller round off errors.
105 # Look at Matrix Computation p.97
106 # for a better justification.
107 aux = b + math.copysign(sqrt_discriminant, b)
108 ta = -aux / (2*a)
109 tb = -2*c / aux
110 return sorted([ta, tb])
112 def solve(self, trust_radius):
113 raise NotImplementedError('The solve method should be implemented by '
114 'the child class')
117def _minimize_trust_region(fun, x0, args=(), jac=None, hess=None, hessp=None,
118 subproblem=None, initial_trust_radius=1.0,
119 max_trust_radius=1000.0, eta=0.15, gtol=1e-4,
120 maxiter=None, disp=False, return_all=False,
121 callback=None, inexact=True, **unknown_options):
122 """
123 Minimization of scalar function of one or more variables using a
124 trust-region algorithm.
126 Options for the trust-region algorithm are:
127 initial_trust_radius : float
128 Initial trust radius.
129 max_trust_radius : float
130 Never propose steps that are longer than this value.
131 eta : float
132 Trust region related acceptance stringency for proposed steps.
133 gtol : float
134 Gradient norm must be less than `gtol`
135 before successful termination.
136 maxiter : int
137 Maximum number of iterations to perform.
138 disp : bool
139 If True, print convergence message.
140 inexact : bool
141 Accuracy to solve subproblems. If True requires less nonlinear
142 iterations, but more vector products. Only effective for method
143 trust-krylov.
145 This function is called by the `minimize` function.
146 It is not supposed to be called directly.
147 """
148 _check_unknown_options(unknown_options)
150 if jac is None:
151 raise ValueError('Jacobian is currently required for trust-region '
152 'methods')
153 if hess is None and hessp is None:
154 raise ValueError('Either the Hessian or the Hessian-vector product '
155 'is currently required for trust-region methods')
156 if subproblem is None:
157 raise ValueError('A subproblem solving strategy is required for '
158 'trust-region methods')
159 if not (0 <= eta < 0.25):
160 raise Exception('invalid acceptance stringency')
161 if max_trust_radius <= 0:
162 raise Exception('the max trust radius must be positive')
163 if initial_trust_radius <= 0:
164 raise ValueError('the initial trust radius must be positive')
165 if initial_trust_radius >= max_trust_radius:
166 raise ValueError('the initial trust radius must be less than the '
167 'max trust radius')
169 # force the initial guess into a nice format
170 x0 = np.asarray(x0).flatten()
172 # A ScalarFunction representing the problem. This caches calls to fun, jac,
173 # hess.
174 sf = _prepare_scalar_function(fun, x0, jac=jac, hess=hess, args=args)
175 fun = sf.fun
176 jac = sf.grad
177 if callable(hess):
178 hess = sf.hess
179 elif callable(hessp):
180 # this elif statement must come before examining whether hess
181 # is estimated by FD methods or a HessianUpdateStrategy
182 pass
183 elif (hess in FD_METHODS or isinstance(hess, HessianUpdateStrategy)):
184 # If the Hessian is being estimated by finite differences or a
185 # Hessian update strategy then ScalarFunction.hess returns a
186 # LinearOperator or a HessianUpdateStrategy. This enables the
187 # calculation/creation of a hessp. BUT you only want to do this
188 # if the user *hasn't* provided a callable(hessp) function.
189 hess = None
191 def hessp(x, p, *args):
192 return sf.hess(x).dot(p)
193 else:
194 raise ValueError('Either the Hessian or the Hessian-vector product '
195 'is currently required for trust-region methods')
197 # ScalarFunction doesn't represent hessp
198 nhessp, hessp = _wrap_function(hessp, args)
200 # limit the number of iterations
201 if maxiter is None:
202 maxiter = len(x0)*200
204 # init the search status
205 warnflag = 0
207 # initialize the search
208 trust_radius = initial_trust_radius
209 x = x0
210 if return_all:
211 allvecs = [x]
212 m = subproblem(x, fun, jac, hess, hessp)
213 k = 0
215 # search for the function min
216 # do not even start if the gradient is small enough
217 while m.jac_mag >= gtol:
219 # Solve the sub-problem.
220 # This gives us the proposed step relative to the current position
221 # and it tells us whether the proposed step
222 # has reached the trust region boundary or not.
223 try:
224 p, hits_boundary = m.solve(trust_radius)
225 except np.linalg.LinAlgError:
226 warnflag = 3
227 break
229 # calculate the predicted value at the proposed point
230 predicted_value = m(p)
232 # define the local approximation at the proposed point
233 x_proposed = x + p
234 m_proposed = subproblem(x_proposed, fun, jac, hess, hessp)
236 # evaluate the ratio defined in equation (4.4)
237 actual_reduction = m.fun - m_proposed.fun
238 predicted_reduction = m.fun - predicted_value
239 if predicted_reduction <= 0:
240 warnflag = 2
241 break
242 rho = actual_reduction / predicted_reduction
244 # update the trust radius according to the actual/predicted ratio
245 if rho < 0.25:
246 trust_radius *= 0.25
247 elif rho > 0.75 and hits_boundary:
248 trust_radius = min(2*trust_radius, max_trust_radius)
250 # if the ratio is high enough then accept the proposed step
251 if rho > eta:
252 x = x_proposed
253 m = m_proposed
255 # append the best guess, call back, increment the iteration count
256 if return_all:
257 allvecs.append(np.copy(x))
258 if callback is not None:
259 callback(np.copy(x))
260 k += 1
262 # check if the gradient is small enough to stop
263 if m.jac_mag < gtol:
264 warnflag = 0
265 break
267 # check if we have looked at enough iterations
268 if k >= maxiter:
269 warnflag = 1
270 break
272 # print some stuff if requested
273 status_messages = (
274 _status_message['success'],
275 _status_message['maxiter'],
276 'A bad approximation caused failure to predict improvement.',
277 'A linalg error occurred, such as a non-psd Hessian.',
278 )
279 if disp:
280 if warnflag == 0:
281 print(status_messages[warnflag])
282 else:
283 warnings.warn(status_messages[warnflag], RuntimeWarning, 3)
284 print(" Current function value: %f" % m.fun)
285 print(" Iterations: %d" % k)
286 print(" Function evaluations: %d" % sf.nfev)
287 print(" Gradient evaluations: %d" % sf.ngev)
288 print(" Hessian evaluations: %d" % (sf.nhev + nhessp[0]))
290 result = OptimizeResult(x=x, success=(warnflag == 0), status=warnflag,
291 fun=m.fun, jac=m.jac, nfev=sf.nfev, njev=sf.ngev,
292 nhev=sf.nhev + nhessp[0], nit=k,
293 message=status_messages[warnflag])
295 if hess is not None:
296 result['hess'] = m.hess
298 if return_all:
299 result['allvecs'] = allvecs
301 return result