naginterfaces.library.opt.lsq_​lincon_​solve

naginterfaces.library.opt.lsq_lincon_solve(bl, bu, x, comm, c=None, cvec=None, istate=None, kx=None, a=None, b=None, io_manager=None)

lsq_lincon_solve solves linearly constrained linear least squares problems and convex quadratic programming problems. It is not intended for large sparse problems.

Note: this function uses optional algorithmic parameters, see also: lsq_lincon_option_file(), lsq_lincon_option_string(), nlp1_init().

For full information please refer to the NAG Library document for e04nc

https://support.nag.com/numeric/nl/nagdoc_30.3/flhtml/e04/e04ncf.html

Parameters

blfloat, array-like, shape $(n+\text{nclin})$

$bl$ must contain the lower bounds for all the constraints

bufloat, array-like, shape $(n+\text{nclin})$

$bu$ must contain the upper bounds for all the constraints

xfloat, array-like, shape $\left(n\right)$

An initial estimate of the solution.

Note: that it may be best to avoid the choice $x=0.0$.

commdict, communication object, modified in place

Communication structure.

This argument must have been initialized by a prior call to nlp1_init().

cNone or float, array-like, shape $(\text{nclin},:)$, optional

Note: the required extent for this argument in dimension 2 is determined as follows: if $\text{nclin}>0$: $n$; otherwise: $1$.

The $\text{i}$th row of $c$ must contain the coefficients of the $\text{i}$th general constraint, for $\text{i}=1,2,\dots ,\text{nclin}$.

cvecNone or float, array-like, shape $\left(n\right)$, optional

The coefficients of the explicit linear term of the objective function.

If the problem is of type FP, QP1, QP3, LS1 (the default) or LS3, $cvec$ is not referenced.

istateNone or int, array-like, shape $(n+\text{nclin})$, optional

Need not be set if the (default) option ‘Cold Start’ is used.

If the option ‘Warm Start’ has been chosen, $istate$ specifies the desired status of the constraints at the start of the feasibility phase.

More precisely, the first $n$ elements of $istate$ refer to the upper and lower bounds on the variables, and the next $n_L$ elements refer to the general linear constraints (if any).

Possible values for $istate[j-1]$ are as follows:

$istate[j-1]$	Meaning
0	The constraint should not be in the initial working set.
1	The constraint should be in the initial working set at its lower bound.
2	The constraint should be in the initial working set at its upper bound.
3	The constraint should be in the initial working set as an equality. This value must not be specified unless $bl[j-1]=bu[j-1]$.

The values $-2$, $-1$ and $4$ are also acceptable but will be reset to zero by the function.

If lsq_lincon_solve has been called previously with the same values of $\text{n}$ and $\text{nclin}$, $istate$ already contains satisfactory information. (See also the description of the option ‘Warm Start’.) The function also adjusts (if necessary) the values supplied in $x$ to be consistent with $istate$.

kxNone or int, array-like, shape $\left(n\right)$, optional

Need not be initialized for problems of type FP, LP, QP1, QP2, LS1 (the default) or LS2.

For problems QP3, QP4, LS3 or LS4, $kx$ must specify the order of the columns of the matrix $A$ with respect to the ordering of $x$.

Thus if column $j$ of $A$ is the column associated with the variable $x_i$ then $kx[j-1]=i$.

aNone or float, array-like, shape $(m,n)$, optional

The matrix $A$.

bNone or float, array-like, shape $\left(m\right)$, optional

The $m$ elements of the vector of observations.

io_managerFileObjManager, optional

Manager for I/O in this routine.

Returns

istateint, ndarray, shape $(n+\text{nclin})$

The status of the constraints in the working set at the point returned in $x$. The significance of each possible value of $istate[j-1]$ is as follows:

$istate[j-1]$	Meaning
$-2$	The constraint violates its lower bound by more than the feasibility tolerance.
$-1$	The constraint violates its upper bound by more than the feasibility tolerance.
$0$	The constraint is satisfied to within the feasibility tolerance, but is not in the working set.
$1$	This inequality constraint is included in the working set at its lower bound.
$2$	This inequality constraint is included in the working set at its upper bound.
$3$	The constraint is included in the working set as an equality. This value of $istate$ can occur only when $bl[j-1]=bu[j-1]$.
$4$	This corresponds to optimality being declared with $x[j-1]$ being temporarily fixed at its current value.

