Coverage for /home/martinb/.local/share/virtualenvs/camcops/lib/python3.6/site-packages/scipy/optimize/_lsq/lsq_linear.py : 14%

Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1"""Linear least squares with bound constraints on independent variables."""
2import numpy as np
3from numpy.linalg import norm
4from scipy.sparse import issparse, csr_matrix
5from scipy.sparse.linalg import LinearOperator, lsmr
6from scipy.optimize import OptimizeResult
8from .common import in_bounds, compute_grad
9from .trf_linear import trf_linear
10from .bvls import bvls
13def prepare_bounds(bounds, n):
14 lb, ub = [np.asarray(b, dtype=float) for b in bounds]
16 if lb.ndim == 0:
17 lb = np.resize(lb, n)
19 if ub.ndim == 0:
20 ub = np.resize(ub, n)
22 return lb, ub
25TERMINATION_MESSAGES = {
26 -1: "The algorithm was not able to make progress on the last iteration.",
27 0: "The maximum number of iterations is exceeded.",
28 1: "The first-order optimality measure is less than `tol`.",
29 2: "The relative change of the cost function is less than `tol`.",
30 3: "The unconstrained solution is optimal."
31}
34def lsq_linear(A, b, bounds=(-np.inf, np.inf), method='trf', tol=1e-10,
35 lsq_solver=None, lsmr_tol=None, max_iter=None, verbose=0):
36 r"""Solve a linear least-squares problem with bounds on the variables.
38 Given a m-by-n design matrix A and a target vector b with m elements,
39 `lsq_linear` solves the following optimization problem::
41 minimize 0.5 * ||A x - b||**2
42 subject to lb <= x <= ub
44 This optimization problem is convex, hence a found minimum (if iterations
45 have converged) is guaranteed to be global.
47 Parameters
48 ----------
49 A : array_like, sparse matrix of LinearOperator, shape (m, n)
50 Design matrix. Can be `scipy.sparse.linalg.LinearOperator`.
51 b : array_like, shape (m,)
52 Target vector.
53 bounds : 2-tuple of array_like, optional
54 Lower and upper bounds on independent variables. Defaults to no bounds.
55 Each array must have shape (n,) or be a scalar, in the latter
56 case a bound will be the same for all variables. Use ``np.inf`` with
57 an appropriate sign to disable bounds on all or some variables.
58 method : 'trf' or 'bvls', optional
59 Method to perform minimization.
61 * 'trf' : Trust Region Reflective algorithm adapted for a linear
62 least-squares problem. This is an interior-point-like method
63 and the required number of iterations is weakly correlated with
64 the number of variables.
65 * 'bvls' : Bounded-variable least-squares algorithm. This is
66 an active set method, which requires the number of iterations
67 comparable to the number of variables. Can't be used when `A` is
68 sparse or LinearOperator.
70 Default is 'trf'.
71 tol : float, optional
72 Tolerance parameter. The algorithm terminates if a relative change
73 of the cost function is less than `tol` on the last iteration.
74 Additionally, the first-order optimality measure is considered:
76 * ``method='trf'`` terminates if the uniform norm of the gradient,
77 scaled to account for the presence of the bounds, is less than
78 `tol`.
79 * ``method='bvls'`` terminates if Karush-Kuhn-Tucker conditions
80 are satisfied within `tol` tolerance.
82 lsq_solver : {None, 'exact', 'lsmr'}, optional
83 Method of solving unbounded least-squares problems throughout
84 iterations:
86 * 'exact' : Use dense QR or SVD decomposition approach. Can't be
87 used when `A` is sparse or LinearOperator.
88 * 'lsmr' : Use `scipy.sparse.linalg.lsmr` iterative procedure
89 which requires only matrix-vector product evaluations. Can't
90 be used with ``method='bvls'``.
92 If None (default), the solver is chosen based on type of `A`.
93 lsmr_tol : None, float or 'auto', optional
94 Tolerance parameters 'atol' and 'btol' for `scipy.sparse.linalg.lsmr`
95 If None (default), it is set to ``1e-2 * tol``. If 'auto', the
96 tolerance will be adjusted based on the optimality of the current
97 iterate, which can speed up the optimization process, but is not always
98 reliable.
99 max_iter : None or int, optional
100 Maximum number of iterations before termination. If None (default), it
101 is set to 100 for ``method='trf'`` or to the number of variables for
102 ``method='bvls'`` (not counting iterations for 'bvls' initialization).
103 verbose : {0, 1, 2}, optional
104 Level of algorithm's verbosity:
106 * 0 : work silently (default).
107 * 1 : display a termination report.
108 * 2 : display progress during iterations.
110 Returns
111 -------
112 OptimizeResult with the following fields defined:
113 x : ndarray, shape (n,)
114 Solution found.
115 cost : float
116 Value of the cost function at the solution.
117 fun : ndarray, shape (m,)
118 Vector of residuals at the solution.
119 optimality : float
120 First-order optimality measure. The exact meaning depends on `method`,
121 refer to the description of `tol` parameter.
122 active_mask : ndarray of int, shape (n,)
123 Each component shows whether a corresponding constraint is active
124 (that is, whether a variable is at the bound):
126 * 0 : a constraint is not active.
127 * -1 : a lower bound is active.
128 * 1 : an upper bound is active.
130 Might be somewhat arbitrary for the `trf` method as it generates a
131 sequence of strictly feasible iterates and active_mask is determined
132 within a tolerance threshold.
133 nit : int
134 Number of iterations. Zero if the unconstrained solution is optimal.
135 status : int
136 Reason for algorithm termination:
138 * -1 : the algorithm was not able to make progress on the last
139 iteration.
140 * 0 : the maximum number of iterations is exceeded.
141 * 1 : the first-order optimality measure is less than `tol`.
142 * 2 : the relative change of the cost function is less than `tol`.
143 * 3 : the unconstrained solution is optimal.
145 message : str
146 Verbal description of the termination reason.
147 success : bool
148 True if one of the convergence criteria is satisfied (`status` > 0).
150 See Also
151 --------
152 nnls : Linear least squares with non-negativity constraint.
153 least_squares : Nonlinear least squares with bounds on the variables.
155 Notes
156 -----
157 The algorithm first computes the unconstrained least-squares solution by
158 `numpy.linalg.lstsq` or `scipy.sparse.linalg.lsmr` depending on
159 `lsq_solver`. This solution is returned as optimal if it lies within the
160 bounds.
162 Method 'trf' runs the adaptation of the algorithm described in [STIR]_ for
163 a linear least-squares problem. The iterations are essentially the same as
164 in the nonlinear least-squares algorithm, but as the quadratic function
165 model is always accurate, we don't need to track or modify the radius of
166 a trust region. The line search (backtracking) is used as a safety net
167 when a selected step does not decrease the cost function. Read more
168 detailed description of the algorithm in `scipy.optimize.least_squares`.
170 Method 'bvls' runs a Python implementation of the algorithm described in
171 [BVLS]_. The algorithm maintains active and free sets of variables, on
172 each iteration chooses a new variable to move from the active set to the
173 free set and then solves the unconstrained least-squares problem on free
174 variables. This algorithm is guaranteed to give an accurate solution
175 eventually, but may require up to n iterations for a problem with n
176 variables. Additionally, an ad-hoc initialization procedure is
177 implemented, that determines which variables to set free or active
178 initially. It takes some number of iterations before actual BVLS starts,
179 but can significantly reduce the number of further iterations.
181 References
182 ----------
183 .. [STIR] M. A. Branch, T. F. Coleman, and Y. Li, "A Subspace, Interior,
184 and Conjugate Gradient Method for Large-Scale Bound-Constrained
185 Minimization Problems," SIAM Journal on Scientific Computing,
186 Vol. 21, Number 1, pp 1-23, 1999.
187 .. [BVLS] P. B. Start and R. L. Parker, "Bounded-Variable Least-Squares:
188 an Algorithm and Applications", Computational Statistics, 10,
189 129-141, 1995.
191 Examples
192 --------
193 In this example, a problem with a large sparse matrix and bounds on the
194 variables is solved.
196 >>> from scipy.sparse import rand
197 >>> from scipy.optimize import lsq_linear
198 ...
199 >>> np.random.seed(0)
200 ...
201 >>> m = 20000
202 >>> n = 10000
203 ...
204 >>> A = rand(m, n, density=1e-4)
205 >>> b = np.random.randn(m)
206 ...
207 >>> lb = np.random.randn(n)
208 >>> ub = lb + 1
209 ...
210 >>> res = lsq_linear(A, b, bounds=(lb, ub), lsmr_tol='auto', verbose=1)
211 # may vary
212 The relative change of the cost function is less than `tol`.
213 Number of iterations 16, initial cost 1.5039e+04, final cost 1.1112e+04,
214 first-order optimality 4.66e-08.
215 """
216 if method not in ['trf', 'bvls']:
217 raise ValueError("`method` must be 'trf' or 'bvls'")
219 if lsq_solver not in [None, 'exact', 'lsmr']:
220 raise ValueError("`solver` must be None, 'exact' or 'lsmr'.")
222 if verbose not in [0, 1, 2]:
223 raise ValueError("`verbose` must be in [0, 1, 2].")
225 if issparse(A):
226 A = csr_matrix(A)
227 elif not isinstance(A, LinearOperator):
228 A = np.atleast_2d(A)
230 if method == 'bvls':
231 if lsq_solver == 'lsmr':
232 raise ValueError("method='bvls' can't be used with "
233 "lsq_solver='lsmr'")
235 if not isinstance(A, np.ndarray):
236 raise ValueError("method='bvls' can't be used with `A` being "
237 "sparse or LinearOperator.")
239 if lsq_solver is None:
240 if isinstance(A, np.ndarray):
241 lsq_solver = 'exact'
242 else:
243 lsq_solver = 'lsmr'
244 elif lsq_solver == 'exact' and not isinstance(A, np.ndarray):
245 raise ValueError("`exact` solver can't be used when `A` is "
246 "sparse or LinearOperator.")
248 if len(A.shape) != 2: # No ndim for LinearOperator.
249 raise ValueError("`A` must have at most 2 dimensions.")
251 if len(bounds) != 2:
252 raise ValueError("`bounds` must contain 2 elements.")
254 if max_iter is not None and max_iter <= 0:
255 raise ValueError("`max_iter` must be None or positive integer.")
257 m, n = A.shape
259 b = np.atleast_1d(b)
260 if b.ndim != 1:
261 raise ValueError("`b` must have at most 1 dimension.")
263 if b.size != m:
264 raise ValueError("Inconsistent shapes between `A` and `b`.")
266 lb, ub = prepare_bounds(bounds, n)
268 if lb.shape != (n,) and ub.shape != (n,):
269 raise ValueError("Bounds have wrong shape.")
271 if np.any(lb >= ub):
272 raise ValueError("Each lower bound must be strictly less than each "
273 "upper bound.")
275 if lsq_solver == 'exact':
276 x_lsq = np.linalg.lstsq(A, b, rcond=-1)[0]
277 elif lsq_solver == 'lsmr':
278 x_lsq = lsmr(A, b, atol=tol, btol=tol)[0]
280 if in_bounds(x_lsq, lb, ub):
281 r = A.dot(x_lsq) - b
282 cost = 0.5 * np.dot(r, r)
283 termination_status = 3
284 termination_message = TERMINATION_MESSAGES[termination_status]
285 g = compute_grad(A, r)
286 g_norm = norm(g, ord=np.inf)
288 if verbose > 0:
289 print(termination_message)
290 print("Final cost {0:.4e}, first-order optimality {1:.2e}"
291 .format(cost, g_norm))
293 return OptimizeResult(
294 x=x_lsq, fun=r, cost=cost, optimality=g_norm,
295 active_mask=np.zeros(n), nit=0, status=termination_status,
296 message=termination_message, success=True)
298 if method == 'trf':
299 res = trf_linear(A, b, x_lsq, lb, ub, tol, lsq_solver, lsmr_tol,
300 max_iter, verbose)
301 elif method == 'bvls':
302 res = bvls(A, b, x_lsq, lb, ub, tol, max_iter, verbose)
304 res.message = TERMINATION_MESSAGES[res.status]
305 res.success = res.status > 0
307 if verbose > 0:
308 print(res.message)
309 print("Number of iterations {0}, initial cost {1:.4e}, "
310 "final cost {2:.4e}, first-order optimality {3:.2e}."
311 .format(res.nit, res.initial_cost, res.cost, res.optimality))
313 del res.initial_cost
315 return res