Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/optimize/_trustregion_ncg.py: 16%
50 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
1"""Newton-CG trust-region optimization."""
2import math
4import numpy as np
5import scipy.linalg
6from ._trustregion import (_minimize_trust_region, BaseQuadraticSubproblem)
8__all__ = []
11def _minimize_trust_ncg(fun, x0, args=(), jac=None, hess=None, hessp=None,
12 **trust_region_options):
13 """
14 Minimization of scalar function of one or more variables using
15 the Newton conjugate gradient trust-region algorithm.
17 Options
18 -------
19 initial_trust_radius : float
20 Initial trust-region radius.
21 max_trust_radius : float
22 Maximum value of the trust-region radius. No steps that are longer
23 than this value will be proposed.
24 eta : float
25 Trust region related acceptance stringency for proposed steps.
26 gtol : float
27 Gradient norm must be less than `gtol` before successful
28 termination.
30 """
31 if jac is None:
32 raise ValueError('Jacobian is required for Newton-CG trust-region '
33 'minimization')
34 if hess is None and hessp is None:
35 raise ValueError('Either the Hessian or the Hessian-vector product '
36 'is required for Newton-CG trust-region minimization')
37 return _minimize_trust_region(fun, x0, args=args, jac=jac, hess=hess,
38 hessp=hessp, subproblem=CGSteihaugSubproblem,
39 **trust_region_options)
42class CGSteihaugSubproblem(BaseQuadraticSubproblem):
43 """Quadratic subproblem solved by a conjugate gradient method"""
44 def solve(self, trust_radius):
45 """
46 Solve the subproblem using a conjugate gradient method.
48 Parameters
49 ----------
50 trust_radius : float
51 We are allowed to wander only this far away from the origin.
53 Returns
54 -------
55 p : ndarray
56 The proposed step.
57 hits_boundary : bool
58 True if the proposed step is on the boundary of the trust region.
60 Notes
61 -----
62 This is algorithm (7.2) of Nocedal and Wright 2nd edition.
63 Only the function that computes the Hessian-vector product is required.
64 The Hessian itself is not required, and the Hessian does
65 not need to be positive semidefinite.
66 """
68 # get the norm of jacobian and define the origin
69 p_origin = np.zeros_like(self.jac)
71 # define a default tolerance
72 tolerance = min(0.5, math.sqrt(self.jac_mag)) * self.jac_mag
74 # Stop the method if the search direction
75 # is a direction of nonpositive curvature.
76 if self.jac_mag < tolerance:
77 hits_boundary = False
78 return p_origin, hits_boundary
80 # init the state for the first iteration
81 z = p_origin
82 r = self.jac
83 d = -r
85 # Search for the min of the approximation of the objective function.
86 while True:
88 # do an iteration
89 Bd = self.hessp(d)
90 dBd = np.dot(d, Bd)
91 if dBd <= 0:
92 # Look at the two boundary points.
93 # Find both values of t to get the boundary points such that
94 # ||z + t d|| == trust_radius
95 # and then choose the one with the predicted min value.
96 ta, tb = self.get_boundaries_intersections(z, d, trust_radius)
97 pa = z + ta * d
98 pb = z + tb * d
99 if self(pa) < self(pb):
100 p_boundary = pa
101 else:
102 p_boundary = pb
103 hits_boundary = True
104 return p_boundary, hits_boundary
105 r_squared = np.dot(r, r)
106 alpha = r_squared / dBd
107 z_next = z + alpha * d
108 if scipy.linalg.norm(z_next) >= trust_radius:
109 # Find t >= 0 to get the boundary point such that
110 # ||z + t d|| == trust_radius
111 ta, tb = self.get_boundaries_intersections(z, d, trust_radius)
112 p_boundary = z + tb * d
113 hits_boundary = True
114 return p_boundary, hits_boundary
115 r_next = r + alpha * Bd
116 r_next_squared = np.dot(r_next, r_next)
117 if math.sqrt(r_next_squared) < tolerance:
118 hits_boundary = False
119 return z_next, hits_boundary
120 beta_next = r_next_squared / r_squared
121 d_next = -r_next + beta_next * d
123 # update the state for the next iteration
124 z = z_next
125 r = r_next
126 d = d_next