Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/optimize/_trustregion_constr/qp_subproblem.py: 6%
214 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
1"""Equality-constrained quadratic programming solvers."""
3from scipy.sparse import (linalg, bmat, csc_matrix)
4from math import copysign
5import numpy as np
6from numpy.linalg import norm
8__all__ = [
9 'eqp_kktfact',
10 'sphere_intersections',
11 'box_intersections',
12 'box_sphere_intersections',
13 'inside_box_boundaries',
14 'modified_dogleg',
15 'projected_cg'
16]
19# For comparison with the projected CG
20def eqp_kktfact(H, c, A, b):
21 """Solve equality-constrained quadratic programming (EQP) problem.
23 Solve ``min 1/2 x.T H x + x.t c`` subject to ``A x + b = 0``
24 using direct factorization of the KKT system.
26 Parameters
27 ----------
28 H : sparse matrix, shape (n, n)
29 Hessian matrix of the EQP problem.
30 c : array_like, shape (n,)
31 Gradient of the quadratic objective function.
32 A : sparse matrix
33 Jacobian matrix of the EQP problem.
34 b : array_like, shape (m,)
35 Right-hand side of the constraint equation.
37 Returns
38 -------
39 x : array_like, shape (n,)
40 Solution of the KKT problem.
41 lagrange_multipliers : ndarray, shape (m,)
42 Lagrange multipliers of the KKT problem.
43 """
44 n, = np.shape(c) # Number of parameters
45 m, = np.shape(b) # Number of constraints
47 # Karush-Kuhn-Tucker matrix of coefficients.
48 # Defined as in Nocedal/Wright "Numerical
49 # Optimization" p.452 in Eq. (16.4).
50 kkt_matrix = csc_matrix(bmat([[H, A.T], [A, None]]))
51 # Vector of coefficients.
52 kkt_vec = np.hstack([-c, -b])
54 # TODO: Use a symmetric indefinite factorization
55 # to solve the system twice as fast (because
56 # of the symmetry).
57 lu = linalg.splu(kkt_matrix)
58 kkt_sol = lu.solve(kkt_vec)
59 x = kkt_sol[:n]
60 lagrange_multipliers = -kkt_sol[n:n+m]
62 return x, lagrange_multipliers
65def sphere_intersections(z, d, trust_radius,
66 entire_line=False):
67 """Find the intersection between segment (or line) and spherical constraints.
69 Find the intersection between the segment (or line) defined by the
70 parametric equation ``x(t) = z + t*d`` and the ball
71 ``||x|| <= trust_radius``.
73 Parameters
74 ----------
75 z : array_like, shape (n,)
76 Initial point.
77 d : array_like, shape (n,)
78 Direction.
79 trust_radius : float
80 Ball radius.
81 entire_line : bool, optional
82 When ``True``, the function returns the intersection between the line
83 ``x(t) = z + t*d`` (``t`` can assume any value) and the ball
84 ``||x|| <= trust_radius``. When ``False``, the function returns the intersection
85 between the segment ``x(t) = z + t*d``, ``0 <= t <= 1``, and the ball.
87 Returns
88 -------
89 ta, tb : float
90 The line/segment ``x(t) = z + t*d`` is inside the ball for
91 for ``ta <= t <= tb``.
92 intersect : bool
93 When ``True``, there is a intersection between the line/segment
94 and the sphere. On the other hand, when ``False``, there is no
95 intersection.
96 """
97 # Special case when d=0
98 if norm(d) == 0:
99 return 0, 0, False
100 # Check for inf trust_radius
101 if np.isinf(trust_radius):
102 if entire_line:
103 ta = -np.inf
104 tb = np.inf
105 else:
106 ta = 0
107 tb = 1
108 intersect = True
109 return ta, tb, intersect
111 a = np.dot(d, d)
112 b = 2 * np.dot(z, d)
113 c = np.dot(z, z) - trust_radius**2
114 discriminant = b*b - 4*a*c
115 if discriminant < 0:
116 intersect = False
117 return 0, 0, intersect
118 sqrt_discriminant = np.sqrt(discriminant)
120 # The following calculation is mathematically
121 # equivalent to:
122 # ta = (-b - sqrt_discriminant) / (2*a)
123 # tb = (-b + sqrt_discriminant) / (2*a)
124 # but produce smaller round off errors.
125 # Look at Matrix Computation p.97
126 # for a better justification.
127 aux = b + copysign(sqrt_discriminant, b)
128 ta = -aux / (2*a)
129 tb = -2*c / aux
130 ta, tb = sorted([ta, tb])
132 if entire_line:
133 intersect = True
134 else:
135 # Checks to see if intersection happens
136 # within vectors length.
137 if tb < 0 or ta > 1:
138 intersect = False
139 ta = 0
140 tb = 0
141 else:
142 intersect = True
143 # Restrict intersection interval
144 # between 0 and 1.
145 ta = max(0, ta)
146 tb = min(1, tb)
148 return ta, tb, intersect
151def box_intersections(z, d, lb, ub,
152 entire_line=False):
153 """Find the intersection between segment (or line) and box constraints.
155 Find the intersection between the segment (or line) defined by the
156 parametric equation ``x(t) = z + t*d`` and the rectangular box
157 ``lb <= x <= ub``.
159 Parameters
160 ----------
161 z : array_like, shape (n,)
162 Initial point.
163 d : array_like, shape (n,)
164 Direction.
165 lb : array_like, shape (n,)
166 Lower bounds to each one of the components of ``x``. Used
167 to delimit the rectangular box.
168 ub : array_like, shape (n, )
169 Upper bounds to each one of the components of ``x``. Used
170 to delimit the rectangular box.
171 entire_line : bool, optional
172 When ``True``, the function returns the intersection between the line
173 ``x(t) = z + t*d`` (``t`` can assume any value) and the rectangular
174 box. When ``False``, the function returns the intersection between the segment
175 ``x(t) = z + t*d``, ``0 <= t <= 1``, and the rectangular box.
177 Returns
178 -------
179 ta, tb : float
180 The line/segment ``x(t) = z + t*d`` is inside the box for
181 for ``ta <= t <= tb``.
182 intersect : bool
183 When ``True``, there is a intersection between the line (or segment)
184 and the rectangular box. On the other hand, when ``False``, there is no
185 intersection.
186 """
187 # Make sure it is a numpy array
188 z = np.asarray(z)
189 d = np.asarray(d)
190 lb = np.asarray(lb)
191 ub = np.asarray(ub)
192 # Special case when d=0
193 if norm(d) == 0:
194 return 0, 0, False
196 # Get values for which d==0
197 zero_d = (d == 0)
198 # If the boundaries are not satisfied for some coordinate
199 # for which "d" is zero, there is no box-line intersection.
200 if (z[zero_d] < lb[zero_d]).any() or (z[zero_d] > ub[zero_d]).any():
201 intersect = False
202 return 0, 0, intersect
203 # Remove values for which d is zero
204 not_zero_d = np.logical_not(zero_d)
205 z = z[not_zero_d]
206 d = d[not_zero_d]
207 lb = lb[not_zero_d]
208 ub = ub[not_zero_d]
210 # Find a series of intervals (t_lb[i], t_ub[i]).
211 t_lb = (lb-z) / d
212 t_ub = (ub-z) / d
213 # Get the intersection of all those intervals.
214 ta = max(np.minimum(t_lb, t_ub))
215 tb = min(np.maximum(t_lb, t_ub))
217 # Check if intersection is feasible
218 if ta <= tb:
219 intersect = True
220 else:
221 intersect = False
222 # Checks to see if intersection happens within vectors length.
223 if not entire_line:
224 if tb < 0 or ta > 1:
225 intersect = False
226 ta = 0
227 tb = 0
228 else:
229 # Restrict intersection interval between 0 and 1.
230 ta = max(0, ta)
231 tb = min(1, tb)
233 return ta, tb, intersect
236def box_sphere_intersections(z, d, lb, ub, trust_radius,
237 entire_line=False,
238 extra_info=False):
239 """Find the intersection between segment (or line) and box/sphere constraints.
241 Find the intersection between the segment (or line) defined by the
242 parametric equation ``x(t) = z + t*d``, the rectangular box
243 ``lb <= x <= ub`` and the ball ``||x|| <= trust_radius``.
245 Parameters
246 ----------
247 z : array_like, shape (n,)
248 Initial point.
249 d : array_like, shape (n,)
250 Direction.
251 lb : array_like, shape (n,)
252 Lower bounds to each one of the components of ``x``. Used
253 to delimit the rectangular box.
254 ub : array_like, shape (n, )
255 Upper bounds to each one of the components of ``x``. Used
256 to delimit the rectangular box.
257 trust_radius : float
258 Ball radius.
259 entire_line : bool, optional
260 When ``True``, the function returns the intersection between the line
261 ``x(t) = z + t*d`` (``t`` can assume any value) and the constraints.
262 When ``False``, the function returns the intersection between the segment
263 ``x(t) = z + t*d``, ``0 <= t <= 1`` and the constraints.
264 extra_info : bool, optional
265 When ``True``, the function returns ``intersect_sphere`` and ``intersect_box``.
267 Returns
268 -------
269 ta, tb : float
270 The line/segment ``x(t) = z + t*d`` is inside the rectangular box and
271 inside the ball for ``ta <= t <= tb``.
272 intersect : bool
273 When ``True``, there is a intersection between the line (or segment)
274 and both constraints. On the other hand, when ``False``, there is no
275 intersection.
276 sphere_info : dict, optional
277 Dictionary ``{ta, tb, intersect}`` containing the interval ``[ta, tb]``
278 for which the line intercepts the ball. And a boolean value indicating
279 whether the sphere is intersected by the line.
280 box_info : dict, optional
281 Dictionary ``{ta, tb, intersect}`` containing the interval ``[ta, tb]``
282 for which the line intercepts the box. And a boolean value indicating
283 whether the box is intersected by the line.
284 """
285 ta_b, tb_b, intersect_b = box_intersections(z, d, lb, ub,
286 entire_line)
287 ta_s, tb_s, intersect_s = sphere_intersections(z, d,
288 trust_radius,
289 entire_line)
290 ta = np.maximum(ta_b, ta_s)
291 tb = np.minimum(tb_b, tb_s)
292 if intersect_b and intersect_s and ta <= tb:
293 intersect = True
294 else:
295 intersect = False
297 if extra_info:
298 sphere_info = {'ta': ta_s, 'tb': tb_s, 'intersect': intersect_s}
299 box_info = {'ta': ta_b, 'tb': tb_b, 'intersect': intersect_b}
300 return ta, tb, intersect, sphere_info, box_info
301 else:
302 return ta, tb, intersect
305def inside_box_boundaries(x, lb, ub):
306 """Check if lb <= x <= ub."""
307 return (lb <= x).all() and (x <= ub).all()
310def reinforce_box_boundaries(x, lb, ub):
311 """Return clipped value of x"""
312 return np.minimum(np.maximum(x, lb), ub)
315def modified_dogleg(A, Y, b, trust_radius, lb, ub):
316 """Approximately minimize ``1/2*|| A x + b ||^2`` inside trust-region.
318 Approximately solve the problem of minimizing ``1/2*|| A x + b ||^2``
319 subject to ``||x|| < Delta`` and ``lb <= x <= ub`` using a modification
320 of the classical dogleg approach.
322 Parameters
323 ----------
324 A : LinearOperator (or sparse matrix or ndarray), shape (m, n)
325 Matrix ``A`` in the minimization problem. It should have
326 dimension ``(m, n)`` such that ``m < n``.
327 Y : LinearOperator (or sparse matrix or ndarray), shape (n, m)
328 LinearOperator that apply the projection matrix
329 ``Q = A.T inv(A A.T)`` to the vector. The obtained vector
330 ``y = Q x`` being the minimum norm solution of ``A y = x``.
331 b : array_like, shape (m,)
332 Vector ``b``in the minimization problem.
333 trust_radius: float
334 Trust radius to be considered. Delimits a sphere boundary
335 to the problem.
336 lb : array_like, shape (n,)
337 Lower bounds to each one of the components of ``x``.
338 It is expected that ``lb <= 0``, otherwise the algorithm
339 may fail. If ``lb[i] = -Inf``, the lower
340 bound for the ith component is just ignored.
341 ub : array_like, shape (n, )
342 Upper bounds to each one of the components of ``x``.
343 It is expected that ``ub >= 0``, otherwise the algorithm
344 may fail. If ``ub[i] = Inf``, the upper bound for the ith
345 component is just ignored.
347 Returns
348 -------
349 x : array_like, shape (n,)
350 Solution to the problem.
352 Notes
353 -----
354 Based on implementations described in pp. 885-886 from [1]_.
356 References
357 ----------
358 .. [1] Byrd, Richard H., Mary E. Hribar, and Jorge Nocedal.
359 "An interior point algorithm for large-scale nonlinear
360 programming." SIAM Journal on Optimization 9.4 (1999): 877-900.
361 """
362 # Compute minimum norm minimizer of 1/2*|| A x + b ||^2.
363 newton_point = -Y.dot(b)
364 # Check for interior point
365 if inside_box_boundaries(newton_point, lb, ub) \
366 and norm(newton_point) <= trust_radius:
367 x = newton_point
368 return x
370 # Compute gradient vector ``g = A.T b``
371 g = A.T.dot(b)
372 # Compute Cauchy point
373 # `cauchy_point = g.T g / (g.T A.T A g)``.
374 A_g = A.dot(g)
375 cauchy_point = -np.dot(g, g) / np.dot(A_g, A_g) * g
376 # Origin
377 origin_point = np.zeros_like(cauchy_point)
379 # Check the segment between cauchy_point and newton_point
380 # for a possible solution.
381 z = cauchy_point
382 p = newton_point - cauchy_point
383 _, alpha, intersect = box_sphere_intersections(z, p, lb, ub,
384 trust_radius)
385 if intersect:
386 x1 = z + alpha*p
387 else:
388 # Check the segment between the origin and cauchy_point
389 # for a possible solution.
390 z = origin_point
391 p = cauchy_point
392 _, alpha, _ = box_sphere_intersections(z, p, lb, ub,
393 trust_radius)
394 x1 = z + alpha*p
396 # Check the segment between origin and newton_point
397 # for a possible solution.
398 z = origin_point
399 p = newton_point
400 _, alpha, _ = box_sphere_intersections(z, p, lb, ub,
401 trust_radius)
402 x2 = z + alpha*p
404 # Return the best solution among x1 and x2.
405 if norm(A.dot(x1) + b) < norm(A.dot(x2) + b):
406 return x1
407 else:
408 return x2
411def projected_cg(H, c, Z, Y, b, trust_radius=np.inf,
412 lb=None, ub=None, tol=None,
413 max_iter=None, max_infeasible_iter=None,
414 return_all=False):
415 """Solve EQP problem with projected CG method.
417 Solve equality-constrained quadratic programming problem
418 ``min 1/2 x.T H x + x.t c`` subject to ``A x + b = 0`` and,
419 possibly, to trust region constraints ``||x|| < trust_radius``
420 and box constraints ``lb <= x <= ub``.
422 Parameters
423 ----------
424 H : LinearOperator (or sparse matrix or ndarray), shape (n, n)
425 Operator for computing ``H v``.
426 c : array_like, shape (n,)
427 Gradient of the quadratic objective function.
428 Z : LinearOperator (or sparse matrix or ndarray), shape (n, n)
429 Operator for projecting ``x`` into the null space of A.
430 Y : LinearOperator, sparse matrix, ndarray, shape (n, m)
431 Operator that, for a given a vector ``b``, compute smallest
432 norm solution of ``A x + b = 0``.
433 b : array_like, shape (m,)
434 Right-hand side of the constraint equation.
435 trust_radius : float, optional
436 Trust radius to be considered. By default, uses ``trust_radius=inf``,
437 which means no trust radius at all.
438 lb : array_like, shape (n,), optional
439 Lower bounds to each one of the components of ``x``.
440 If ``lb[i] = -Inf`` the lower bound for the i-th
441 component is just ignored (default).
442 ub : array_like, shape (n, ), optional
443 Upper bounds to each one of the components of ``x``.
444 If ``ub[i] = Inf`` the upper bound for the i-th
445 component is just ignored (default).
446 tol : float, optional
447 Tolerance used to interrupt the algorithm.
448 max_iter : int, optional
449 Maximum algorithm iterations. Where ``max_inter <= n-m``.
450 By default, uses ``max_iter = n-m``.
451 max_infeasible_iter : int, optional
452 Maximum infeasible (regarding box constraints) iterations the
453 algorithm is allowed to take.
454 By default, uses ``max_infeasible_iter = n-m``.
455 return_all : bool, optional
456 When ``true``, return the list of all vectors through the iterations.
458 Returns
459 -------
460 x : array_like, shape (n,)
461 Solution of the EQP problem.
462 info : Dict
463 Dictionary containing the following:
465 - niter : Number of iterations.
466 - stop_cond : Reason for algorithm termination:
467 1. Iteration limit was reached;
468 2. Reached the trust-region boundary;
469 3. Negative curvature detected;
470 4. Tolerance was satisfied.
471 - allvecs : List containing all intermediary vectors (optional).
472 - hits_boundary : True if the proposed step is on the boundary
473 of the trust region.
475 Notes
476 -----
477 Implementation of Algorithm 6.2 on [1]_.
479 In the absence of spherical and box constraints, for sufficient
480 iterations, the method returns a truly optimal result.
481 In the presence of those constraints, the value returned is only
482 a inexpensive approximation of the optimal value.
484 References
485 ----------
486 .. [1] Gould, Nicholas IM, Mary E. Hribar, and Jorge Nocedal.
487 "On the solution of equality constrained quadratic
488 programming problems arising in optimization."
489 SIAM Journal on Scientific Computing 23.4 (2001): 1376-1395.
490 """
491 CLOSE_TO_ZERO = 1e-25
493 n, = np.shape(c) # Number of parameters
494 m, = np.shape(b) # Number of constraints
496 # Initial Values
497 x = Y.dot(-b)
498 r = Z.dot(H.dot(x) + c)
499 g = Z.dot(r)
500 p = -g
502 # Store ``x`` value
503 if return_all:
504 allvecs = [x]
505 # Values for the first iteration
506 H_p = H.dot(p)
507 rt_g = norm(g)**2 # g.T g = r.T Z g = r.T g (ref [1]_ p.1389)
509 # If x > trust-region the problem does not have a solution.
510 tr_distance = trust_radius - norm(x)
511 if tr_distance < 0:
512 raise ValueError("Trust region problem does not have a solution.")
513 # If x == trust_radius, then x is the solution
514 # to the optimization problem, since x is the
515 # minimum norm solution to Ax=b.
516 elif tr_distance < CLOSE_TO_ZERO:
517 info = {'niter': 0, 'stop_cond': 2, 'hits_boundary': True}
518 if return_all:
519 allvecs.append(x)
520 info['allvecs'] = allvecs
521 return x, info
523 # Set default tolerance
524 if tol is None:
525 tol = max(min(0.01 * np.sqrt(rt_g), 0.1 * rt_g), CLOSE_TO_ZERO)
526 # Set default lower and upper bounds
527 if lb is None:
528 lb = np.full(n, -np.inf)
529 if ub is None:
530 ub = np.full(n, np.inf)
531 # Set maximum iterations
532 if max_iter is None:
533 max_iter = n-m
534 max_iter = min(max_iter, n-m)
535 # Set maximum infeasible iterations
536 if max_infeasible_iter is None:
537 max_infeasible_iter = n-m
539 hits_boundary = False
540 stop_cond = 1
541 counter = 0
542 last_feasible_x = np.zeros_like(x)
543 k = 0
544 for i in range(max_iter):
545 # Stop criteria - Tolerance : r.T g < tol
546 if rt_g < tol:
547 stop_cond = 4
548 break
549 k += 1
550 # Compute curvature
551 pt_H_p = H_p.dot(p)
552 # Stop criteria - Negative curvature
553 if pt_H_p <= 0:
554 if np.isinf(trust_radius):
555 raise ValueError("Negative curvature not allowed "
556 "for unrestricted problems.")
557 else:
558 # Find intersection with constraints
559 _, alpha, intersect = box_sphere_intersections(
560 x, p, lb, ub, trust_radius, entire_line=True)
561 # Update solution
562 if intersect:
563 x = x + alpha*p
564 # Reinforce variables are inside box constraints.
565 # This is only necessary because of roundoff errors.
566 x = reinforce_box_boundaries(x, lb, ub)
567 # Attribute information
568 stop_cond = 3
569 hits_boundary = True
570 break
572 # Get next step
573 alpha = rt_g / pt_H_p
574 x_next = x + alpha*p
576 # Stop criteria - Hits boundary
577 if np.linalg.norm(x_next) >= trust_radius:
578 # Find intersection with box constraints
579 _, theta, intersect = box_sphere_intersections(x, alpha*p, lb, ub,
580 trust_radius)
581 # Update solution
582 if intersect:
583 x = x + theta*alpha*p
584 # Reinforce variables are inside box constraints.
585 # This is only necessary because of roundoff errors.
586 x = reinforce_box_boundaries(x, lb, ub)
587 # Attribute information
588 stop_cond = 2
589 hits_boundary = True
590 break
592 # Check if ``x`` is inside the box and start counter if it is not.
593 if inside_box_boundaries(x_next, lb, ub):
594 counter = 0
595 else:
596 counter += 1
597 # Whenever outside box constraints keep looking for intersections.
598 if counter > 0:
599 _, theta, intersect = box_sphere_intersections(x, alpha*p, lb, ub,
600 trust_radius)
601 if intersect:
602 last_feasible_x = x + theta*alpha*p
603 # Reinforce variables are inside box constraints.
604 # This is only necessary because of roundoff errors.
605 last_feasible_x = reinforce_box_boundaries(last_feasible_x,
606 lb, ub)
607 counter = 0
608 # Stop after too many infeasible (regarding box constraints) iteration.
609 if counter > max_infeasible_iter:
610 break
611 # Store ``x_next`` value
612 if return_all:
613 allvecs.append(x_next)
615 # Update residual
616 r_next = r + alpha*H_p
617 # Project residual g+ = Z r+
618 g_next = Z.dot(r_next)
619 # Compute conjugate direction step d
620 rt_g_next = norm(g_next)**2 # g.T g = r.T g (ref [1]_ p.1389)
621 beta = rt_g_next / rt_g
622 p = - g_next + beta*p
623 # Prepare for next iteration
624 x = x_next
625 g = g_next
626 r = g_next
627 rt_g = norm(g)**2 # g.T g = r.T Z g = r.T g (ref [1]_ p.1389)
628 H_p = H.dot(p)
630 if not inside_box_boundaries(x, lb, ub):
631 x = last_feasible_x
632 hits_boundary = True
633 info = {'niter': k, 'stop_cond': stop_cond,
634 'hits_boundary': hits_boundary}
635 if return_all:
636 info['allvecs'] = allvecs
637 return x, info