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