Coverage for /usr/lib/python3/dist-packages/scipy/optimize/_trustregion_constr/minimize_trustregion_constr.py: 11%
193 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1import time
2import numpy as np
3from scipy.sparse.linalg import LinearOperator
4from .._differentiable_functions import VectorFunction
5from .._constraints import (
6 NonlinearConstraint, LinearConstraint, PreparedConstraint, Bounds, strict_bounds)
7from .._hessian_update_strategy import BFGS
8from .._optimize import OptimizeResult
9from .._differentiable_functions import ScalarFunction
10from .equality_constrained_sqp import equality_constrained_sqp
11from .canonical_constraint import (CanonicalConstraint,
12 initial_constraints_as_canonical)
13from .tr_interior_point import tr_interior_point
14from .report import BasicReport, SQPReport, IPReport
17TERMINATION_MESSAGES = {
18 0: "The maximum number of function evaluations is exceeded.",
19 1: "`gtol` termination condition is satisfied.",
20 2: "`xtol` termination condition is satisfied.",
21 3: "`callback` function requested termination."
22}
25class HessianLinearOperator:
26 """Build LinearOperator from hessp"""
27 def __init__(self, hessp, n):
28 self.hessp = hessp
29 self.n = n
31 def __call__(self, x, *args):
32 def matvec(p):
33 return self.hessp(x, p, *args)
35 return LinearOperator((self.n, self.n), matvec=matvec)
38class LagrangianHessian:
39 """The Hessian of the Lagrangian as LinearOperator.
41 The Lagrangian is computed as the objective function plus all the
42 constraints multiplied with some numbers (Lagrange multipliers).
43 """
44 def __init__(self, n, objective_hess, constraints_hess):
45 self.n = n
46 self.objective_hess = objective_hess
47 self.constraints_hess = constraints_hess
49 def __call__(self, x, v_eq=np.empty(0), v_ineq=np.empty(0)):
50 H_objective = self.objective_hess(x)
51 H_constraints = self.constraints_hess(x, v_eq, v_ineq)
53 def matvec(p):
54 return H_objective.dot(p) + H_constraints.dot(p)
56 return LinearOperator((self.n, self.n), matvec)
59def update_state_sqp(state, x, last_iteration_failed, objective, prepared_constraints,
60 start_time, tr_radius, constr_penalty, cg_info):
61 state.nit += 1
62 state.nfev = objective.nfev
63 state.njev = objective.ngev
64 state.nhev = objective.nhev
65 state.constr_nfev = [c.fun.nfev if isinstance(c.fun, VectorFunction) else 0
66 for c in prepared_constraints]
67 state.constr_njev = [c.fun.njev if isinstance(c.fun, VectorFunction) else 0
68 for c in prepared_constraints]
69 state.constr_nhev = [c.fun.nhev if isinstance(c.fun, VectorFunction) else 0
70 for c in prepared_constraints]
72 if not last_iteration_failed:
73 state.x = x
74 state.fun = objective.f
75 state.grad = objective.g
76 state.v = [c.fun.v for c in prepared_constraints]
77 state.constr = [c.fun.f for c in prepared_constraints]
78 state.jac = [c.fun.J for c in prepared_constraints]
79 # Compute Lagrangian Gradient
80 state.lagrangian_grad = np.copy(state.grad)
81 for c in prepared_constraints:
82 state.lagrangian_grad += c.fun.J.T.dot(c.fun.v)
83 state.optimality = np.linalg.norm(state.lagrangian_grad, np.inf)
84 # Compute maximum constraint violation
85 state.constr_violation = 0
86 for i in range(len(prepared_constraints)):
87 lb, ub = prepared_constraints[i].bounds
88 c = state.constr[i]
89 state.constr_violation = np.max([state.constr_violation,
90 np.max(lb - c),
91 np.max(c - ub)])
93 state.execution_time = time.time() - start_time
94 state.tr_radius = tr_radius
95 state.constr_penalty = constr_penalty
96 state.cg_niter += cg_info["niter"]
97 state.cg_stop_cond = cg_info["stop_cond"]
99 return state
102def update_state_ip(state, x, last_iteration_failed, objective,
103 prepared_constraints, start_time,
104 tr_radius, constr_penalty, cg_info,
105 barrier_parameter, barrier_tolerance):
106 state = update_state_sqp(state, x, last_iteration_failed, objective,
107 prepared_constraints, start_time, tr_radius,
108 constr_penalty, cg_info)
109 state.barrier_parameter = barrier_parameter
110 state.barrier_tolerance = barrier_tolerance
111 return state
114def _minimize_trustregion_constr(fun, x0, args, grad,
115 hess, hessp, bounds, constraints,
116 xtol=1e-8, gtol=1e-8,
117 barrier_tol=1e-8,
118 sparse_jacobian=None,
119 callback=None, maxiter=1000,
120 verbose=0, finite_diff_rel_step=None,
121 initial_constr_penalty=1.0, initial_tr_radius=1.0,
122 initial_barrier_parameter=0.1,
123 initial_barrier_tolerance=0.1,
124 factorization_method=None,
125 disp=False):
126 """Minimize a scalar function subject to constraints.
128 Parameters
129 ----------
130 gtol : float, optional
131 Tolerance for termination by the norm of the Lagrangian gradient.
132 The algorithm will terminate when both the infinity norm (i.e., max
133 abs value) of the Lagrangian gradient and the constraint violation
134 are smaller than ``gtol``. Default is 1e-8.
135 xtol : float, optional
136 Tolerance for termination by the change of the independent variable.
137 The algorithm will terminate when ``tr_radius < xtol``, where
138 ``tr_radius`` is the radius of the trust region used in the algorithm.
139 Default is 1e-8.
140 barrier_tol : float, optional
141 Threshold on the barrier parameter for the algorithm termination.
142 When inequality constraints are present, the algorithm will terminate
143 only when the barrier parameter is less than `barrier_tol`.
144 Default is 1e-8.
145 sparse_jacobian : {bool, None}, optional
146 Determines how to represent Jacobians of the constraints. If bool,
147 then Jacobians of all the constraints will be converted to the
148 corresponding format. If None (default), then Jacobians won't be
149 converted, but the algorithm can proceed only if they all have the
150 same format.
151 initial_tr_radius: float, optional
152 Initial trust radius. The trust radius gives the maximum distance
153 between solution points in consecutive iterations. It reflects the
154 trust the algorithm puts in the local approximation of the optimization
155 problem. For an accurate local approximation the trust-region should be
156 large and for an approximation valid only close to the current point it
157 should be a small one. The trust radius is automatically updated throughout
158 the optimization process, with ``initial_tr_radius`` being its initial value.
159 Default is 1 (recommended in [1]_, p. 19).
160 initial_constr_penalty : float, optional
161 Initial constraints penalty parameter. The penalty parameter is used for
162 balancing the requirements of decreasing the objective function
163 and satisfying the constraints. It is used for defining the merit function:
164 ``merit_function(x) = fun(x) + constr_penalty * constr_norm_l2(x)``,
165 where ``constr_norm_l2(x)`` is the l2 norm of a vector containing all
166 the constraints. The merit function is used for accepting or rejecting
167 trial points and ``constr_penalty`` weights the two conflicting goals
168 of reducing objective function and constraints. The penalty is automatically
169 updated throughout the optimization process, with
170 ``initial_constr_penalty`` being its initial value. Default is 1
171 (recommended in [1]_, p 19).
172 initial_barrier_parameter, initial_barrier_tolerance: float, optional
173 Initial barrier parameter and initial tolerance for the barrier subproblem.
174 Both are used only when inequality constraints are present. For dealing with
175 optimization problems ``min_x f(x)`` subject to inequality constraints
176 ``c(x) <= 0`` the algorithm introduces slack variables, solving the problem
177 ``min_(x,s) f(x) + barrier_parameter*sum(ln(s))`` subject to the equality
178 constraints ``c(x) + s = 0`` instead of the original problem. This subproblem
179 is solved for decreasing values of ``barrier_parameter`` and with decreasing
180 tolerances for the termination, starting with ``initial_barrier_parameter``
181 for the barrier parameter and ``initial_barrier_tolerance`` for the
182 barrier tolerance. Default is 0.1 for both values (recommended in [1]_ p. 19).
183 Also note that ``barrier_parameter`` and ``barrier_tolerance`` are updated
184 with the same prefactor.
185 factorization_method : string or None, optional
186 Method to factorize the Jacobian of the constraints. Use None (default)
187 for the auto selection or one of:
189 - 'NormalEquation' (requires scikit-sparse)
190 - 'AugmentedSystem'
191 - 'QRFactorization'
192 - 'SVDFactorization'
194 The methods 'NormalEquation' and 'AugmentedSystem' can be used only
195 with sparse constraints. The projections required by the algorithm
196 will be computed using, respectively, the normal equation and the
197 augmented system approaches explained in [1]_. 'NormalEquation'
198 computes the Cholesky factorization of ``A A.T`` and 'AugmentedSystem'
199 performs the LU factorization of an augmented system. They usually
200 provide similar results. 'AugmentedSystem' is used by default for
201 sparse matrices.
203 The methods 'QRFactorization' and 'SVDFactorization' can be used
204 only with dense constraints. They compute the required projections
205 using, respectively, QR and SVD factorizations. The 'SVDFactorization'
206 method can cope with Jacobian matrices with deficient row rank and will
207 be used whenever other factorization methods fail (which may imply the
208 conversion of sparse matrices to a dense format when required).
209 By default, 'QRFactorization' is used for dense matrices.
210 finite_diff_rel_step : None or array_like, optional
211 Relative step size for the finite difference approximation.
212 maxiter : int, optional
213 Maximum number of algorithm iterations. Default is 1000.
214 verbose : {0, 1, 2}, optional
215 Level of algorithm's verbosity:
217 * 0 (default) : work silently.
218 * 1 : display a termination report.
219 * 2 : display progress during iterations.
220 * 3 : display progress during iterations (more complete report).
222 disp : bool, optional
223 If True (default), then `verbose` will be set to 1 if it was 0.
225 Returns
226 -------
227 `OptimizeResult` with the fields documented below. Note the following:
229 1. All values corresponding to the constraints are ordered as they
230 were passed to the solver. And values corresponding to `bounds`
231 constraints are put *after* other constraints.
232 2. All numbers of function, Jacobian or Hessian evaluations correspond
233 to numbers of actual Python function calls. It means, for example,
234 that if a Jacobian is estimated by finite differences, then the
235 number of Jacobian evaluations will be zero and the number of
236 function evaluations will be incremented by all calls during the
237 finite difference estimation.
239 x : ndarray, shape (n,)
240 Solution found.
241 optimality : float
242 Infinity norm of the Lagrangian gradient at the solution.
243 constr_violation : float
244 Maximum constraint violation at the solution.
245 fun : float
246 Objective function at the solution.
247 grad : ndarray, shape (n,)
248 Gradient of the objective function at the solution.
249 lagrangian_grad : ndarray, shape (n,)
250 Gradient of the Lagrangian function at the solution.
251 nit : int
252 Total number of iterations.
253 nfev : integer
254 Number of the objective function evaluations.
255 njev : integer
256 Number of the objective function gradient evaluations.
257 nhev : integer
258 Number of the objective function Hessian evaluations.
259 cg_niter : int
260 Total number of the conjugate gradient method iterations.
261 method : {'equality_constrained_sqp', 'tr_interior_point'}
262 Optimization method used.
263 constr : list of ndarray
264 List of constraint values at the solution.
265 jac : list of {ndarray, sparse matrix}
266 List of the Jacobian matrices of the constraints at the solution.
267 v : list of ndarray
268 List of the Lagrange multipliers for the constraints at the solution.
269 For an inequality constraint a positive multiplier means that the upper
270 bound is active, a negative multiplier means that the lower bound is
271 active and if a multiplier is zero it means the constraint is not
272 active.
273 constr_nfev : list of int
274 Number of constraint evaluations for each of the constraints.
275 constr_njev : list of int
276 Number of Jacobian matrix evaluations for each of the constraints.
277 constr_nhev : list of int
278 Number of Hessian evaluations for each of the constraints.
279 tr_radius : float
280 Radius of the trust region at the last iteration.
281 constr_penalty : float
282 Penalty parameter at the last iteration, see `initial_constr_penalty`.
283 barrier_tolerance : float
284 Tolerance for the barrier subproblem at the last iteration.
285 Only for problems with inequality constraints.
286 barrier_parameter : float
287 Barrier parameter at the last iteration. Only for problems
288 with inequality constraints.
289 execution_time : float
290 Total execution time.
291 message : str
292 Termination message.
293 status : {0, 1, 2, 3}
294 Termination status:
296 * 0 : The maximum number of function evaluations is exceeded.
297 * 1 : `gtol` termination condition is satisfied.
298 * 2 : `xtol` termination condition is satisfied.
299 * 3 : `callback` function requested termination.
301 cg_stop_cond : int
302 Reason for CG subproblem termination at the last iteration:
304 * 0 : CG subproblem not evaluated.
305 * 1 : Iteration limit was reached.
306 * 2 : Reached the trust-region boundary.
307 * 3 : Negative curvature detected.
308 * 4 : Tolerance was satisfied.
310 References
311 ----------
312 .. [1] Conn, A. R., Gould, N. I., & Toint, P. L.
313 Trust region methods. 2000. Siam. pp. 19.
314 """
315 x0 = np.atleast_1d(x0).astype(float)
316 n_vars = np.size(x0)
317 if hess is None:
318 if callable(hessp):
319 hess = HessianLinearOperator(hessp, n_vars)
320 else:
321 hess = BFGS()
322 if disp and verbose == 0:
323 verbose = 1
325 if bounds is not None:
326 modified_lb = np.nextafter(bounds.lb, -np.inf, where=bounds.lb > -np.inf)
327 modified_ub = np.nextafter(bounds.ub, np.inf, where=bounds.ub < np.inf)
328 modified_lb = np.where(np.isfinite(bounds.lb), modified_lb, bounds.lb)
329 modified_ub = np.where(np.isfinite(bounds.ub), modified_ub, bounds.ub)
330 bounds = Bounds(modified_lb, modified_ub, keep_feasible=bounds.keep_feasible)
331 finite_diff_bounds = strict_bounds(bounds.lb, bounds.ub,
332 bounds.keep_feasible, n_vars)
333 else:
334 finite_diff_bounds = (-np.inf, np.inf)
336 # Define Objective Function
337 objective = ScalarFunction(fun, x0, args, grad, hess,
338 finite_diff_rel_step, finite_diff_bounds)
340 # Put constraints in list format when needed.
341 if isinstance(constraints, (NonlinearConstraint, LinearConstraint)):
342 constraints = [constraints]
344 # Prepare constraints.
345 prepared_constraints = [
346 PreparedConstraint(c, x0, sparse_jacobian, finite_diff_bounds)
347 for c in constraints]
349 # Check that all constraints are either sparse or dense.
350 n_sparse = sum(c.fun.sparse_jacobian for c in prepared_constraints)
351 if 0 < n_sparse < len(prepared_constraints):
352 raise ValueError("All constraints must have the same kind of the "
353 "Jacobian --- either all sparse or all dense. "
354 "You can set the sparsity globally by setting "
355 "`sparse_jacobian` to either True of False.")
356 if prepared_constraints:
357 sparse_jacobian = n_sparse > 0
359 if bounds is not None:
360 if sparse_jacobian is None:
361 sparse_jacobian = True
362 prepared_constraints.append(PreparedConstraint(bounds, x0,
363 sparse_jacobian))
365 # Concatenate initial constraints to the canonical form.
366 c_eq0, c_ineq0, J_eq0, J_ineq0 = initial_constraints_as_canonical(
367 n_vars, prepared_constraints, sparse_jacobian)
369 # Prepare all canonical constraints and concatenate it into one.
370 canonical_all = [CanonicalConstraint.from_PreparedConstraint(c)
371 for c in prepared_constraints]
373 if len(canonical_all) == 0:
374 canonical = CanonicalConstraint.empty(n_vars)
375 elif len(canonical_all) == 1:
376 canonical = canonical_all[0]
377 else:
378 canonical = CanonicalConstraint.concatenate(canonical_all,
379 sparse_jacobian)
381 # Generate the Hessian of the Lagrangian.
382 lagrangian_hess = LagrangianHessian(n_vars, objective.hess, canonical.hess)
384 # Choose appropriate method
385 if canonical.n_ineq == 0:
386 method = 'equality_constrained_sqp'
387 else:
388 method = 'tr_interior_point'
390 # Construct OptimizeResult
391 state = OptimizeResult(
392 nit=0, nfev=0, njev=0, nhev=0,
393 cg_niter=0, cg_stop_cond=0,
394 fun=objective.f, grad=objective.g,
395 lagrangian_grad=np.copy(objective.g),
396 constr=[c.fun.f for c in prepared_constraints],
397 jac=[c.fun.J for c in prepared_constraints],
398 constr_nfev=[0 for c in prepared_constraints],
399 constr_njev=[0 for c in prepared_constraints],
400 constr_nhev=[0 for c in prepared_constraints],
401 v=[c.fun.v for c in prepared_constraints],
402 method=method)
404 # Start counting
405 start_time = time.time()
407 # Define stop criteria
408 if method == 'equality_constrained_sqp':
409 def stop_criteria(state, x, last_iteration_failed,
410 optimality, constr_violation,
411 tr_radius, constr_penalty, cg_info):
412 state = update_state_sqp(state, x, last_iteration_failed,
413 objective, prepared_constraints,
414 start_time, tr_radius, constr_penalty,
415 cg_info)
416 if verbose == 2:
417 BasicReport.print_iteration(state.nit,
418 state.nfev,
419 state.cg_niter,
420 state.fun,
421 state.tr_radius,
422 state.optimality,
423 state.constr_violation)
424 elif verbose > 2:
425 SQPReport.print_iteration(state.nit,
426 state.nfev,
427 state.cg_niter,
428 state.fun,
429 state.tr_radius,
430 state.optimality,
431 state.constr_violation,
432 state.constr_penalty,
433 state.cg_stop_cond)
434 state.status = None
435 state.niter = state.nit # Alias for callback (backward-compatibility)
436 if callback is not None:
437 callback_stop = False
438 try:
439 callback_stop = callback(state)
440 except StopIteration:
441 callback_stop = True
442 if callback_stop:
443 state.status = 3
444 return True
445 if state.optimality < gtol and state.constr_violation < gtol:
446 state.status = 1
447 elif state.tr_radius < xtol:
448 state.status = 2
449 elif state.nit >= maxiter:
450 state.status = 0
451 return state.status in (0, 1, 2, 3)
452 elif method == 'tr_interior_point':
453 def stop_criteria(state, x, last_iteration_failed, tr_radius,
454 constr_penalty, cg_info, barrier_parameter,
455 barrier_tolerance):
456 state = update_state_ip(state, x, last_iteration_failed,
457 objective, prepared_constraints,
458 start_time, tr_radius, constr_penalty,
459 cg_info, barrier_parameter, barrier_tolerance)
460 if verbose == 2:
461 BasicReport.print_iteration(state.nit,
462 state.nfev,
463 state.cg_niter,
464 state.fun,
465 state.tr_radius,
466 state.optimality,
467 state.constr_violation)
468 elif verbose > 2:
469 IPReport.print_iteration(state.nit,
470 state.nfev,
471 state.cg_niter,
472 state.fun,
473 state.tr_radius,
474 state.optimality,
475 state.constr_violation,
476 state.constr_penalty,
477 state.barrier_parameter,
478 state.cg_stop_cond)
479 state.status = None
480 state.niter = state.nit # Alias for callback (backward compatibility)
481 if callback is not None:
482 callback_stop = False
483 try:
484 callback_stop = callback(state)
485 except StopIteration:
486 callback_stop = True
487 if callback_stop:
488 state.status = 3
489 return True
490 if state.optimality < gtol and state.constr_violation < gtol:
491 state.status = 1
492 elif (state.tr_radius < xtol
493 and state.barrier_parameter < barrier_tol):
494 state.status = 2
495 elif state.nit >= maxiter:
496 state.status = 0
497 return state.status in (0, 1, 2, 3)
499 if verbose == 2:
500 BasicReport.print_header()
501 elif verbose > 2:
502 if method == 'equality_constrained_sqp':
503 SQPReport.print_header()
504 elif method == 'tr_interior_point':
505 IPReport.print_header()
507 # Call inferior function to do the optimization
508 if method == 'equality_constrained_sqp':
509 def fun_and_constr(x):
510 f = objective.fun(x)
511 c_eq, _ = canonical.fun(x)
512 return f, c_eq
514 def grad_and_jac(x):
515 g = objective.grad(x)
516 J_eq, _ = canonical.jac(x)
517 return g, J_eq
519 _, result = equality_constrained_sqp(
520 fun_and_constr, grad_and_jac, lagrangian_hess,
521 x0, objective.f, objective.g,
522 c_eq0, J_eq0,
523 stop_criteria, state,
524 initial_constr_penalty, initial_tr_radius,
525 factorization_method)
527 elif method == 'tr_interior_point':
528 _, result = tr_interior_point(
529 objective.fun, objective.grad, lagrangian_hess,
530 n_vars, canonical.n_ineq, canonical.n_eq,
531 canonical.fun, canonical.jac,
532 x0, objective.f, objective.g,
533 c_ineq0, J_ineq0, c_eq0, J_eq0,
534 stop_criteria,
535 canonical.keep_feasible,
536 xtol, state, initial_barrier_parameter,
537 initial_barrier_tolerance,
538 initial_constr_penalty, initial_tr_radius,
539 factorization_method)
541 # Status 3 occurs when the callback function requests termination,
542 # this is assumed to not be a success.
543 result.success = True if result.status in (1, 2) else False
544 result.message = TERMINATION_MESSAGES[result.status]
546 # Alias (for backward compatibility with 1.1.0)
547 result.niter = result.nit
549 if verbose == 2:
550 BasicReport.print_footer()
551 elif verbose > 2:
552 if method == 'equality_constrained_sqp':
553 SQPReport.print_footer()
554 elif method == 'tr_interior_point':
555 IPReport.print_footer()
556 if verbose >= 1:
557 print(result.message)
558 print("Number of iterations: {}, function evaluations: {}, "
559 "CG iterations: {}, optimality: {:.2e}, "
560 "constraint violation: {:.2e}, execution time: {:4.2} s."
561 .format(result.nit, result.nfev, result.cg_niter,
562 result.optimality, result.constr_violation,
563 result.execution_time))
564 return result