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