Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/optimize/_lsq/lsq_linear.py: 12%
90 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
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,
36 verbose=0, *, lsmr_maxiter=None,):
37 r"""Solve a linear least-squares problem with bounds on the variables.
39 Given a m-by-n design matrix A and a target vector b with m elements,
40 `lsq_linear` solves the following optimization problem::
42 minimize 0.5 * ||A x - b||**2
43 subject to lb <= x <= ub
45 This optimization problem is convex, hence a found minimum (if iterations
46 have converged) is guaranteed to be global.
48 Parameters
49 ----------
50 A : array_like, sparse matrix of LinearOperator, shape (m, n)
51 Design matrix. Can be `scipy.sparse.linalg.LinearOperator`.
52 b : array_like, shape (m,)
53 Target vector.
54 bounds : 2-tuple of array_like, optional
55 Lower and upper bounds on independent variables. Defaults to no bounds.
56 Each array must have shape (n,) or be a scalar, in the latter
57 case a bound will be the same for all variables. Use ``np.inf`` with
58 an appropriate sign to disable bounds on all or some variables.
59 method : 'trf' or 'bvls', optional
60 Method to perform minimization.
62 * 'trf' : Trust Region Reflective algorithm adapted for a linear
63 least-squares problem. This is an interior-point-like method
64 and the required number of iterations is weakly correlated with
65 the number of variables.
66 * 'bvls' : Bounded-variable least-squares algorithm. This is
67 an active set method, which requires the number of iterations
68 comparable to the number of variables. Can't be used when `A` is
69 sparse or LinearOperator.
71 Default is 'trf'.
72 tol : float, optional
73 Tolerance parameter. The algorithm terminates if a relative change
74 of the cost function is less than `tol` on the last iteration.
75 Additionally, the first-order optimality measure is considered:
77 * ``method='trf'`` terminates if the uniform norm of the gradient,
78 scaled to account for the presence of the bounds, is less than
79 `tol`.
80 * ``method='bvls'`` terminates if Karush-Kuhn-Tucker conditions
81 are satisfied within `tol` tolerance.
83 lsq_solver : {None, 'exact', 'lsmr'}, optional
84 Method of solving unbounded least-squares problems throughout
85 iterations:
87 * 'exact' : Use dense QR or SVD decomposition approach. Can't be
88 used when `A` is sparse or LinearOperator.
89 * 'lsmr' : Use `scipy.sparse.linalg.lsmr` iterative procedure
90 which requires only matrix-vector product evaluations. Can't
91 be used with ``method='bvls'``.
93 If None (default), the solver is chosen based on type of `A`.
94 lsmr_tol : None, float or 'auto', optional
95 Tolerance parameters 'atol' and 'btol' for `scipy.sparse.linalg.lsmr`
96 If None (default), it is set to ``1e-2 * tol``. If 'auto', the
97 tolerance will be adjusted based on the optimality of the current
98 iterate, which can speed up the optimization process, but is not always
99 reliable.
100 max_iter : None or int, optional
101 Maximum number of iterations before termination. If None (default), it
102 is set to 100 for ``method='trf'`` or to the number of variables for
103 ``method='bvls'`` (not counting iterations for 'bvls' initialization).
104 verbose : {0, 1, 2}, optional
105 Level of algorithm's verbosity:
107 * 0 : work silently (default).
108 * 1 : display a termination report.
109 * 2 : display progress during iterations.
110 lsmr_maxiter : None or int, optional
111 Maximum number of iterations for the lsmr least squares solver,
112 if it is used (by setting ``lsq_solver='lsmr'``). If None (default), it
113 uses lsmr's default of ``min(m, n)`` where ``m`` and ``n`` are the
114 number of rows and columns of `A`, respectively. Has no effect if
115 ``lsq_solver='exact'``.
117 Returns
118 -------
119 OptimizeResult with the following fields defined:
120 x : ndarray, shape (n,)
121 Solution found.
122 cost : float
123 Value of the cost function at the solution.
124 fun : ndarray, shape (m,)
125 Vector of residuals at the solution.
126 optimality : float
127 First-order optimality measure. The exact meaning depends on `method`,
128 refer to the description of `tol` parameter.
129 active_mask : ndarray of int, shape (n,)
130 Each component shows whether a corresponding constraint is active
131 (that is, whether a variable is at the bound):
133 * 0 : a constraint is not active.
134 * -1 : a lower bound is active.
135 * 1 : an upper bound is active.
137 Might be somewhat arbitrary for the `trf` method as it generates a
138 sequence of strictly feasible iterates and active_mask is determined
139 within a tolerance threshold.
140 unbounded_sol : tuple
141 Unbounded least squares solution tuple returned by the least squares
142 solver (set with `lsq_solver` option). If `lsq_solver` is not set or is
143 set to ``'exact'``, the tuple contains an ndarray of shape (n,) with
144 the unbounded solution, an ndarray with the sum of squared residuals,
145 an int with the rank of `A`, and an ndarray with the singular values
146 of `A` (see NumPy's ``linalg.lstsq`` for more information). If
147 `lsq_solver` is set to ``'lsmr'``, the tuple contains an ndarray of
148 shape (n,) with the unbounded solution, an int with the exit code,
149 an int with the number of iterations, and five floats with
150 various norms and the condition number of `A` (see SciPy's
151 ``sparse.linalg.lsmr`` for more information). This output can be
152 useful for determining the convergence of the least squares solver,
153 particularly the iterative ``'lsmr'`` solver. The unbounded least
154 squares problem is to minimize ``0.5 * ||A x - b||**2``.
155 nit : int
156 Number of iterations. Zero if the unconstrained solution is optimal.
157 status : int
158 Reason for algorithm termination:
160 * -1 : the algorithm was not able to make progress on the last
161 iteration.
162 * 0 : the maximum number of iterations is exceeded.
163 * 1 : the first-order optimality measure is less than `tol`.
164 * 2 : the relative change of the cost function is less than `tol`.
165 * 3 : the unconstrained solution is optimal.
167 message : str
168 Verbal description of the termination reason.
169 success : bool
170 True if one of the convergence criteria is satisfied (`status` > 0).
172 See Also
173 --------
174 nnls : Linear least squares with non-negativity constraint.
175 least_squares : Nonlinear least squares with bounds on the variables.
177 Notes
178 -----
179 The algorithm first computes the unconstrained least-squares solution by
180 `numpy.linalg.lstsq` or `scipy.sparse.linalg.lsmr` depending on
181 `lsq_solver`. This solution is returned as optimal if it lies within the
182 bounds.
184 Method 'trf' runs the adaptation of the algorithm described in [STIR]_ for
185 a linear least-squares problem. The iterations are essentially the same as
186 in the nonlinear least-squares algorithm, but as the quadratic function
187 model is always accurate, we don't need to track or modify the radius of
188 a trust region. The line search (backtracking) is used as a safety net
189 when a selected step does not decrease the cost function. Read more
190 detailed description of the algorithm in `scipy.optimize.least_squares`.
192 Method 'bvls' runs a Python implementation of the algorithm described in
193 [BVLS]_. The algorithm maintains active and free sets of variables, on
194 each iteration chooses a new variable to move from the active set to the
195 free set and then solves the unconstrained least-squares problem on free
196 variables. This algorithm is guaranteed to give an accurate solution
197 eventually, but may require up to n iterations for a problem with n
198 variables. Additionally, an ad-hoc initialization procedure is
199 implemented, that determines which variables to set free or active
200 initially. It takes some number of iterations before actual BVLS starts,
201 but can significantly reduce the number of further iterations.
203 References
204 ----------
205 .. [STIR] M. A. Branch, T. F. Coleman, and Y. Li, "A Subspace, Interior,
206 and Conjugate Gradient Method for Large-Scale Bound-Constrained
207 Minimization Problems," SIAM Journal on Scientific Computing,
208 Vol. 21, Number 1, pp 1-23, 1999.
209 .. [BVLS] P. B. Start and R. L. Parker, "Bounded-Variable Least-Squares:
210 an Algorithm and Applications", Computational Statistics, 10,
211 129-141, 1995.
213 Examples
214 --------
215 In this example, a problem with a large sparse matrix and bounds on the
216 variables is solved.
218 >>> import numpy as np
219 >>> from scipy.sparse import rand
220 >>> from scipy.optimize import lsq_linear
221 >>> rng = np.random.default_rng()
222 ...
223 >>> m = 20000
224 >>> n = 10000
225 ...
226 >>> A = rand(m, n, density=1e-4, random_state=rng)
227 >>> b = rng.standard_normal(m)
228 ...
229 >>> lb = rng.standard_normal(n)
230 >>> ub = lb + 1
231 ...
232 >>> res = lsq_linear(A, b, bounds=(lb, ub), lsmr_tol='auto', verbose=1)
233 # may vary
234 The relative change of the cost function is less than `tol`.
235 Number of iterations 16, initial cost 1.5039e+04, final cost 1.1112e+04,
236 first-order optimality 4.66e-08.
237 """
238 if method not in ['trf', 'bvls']:
239 raise ValueError("`method` must be 'trf' or 'bvls'")
241 if lsq_solver not in [None, 'exact', 'lsmr']:
242 raise ValueError("`solver` must be None, 'exact' or 'lsmr'.")
244 if verbose not in [0, 1, 2]:
245 raise ValueError("`verbose` must be in [0, 1, 2].")
247 if issparse(A):
248 A = csr_matrix(A)
249 elif not isinstance(A, LinearOperator):
250 A = np.atleast_2d(np.asarray(A))
252 if method == 'bvls':
253 if lsq_solver == 'lsmr':
254 raise ValueError("method='bvls' can't be used with "
255 "lsq_solver='lsmr'")
257 if not isinstance(A, np.ndarray):
258 raise ValueError("method='bvls' can't be used with `A` being "
259 "sparse or LinearOperator.")
261 if lsq_solver is None:
262 if isinstance(A, np.ndarray):
263 lsq_solver = 'exact'
264 else:
265 lsq_solver = 'lsmr'
266 elif lsq_solver == 'exact' and not isinstance(A, np.ndarray):
267 raise ValueError("`exact` solver can't be used when `A` is "
268 "sparse or LinearOperator.")
270 if len(A.shape) != 2: # No ndim for LinearOperator.
271 raise ValueError("`A` must have at most 2 dimensions.")
273 if len(bounds) != 2:
274 raise ValueError("`bounds` must contain 2 elements.")
276 if max_iter is not None and max_iter <= 0:
277 raise ValueError("`max_iter` must be None or positive integer.")
279 m, n = A.shape
281 b = np.atleast_1d(b)
282 if b.ndim != 1:
283 raise ValueError("`b` must have at most 1 dimension.")
285 if b.size != m:
286 raise ValueError("Inconsistent shapes between `A` and `b`.")
288 lb, ub = prepare_bounds(bounds, n)
290 if lb.shape != (n,) and ub.shape != (n,):
291 raise ValueError("Bounds have wrong shape.")
293 if np.any(lb >= ub):
294 raise ValueError("Each lower bound must be strictly less than each "
295 "upper bound.")
297 if lsmr_maxiter is not None and lsmr_maxiter < 1:
298 raise ValueError("`lsmr_maxiter` must be None or positive integer.")
300 if not ((isinstance(lsmr_tol, float) and lsmr_tol > 0) or
301 lsmr_tol in ('auto', None)):
302 raise ValueError("`lsmr_tol` must be None, 'auto', or positive float.")
304 if lsq_solver == 'exact':
305 unbd_lsq = np.linalg.lstsq(A, b, rcond=-1)
306 elif lsq_solver == 'lsmr':
307 first_lsmr_tol = lsmr_tol # tol of first call to lsmr
308 if lsmr_tol is None or lsmr_tol == 'auto':
309 first_lsmr_tol = 1e-2 * tol # default if lsmr_tol not defined
310 unbd_lsq = lsmr(A, b, maxiter=lsmr_maxiter,
311 atol=first_lsmr_tol, btol=first_lsmr_tol)
312 x_lsq = unbd_lsq[0] # extract the solution from the least squares solver
314 if in_bounds(x_lsq, lb, ub):
315 r = A @ x_lsq - b
316 cost = 0.5 * np.dot(r, r)
317 termination_status = 3
318 termination_message = TERMINATION_MESSAGES[termination_status]
319 g = compute_grad(A, r)
320 g_norm = norm(g, ord=np.inf)
322 if verbose > 0:
323 print(termination_message)
324 print("Final cost {0:.4e}, first-order optimality {1:.2e}"
325 .format(cost, g_norm))
327 return OptimizeResult(
328 x=x_lsq, fun=r, cost=cost, optimality=g_norm,
329 active_mask=np.zeros(n), unbounded_sol=unbd_lsq,
330 nit=0, status=termination_status,
331 message=termination_message, success=True)
333 if method == 'trf':
334 res = trf_linear(A, b, x_lsq, lb, ub, tol, lsq_solver, lsmr_tol,
335 max_iter, verbose, lsmr_maxiter=lsmr_maxiter)
336 elif method == 'bvls':
337 res = bvls(A, b, x_lsq, lb, ub, tol, max_iter, verbose)
339 res.unbounded_sol = unbd_lsq
340 res.message = TERMINATION_MESSAGES[res.status]
341 res.success = res.status > 0
343 if verbose > 0:
344 print(res.message)
345 print("Number of iterations {0}, initial cost {1:.4e}, "
346 "final cost {2:.4e}, first-order optimality {3:.2e}."
347 .format(res.nit, res.initial_cost, res.cost, res.optimality))
349 del res.initial_cost
351 return res