scijit.optimize.lsq_linear¶

scijit.optimize.lsq_linear(A, b, bounds=(-inf, inf), method='trf', tol=1e-10, lsq_solver=None, lsmr_tol=None, max_iter=None, verbose=0, *, lsmr_maxiter=None)¶

Bounded-variable linear least squares.

Solves min 0.5 * ||A x - b||^2 subject to lb <= x <= ub, by the trust-region-reflective method by default and by Stark and Parker’s BVLS active-set method on request. Callable from Python and from inside @njit; both entries run the same compiled cores.

Parameters:

Aarray_like, shape (m, n): Design matrix. Dense only. A one-dimensional A is read as a single row.
barray_like, shape (m,): Right-hand side.
bounds2-tuple, optional: (lb, ub), exactly two elements. Each is a scalar or an array of length n. Default (-inf, inf), unbounded.
method{‘trf’, ‘bvls’}, optional: 'trf', the default, is the trust-region-reflective solver, whose iterates stay strictly inside the box. 'bvls' is the active-set solver, which reaches the bounds exactly. Inside @njit this must be a literal string, since the two report nit and active_mask differently.
tolfloat, optional: Termination tolerance on the first-order optimality measure and on the relative change in cost. Default 1e-10.
lsq_solver{None, ‘exact’}, optional: The inner solve. None and 'exact' both run the dense pivoted QR. 'lsmr' raises: on method='bvls' it is not a valid pairing, and on method='trf' it selects an iterative solver this package does not have. Any other string raises.
lsmr_tolNone, ‘auto’ or float, optional: Tunes the LSMR inner solver, which the 'lsmr' value of lsq_solver selects, so it is not read here. Validated to None, 'auto' or a positive float; an integer is refused whatever its value. Default None.
max_iterint or None, optional: Iteration cap. None (default) resolves to 100 on 'trf' and to n on 'bvls'. Anything at or below zero raises.
verbose{0, 1, 2}, optional: Progress reporting on stdout. 0, the default, is silent. 1 prints the termination message and a summary line after the solve. 2 adds a five-column table with one row per iteration. See Notes.
lsmr_maxiterint or None, keyword-only, optional: Tunes the LSMR inner solver, so it is not read here. Validated to None or an integer at least 1. Default None.

Returns:

resLsqLinearResult: x, fun, cost, optimality, active_mask, nit, status, unbounded_sol, message, success, reached as attributes. status is -1 no progress on the last iteration, reachable on 'trf' only, 0 iteration limit, 1 optimality below tol, 2 cost change below tol, 3 the unconstrained solution was already optimal. optimality is the uniform norm of the scaled gradient on 'trf' and the KKT measure on 'bvls'. active_mask is -1 on a lower bound, +1 on an upper one and 0 free, and is all zero whenever status is 3. On 'trf' every iterate is strictly interior, so the mask is decided within tol of a bound rather than on equality. unbounded_sol is the whole numpy.linalg.lstsq tuple, (solution, residuals, rank, singular_values), taken at rcond=-1.

Raises:

ValueError: If method is neither 'trf' nor 'bvls'; if lsq_solver is not None, 'exact' or 'lsmr', or is 'lsmr'; if verbose is outside 0, 1, 2; if A has more than two dimensions; if max_iter is at or below zero; if b has more than one dimension or a length other than A.shape[0]; if bounds does not hold exactly two elements, a bound has the wrong length, or a lower bound is not strictly below its upper one; if lsmr_maxiter is below 1 or lsmr_tol is neither None, 'auto' nor a positive float; or if either dimension of A is zero.
TypeError: If bounds has no length at all, which is CPython’s own message.
numpy.linalg.LinAlgError: On a non-finite entry in A or b, from the unconstrained solve.

See also

scipy.optimize.lsq_linear: The scipy routine this mirrors.
scijit.optimize.nnls: The x >= 0 special case.
scijit.optimize.leastsq: Nonlinear least squares.
scijit.optimize.curve_fit: Fit a model to data.

Notes

The two methods solve the same problem and stop on different criteria, so they land in different places and report different accounting. On the 3x2 problem below, 'trf' takes 16 iterations to x[0] = 1.4999954223632812 and 'bvls' takes 2 to an exact 1.5, the two 4.578e-06 apart. scipy behaves the same way, and the default is 'trf' on both sides.

res carries the ten fields in ONE order on every exit path: unbounded_sol sits between active_mask and nit. scipy has two orders, that one on the status = 3 short circuit and a second elsewhere with unbounded_sol after status, because its two paths build the object differently. Reaching a field by name or by attribute is unaffected on both sides; what differs is list(res.keys()) on the paths where status is not 3.

unbounded_sol’s third member, the rank, is an int64 here and an int32 in scipy.

A sparse or LinearOperator A is not accepted, and neither is a bounds given as a scipy.optimize.Bounds object; pass the pair, nor lsq_solver='lsmr', which selects an iterative inner solve.

verbose above 0 prints, and printing takes the GIL, so a prange loop over solves is serialized for the duration of each line. verbose=0 reaches no printing code.

An A with a zero dimension raises. scipy 1.18 answers instead: shape (12, 0) returns [] while LAPACK prints On entry to DLASCL parameter number 4 had an illegal value. That answer is produced through an invalid LAPACK call.

No callback of either style and no module state, so it is safe to call from a numba.prange loop; 'trf' reaches LAPACK dgeqp3 through scijit.linalg.qr_pivot, which is prange-safe too (verified: 24 concurrent solves reproduce the serial answer to 2.78e-16).

https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.lsq_linear.html

Examples

>>> import numpy as np
>>> from numba import njit
>>> from scijit.optimize import lsq_linear
>>> A = np.array([[1.0, 0.0], [1.0, 0.0], [0.0, 1.0]])
>>> b = np.array([2.0, 1.0, -1.0])
>>> @njit
... def run():
...     return lsq_linear(A, b, (np.array([0.0, 0.0]),
...                              np.array([1.5, 1.0])))
>>> res = run()
>>> res.x
array([1.49999542e+000, 4.94065646e-324])
>>> res.status, res.success
(1, True)

The active-set method reaches the bound exactly, in two iterations:

>>> bvls = lsq_linear(A, b, (np.array([0.0, 0.0]),
...                          np.array([1.5, 1.0])), 'bvls')
>>> bvls.x, bvls.nit
(array([1.5, 0. ]), 2)