imsl.linalg.lin_lsq_lin_constraints

lin_lsq_lin_constraints(a, b, c, bl, bu, con_type, xlb=None, xub=None, abs_tol=None, rel_tol=None, max_iter=None)

Solve a linear least-squares problem with linear constraints.

Parameters:

a : (M, N) array_like

Array containing the coefficients of the M least-squares equations.

b : (M,) array_like

Array containing the right-hand sides of the least-squares equations.

c : (K, N) array_like

Array containing the coefficients of the K general constraints. If the problem contains no general constraints, set c = None.

bl : (K,) array_like

Array containing the lower limits of the general constraints. If there is no lower limit on the i-th constraint, then bl[i] will not be referenced.

bu : (K,) array_like

Array containing the upper limits of the general constraints. If there is no upper limit on the i-th constraint, then bl[i] will not be referenced.

con_type : (K,) array_like

Array indicating the type of the general constraints, where con_type[i] = 0, 1, 2, 3 indicates =, <=, >= and range constraints, respectively.

xlb : (N,) array_like, optional

Array containing the lower bounds on the variables. If there is no lower bound on the i-th variable, then xlb[i] should be set to constant UNBOUNDED_ABOVE in module imsl.constants.

Default: The variables have no lower bounds.

xub : (N,) array_like, optional

Array containing the upper bounds on the variables. If there is no upper bound on the i-th variable, then xlb[i] should be set to constant UNBOUNDED_BELOW in module imsl.constants.

Default: The variables have no upper bounds.

abs_tol : float, optional

Absolute rank determination tolerance to be used.

Default: abs_tol = math.sqrt(numpy.finfo(numpy.float64).eps).

rel_tol : float, optional

Relative rank determination tolerance to be used.

Default: rel_tol = math.sqrt(numpy.finfo(numpy.float64).eps).

max_iter : int, optional

The maximum number of add/drop iterations.

Default: max_iter = 5 * max(M,N).

Returns:

A named tuple with the following fields:

solution : (N,) ndarray

The approximate solution of the least-squares problem.

residual : (M,) ndarray

The residuals b-Ax of the least-squares equations at the approximate solution.

Notes

Function lin_lsq_lin_constraints solves linear least-squares problems with linear constraints. These are systems of least-squares equations of the form

\[Ax \cong b\]

subject to

\[\begin{split}b_l \le Cx \le b_u \\ x_l \le x \le x_u\end{split}\]

Here A is the coefficient matrix of the least-squares equations, b is the right-hand side, and C is the coefficient matrix of the constraints. The vectors \(b_l, b_u, x_l\) and \(x_u\) are the lower and upper bounds on the constraints and the variables, respectively. The system is solved by defining dependent variables \(y \equiv Cx\) and then solving the least-squares system with the lower and upper bounds on x and y. The equation \(Cx-y=0\) is a set of equality constraints. These constraints are realized by heavy weighting, i.e., a penalty method ([R14]).

References

[R14]

(1, 2) Hanson, Richard J. (1986), Linear Least Squares with Bounds and Linear Constraints, SIAM J. Sci. and Stat. Computing, 7(3), 826-834.

Examples

The following problem is solved in the least-squares sense:

\[\begin{split}3x_1+2x_2+x_3=3.3 \\ 4x_1+2x_2+x_3=2.3 \\ 2x_1+2x_2+x_3=1.3 \\ x_1+x_2+x_3=1.0\end{split}\]

subject to

\[\begin{split}x_1=x_2+x_3 \le 1 \\ 0 \le x_1 \le 0.5 \\ 0 \le x_2 \le 0.5 \\ 0 \le x_3 \le 0.5\end{split}\]

The approximate solution of the least-squares problem, the residuals of the least-squares equations at the solution and the norm of the residual vector are printed.

>>> import imsl.linalg as la
>>> import numpy as np
>>> a = [[3.0, 2.0, 1.0], [4.0, 2.0, 1.0], [2.0, 2.0, 1.0],
...      [1.0, 1.0, 1.0]]
>>> b = [3.3, 2.3, 1.3, 1.0]
>>> c = [[1.0, 1.0, 1.0]]
>>> xlb = [0.0, 0.0, 0.0]
>>> xub = [0.5, 0.5, 0.5]
>>> con_type = [1]
>>> bc = [1.0]
>>> x = la.lin_lsq_lin_constraints(a, b, c, bc, bc, con_type, xlb=xlb,
...                                xub=xub)
>>> print(x.solution) 
[ 0.5  0.3  0.2]
>>> print(x.residual) 
[-1.   0.5  0.5  0. ]
>>> print("{0:8.6f}".format(np.linalg.norm(x.residual)))
1.224745