kxint, ndarray, shape $\left(n\right)$

Defines the order of the columns of $a$ with respect to the ordering of $x$, as described above.

xfloat, ndarray, shape $\left(n\right)$

The point at which lsq_lincon_solve terminated. If the function exits successfully or $errno$ = 1 or 4, $x$ contains an estimate of the solution.

afloat, ndarray, shape $(m,n)$

The matrix $A$.

bNone or float, ndarray, shape $\left(m\right)$

The transformed residual vector of (0) (see Main Iteration).

If the problem is of type FP, LP, QP1, QP2, QP3 or QP4, $b$ is not referenced.

iteraint

The total number of iterations performed.

objfloat

The value of the objective function at $x$ if $x$ is feasible, or the sum of infeasibiliites at $x$ otherwise. If the problem is of type FP and $x$ is feasible, $obj$ is set to zero.

clamdafloat, ndarray, shape $(n+\text{nclin})$

The values of the Lagrange multipliers for each constraint with respect to the current working set. The first $n$ elements contain the multipliers for the bound constraints on the variables, and the next $n_L$ elements contain the multipliers for the general linear constraints (if any). If $istate[j-1]=0$ (i.e., constraint $j$ is not in the working set), $clamda[j-1]$ is zero. If $x$ is optimal, $clamda[j-1]$ should be non-negative if $istate[j-1]=1$, non-positive if $istate[j-1]=2$ and zero if $istate[j-1]=4$.

Other Parameters

‘Cold Start’valueless

Default

This option specifies how the initial working set is chosen. With a ‘Cold Start’, lsq_lincon_solve chooses the initial working set based on the values of the variables and constraints at the initial point. Broadly speaking, the initial working set will include equality constraints and bounds or inequality constraints that violate or ‘nearly’ satisfy their bounds (to within ‘Crash Tolerance’).

With a ‘Warm Start’, you must provide a valid definition of every element of the array $istate$. lsq_lincon_solve will override your specification of $istate$ if necessary, so that a poor choice of the working set will not cause a fatal error. For instance, any elements of $istate$ which are set to $-2$, $-1$ or $4$ will be reset to zero, as will any elements which are set to $3$ when the corresponding elements of $bl$ and $bu$ are not equal. A warm start will be advantageous if a good estimate of the initial working set is available – for example, when lsq_lincon_solve is called repeatedly to solve related problems.

‘Warm Start’valueless

This option specifies how the initial working set is chosen. With a ‘Cold Start’, lsq_lincon_solve chooses the initial working set based on the values of the variables and constraints at the initial point. Broadly speaking, the initial working set will include equality constraints and bounds or inequality constraints that violate or ‘nearly’ satisfy their bounds (to within ‘Crash Tolerance’).

With a ‘Warm Start’, you must provide a valid definition of every element of the array $istate$. lsq_lincon_solve will override your specification of $istate$ if necessary, so that a poor choice of the working set will not cause a fatal error. For instance, any elements of $istate$ which are set to $-2$, $-1$ or $4$ will be reset to zero, as will any elements which are set to $3$ when the corresponding elements of $bl$ and $bu$ are not equal. A warm start will be advantageous if a good estimate of the initial working set is available – for example, when lsq_lincon_solve is called repeatedly to solve related problems.

‘Crash Tolerance’float

Default $=0.01$

This value is used in conjunction with the option ‘Cold Start’ (the default value) when lsq_lincon_solve selects an initial working set. If $0\le r\le 1$, the initial working set will include (if possible) bounds or general inequality constraints that lie within $r$ of their bounds. In particular, a constraint of the form $c_j^{T}x\ge l$ will be included in the initial working set if $|c_j^{T}x-l|\le r(1+\left|l\right|)$. If $r<0$ or $r>1$, the default value is used.

‘Defaults’valueless

