Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/optimize/_trustregion_constr/tr_interior_point.py: 16%
147 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
1"""Trust-region interior point method.
3References
4----------
5.. [1] Byrd, Richard H., Mary E. Hribar, and Jorge Nocedal.
6 "An interior point algorithm for large-scale nonlinear
7 programming." SIAM Journal on Optimization 9.4 (1999): 877-900.
8.. [2] Byrd, Richard H., Guanghui Liu, and Jorge Nocedal.
9 "On the local behavior of an interior point method for
10 nonlinear programming." Numerical analysis 1997 (1997): 37-56.
11.. [3] Nocedal, Jorge, and Stephen J. Wright. "Numerical optimization"
12 Second Edition (2006).
13"""
15import scipy.sparse as sps
16import numpy as np
17from .equality_constrained_sqp import equality_constrained_sqp
18from scipy.sparse.linalg import LinearOperator
20__all__ = ['tr_interior_point']
23class BarrierSubproblem:
24 """
25 Barrier optimization problem:
26 minimize fun(x) - barrier_parameter*sum(log(s))
27 subject to: constr_eq(x) = 0
28 constr_ineq(x) + s = 0
29 """
31 def __init__(self, x0, s0, fun, grad, lagr_hess, n_vars, n_ineq, n_eq,
32 constr, jac, barrier_parameter, tolerance,
33 enforce_feasibility, global_stop_criteria,
34 xtol, fun0, grad0, constr_ineq0, jac_ineq0, constr_eq0,
35 jac_eq0):
36 # Store parameters
37 self.n_vars = n_vars
38 self.x0 = x0
39 self.s0 = s0
40 self.fun = fun
41 self.grad = grad
42 self.lagr_hess = lagr_hess
43 self.constr = constr
44 self.jac = jac
45 self.barrier_parameter = barrier_parameter
46 self.tolerance = tolerance
47 self.n_eq = n_eq
48 self.n_ineq = n_ineq
49 self.enforce_feasibility = enforce_feasibility
50 self.global_stop_criteria = global_stop_criteria
51 self.xtol = xtol
52 self.fun0 = self._compute_function(fun0, constr_ineq0, s0)
53 self.grad0 = self._compute_gradient(grad0)
54 self.constr0 = self._compute_constr(constr_ineq0, constr_eq0, s0)
55 self.jac0 = self._compute_jacobian(jac_eq0, jac_ineq0, s0)
56 self.terminate = False
58 def update(self, barrier_parameter, tolerance):
59 self.barrier_parameter = barrier_parameter
60 self.tolerance = tolerance
62 def get_slack(self, z):
63 return z[self.n_vars:self.n_vars+self.n_ineq]
65 def get_variables(self, z):
66 return z[:self.n_vars]
68 def function_and_constraints(self, z):
69 """Returns barrier function and constraints at given point.
71 For z = [x, s], returns barrier function:
72 function(z) = fun(x) - barrier_parameter*sum(log(s))
73 and barrier constraints:
74 constraints(z) = [ constr_eq(x) ]
75 [ constr_ineq(x) + s ]
77 """
78 # Get variables and slack variables
79 x = self.get_variables(z)
80 s = self.get_slack(z)
81 # Compute function and constraints
82 f = self.fun(x)
83 c_eq, c_ineq = self.constr(x)
84 # Return objective function and constraints
85 return (self._compute_function(f, c_ineq, s),
86 self._compute_constr(c_ineq, c_eq, s))
88 def _compute_function(self, f, c_ineq, s):
89 # Use technique from Nocedal and Wright book, ref [3]_, p.576,
90 # to guarantee constraints from `enforce_feasibility`
91 # stay feasible along iterations.
92 s[self.enforce_feasibility] = -c_ineq[self.enforce_feasibility]
93 log_s = [np.log(s_i) if s_i > 0 else -np.inf for s_i in s]
94 # Compute barrier objective function
95 return f - self.barrier_parameter*np.sum(log_s)
97 def _compute_constr(self, c_ineq, c_eq, s):
98 # Compute barrier constraint
99 return np.hstack((c_eq,
100 c_ineq + s))
102 def scaling(self, z):
103 """Returns scaling vector.
104 Given by:
105 scaling = [ones(n_vars), s]
106 """
107 s = self.get_slack(z)
108 diag_elements = np.hstack((np.ones(self.n_vars), s))
110 # Diagonal matrix
111 def matvec(vec):
112 return diag_elements*vec
113 return LinearOperator((self.n_vars+self.n_ineq,
114 self.n_vars+self.n_ineq),
115 matvec)
117 def gradient_and_jacobian(self, z):
118 """Returns scaled gradient.
120 Return scaled gradient:
121 gradient = [ grad(x) ]
122 [ -barrier_parameter*ones(n_ineq) ]
123 and scaled Jacobian matrix:
124 jacobian = [ jac_eq(x) 0 ]
125 [ jac_ineq(x) S ]
126 Both of them scaled by the previously defined scaling factor.
127 """
128 # Get variables and slack variables
129 x = self.get_variables(z)
130 s = self.get_slack(z)
131 # Compute first derivatives
132 g = self.grad(x)
133 J_eq, J_ineq = self.jac(x)
134 # Return gradient and Jacobian
135 return (self._compute_gradient(g),
136 self._compute_jacobian(J_eq, J_ineq, s))
138 def _compute_gradient(self, g):
139 return np.hstack((g, -self.barrier_parameter*np.ones(self.n_ineq)))
141 def _compute_jacobian(self, J_eq, J_ineq, s):
142 if self.n_ineq == 0:
143 return J_eq
144 else:
145 if sps.issparse(J_eq) or sps.issparse(J_ineq):
146 # It is expected that J_eq and J_ineq
147 # are already `csr_matrix` because of
148 # the way ``BoxConstraint``, ``NonlinearConstraint``
149 # and ``LinearConstraint`` are defined.
150 J_eq = sps.csr_matrix(J_eq)
151 J_ineq = sps.csr_matrix(J_ineq)
152 return self._assemble_sparse_jacobian(J_eq, J_ineq, s)
153 else:
154 S = np.diag(s)
155 zeros = np.zeros((self.n_eq, self.n_ineq))
156 # Convert to matrix
157 if sps.issparse(J_ineq):
158 J_ineq = J_ineq.toarray()
159 if sps.issparse(J_eq):
160 J_eq = J_eq.toarray()
161 # Concatenate matrices
162 return np.block([[J_eq, zeros],
163 [J_ineq, S]])
165 def _assemble_sparse_jacobian(self, J_eq, J_ineq, s):
166 """Assemble sparse Jacobian given its components.
168 Given ``J_eq``, ``J_ineq`` and ``s`` returns:
169 jacobian = [ J_eq, 0 ]
170 [ J_ineq, diag(s) ]
172 It is equivalent to:
173 sps.bmat([[ J_eq, None ],
174 [ J_ineq, diag(s) ]], "csr")
175 but significantly more efficient for this
176 given structure.
177 """
178 n_vars, n_ineq, n_eq = self.n_vars, self.n_ineq, self.n_eq
179 J_aux = sps.vstack([J_eq, J_ineq], "csr")
180 indptr, indices, data = J_aux.indptr, J_aux.indices, J_aux.data
181 new_indptr = indptr + np.hstack((np.zeros(n_eq, dtype=int),
182 np.arange(n_ineq+1, dtype=int)))
183 size = indices.size+n_ineq
184 new_indices = np.empty(size)
185 new_data = np.empty(size)
186 mask = np.full(size, False, bool)
187 mask[new_indptr[-n_ineq:]-1] = True
188 new_indices[mask] = n_vars+np.arange(n_ineq)
189 new_indices[~mask] = indices
190 new_data[mask] = s
191 new_data[~mask] = data
192 J = sps.csr_matrix((new_data, new_indices, new_indptr),
193 (n_eq + n_ineq, n_vars + n_ineq))
194 return J
196 def lagrangian_hessian_x(self, z, v):
197 """Returns Lagrangian Hessian (in relation to `x`) -> Hx"""
198 x = self.get_variables(z)
199 # Get lagrange multipliers relatated to nonlinear equality constraints
200 v_eq = v[:self.n_eq]
201 # Get lagrange multipliers relatated to nonlinear ineq. constraints
202 v_ineq = v[self.n_eq:self.n_eq+self.n_ineq]
203 lagr_hess = self.lagr_hess
204 return lagr_hess(x, v_eq, v_ineq)
206 def lagrangian_hessian_s(self, z, v):
207 """Returns scaled Lagrangian Hessian (in relation to`s`) -> S Hs S"""
208 s = self.get_slack(z)
209 # Using the primal formulation:
210 # S Hs S = diag(s)*diag(barrier_parameter/s**2)*diag(s).
211 # Reference [1]_ p. 882, formula (3.1)
212 primal = self.barrier_parameter
213 # Using the primal-dual formulation
214 # S Hs S = diag(s)*diag(v/s)*diag(s)
215 # Reference [1]_ p. 883, formula (3.11)
216 primal_dual = v[-self.n_ineq:]*s
217 # Uses the primal-dual formulation for
218 # positives values of v_ineq, and primal
219 # formulation for the remaining ones.
220 return np.where(v[-self.n_ineq:] > 0, primal_dual, primal)
222 def lagrangian_hessian(self, z, v):
223 """Returns scaled Lagrangian Hessian"""
224 # Compute Hessian in relation to x and s
225 Hx = self.lagrangian_hessian_x(z, v)
226 if self.n_ineq > 0:
227 S_Hs_S = self.lagrangian_hessian_s(z, v)
229 # The scaled Lagragian Hessian is:
230 # [ Hx 0 ]
231 # [ 0 S Hs S ]
232 def matvec(vec):
233 vec_x = self.get_variables(vec)
234 vec_s = self.get_slack(vec)
235 if self.n_ineq > 0:
236 return np.hstack((Hx.dot(vec_x), S_Hs_S*vec_s))
237 else:
238 return Hx.dot(vec_x)
239 return LinearOperator((self.n_vars+self.n_ineq,
240 self.n_vars+self.n_ineq),
241 matvec)
243 def stop_criteria(self, state, z, last_iteration_failed,
244 optimality, constr_violation,
245 trust_radius, penalty, cg_info):
246 """Stop criteria to the barrier problem.
247 The criteria here proposed is similar to formula (2.3)
248 from [1]_, p.879.
249 """
250 x = self.get_variables(z)
251 if self.global_stop_criteria(state, x,
252 last_iteration_failed,
253 trust_radius, penalty,
254 cg_info,
255 self.barrier_parameter,
256 self.tolerance):
257 self.terminate = True
258 return True
259 else:
260 g_cond = (optimality < self.tolerance and
261 constr_violation < self.tolerance)
262 x_cond = trust_radius < self.xtol
263 return g_cond or x_cond
266def tr_interior_point(fun, grad, lagr_hess, n_vars, n_ineq, n_eq,
267 constr, jac, x0, fun0, grad0,
268 constr_ineq0, jac_ineq0, constr_eq0,
269 jac_eq0, stop_criteria,
270 enforce_feasibility, xtol, state,
271 initial_barrier_parameter,
272 initial_tolerance,
273 initial_penalty,
274 initial_trust_radius,
275 factorization_method):
276 """Trust-region interior points method.
278 Solve problem:
279 minimize fun(x)
280 subject to: constr_ineq(x) <= 0
281 constr_eq(x) = 0
282 using trust-region interior point method described in [1]_.
283 """
284 # BOUNDARY_PARAMETER controls the decrease on the slack
285 # variables. Represents ``tau`` from [1]_ p.885, formula (3.18).
286 BOUNDARY_PARAMETER = 0.995
287 # BARRIER_DECAY_RATIO controls the decay of the barrier parameter
288 # and of the subproblem toloerance. Represents ``theta`` from [1]_ p.879.
289 BARRIER_DECAY_RATIO = 0.2
290 # TRUST_ENLARGEMENT controls the enlargement on trust radius
291 # after each iteration
292 TRUST_ENLARGEMENT = 5
294 # Default enforce_feasibility
295 if enforce_feasibility is None:
296 enforce_feasibility = np.zeros(n_ineq, bool)
297 # Initial Values
298 barrier_parameter = initial_barrier_parameter
299 tolerance = initial_tolerance
300 trust_radius = initial_trust_radius
301 # Define initial value for the slack variables
302 s0 = np.maximum(-1.5*constr_ineq0, np.ones(n_ineq))
303 # Define barrier subproblem
304 subprob = BarrierSubproblem(
305 x0, s0, fun, grad, lagr_hess, n_vars, n_ineq, n_eq, constr, jac,
306 barrier_parameter, tolerance, enforce_feasibility,
307 stop_criteria, xtol, fun0, grad0, constr_ineq0, jac_ineq0,
308 constr_eq0, jac_eq0)
309 # Define initial parameter for the first iteration.
310 z = np.hstack((x0, s0))
311 fun0_subprob, constr0_subprob = subprob.fun0, subprob.constr0
312 grad0_subprob, jac0_subprob = subprob.grad0, subprob.jac0
313 # Define trust region bounds
314 trust_lb = np.hstack((np.full(subprob.n_vars, -np.inf),
315 np.full(subprob.n_ineq, -BOUNDARY_PARAMETER)))
316 trust_ub = np.full(subprob.n_vars+subprob.n_ineq, np.inf)
318 # Solves a sequence of barrier problems
319 while True:
320 # Solve SQP subproblem
321 z, state = equality_constrained_sqp(
322 subprob.function_and_constraints,
323 subprob.gradient_and_jacobian,
324 subprob.lagrangian_hessian,
325 z, fun0_subprob, grad0_subprob,
326 constr0_subprob, jac0_subprob, subprob.stop_criteria,
327 state, initial_penalty, trust_radius,
328 factorization_method, trust_lb, trust_ub, subprob.scaling)
329 if subprob.terminate:
330 break
331 # Update parameters
332 trust_radius = max(initial_trust_radius,
333 TRUST_ENLARGEMENT*state.tr_radius)
334 # TODO: Use more advanced strategies from [2]_
335 # to update this parameters.
336 barrier_parameter *= BARRIER_DECAY_RATIO
337 tolerance *= BARRIER_DECAY_RATIO
338 # Update Barrier Problem
339 subprob.update(barrier_parameter, tolerance)
340 # Compute initial values for next iteration
341 fun0_subprob, constr0_subprob = subprob.function_and_constraints(z)
342 grad0_subprob, jac0_subprob = subprob.gradient_and_jacobian(z)
344 # Get x and s
345 x = subprob.get_variables(z)
346 return x, state