This special keyword may be used to reset all options to their default values.

‘Feasibility Phase Iteration Limit’int

Default $=max(50,5(n+n_L))$

The scalars $i_1$ and $i_2$ specify the maximum number of iterations allowed in the feasibility and optimality phases. Option ‘Optimality Phase Iteration Limit’ is equivalent to option ‘Iteration Limit’. Setting $i_2=0$ and $\text{‘Print Level'}>0$ means that the workspace needed will be computed and printed, but no iterations will be performed. If $i_1<0$ or $i_2<0$, the default value is used.

‘Optimality Phase Iteration Limit’int

Default $=max(50,5(n+n_L))$

The scalars $i_1$ and $i_2$ specify the maximum number of iterations allowed in the feasibility and optimality phases. Option ‘Optimality Phase Iteration Limit’ is equivalent to option ‘Iteration Limit’. Setting $i_2=0$ and $\text{‘Print Level'}>0$ means that the workspace needed will be computed and printed, but no iterations will be performed. If $i_1<0$ or $i_2<0$, the default value is used.

‘Feasibility Tolerance’float

Default $=\sqrt{ϵ}$

If $r>ϵ$, $r$ defines the maximum acceptable absolute violation in each constraint at a ‘feasible’ point. For example, if the variables and the coefficients in the general constraints are of order unity, and the latter are correct to about $6$ decimal digits, it would be appropriate to specify $r$ as ${10}^{-6}$. If $0\le r<ϵ$, the default value is used.

Note that a ‘feasible solution’ is a solution that satisfies the current constraints to within the tolerance $r$.

‘Hessian’str

Default $=NO$

This option controls the contents of the upper triangular matrix $R$ (see the description of $a$ in Parameters). lsq_lincon_solve works exclusively with the transformed and reordered matrix $H_Q$ (8), and hence extra computation is required to form the Hessian itself. If $\text{‘Hessian'}=\text{'NO'}$, $a$ contains the Cholesky factor of the matrix $H_Q$ with columns ordered as indicated by $kx$ (see Parameters). If $\text{‘Hessian'}=\text{'YES'}$, $a$ contains the Cholesky factor of the matrix $H$, with columns ordered as indicated by $kx$.

‘Infinite Bound Size’float

Default $={10}^{20}$

If $r>0$, $r$ defines the ‘infinite’ bound $\text{bigbnd}$ in the definition of the problem constraints. Any upper bound greater than or equal to $\text{bigbnd}$ will be regarded as $+\infty$ (and similarly any lower bound less than or equal to $-\text{bigbnd}$ will be regarded as $-\infty$). If $r<0$, the default value is used.

‘Infinite Step Size’float

Default $=max(\text{bigbnd},{10}^{20})$

If $r>0$, $r$ specifies the magnitude of the change in variables that will be considered a step to an unbounded solution. (Note that an unbounded solution can occur only when the Hessian is singular and the objective contains an explicit linear term.) If the change in $x$ during an iteration would exceed the value of $r$, the objective function is considered to be unbounded below in the feasible region. If $r\le 0$, the default value is used.

‘Iteration Limit’int

Default $=max(50,5(n+n_L))$

See option ‘Feasibility Phase Iteration Limit’.

‘Iters’int

Default $=max(50,5(n+n_L))$

See option ‘Feasibility Phase Iteration Limit’.

‘Itns’int

Default $=max(50,5(n+n_L))$

See option ‘Feasibility Phase Iteration Limit’.

‘List’valueless

Option ‘List’ enables printing of each option specification as it is supplied. ‘Nolist’ suppresses this printing.

‘Nolist’valueless

Default $=\text{‘Nolist'}$

Option ‘List’ enables printing of each option specification as it is supplied. ‘Nolist’ suppresses this printing.

‘Monitoring File’int

Default $=-1$

If $i\ge 0$ and $\text{‘Print Level'}\ge 5$, monitoring information produced by lsq_lincon_solve at every iteration is sent to a file with logical unit number $i$. If $i<0$ and/or $\text{‘Print Level'}<5$, no monitoring information is produced.

‘Print Level’int

Default $=0$

The value of $i$ controls the amount of printout produced by lsq_lincon_solve, as indicated below. A detailed description of the printed output is given in Further Comments (summary output at each iteration and the final solution) and Monitoring Information (monitoring information at each iteration).

The following printout is sent to the file object associated with the advisory I/O unit (see FileObjManager):

$i$	Output
$0$	No output.
$1$	The final solution only.
$5$	One line of summary output ($<80$ characters; see Further Comments) for each iteration (no printout of the final solution).
$\ge 10$	The final solution and one line of summary output for each iteration.

The following printout is sent to the unit number given by the option ‘Monitoring File’:

$i$	Output
$<5$	No output.
$\ge 5$	One long line of output ($>80$ characters; see Monitoring Information) for each iteration (no printout of the final solution).
$\ge 20$	At each iteration, the Lagrange multipliers, the variables $x$, the constraint values $Cx$ and the constraint status.
$\ge 30$	At each iteration, the diagonal elements of the matrix $T$ associated with the $TQ$ factorization (4) (see Definition of Search Direction) of the working set, and the diagonal elements of the upper triangular matrix $R$.

If $\text{‘Print Level'}\ge 5$ and the unit number defined by the option ‘Monitoring File’ is the advisory unit number, the summary output for each major iteration is suppressed.

‘Problem Type’str

Default $=$ LS1

This option specifies the type of objective function to be minimized during the optimality phase. The following are the nine optional keywords and the dimensions of the arrays that must be specified in order to define the objective function:

LP	$a$ and $b$ not referenced, length-$\text{n}$ $cvec$;
QP1	shape-$(m,n)$ $a$ symmetric, $b$ and $cvec$ not referenced;
QP2	shape-$(m,n)$ $a$ symmetric, $b$ not referenced, length-$\text{n}$ $cvec$;
QP3	shape-$(m,n)$ $a$ upper trapezoidal, length-$\text{n}$ $kx$, $b$ and $cvec$ not referenced;
QP4	shape-$(m,n)$ $a$ upper trapezoidal, length-$\text{n}$ $kx$, $b$ not referenced, length-$\text{n}$ $cvec$;
LS1	shape-$(m,n)$ $a$, length-$\text{m}$ $b$, $cvec$ not referenced;
LS2	shape-$(m,n)$ $a$, length-$\text{m}$ $b$, length-$\text{n}$ $cvec$;
LS3	shape-$(m,n)$ $a$ upper trapezoidal, length-$\text{n}$ $kx$, length-$\text{m}$ $b$, $cvec$ not referenced;
LS4	shape-$(m,n)$ $a$ upper trapezoidal, length-$\text{n}$ $kx$, length-$\text{m}$ $b$, length-$\text{n}$ $cvec$.

For problems of type FP, the objective function is omitted and $a$, $b$ and $cvec$ are not referenced.

The following keywords are also acceptable.

$a$	Option
Least	LS1
Quadratic	QP2
Linear	LP

In addition, the keywords LS and LSQ are equivalent to the default option LS1, and the keyword QP is equivalent to the option QP2.

If $A=0$, i.e., the objective function is purely linear, the efficiency of lsq_lincon_solve may be increased by specifying $a$ as LP.

‘Rank Tolerance’float

Default $=100ϵ$ or $10\sqrt{ϵ}$ (see below)

Note that this option does not apply to problems of type FP or LP.

The default value of $r$ depends on the problem type. If $A$ occurs as a least squares matrix, as it does in problem types QP1, LS1 and LS3, then the default value of $r$ is $100ϵ$. In all other cases, $A$ is treated as the ‘square root’ of the Hessian matrix $H$ and $r$ has the default value $10\sqrt{ϵ}$.

This argument enables you to control the estimate of the triangular factor $R_1$ (see Main Iteration). If ${\rho }_i$ denotes the function ${\rho }_i=max\{\left|R_{11}\right|,\left|R_{22}\right|,\dots ,\left|R_{ii}\right|\}$, the rank of $R$ is defined to be smallest index i such that $\left|R_{i+1,i+1}\right|\le r\left|{\rho }_{i+1}\right|$. If $r\le 0$, the default value is used.

Raises

NagValueError

(errno $4$)

Too many iterations.

(errno $5$)

Too many iterations without changing $x$.

(errno $6$)

On entry, $n=⟨value⟩$.

Constraint: $n>0$.

(errno $6$)

On entry, $\text{nclin}=⟨value⟩$.

Constraint: $\text{nclin}\ge 0$.

(errno $6$)

On entry, $m=⟨value⟩$.

Constraint: $m>0$.

(errno $6$)

On entry, $kx$ has not been supplied as a valid permutation.

(errno $6$)

On entry, the equal bounds on $⟨value⟩$ are infinite, because $bl[⟨value⟩]=\text{beta}$ and $bu[⟨value⟩]=\text{beta}$, but $\left|\text{beta}\right|\ge \text{bigbnd}$: $\text{beta}=⟨value⟩$ and $\text{bigbnd}=⟨value⟩$.

(errno $6$)

On entry, the bounds on $⟨value⟩$ are inconsistent: $bl[⟨value⟩]=⟨value⟩$ and $bu[⟨value⟩]=⟨value⟩$.

(errno $6$)

On entry with a Warm Start, $istate[⟨value⟩]=⟨value⟩$.

(errno $6$)

On entry, the equal bounds on variable $⟨value⟩$ are infinite, because $bl[⟨value⟩]=\text{beta}$ and $bu[⟨value⟩]=\text{beta}$, but $\left|\text{beta}\right|\ge \text{bigbnd}$: $\text{beta}=⟨value⟩$ and $\text{bigbnd}=⟨value⟩$.

(errno $6$)

On entry, the equal bounds on linear constraint $⟨value⟩$ are infinite, because $bl[⟨value⟩]=\text{beta}$ and $bu[⟨value⟩]=\text{beta}$, but $\left|\text{beta}\right|\ge \text{bigbnd}$: $\text{beta}=⟨value⟩$ and $\text{bigbnd}=⟨value⟩$.

(errno $6$)

On entry, the equal bounds on nonlinear constraint $⟨value⟩$ are infinite, because $bl[⟨value⟩]=\text{beta}$ and $bu[⟨value⟩]=\text{beta}$, but $\left|\text{beta}\right|\ge \text{bigbnd}$: $\text{beta}=⟨value⟩$ and $\text{bigbnd}=⟨value⟩$.

(errno $6$)

On entry, the bounds on variable $⟨value⟩$ are inconsistent: $bl[⟨value⟩]=⟨value⟩$ and $bu[⟨value⟩]=⟨value⟩$.

(errno $6$)

On entry, the bounds on linear constraint $⟨value⟩$ are inconsistent: $bl[⟨value⟩]=⟨value⟩$ and $bu[⟨value⟩]=⟨value⟩$.

(errno $6$)

On entry, the bounds on nonlinear constraint $⟨value⟩$ are inconsistent: $bl[⟨value⟩]=⟨value⟩$ and $bu[⟨value⟩]=⟨value⟩$.

(errno $7$)

The problem to be solved is of type QP1 or QP2, but the Hessian matrix supplied in $a$ is not positive semidefinite.

Warns

NagAlgorithmicWarning

(errno $1$)

Weak $⟨value⟩$ solution.

(errno $2$)

$⟨value⟩$ solution is unbounded.

(errno $3$)

Cannot satisfy the linear constraints.

Notes

In the NAG Library the traditional C interface for this routine uses a different algorithmic base. Please contact NAG if you have any questions about compatibility.

lsq_lincon_solve is designed to solve a class of quadratic programming problems of the following general form:

$$\begin{array}{c}{minimize}_{x\in R^n}F\left(x\right)\text{subject to}l\le \left\{\begin{array}{c}x\\ Cx\end{array}\right\}\le u\end{array}$$

where $c$ is an $n_L\times n$ matrix and the objective function $F\left(x\right)$ may be specified in a variety of ways depending upon the particular problem to be solved. The available forms for $F\left(x\right)$ are listed in Table [label omitted], in which the prefixes FP, LP, QP and LS stand for ‘feasible point’, ‘linear programming’, ‘quadratic programming’ and ‘least squares’ respectively, $c$ is an $n$-element vector, $b$ is an $m$ element vector and $\parallel z\parallel$ denotes the Euclidean length of $z$.

Problem type	$F\left(x\right)$	Matrix $A$
FP	None	Not applicable
LP	$c^{T}x$	Not applicable
QP1	$\frac{1}{2}{x}^{T}Ax$	$n\times n$ symmetric positive semidefinite
QP2	$c^{T}x+\frac{1}{2}{x}^{T}Ax$	$n\times n$ symmetric positive semidefinite
QP3	$\frac{1}{2}{x}^T{A}^{T}Ax$	$m\times n$ upper trapezoidal
QP4	$c^{T}x+\frac{1}{2}{x}^T{A}^{T}Ax$	$m\times n$ upper trapezoidal
LS1	$\frac{1}{2}{\parallel b-Ax\parallel }^2$	$m\times n$
LS2	$c^{T}x+\frac{1}{2}{\parallel b-Ax\parallel }^2$	$m\times n$
LS3	$\frac{1}{2}{\parallel b-Ax\parallel }^2$	$m\times n$ upper trapezoidal
LS4	$c^{T}x+\frac{1}{2}{\parallel b-Ax\parallel }^2$	$m\times n$ upper trapezoidal

In the standard LS problem $F\left(x\right)$ will usually have the form LS1, and in the standard convex QP problem $F\left(x\right)$ will usually have the form QP2. The default problem type is LS1 and other objective functions are selected by using the option ‘Problem Type’.

When $A$ is upper trapezoidal it will usually be the case that $m=n$, so that $A$ is upper triangular, but full generality has been allowed for in the specification of the problem. The upper trapezoidal form is intended for cases where a previous factorization, such as a $QR$ factorization, has been performed.

The constraints involving $c$ are called the general constraints. Note that upper and lower bounds are specified for all the variables and for all the general constraints. An equality constraint can be specified by setting $l_i=u_i$. If certain bounds are not present, the associated elements of $l$ or $u$ can be set to special values that will be treated as $-\infty$ or $+\infty$. (See the description of the option ‘Infinite Bound Size’.)

The defining feature of a quadratic function $F\left(x\right)$ is that the second-derivative matrix $H$ (the Hessian matrix) is constant. For the LP case $H=0$; for QP1 and QP2, $H=A$; for QP3 and QP4, $H=A^{T}A$ and for LS1 (the default), LS2, LS3 and LS4, $H=A^{T}A$.

Problems of type QP3 and QP4 for which $A$ is not in upper trapezoidal form should be solved as types LS1 and LS2 respectively, with $b=0$.

For problems of type LS, we refer to $A$ as the least squares matrix, or the matrix of observations and to $b$ as the vector of observations.

You must supply an initial estimate of the solution.

If $H$ is nonsingular then lsq_lincon_solve will obtain the unique (global) minimum. If $H$ is singular then the solution may still be a global minimum if all active constraints have nonzero Lagrange multipliers. Otherwise the solution obtained will be either a weak minimum (i.e., with a unique optimal objective value, but an infinite set of optimal $x$), or else the objective function is unbounded below in the feasible region. The last case can only occur when $F\left(x\right)$ contains an explicit linear term (as in problems LP, QP2, QP4, LS2 and LS4).

The method used by lsq_lincon_solve is described in detail in Algorithmic Details.

References

Gill, P E, Hammarling, S, Murray, W, Saunders, M A and Wright, M H, 1986, Users’ guide for LSSOL (Version 1.0), Report SOL 86-1, Department of Operations Research, Stanford University

Gill, P E, Murray, W, Saunders, M A and Wright, M H, 1984, Procedures for optimization problems with a mixture of bounds and general linear constraints, ACM Trans. Math. Software (10), 282–298

Gill, P E, Murray, W and Wright, M H, 1981, Practical Optimization, Academic Press

Stoer, J, 1971, On the numerical solution of constrained least squares problems, SIAM J. Numer. Anal. (8), 382–411