naginterfaces.library.opt.nlp2_​sparse_​solve

naginterfaces.library.opt.nlp2_sparse_solve(start, objadd, objrow, prob, usrfun, iafun, javar, nea, a, igfun, jgvar, neg, xlow, xupp, xnames, flow, fupp, fnames, x, xstate, f, fstate, fmul, ns, comm, data=None, io_manager=None)

nlp2_sparse_solve solves sparse linear and nonlinear programming problems.

Note: this function uses optional algorithmic parameters, see also: nlp2_sparse_option_file(), nlp2_sparse_option_string(), nlp2_sparse_option_integer_set(), nlp2_sparse_option_double_set(), nlp2_sparse_init(), nlp2_sparse_jacobian().

Deprecated since version 28.3.0.0: nlp2_sparse_solve is deprecated. Please use handle_solve_ssqp() instead. See also the Replacement Calls document.

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

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

Parameters

startint

Indicates how a starting point is to be obtained.

$start=0$

Requests that the Crash procedure be used, unless a Basis file is provided via options ‘Old Basis File’, ‘Insert File’ or ‘Load File’.

$start=1$

Is the same as $start=0$ but is more meaningful when a Basis file is given.

$start=2$

Means that $xstate$ and $fstate$ define a valid starting point (probably from an earlier call, though not necessarily).

objaddfloat

Is a constant that will be added to the objective row $F_{obj}$ for printing purposes. Typically, $objadd=0.0e+0$.

objrowint

Says which row of $F\left(x\right)$ is to act as the objective function. If there is no such row, set $objrow=0$. Then nlp2_sparse_solve will seek a feasible point such that $l_F\le F\left(x\right)\le u_F$ and $l_x\le x\le u_x$.

probstr, length 8

Is an $8$-character name for the problem. $prob$ is used in the printed solution and in some functions that output Basis files. A blank name may be used.

usrfuncallable (status, f, g) = usrfun(status, x, needf, f, needg, g, data=None)

$usrfun$ must define the nonlinear portion $f\left(x\right)$ of the problem functions $F\left(x\right)=f\left(x\right)+Ax$, along with its gradient elements $G_{ij}\left(x\right)=\frac{\partial f_i\left(x\right)}{\partial x_j}$. (A dummy function is needed even if $f\equiv 0$ and all functions are linear.)

In general, $usrfun$ should return all function and gradient values on every entry except perhaps the last.

This provides maximum reliability and corresponds to the default option setting, $\text{‘Derivative Option'}=1$.

The elements of $G\left(x\right)$ are stored in the array $g(1:leng)$ in the order specified by the input arrays $igfun$ and $jgvar$.

In practice it is often convenient not to code gradients. nlp2_sparse_solve is able to estimate them by finite differences, using a call to $usrfun$ for each variable $x_j$ for which some $\frac{\partial f_i\left(x\right)}{\partial x_j}$ needs to be estimated.

However, this reduces the reliability of the optimization algorithm, and it can be very expensive if there are many such variables $x_j$.

As a compromise, nlp2_sparse_solve allows you to code as many gradients as you like.

This option is implemented as follows.

Just before $usrfun$ is called, each element of the derivative array $g$ is initialized to a specific value.

On exit, any element retaining that value must be estimated by finite differences.

Some rules of thumb follow:

1. for maximum reliability, compute all gradients;

2. if the gradients are expensive to compute, specify option ‘Nonderivative Linesearch’ and use the value of the input argument $needg$ to avoid computing them on certain entries. (There is no need to compute gradients if $needg=0$ on entry to $usrfun$.);

3. if not all gradients are known, you must specify $\text{‘Derivative Option'}=0$. You should still compute as many gradients as you can. (It often happens that some of them are constant or zero.);

4. again, if the known gradients are expensive, don’t compute them if $needg=0$ on entry to $usrfun$;

5. use the input argument $status$ to test for special actions on the first or last entries;

6. while $usrfun$ is being developed, use the option ‘Verify Level’ to check the computation of gradients that are supposedly known;

7. $usrfun$ is not called until the linear constraints and bounds on $x$ are satisfied. This helps confine $x$ to regions where the functions $f_i\left(x\right)$ are likely to be defined. However, be aware of the option ‘Minor Feasibility Tolerance’ if the functions have singularities on the constraint boundaries;

8. set $status=-1$ if some of the functions are undefined. The linesearch will shorten the step and try again;

9. set $status\le -2$ if you want nlp2_sparse_solve to stop.

Parameters

statusint

Indicates the first and last calls to $usrfun$.

$status=0$

There is nothing special about the current call to $usrfun$.

$status=1$

nlp2_sparse_solve is calling your function for the first time. You may wish to do something special such as read data from a file.

$status\ge 2$

nlp2_sparse_solve is calling your function for the last time. This argument setting allows you to perform some additional computation on the final solution.

$status=2$

The current $x$ is optimal.

$status=3$

The problem appears to be infeasible.

$status=4$

The problem appears to be unbounded.

$status=5$

An iterations limit was reached.

If the functions are expensive to evaluate, it may be desirable to do nothing on the last call.

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

The variables $x$ at which the problem functions are to be calculated. The array $x$ must not be altered.

needfint

Indicates whether $f$ must be assigned during this call of $usrfun$.

$needf=0$

$f$ is not required and is ignored.

$needf>0$

The components of $f\left(x\right)$ corresponding to the nonlinear part of $F\left(x\right)$ must be calculated and assigned to $f$.

If $F_i\left(x\right)$ is linear and completely defined by the $i$th row of $A$, $A_i^{\prime }$, the associated value $f_i\left(x\right)$ is ignored and need not be assigned.

However, if $F_i\left(x\right)$ has a nonlinear portion $f_i\left(x\right)$ that happens to be zero at $x$, it is still necessary to set $f_i\left(x\right)=0$.

If the linear part $A_i^{\prime }$ of a nonlinear $F_i\left(x\right)$ is provided using the arrays $iafun$, $javar$ and $a$, it must not be computed again as part of $f_i\left(x\right)$.

To simplify the code, you may ignore the value of $needf$ and compute $f\left(x\right)$ on every entry to $usrfun$.

$needf$ may also be ignored with ‘Derivative Linesearch’ and $\text{‘Derivative Option'}=1$.

In this case, $needf$ is always $1$, and $f$ must always be assigned.

ffloat, ndarray, shape $\left(\text{nf}\right)$

Concerns the calculation of $f\left(x\right)$.

needgint

Indicates whether $g$ must be assigned during this call of $usrfun$.

$needg=0$

$g$ is not required and is ignored.

$needg>0$

The partial derivatives of $f\left(x\right)$ must be calculated and assigned to $g$. The value of $g[k-1]$ should be $\frac{\partial f_i\left(x\right)}{\partial x_j}$, where $i=igfun[k-1]$, $j=jgvar[k-1]$ and $k=1,2,\dots ,\text{leng}$.

gfloat, ndarray, shape $\left(\text{leng}\right)$

Concerns the calculations of the derivatives of the function $f\left(x\right)$.

dataarbitrary, optional, modifiable in place

User-communication data for callback functions.

Returns

statusint

May be used to indicate that you are unable to evaluate $f$ or its gradients at the current $x$. (For example, the problem functions may not be defined there.)

During the linesearch, $f\left(x\right)$ is evaluated at points $x=x_k+\alpha p_k$ for various step lengths $\alpha$, where $f\left(x_k\right)$ has already been evaluated satisfactorily.

For any such $x$, if you set $status=-1$, nlp2_sparse_solve will reduce $\alpha$ and evaluate $f$ again (closer to $x_k$, where $f\left(x_k\right)$ is more likely to be defined).

If for some reason you wish to terminate the current problem, set $status\le -2$.

ffloat, array-like, shape $\left(\text{nf}\right)$

$f$ contains the computed functions $f\left(x\right)$ (except perhaps if $needf=0$).

gfloat, array-like, shape $\left(\text{leng}\right)$

Contains the computed derivatives $G\left(x\right)$ (unless $needg=0$).

These derivative elements must be stored in $g$ in exactly the same positions as implied by the definitions of arrays $igfun$ and $jgvar$.

There is no internal check for consistency (except indirectly via the option ‘Verify Level’), so great care is essential.

iafunint, array-like, shape $\left(\text{lena}\right)$

$iafun$ and $javar$ define the coordinates $(i,j)$ of the nonzero elements of the linear part $A$ of the function $F\left(x\right)=f\left(x\right)+Ax$.

In particular, $nea$ triples $\left(iafun\right[k-1],javar[k-1],a[k-1\left]\right)$ define the row and column indices $i=iafun[k-1]$ and $j=javar[k-1]$ of the element $A_{ij}=a[k-1]$.

The coordinates may define the elements of $A$ in any order.

javarint, array-like, shape $\left(\text{lena}\right)$

$iafun$ and $javar$ define the coordinates $(i,j)$ of the nonzero elements of the linear part $A$ of the function $F\left(x\right)=f\left(x\right)+Ax$.

In particular, $nea$ triples $\left(iafun\right[k-1],javar[k-1],a[k-1\left]\right)$ define the row and column indices $i=iafun[k-1]$ and $j=javar[k-1]$ of the element $A_{ij}=a[k-1]$.

The coordinates may define the elements of $A$ in any order.

neaint

Is the number of nonzero entries in $A$ such that $F\left(x\right)=f\left(x\right)+Ax$.

afloat, array-like, shape $\left(\text{lena}\right)$

Defines the values $A_{ij}$ of the nonzero elements of the linear part $A$ of the function $F\left(x\right)=f\left(x\right)+Ax$.

In particular, $nea$ triples $\left(iafun\right[k-1],javar[k-1],a[k-1\left]\right)$ define the row and column indices $i=iafun[k-1]$ and $j=javar[k-1]$ of the element $A_{ij}=a[k-1]$.

The coordinates may define the elements of $A$ in any order.

igfunint, array-like, shape $\left(\text{leng}\right)$

$igfun$ and $jgvar$ define the coordinates $(i,j)$ of the nonzero elements of $G$, the nonlinear part of the derivative $J\left(x\right)=G\left(x\right)+A$ of the function $F\left(x\right)=f\left(x\right)+Ax$. nlp2_sparse_jacobian() may be used to define these two arrays.

The coordinates can define the elements of $G$ in any order.

However, $usrfun$ must define the actual elements of $g$ in exactly the same order as defined by the coordinates $(igfun,jgvar)$.

jgvarint, array-like, shape $\left(\text{leng}\right)$

$igfun$ and $jgvar$ define the coordinates $(i,j)$ of the nonzero elements of $G$, the nonlinear part of the derivative $J\left(x\right)=G\left(x\right)+A$ of the function $F\left(x\right)=f\left(x\right)+Ax$. nlp2_sparse_jacobian() may be used to define these two arrays.

The coordinates can define the elements of $G$ in any order.

However, $usrfun$ must define the actual elements of $g$ in exactly the same order as defined by the coordinates $(igfun,jgvar)$.

negint

The number of nonzero entries in $G$.

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

Contain the lower and upper bounds $l_x$ and $u_x$ on the variables $x$.

To specify a nonexistent lower bound ${\left[l_x\right]}_j=-\infty$, set $xlow[j-1]\le -\text{bigbnd}$, where $\text{bigbnd}$ is the option ‘Infinite Bound Size’.

To specify a nonexistent upper bound ${\left[u_x\right]}_j=\infty$, set $xupp[j-1]\ge \text{bigbnd}$.

To fix the $j$th variable at $x_j=\beta$, where $\left|\beta \right|<\text{bigbnd}$, set $xlow[j-1]=xupp[j-1]=\beta$.

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

Contain the lower and upper bounds $l_x$ and $u_x$ on the variables $x$.

To specify a nonexistent lower bound ${\left[l_x\right]}_j=-\infty$, set $xlow[j-1]\le -\text{bigbnd}$, where $\text{bigbnd}$ is the option ‘Infinite Bound Size’.

To specify a nonexistent upper bound ${\left[u_x\right]}_j=\infty$, set $xupp[j-1]\ge \text{bigbnd}$.

To fix the $j$th variable at $x_j=\beta$, where $\left|\beta \right|<\text{bigbnd}$, set $xlow[j-1]=xupp[j-1]=\beta$.

xnamesstr, length 8, array-like, shape $\left(\text{nxname}\right)$

The optional names for the variables.

If $\text{nxname}=1$, $xnames$ is not referenced and default names will be used for output.

If $\text{nxname}=n$, $xnames[j-1]$ should contain the $8$-character name of the $j$th variable.

flowfloat, array-like, shape $\left(\text{nf}\right)$

Contain the lower and upper bounds $l_F$ and $u_F$ on $F\left(x\right)$.

To specify a nonexistent lower bound ${\left[l_F\right]}_i=-\infty$, set $flow[i-1]\le -\text{bigbnd}$.

For a nonexistent upper bound ${\left[u_F\right]}_i=\infty$, set $fupp[i-1]\ge \text{bigbnd}$.

To make the $i$th constraint an equality at $F_i=\beta$, where $\left|\beta \right|<\text{bigbnd}$, set $flow[i-1]=fupp[i-1]=\beta$.

fuppfloat, array-like, shape $\left(\text{nf}\right)$

Contain the lower and upper bounds $l_F$ and $u_F$ on $F\left(x\right)$.

To specify a nonexistent lower bound ${\left[l_F\right]}_i=-\infty$, set $flow[i-1]\le -\text{bigbnd}$.

For a nonexistent upper bound ${\left[u_F\right]}_i=\infty$, set $fupp[i-1]\ge \text{bigbnd}$.

To make the $i$th constraint an equality at $F_i=\beta$, where $\left|\beta \right|<\text{bigbnd}$, set $flow[i-1]=fupp[i-1]=\beta$.

fnamesstr, length 8, array-like, shape $\left(\text{nfname}\right)$

The optional names for the problem functions.

If $\text{nfname}=1$, $fnames$ is not referenced and default names will be used for output.

If $\text{nfname}=\text{nf}$, $fnames[i-1]$ should contain the $8$-character name of the $i$th row of $F$.

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

An initial estimate of the variables $x$. See the following description of $xstate$.

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

The initial state for each variable $x$.

If $start=0$ or $1$ and no basis information is provided (the options ‘Old Basis File’, ‘Insert File’ and ‘Load File’ are all set to $0$; the default) $x$ and $xstate$ must be defined.

If nothing special is known about the problem, or if there is no wish to provide special information, you may set $x[j-1]=0.0$, $xstate[j-1]=0$, for all $j=1,2,\dots ,n$.

If you set $x[j-1]=xlow[j-1]$ set $xstate[j-1]=4$; if you set $x[j-1]=xupp[j-1]$, then set $xstate[j-1]=5$.

In this case a Crash procedure is used to select an initial basis.

If $start=0$ or $1$ and basis information is provided (at least one of the options ‘Old Basis File’, ‘Insert File’ and ‘Load File’ is nonzero) $x$ and $xstate$ need not be set.

If $start=2$ (Warm Start), $x$ and $xstate$ must be set (probably from a previous call).

In this case $xstate[\text{j}-1]$ must be $0$, $1$, $2$ or $3$, for $\text{j}=1,2,\dots ,n$.

ffloat, array-like, shape $\left(\text{nf}\right)$

An initial value for the problem functions $F$. See the following description of $fstate$.

fstateint, array-like, shape $\left(\text{nf}\right)$

The initial state for the problem functions $F$.

If $start=0$ or $1$ and no basis information is provided (the options ‘Old Basis File’, ‘Insert File’ and ‘Load File’ are all set to $0$; the default, $f$ and $fstate$ must be defined.

If nothing special is known about the problem, or if there is no wish to provide special information, you may set $f[i-1]=0.0$, $fstate[i-1]=0$, for all $i=1,2,\dots ,\text{nf}$.

Less trivially, to say that the optimal value of function $f[i-1]$ will probably be equal to one of its bounds, set $f[i-1]=flow[i-1]$ and $fstate[i-1]=4$ or $f[i-1]=fupp[i-1]$ and $fstate[i-1]=5$ as appropriate.

In this case a Crash procedure is used to select an initial basis.

If $start=0$ or $1$ and basis information is provided (at least one of the options ‘Old Basis File’, ‘Insert File’ and ‘Load File’ is nonzero), $f$ and $fstate$ need not be set.

If $start=2$ (Warm Start), $f$ and $fstate$ must be set (probably from a previous call).

In this case $fstate[\text{i}-1]$ must be $0$, $1$, $2$ or $3$, for $\text{i}=1,2,\dots ,\text{nf}$.

fmulfloat, array-like, shape $\left(\text{nf}\right)$

An estimate of $\gamma$, the vector of Lagrange multipliers (shadow prices) for the constraints $l_F\le F\left(x\right)\le u_F$. All $\text{nf}$ components must be defined. If nothing is known about $\gamma$, set $fmul[\text{i}-1]=0.0$, for $\text{i}=1,2,\dots ,\text{nf}$. For warm start use the values from a previous call.

nsint

The number of superbasic variables. $ns$ need not be specified for cold starts, but should retain its value from a previous call when warm start is used.

commdict, communication object, modified in place

Communication structure.

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

dataarbitrary, optional

User-communication data for callback functions.

io_managerFileObjManager, optional

Manager for I/O in this routine.

Returns

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

The final values of the variable $x$.

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

The final state of the variables.

$xstate[j-1]$	State of variable $j$	Usual value of $x[j-1]$
0	nonbasic	$xlow[j-1]$
1	nonbasic	$xupp[j-1]$
2	superbasic	Between $xlow[j-1]$ and $xupp[j-1]$
3	basic	Between $xlow[j-1]$ and $xupp[j-1]$

Basic and superbasic variables may be outside their bounds by as much as the option ‘Minor Feasibility Tolerance’.

Note that if scaling is specified, the feasibility tolerance applies to the variables of the scaled problem.

In this case, the variables of the original problem may be as much as $0.1$ outside their bounds, but this is unlikely unless the problem is very badly scaled.

Check the value of Primal infeasibility output to the unit number (see unit_from_fileobj()) associated with the option ‘Print File’.

Very occasionally some nonbasic variables may be outside their bounds by as much as the option ‘Minor Feasibility Tolerance’, and there may be some nonbasics for which $x[j-1]$ lies strictly between its bounds.

If $ninf>0$, some basic and superbasic variables may be outside their bounds by an arbitrary amount (bounded by $sinf$ if scaling was not used).

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

The vector of the dual variables (Lagrange multipliers) for the simple bounds $l_x\le x\le u_x$.

ffloat, ndarray, shape $\left(\text{nf}\right)$

The final values for the problem functions $F$ (the values $F$ at the final point $x$).

fstateint, ndarray, shape $\left(\text{nf}\right)$

The final state of the variables. The elements of $fstate$ have the following meaning:

$fstate[i-1]$	State of the corresponding slack variable	Usual value of $f[i-1]$
0	nonbasic	$flow[i-1]$
1	nonbasic	$fupp[i-1]$
2	superbasic	Between $flow[i-1]$ and $fupp[i-1]$
3	basic	Between $flow[i-1]$ and $fupp[i-1]$

Basic and superbasic slack variables may lead to the corresponding functions being outside their bounds by as much as the option ‘Minor Feasibility Tolerance’.

Very occasionally some functions may be outside their bounds by as much as the option ‘Minor Feasibility Tolerance’, and there may be some nonbasics for which $f[i-1]$ lies strictly between its bounds.

If $ninf>0$, some basic and superbasic variables may be outside their bounds by an arbitrary amount (bounded by $sinf$ if scaling was not used).

fmulfloat, ndarray, shape $\left(\text{nf}\right)$

The vector of the dual variables (Lagrange multipliers) for the general constraints $l_F\le F\left(x\right)\le u_F$

nsint

The final number of superbasic variables.

ninfint

Are the number and the sum of the infeasibilities of constraints that lie outside one of their bounds by more than the option ‘Minor Feasibility Tolerance’ before the solution is unscaled.

If any linear constraints are infeasible, $x$ minimizes the sum of the infeasibilities of the linear constraints subject to the upper and lower bounds being satisfied.

In this case $ninf$ gives the number of variables and linear constraints lying outside their upper or lower bounds.

The nonlinear constraints are not evaluated.

Otherwise, $x$ minimizes the sum of infeasibilities of the nonlinear constraints subject to the linear constraints and upper and lower bounds being satisfied.

In this case $ninf$ gives the number of components of $F\left(x\right)$ lying outside their bounds by more than the option ‘Minor Feasibility Tolerance’.

Again this is before the solution is unscaled.

sinffloat

Are the number and the sum of the infeasibilities of constraints that lie outside one of their bounds by more than the option ‘Minor Feasibility Tolerance’ before the solution is unscaled.

If any linear constraints are infeasible, $x$ minimizes the sum of the infeasibilities of the linear constraints subject to the upper and lower bounds being satisfied.

In this case $ninf$ gives the number of variables and linear constraints lying outside their upper or lower bounds.

The nonlinear constraints are not evaluated.

Otherwise, $x$ minimizes the sum of infeasibilities of the nonlinear constraints subject to the linear constraints and upper and lower bounds being satisfied.

In this case $ninf$ gives the number of components of $F\left(x\right)$ lying outside their bounds by more than the option ‘Minor Feasibility Tolerance’.

Again this is before the solution is unscaled.

Other Parameters

‘Central Difference Interval’float

Default $={ϵ}_r^{\frac{1}{3}}$

When $\text{‘Derivative Option'}=0$, the central-difference interval $r$ is used near an optimal solution to obtain more accurate (but more expensive) estimates of gradients. Twice as many function evaluations are required compared to forward differencing. The interval used for the $j$th variable is $h_j=r(1+\left|x_j\right|)$. The resulting derivative estimates should be accurate to $O\left(r^2\right)$, unless the functions are badly scaled.

If you supply a value for this option, a small value between $0.0$ and $1.0$ is appropriate.

‘Check Frequency’int

Default $=60$

Every $i$th minor iteration after the most recent basis factorization, a numerical test is made to see if the current solution $x$ satisfies the general linear constraints (the linear constraints and the linearized nonlinear constraints, if any). The constraints are of the form $Ax-s=b$, where $s$ is the set of slack variables. To perform the numerical test, the residual vector $r=b-Ax+s$ is computed. If the largest component of $r$ is judged to be too large, the current basis is refactorized and the basic variables are recomputed to satisfy the general constraints more accurately. If $i\le 0$, the value of $i=99999999$ is used and effectively no checks are made.

$\text{‘Check Frequency'}=1$ is useful for debugging purposes, but otherwise this option should not be needed.

‘Crash Option’int

Default $=3$

Except on restarts, an internal Crash procedure is used to select an initial basis from certain rows and columns of the constraint matrix $\left(\begin{array}{cc}A& -I\end{array}\right)$. The ‘Crash Option’ $i$ determines which rows and columns of $A$ are eligible initially, and how many times the Crash procedure is called. Columns of $-I$ are used to pad the basis where necessary.

$i$	Meaning
$0$	The initial basis contains only slack variables: $B=I$.
$1$	The Crash procedure is called once, looking for a triangular basis in all rows and columns of $A$.
$2$	The Crash procedure is called twice (if there are nonlinear constraints). The first call looks for a triangular basis in linear rows, and the iteration proceeds with simplex iterations until the linear constraints are satisfied. The Jacobian is then evaluated for the first major iteration and the Crash procedure is called again to find a triangular basis in the nonlinear rows (retaining the current basis for linear rows).
$3$	The Crash procedure is called up to three times (if there are nonlinear constraints). The first two calls treat linear equalities and linear inequalities separately. As before, the last call treats nonlinear rows before the first major iteration.

If $i\ge 1$, certain slacks on inequality rows are selected for the basis first. (If $i\ge 2$, numerical values are used to exclude slacks that are close to a bound). The Crash procedure then makes several passes through the columns of $A$, searching for a basis matrix that is essentially triangular. A column is assigned to ‘pivot’ on a particular row if the column contains a suitably large element in a row that has not yet been assigned. (The pivot elements ultimately form the diagonals of the triangular basis.) For remaining unassigned rows, slack variables are inserted to complete the basis.

The ‘Crash Tolerance’ $r$ allows the starting Crash procedure to ignore certain ‘small’ nonzeros in each column of $A$. If $a_{max}$ is the largest element in column $j$, other nonzeros of $a_{ij}$ in the columns are ignored if $\left|a_{ij}\right|\le a_{max}\times r$. (To be meaningful, $r$ should be in the range $0\le r<1$.)

When $r>0.0$, the basis obtained by the Crash procedure may not be strictly triangular, but it is likely to be nonsingular and almost triangular. The intention is to obtain a starting basis containing more columns of $A$ and fewer (arbitrary) slacks. A feasible solution may be reached sooner on some problems.

For example, suppose the first $m$ columns of $A$ form the matrix shown under ‘LU Factor Tolerance’; i.e., a tridiagonal matrix with entries $-1$, $4$, $-1$. To help the Crash procedure choose all $m$ columns for the initial basis, we would specify a ‘Crash Tolerance’ of $r$ for some value of $r>0.5$.

‘Crash Tolerance’float

Default $=0.1$

Except on restarts, an internal Crash procedure is used to select an initial basis from certain rows and columns of the constraint matrix $\left(\begin{array}{cc}A& -I\end{array}\right)$. The ‘Crash Option’ $i$ determines which rows and columns of $A$ are eligible initially, and how many times the Crash procedure is called. Columns of $-I$ are used to pad the basis where necessary.

$i$	Meaning
$0$	The initial basis contains only slack variables: $B=I$.
$1$	The Crash procedure is called once, looking for a triangular basis in all rows and columns of $A$.
$2$	The Crash procedure is called twice (if there are nonlinear constraints). The first call looks for a triangular basis in linear rows, and the iteration proceeds with simplex iterations until the linear constraints are satisfied. The Jacobian is then evaluated for the first major iteration and the Crash procedure is called again to find a triangular basis in the nonlinear rows (retaining the current basis for linear rows).
$3$	The Crash procedure is called up to three times (if there are nonlinear constraints). The first two calls treat linear equalities and linear inequalities separately. As before, the last call treats nonlinear rows before the first major iteration.

If $i\ge 1$, certain slacks on inequality rows are selected for the basis first. (If $i\ge 2$, numerical values are used to exclude slacks that are close to a bound). The Crash procedure then makes several passes through the columns of $A$, searching for a basis matrix that is essentially triangular. A column is assigned to ‘pivot’ on a particular row if the column contains a suitably large element in a row that has not yet been assigned. (The pivot elements ultimately form the diagonals of the triangular basis.) For remaining unassigned rows, slack variables are inserted to complete the basis.

The ‘Crash Tolerance’ $r$ allows the starting Crash procedure to ignore certain ‘small’ nonzeros in each column of $A$. If $a_{max}$ is the largest element in column $j$, other nonzeros of $a_{ij}$ in the columns are ignored if $\left|a_{ij}\right|\le a_{max}\times r$. (To be meaningful, $r$ should be in the range $0\le r<1$.)

When $r>0.0$, the basis obtained by the Crash procedure may not be strictly triangular, but it is likely to be nonsingular and almost triangular. The intention is to obtain a starting basis containing more columns of $A$ and fewer (arbitrary) slacks. A feasible solution may be reached sooner on some problems.

For example, suppose the first $m$ columns of $A$ form the matrix shown under ‘LU Factor Tolerance’; i.e., a tridiagonal matrix with entries $-1$, $4$, $-1$. To help the Crash procedure choose all $m$ columns for the initial basis, we would specify a ‘Crash Tolerance’ of $r$ for some value of $r>0.5$.

‘Defaults’valueless

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

‘Derivative Option’int

Default $=1$

Option ‘Derivative Option’ specifies which nonlinear function gradients are known analytically and will be supplied to nlp2_sparse_solve by $usrfun$.

$i$	Meaning
$0$	Some problem derivatives are unknown.
$1$	All problem derivatives are known.

The value $i=1$ should be used whenever possible. It is the most reliable and will usually be the most efficient.

If $i=0$, nlp2_sparse_solve will estimate the missing components of $G\left(x\right)$ using finite differences. This may simplify the coding of $usrfun$. However, it could increase the total run-time substantially (since a special call to $usrfun$ is required for each column of the Jacobian that has a missing element), and there is less assurance that an acceptable solution will be located. If the nonlinear variables are not well scaled, it may be necessary to specify a nonstandard option ‘Difference Interval’.

For each column of the Jacobian, one call to $usrfun$ is needed to estimate all missing elements in that column, if any.

At times, central differences are used rather than forward differences. Twice as many calls to $usrfun$ are needed. (This is not under your control.)

‘Derivative Linesearch’valueless

Default

At each major iteration a linesearch is used to improve the merit function. Option ‘Derivative Linesearch’ uses safeguarded cubic interpolation and requires both function and gradient values to compute estimates of the step ${\alpha }_k$. If some analytic derivatives are not provided, or option ‘Nonderivative Linesearch’ is specified, nlp2_sparse_solve employs a linesearch based upon safeguarded quadratic interpolation, which does not require gradient evaluations.

A nonderivative linesearch can be slightly less robust on difficult problems, and it is recommended that the default be used if the functions and derivatives can be computed at approximately the same cost. If the gradients are very expensive relative to the functions, a nonderivative linesearch may give a significant decrease in computation time.

If ‘Nonderivative Linesearch’ is selected, nlp2_sparse_solve signals the evaluation of the linesearch by calling $usrfun$ with $needg=0$. Once the linesearch is completed, the problem functions are called again with $needf=0$ and $needg=0$. If the potential saving provided by a nonderivative linesearch is to be realised, it is essential that $usrfun$ be coded so that derivatives are not computed when $needg=0$.

‘Nonderivative Linesearch’valueless

At each major iteration a linesearch is used to improve the merit function. Option ‘Derivative Linesearch’ uses safeguarded cubic interpolation and requires both function and gradient values to compute estimates of the step ${\alpha }_k$. If some analytic derivatives are not provided, or option ‘Nonderivative Linesearch’ is specified, nlp2_sparse_solve employs a linesearch based upon safeguarded quadratic interpolation, which does not require gradient evaluations.

A nonderivative linesearch can be slightly less robust on difficult problems, and it is recommended that the default be used if the functions and derivatives can be computed at approximately the same cost. If the gradients are very expensive relative to the functions, a nonderivative linesearch may give a significant decrease in computation time.

If ‘Nonderivative Linesearch’ is selected, nlp2_sparse_solve signals the evaluation of the linesearch by calling $usrfun$ with $needg=0$. Once the linesearch is completed, the problem functions are called again with $needf=0$ and $needg=0$. If the potential saving provided by a nonderivative linesearch is to be realised, it is essential that $usrfun$ be coded so that derivatives are not computed when $needg=0$.

‘Difference Interval’float

Default $=\sqrt{{ϵ}_r}$

This alters the interval $r$ used to estimate gradients by forward differences. It does so in the following circumstances:

• in the interval (‘cheap’) phase of verifying the problem derivatives;

• for verifying the problem derivatives;

• for estimating missing derivatives.

In all cases, a derivative with respect to $x_j$ is estimated by perturbing that component of $x$ to the value $x_j+r(1+\left|x_j\right|)$, and then evaluating $F_{obj}\left(x\right)$ or $f\left(x\right)$ at the perturbed point. The resulting gradient estimates should be accurate to $O\left(r\right)$ unless the functions are badly scaled. Judicious alteration of $r$ may sometimes lead to greater accuracy.

If you supply a value for this option, a small value between $0.0$ and $1.0$ is appropriate.

‘Dump File’int

Default $=0$

Options ‘Dump File’ and ‘Load File’ are similar to options ‘Punch File’ and ‘Insert File’, but they record solution information in a manner that is more direct and more easily modified. A full description of information recorded in options ‘Dump File’ and ‘Load File’ is given in Gill et al. (2005a).

If $i_1>0$, the last solution obtained will be output to the file with unit number $i_1$.

If $i_2>0$, the ‘Load File’, containing basis information, will be read. The file will usually have been output previously as a ‘Dump File’. The file will not be accessed if options ‘Old Basis File’ or ‘Insert File’ are specified.

‘Load File’int

Default $=0$

Options ‘Dump File’ and ‘Load File’ are similar to options ‘Punch File’ and ‘Insert File’, but they record solution information in a manner that is more direct and more easily modified. A full description of information recorded in options ‘Dump File’ and ‘Load File’ is given in Gill et al. (2005a).

If $i_1>0$, the last solution obtained will be output to the file with unit number $i_1$.

If $i_2>0$, the ‘Load File’, containing basis information, will be read. The file will usually have been output previously as a ‘Dump File’. The file will not be accessed if options ‘Old Basis File’ or ‘Insert File’ are specified.

‘Elastic Weight’float

Default $={10}^4$

This keyword determines the initial weight $\gamma$ associated with the problem (1) (see Treatment of Constraint Infeasibilities).

At major iteration $k$, if elastic mode has not yet started, a scale factor ${\sigma }_k=1+{\parallel g\left(x_k\right)\parallel }_{\infty }$ is defined from the current objective gradient. Elastic mode is then started if the QP subproblem is infeasible, or the QP dual variables are larger in magnitude than ${\sigma }_{k}r$. The QP is resolved in elastic mode with $\gamma ={\sigma }_{k}r$.

Thereafter, major iterations continue in elastic mode until they converge to a point that is optimal for (1) (see Treatment of Constraint Infeasibilities). If the point is feasible for equation (1) $(v=w=0)$, it is declared locally optimal. Otherwise, $\gamma$ is increased by a factor of $10$ and major iterations continue. If $\gamma$ has already reached a maximum allowable value, equation (1) is declared locally infeasible.

‘Expand Frequency’int

Default $=10000$

This option is part of the anti-cycling procedure designed to make progress even on highly degenerate problems.

For linear models, the strategy is to force a positive step at every iteration, at the expense of violating the bounds on the variables by a small amount. Suppose that the option ‘Minor Feasibility Tolerance’ is $\delta$. Over a period of $i$ iterations, the tolerance actually used by nlp2_sparse_solve increases from $0.5\delta$ to $\delta$ (in steps of $0.5\delta /i$).

For nonlinear models, the same procedure is used for iterations in which there is only one superbasic variable. (Cycling can occur only when the current solution is at a vertex of the feasible region.) Thus, zero steps are allowed if there is more than one superbasic variable, but otherwise positive steps are enforced.

Increasing $i$ helps reduce the number of slightly infeasible nonbasic variables (most of which are eliminated during a resetting procedure). However, it also diminishes the freedom to choose a large pivot element (see option ‘Pivot Tolerance’).

‘Factorization Frequency’int

Default $=50$

At most $i$ basis changes will occur between factorizations of the basis matrix.

With linear programs, the basis factors are usually updated every iteration. The default $i$ is reasonable for typical problems. Higher values up to $i=100$ (say) may be more efficient on well-scaled problems.

When the objective function is nonlinear, fewer basis updates will occur as an optimum is approached. The number of iterations between basis factorizations will, therefore, increase. During these iterations a test is made regularly (according to the option ‘Check Frequency’) to ensure that the general constraints are satisfied. If necessary the basis will be refactorized before the limit of $i$ updates is reached.

‘Function Precision’float

Default $={ϵ}^{0.8}$

The relative function precision ${ϵ}_r$ is intended to be a measure of the relative accuracy with which the nonlinear functions can be computed. For example, if $f\left(x\right)$ is computed as $1000.56789$ for some relevant $x$ and if the first $6$ significant digits are known to be correct, the appropriate value for ${ϵ}_r$ would be $1.0e-6$.

Ideally the functions $f_i\left(x\right)$ should have magnitude of order $1$. If all functions are substantially less than $1$ in magnitude, ${ϵ}_r$ should be the absolute precision. For example, if $f\left(x\right)=1.23456789e-4$ at some point and if the first $6$ significant digits are known to be correct, the appropriate value for ${ϵ}_r$ would be $1.0e-10$.)

The default value of ${ϵ}_r$ is appropriate for simple analytic functions.

In some cases the function values will be the result of extensive computation, possibly involving a costly iterative procedure that can provide few digits of precision. Specifying an appropriate ‘Function Precision’ may lead to savings, by allowing the linesearch procedure to terminate when the difference between function values along the search direction becomes as small as the absolute error in the values.

‘Hessian Full Memory’valueless

Default if $n_1\le 75$

These options select the method for storing and updating the approximate Hessian. (nlp2_sparse_solve uses a quasi-Newton approximation to the Hessian of the Lagrangian. A BFGS update is applied after each major iteration.)

If ‘Hessian Full Memory’ is specified, the approximate Hessian is treated as a dense matrix and the BFGS updates are applied explicitly. This option is most efficient when the number of variables $n$ is not too large (say, less than $75$). In this case, the storage requirement is fixed and one can expect $1$-step Q-superlinear convergence to the solution.

‘Hessian Limited Memory’ should be used on problems where $n$ is very large. In this case a limited-memory procedure is used to update a diagonal Hessian approximation $H_r$ a limited number of times. (Updates are accumulated as a list of vector pairs. They are discarded at regular intervals after $H_r$ has been reset to their diagonal.)

‘Hessian Limited Memory’valueless

Default if $n_1>75$

These options select the method for storing and updating the approximate Hessian. (nlp2_sparse_solve uses a quasi-Newton approximation to the Hessian of the Lagrangian. A BFGS update is applied after each major iteration.)

If ‘Hessian Full Memory’ is specified, the approximate Hessian is treated as a dense matrix and the BFGS updates are applied explicitly. This option is most efficient when the number of variables $n$ is not too large (say, less than $75$). In this case, the storage requirement is fixed and one can expect $1$-step Q-superlinear convergence to the solution.

‘Hessian Limited Memory’ should be used on problems where $n$ is very large. In this case a limited-memory procedure is used to update a diagonal Hessian approximation $H_r$ a limited number of times. (Updates are accumulated as a list of vector pairs. They are discarded at regular intervals after $H_r$ has been reset to their diagonal.)

‘Hessian Frequency’int

Default $=99999999$

If option ‘Hessian Full Memory’ is in effect and $i$ BFGS updates have already been carried out, the Hessian approximation is reset to the identity matrix. (For certain problems, occasional resets may improve convergence, but in general they should not be necessary.)

‘Hessian Full Memory’ and $\text{‘Hessian Frequency'}=10$ have a similar effect to ‘Hessian Limited Memory’ and $\text{‘Hessian Updates'}=10$ (except that the latter retains the current diagonal during resets).

‘Hessian Updates’int

Default $=\text{‘Hessian Frequency'}$ if ‘Hessian Full Memory’, $10$ otherwise

If option ‘Hessian Limited Memory’ is in effect and $i$ BFGS updates have already been carried out, all but the diagonal elements of the accumulated updates are discarded and the updating process starts again.

Broadly speaking, the more updates stored, the better the quality of the approximate Hessian. However, the more vectors stored, the greater the cost of each QP iteration. The default value is likely to give a robust algorithm without significant expense, but faster convergence can sometimes be obtained with significantly fewer updates (e.g., $i=5$).

‘Infinite Bound Size’float

Default $={10}^{20}$

If $r\ge 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.

‘Iterations Limit’int

Default $=max(10000,10max(n,\text{nf}))$

The value of $i$ specifies the maximum number of minor iterations allowed (i.e., iterations of the simplex method or the QP algorithm), summed over all major iterations. (See also the description of the option ‘Minor Iterations Limit’.)

‘Linesearch Tolerance’float

Default $=0.9$

This tolerance, $r$, controls the accuracy with which a step length will be located along the direction of search each iteration. At the start of each linesearch a target directional derivative for the merit function is identified. This argument determines the accuracy to which this target value is approximated, and it must be a value in the range $0.0\le r\le 1.0$.

The default value $r=0.9$ requests just moderate accuracy in the linesearch.

If the nonlinear functions are cheap to evaluate, a more accurate search may be appropriate; try $r=0.1\text{,}0.01\text{or}0.001$.

If the nonlinear functions are expensive to evaluate, a less accurate search may be appropriate. If all gradients are known, try $r=0.99$. (The number of major iterations might increase, but the total number of function evaluations may decrease enough to compensate.)

If not all gradients are known, a moderately accurate search remains appropriate. Each search will require only $1$–5 function values (typically), but many function calls will then be needed to estimate missing gradients for the next iteration.

‘LU Density Tolerance’float

Default $=0.6$

The density tolerance, $r_1$, is used during $LU$ factorization of the basis matrix $B$. Columns of $L$ and rows of $U$ are formed one at a time, and the remaining rows and columns of the basis are altered appropriately. At any stage, if the density of the remaining matrix exceeds $r_1$, the Markowitz strategy for choosing pivots is terminated, and the remaining matrix is factored by a dense $LU$ procedure. Raising the density tolerance towards $1.0$ may give slightly sparser $LU$ factors, with a slight increase in factorization time.

The singularity tolerance, $r_2$, helps guard against ill-conditioned basis matrices. After $B$ is refactorized, the diagonal elements of $U$ are tested as follows: if $\left|u_{jj}\right|\le r_2$ or $\left|u_{jj}\right|<r_2{max}_i\left|u_{ij}\right|$, the $j$th column of the basis is replaced by the corresponding slack variable. (This is most likely to occur after a restart.)

‘LU Singularity Tolerance’float

Default $={ϵ}^{\frac{2}{3}}$

The density tolerance, $r_1$, is used during $LU$ factorization of the basis matrix $B$. Columns of $L$ and rows of $U$ are formed one at a time, and the remaining rows and columns of the basis are altered appropriately. At any stage, if the density of the remaining matrix exceeds $r_1$, the Markowitz strategy for choosing pivots is terminated, and the remaining matrix is factored by a dense $LU$ procedure. Raising the density tolerance towards $1.0$ may give slightly sparser $LU$ factors, with a slight increase in factorization time.

The singularity tolerance, $r_2$, helps guard against ill-conditioned basis matrices. After $B$ is refactorized, the diagonal elements of $U$ are tested as follows: if $\left|u_{jj}\right|\le r_2$ or $\left|u_{jj}\right|<r_2{max}_i\left|u_{ij}\right|$, the $j$th column of the basis is replaced by the corresponding slack variable. (This is most likely to occur after a restart.)

‘LU Factor Tolerance’float

Default $=3.99$

The values of $r_1$ and $r_2$ affect the stability of the basis factorization $B=LU$, during refactorization and updates respectively. The lower triangular matrix $L$ is a product of matrices of the form

$$\begin{array}{c}\left(\begin{array}{cc}1& \\ \mu & 1\end{array}\right)\end{array}$$

where the multipliers $\mu$ will satisfy $\left|\mu \right|\le r_i$. The default values of $r_1$ and $r_2$ usually strike a good compromise between stability and sparsity. They must satisfy $r_1$, $r_2\ge 1.0$.

For large and relatively dense problems, $r_1=10.0$ or $5.0$ (say) may give a useful improvement in stability without impairing sparsity to a serious degree.

For certain very regular structures (e.g., band matrices) it may be necessary to reduce $r_1\text{and/or}{r}_2$ in order to achieve stability. For example, if the columns of $A$ include a sub-matrix of the form

$$\begin{array}{c}\left(\begin{array}{cccccc}4& -1& & & & \\ -1& 4& -1& & & \\ & -1& 4& -1& & \\ & & \dots & \dots & \dots & \\ & & & -1& 4& -1\\ & & & & -1& 4\end{array}\right)\text{,}\end{array}$$

one should set both $r_1$ and $r_2$ to values in the range $1.0\le r_i<4.0$.

‘LU Update Tolerance’float

Default $=3.99$

The values of $r_1$ and $r_2$ affect the stability of the basis factorization $B=LU$, during refactorization and updates respectively. The lower triangular matrix $L$ is a product of matrices of the form

$$\begin{array}{c}\left(\begin{array}{cc}1& \\ \mu & 1\end{array}\right)\end{array}$$

where the multipliers $\mu$ will satisfy $\left|\mu \right|\le r_i$. The default values of $r_1$ and $r_2$ usually strike a good compromise between stability and sparsity. They must satisfy $r_1$, $r_2\ge 1.0$.

For large and relatively dense problems, $r_1=10.0$ or $5.0$ (say) may give a useful improvement in stability without impairing sparsity to a serious degree.

For certain very regular structures (e.g., band matrices) it may be necessary to reduce $r_1\text{and/or}{r}_2$ in order to achieve stability. For example, if the columns of $A$ include a sub-matrix of the form

$$\begin{array}{c}\left(\begin{array}{cccccc}4& -1& & & & \\ -1& 4& -1& & & \\ & -1& 4& -1& & \\ & & \dots & \dots & \dots & \\ & & & -1& 4& -1\\ & & & & -1& 4\end{array}\right)\text{,}\end{array}$$

one should set both $r_1$ and $r_2$ to values in the range $1.0\le r_i<4.0$.

‘LU Partial Pivoting’valueless

Default

The $LU$ factorization implements a Markowitz-type search for pivots that locally minimize the fill-in subject to a threshold pivoting stability criterion. The default option is to use threshhold partial pivoting. The options ‘LU Rook Pivoting’ and ‘LU Complete Pivoting’ are more expensive than partial pivoting but are more stable and better at revealing rank, as long as ‘LU Factor Tolerance’ is not too large (say $<2.0$). When numerical difficulties are encountered, nlp2_sparse_solve automatically reduces the $LU$ tolerance towards $1.0$ and switches (if necessary) to rook or complete pivoting, before reverting to the default or specified options at the next refactorization (with ‘System Information Yes’, relevant messages are output to the ‘Print File’).

‘LU Complete Pivoting’valueless

The $LU$ factorization implements a Markowitz-type search for pivots that locally minimize the fill-in subject to a threshold pivoting stability criterion. The default option is to use threshhold partial pivoting. The options ‘LU Rook Pivoting’ and ‘LU Complete Pivoting’ are more expensive than partial pivoting but are more stable and better at revealing rank, as long as ‘LU Factor Tolerance’ is not too large (say $<2.0$). When numerical difficulties are encountered, nlp2_sparse_solve automatically reduces the $LU$ tolerance towards $1.0$ and switches (if necessary) to rook or complete pivoting, before reverting to the default or specified options at the next refactorization (with ‘System Information Yes’, relevant messages are output to the ‘Print File’).

‘LU Rook Pivoting’valueless

The $LU$ factorization implements a Markowitz-type search for pivots that locally minimize the fill-in subject to a threshold pivoting stability criterion. The default option is to use threshhold partial pivoting. The options ‘LU Rook Pivoting’ and ‘LU Complete Pivoting’ are more expensive than partial pivoting but are more stable and better at revealing rank, as long as ‘LU Factor Tolerance’ is not too large (say $<2.0$). When numerical difficulties are encountered, nlp2_sparse_solve automatically reduces the $LU$ tolerance towards $1.0$ and switches (if necessary) to rook or complete pivoting, before reverting to the default or specified options at the next refactorization (with ‘System Information Yes’, relevant messages are output to the ‘Print File’).

‘Major Feasibility Tolerance’float

Default $=max({10}^{-6},\sqrt{ϵ})$

This tolerance, $r$, specifies how accurately the nonlinear constraints should be satisfied. The default value is appropriate when the linear and nonlinear constraints contain data to about that accuracy.

Let $v_{max}$ be the maximum nonlinear constraint violation, normalized by the size of the solution, which is required to satisfy

$$v_{max}={max}_i{v}_i/\parallel x\parallel \le r\text{,}$$

where $v_{\text{i}}$ is the violation of the $\text{i}$th nonlinear constraint, for $\text{i}=1,2,\dots ,\text{nf}$.

In the major iteration log (see Minor Iteration Log), $v_{max}$ appears as the quantity labelled ‘Feasible’. If some of the problem functions are known to be of low accuracy, a larger ‘Major Feasibility Tolerance’ may be appropriate.

‘Major Optimality Tolerance’float

Default $=2max({10}^{-6},\sqrt{ϵ})$

This tolerance, $r$, specifies the final accuracy of the dual variables. On successful termination, nlp2_sparse_solve will have computed a solution $(x,s,\pi )$ such that

$$c_{max}={max}_j{c}_j/\parallel \pi \parallel \le r\text{,}$$

where $c_{\text{j}}$ is an estimate of the complementarity slackness for variable $\text{j}$, for $\text{j}=1,2,\dots ,n+\text{nf}$. The values $c_i$ are computed from the final QP solution using the reduced gradients $d_j=g_j-{\pi }^T{a}_j$ (where $g_j$ is the $j$th component of the objective gradient, $a_j$ is the associated column of the constraint matrix $\left(\begin{array}{cc}A& -I\end{array}\right)$, and $\pi$ is the set of QP dual variables):

$$\begin{array}{c}{c}_j=\left\{\begin{array}{cc}{d}_{j}min(x_j-l_j,1)& \text{if}{d}_j\ge 0\text{;}\\ -d_{j}min(u_j-x_j,1)& \text{if}{d}_j<0\text{.}\end{array}\right)\end{array}$$

In the ‘Print File’, $c_{max}$ appears as the quantity labelled ‘Optimal’.

‘Major Iterations Limit’int

Default $=max(1000,3max(n,\text{nf}))$

This is the maximum number of major iterations allowed. It is intended to guard against an excessive number of linearizations of the constraints. If $i=0$, optimality and feasibility are checked.

‘Major Print Level’int

Default $=1$

This controls the amount of output to the options ‘Print File’ and ‘Summary File’ at each major iteration. $\text{‘Major Print Level'}=0$ suppresses most output, except for error messages. $\text{‘Major Print Level'}=1$ gives normal output for linear and nonlinear problems, and $\text{‘Major Print Level'}=11$ gives additional details of the Jacobian factorization that commences each major iteration.

In general, the value being specified may be thought of as a binary number of the form

$$\text{‘Major Print Level'}JFDXbs$$

where each letter stands for a digit that is either $0$ or $1$ as follows:

$s$	a single line that gives a summary of each major iteration. (This entry in $JFDXbs$ is not strictly binary since the summary line is printed whenever $JFDXbs\ge 1$);
$b$	basis statistics, i.e., information relating to the basis matrix whenever it is refactorized. (This output is always provided if $JFDXbs\ge 10$);
$X$	$x_k$, the nonlinear variables involved in the objective function or the constraints. These appear under the heading ‘Jacobian variables’;
$D$	${\pi }_k$, the dual variables for the nonlinear constraints. These appear under the heading ‘Multiplier estimates’;
$F$	$f\left(x_k\right)$, the values of the nonlinear constraint functions;
$J$	$J\left(x_k\right)$, the Jacobian matrix. This appears under the heading ‘$x$ and Jacobian’.

To obtain output of any items $JFDXbs$, set the corresponding digit to $1$, otherwise to $0$.

If $J=1$, the Jacobian matrix will be output column-wise at the start of each major iteration. Column $j$ will be preceded by the value of the corresponding variable $x_j$ and a key to indicate whether the variable is basic, superbasic or nonbasic. (Hence if $J=1$, there is no reason to specify $X=1$ unless the objective contains more nonlinear variables than the Jacobian.) A typical line of output is

3 1.250000e+01 BS 1 1.00000e+00 4 2.00000e+00

which would mean that $x_3$ is basic at value $12.5$, and the third column of the Jacobian has elements of $1.0$ and $2.0$ in rows $1$ and $4$.

‘Major Step Limit’float

Default $=2.0$

This argument limits the change in $x$ during a linesearch. It applies to all nonlinear problems, once a ‘feasible solution’ or ‘feasible subproblem’ has been found.

1. A linesearch determines a step $\alpha$ over the range $0<\alpha \le \beta$, where $\beta$ is $1$ if there are nonlinear constraints or is the step to the nearest upper or lower bound on $x$ if all the constraints are linear. Normally, the first step length tried is ${\alpha }_1=min(1,\beta )$.

2. In some cases, such as $f\left(x\right)=a{e}^{bx}$ or $f\left(x\right)=a{x}^b$, even a moderate change in the components of $x$ can lead to floating-point overflow. The argument $r$ is, therefore, used to define a limit $\overline{\beta }=r(1+\parallel x\parallel )/\parallel p\parallel$ (where $p$ is the search direction), and the first evaluation of $f\left(x\right)$ is at the potentially smaller step length ${\alpha }_1=min(1,\overline{\beta },\beta )$.

3. Wherever possible, upper and lower bounds on $x$ should be used to prevent evaluation of nonlinear functions at meaningless points. The option ‘Major Step Limit’ provides an additional safeguard. The default value $r=2.0$ should not affect progress on well behaved problems, but setting $r=0.1\text{or}0.01$ may be helpful when rapidly varying functions are present. A ‘good’ starting point may be required. An important application is to the class of nonlinear least squares problems.

4. In cases where several local optima exist, specifying a small value for $r$ may help locate an optimum near the starting point.

‘Minimize’valueless

Default

The keywords ‘Minimize’ and ‘Maximize’ specify the required direction of optimization. It applies to both linear and nonlinear terms in the objective.

The keyword ‘Feasible Point’ means ‘Ignore the objective function, while finding a feasible point for the linear and nonlinear constraints’. It can be used to check that the nonlinear constraints are feasible without altering the call to nlp2_sparse_solve.

‘Maximize’valueless

The keywords ‘Minimize’ and ‘Maximize’ specify the required direction of optimization. It applies to both linear and nonlinear terms in the objective.

The keyword ‘Feasible Point’ means ‘Ignore the objective function, while finding a feasible point for the linear and nonlinear constraints’. It can be used to check that the nonlinear constraints are feasible without altering the call to nlp2_sparse_solve.

‘Feasible Point’valueless

The keywords ‘Minimize’ and ‘Maximize’ specify the required direction of optimization. It applies to both linear and nonlinear terms in the objective.

The keyword ‘Feasible Point’ means ‘Ignore the objective function, while finding a feasible point for the linear and nonlinear constraints’. It can be used to check that the nonlinear constraints are feasible without altering the call to nlp2_sparse_solve.

‘Minor Feasibility Tolerance’float

Default $=max({10}^{-6},\sqrt{ϵ})$

nlp2_sparse_solve tries to ensure that all variables eventually satisfy their upper and lower bounds to within this tolerance, $r$. This includes slack variables. Hence, general linear constraints should also be satisfied to within $r$.

Feasibility with respect to nonlinear constraints is judged by the option ‘Major Feasibility Tolerance’ (not by $r$).

If the bounds and linear constraints cannot be satisfied to within $r$, the problem is declared infeasible. If $sinf$ is quite small, it may be appropriate to raise $r$ by a factor of $10$ or $100$. Otherwise, some error in the data should be suspected.

Nonlinear functions will be evaluated only at points that satisfy the bounds and linear constraints. If there are regions where a function is undefined, every attempt should be made to eliminate these regions from the problem.

For example, if $f\left(x\right)=\sqrt{x_1}+\log \left(x_2\right)$, it is essential to place lower bounds on both variables. If $r=1.0e-6$, the bounds $x_1\ge {10}^{-5}$ and $x_2\ge {10}^{-4}$ might be appropriate. (The log singularity is more serious. In general, keep $x$ as far away from singularities as possible.)

If $\text{‘Scale Option'}\ge 1$, feasibility is defined in terms of the scaled problem (since it is then more likely to be meaningful).

In reality, nlp2_sparse_solve uses $r$ as a feasibility tolerance for satisfying the bounds on $x$ and $s$ in each QP subproblem. If the sum of infeasibilities cannot be reduced to zero, the QP subproblem is declared infeasible. nlp2_sparse_solve is then in elastic mode thereafter (with only the linearized nonlinear constraints defined to be elastic). See the description of the option ‘Elastic Weight’.

‘Feasibility Tolerance’float

Default $=max\{{10}^{-6},\sqrt{ϵ}\}$

nlp2_sparse_solve tries to ensure that all variables eventually satisfy their upper and lower bounds to within this tolerance, $r$. This includes slack variables. Hence, general linear constraints should also be satisfied to within $r$.

Feasibility with respect to nonlinear constraints is judged by the option ‘Major Feasibility Tolerance’ (not by $r$).

If the bounds and linear constraints cannot be satisfied to within $r$, the problem is declared infeasible. If $sinf$ is quite small, it may be appropriate to raise $r$ by a factor of $10$ or $100$. Otherwise, some error in the data should be suspected.

Nonlinear functions will be evaluated only at points that satisfy the bounds and linear constraints. If there are regions where a function is undefined, every attempt should be made to eliminate these regions from the problem.

For example, if $f\left(x\right)=\sqrt{x_1}+\log \left(x_2\right)$, it is essential to place lower bounds on both variables. If $r=1.0e-6$, the bounds $x_1\ge {10}^{-5}$ and $x_2\ge {10}^{-4}$ might be appropriate. (The log singularity is more serious. In general, keep $x$ as far away from singularities as possible.)

If $\text{‘Scale Option'}\ge 1$, feasibility is defined in terms of the scaled problem (since it is then more likely to be meaningful).

In reality, nlp2_sparse_solve uses $r$ as a feasibility tolerance for satisfying the bounds on $x$ and $s$ in each QP subproblem. If the sum of infeasibilities cannot be reduced to zero, the QP subproblem is declared infeasible. nlp2_sparse_solve is then in elastic mode thereafter (with only the linearized nonlinear constraints defined to be elastic). See the description of the option ‘Elastic Weight’.

‘Minor Iterations Limit’int

Default $=500$

If the number of minor iterations for the optimality phase of the QP subproblem exceeds $i$, then all nonbasic QP variables that have not yet moved are frozen at their current values and the reduced QP is solved to optimality.

Note that more than $i$ minor iterations may be necessary to solve the reduced QP to optimality. These extra iterations are necessary to ensure that the terminated point gives a suitable direction for the linesearch.

In the major iteration log (see Minor Iteration Log) a t at the end of a line indicates that the corresponding QP was artificially terminated using the limit $i$.

Compare with the option ‘Iterations Limit’, which defines an independent absolute limit on the total number of minor iterations (summed over all QP subproblems).

‘Minor Print Level’int

Default $=1$

This controls the amount of output to the ‘Print File’ and ‘Summary File’ during solution of the QP subproblems. The value of $i$ has the following effect:

$i$	Meaning
$0$	No minor iteration output except error messages.
$\ge 1$	A single line of output at each minor iteration (controlled by options ‘Print Frequency’ and ‘Summary Frequency’.
$\ge 10$	Basis factorization statistics generated during the periodic refactorization of the basis (see the option ‘Factorization Frequency’). Statistics for the first factorization each major iteration are controlled by the option ‘Major Print Level’.

‘New Basis File’int

Default $=0$

‘New Basis File’ and ‘Backup Basis File’ are sometimes referred to as basis maps. They contain the most compact representation of the state of each variable. They are intended for restarting the solution of a problem at a point that was reached by an earlier run. For nontrivial problems, it is advisable to save basis maps at the end of a run, in order to restart the run if necessary.

If $i_1>0$, a basis map will be saved in the file associated with unit $i_1$ every $i_3$th iteration. The first record of the file will contain the word PROCEEDING if the run is still in progress. A basis map will also be saved at the end of a run, with some other word indicating the final solution status.

Use of $i_2>0$ is intended as a safeguard against losing the results of a long run. Suppose that a ‘New Basis File’ is being saved every $100$ (‘Save Frequency’) iterations, and that nlp2_sparse_solve is about to save such a basis at iteration $2000$. It is conceivable that the run may be interrupted during the next few milliseconds (in the middle of the save). In this case the Basis file will be corrupted and the run will have been essentially wasted.

To eliminate this risk, both a ‘New Basis File’ and a ‘Backup Basis File’ may be specified. The following would be suitable for the above example:

Backup Basis File 11
New Basis File 12

The current basis will then be saved every $100$ iterations, first in the file associated with unit $12$ and then immediately in the file associated with unit $11$. If the run is interrupted at iteration $2000$ during the save in the file associated with unit $12$, there will still be a usable basis in the file associated with unit $11$ (corresponding to iteration $1900$).

Note that a new basis will be saved in ‘New Basis File’ at the end of a run if it terminates normally, but it will not be saved in ‘Backup Basis File’. In the above example, if an optimum solution is found at iteration $2050$ (or if the iteration limit is $2050$), the final basis in the file associated with unit $12$ will correspond to iteration $2050$, but the last basis saved in the file associated with unit $11$ will be the one for iteration $2000$.

A full description of information recorded in ‘New Basis File’ and ‘Backup Basis File’ is given in Gill et al. (2005a).

‘Backup Basis File’int

Default $=0$

‘New Basis File’ and ‘Backup Basis File’ are sometimes referred to as basis maps. They contain the most compact representation of the state of each variable. They are intended for restarting the solution of a problem at a point that was reached by an earlier run. For nontrivial problems, it is advisable to save basis maps at the end of a run, in order to restart the run if necessary.

If $i_1>0$, a basis map will be saved in the file associated with unit $i_1$ every $i_3$th iteration. The first record of the file will contain the word PROCEEDING if the run is still in progress. A basis map will also be saved at the end of a run, with some other word indicating the final solution status.

Use of $i_2>0$ is intended as a safeguard against losing the results of a long run. Suppose that a ‘New Basis File’ is being saved every $100$ (‘Save Frequency’) iterations, and that nlp2_sparse_solve is about to save such a basis at iteration $2000$. It is conceivable that the run may be interrupted during the next few milliseconds (in the middle of the save). In this case the Basis file will be corrupted and the run will have been essentially wasted.

To eliminate this risk, both a ‘New Basis File’ and a ‘Backup Basis File’ may be specified. The following would be suitable for the above example:

Backup Basis File 11
New Basis File 12

The current basis will then be saved every $100$ iterations, first in the file associated with unit $12$ and then immediately in the file associated with unit $11$. If the run is interrupted at iteration $2000$ during the save in the file associated with unit $12$, there will still be a usable basis in the file associated with unit $11$ (corresponding to iteration $1900$).

Note that a new basis will be saved in ‘New Basis File’ at the end of a run if it terminates normally, but it will not be saved in ‘Backup Basis File’. In the above example, if an optimum solution is found at iteration $2050$ (or if the iteration limit is $2050$), the final basis in the file associated with unit $12$ will correspond to iteration $2050$, but the last basis saved in the file associated with unit $11$ will be the one for iteration $2000$.

A full description of information recorded in ‘New Basis File’ and ‘Backup Basis File’ is given in Gill et al. (2005a).

‘Save Frequency’int

Default $=100$

‘New Basis File’ and ‘Backup Basis File’ are sometimes referred to as basis maps. They contain the most compact representation of the state of each variable. They are intended for restarting the solution of a problem at a point that was reached by an earlier run. For nontrivial problems, it is advisable to save basis maps at the end of a run, in order to restart the run if necessary.

If $i_1>0$, a basis map will be saved in the file associated with unit $i_1$ every $i_3$th iteration. The first record of the file will contain the word PROCEEDING if the run is still in progress. A basis map will also be saved at the end of a run, with some other word indicating the final solution status.

Use of $i_2>0$ is intended as a safeguard against losing the results of a long run. Suppose that a ‘New Basis File’ is being saved every $100$ (‘Save Frequency’) iterations, and that nlp2_sparse_solve is about to save such a basis at iteration $2000$. It is conceivable that the run may be interrupted during the next few milliseconds (in the middle of the save). In this case the Basis file will be corrupted and the run will have been essentially wasted.

To eliminate this risk, both a ‘New Basis File’ and a ‘Backup Basis File’ may be specified. The following would be suitable for the above example:

Backup Basis File 11
New Basis File 12

The current basis will then be saved every $100$ iterations, first in the file associated with unit $12$ and then immediately in the file associated with unit $11$. If the run is interrupted at iteration $2000$ during the save in the file associated with unit $12$, there will still be a usable basis in the file associated with unit $11$ (corresponding to iteration $1900$).

Note that a new basis will be saved in ‘New Basis File’ at the end of a run if it terminates normally, but it will not be saved in ‘Backup Basis File’. In the above example, if an optimum solution is found at iteration $2050$ (or if the iteration limit is $2050$), the final basis in the file associated with unit $12$ will correspond to iteration $2050$, but the last basis saved in the file associated with unit $11$ will be the one for iteration $2000$.

A full description of information recorded in ‘New Basis File’ and ‘Backup Basis File’ is given in Gill et al. (2005a).

‘New Superbasics Limit’int

Default $=99$

This option causes early termination of the QP subproblems if the number of free variables has increased significantly since the first feasible point. If the number of new superbasics is greater than $i$, the nonbasic variables that have not yet moved are frozen and the resulting smaller QP is solved to optimality.

In the major iteration log (see Major Iteration Log), a t at the end of a line indicates that the QP was terminated early in this way.

‘Nolist’valueless

Default

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

‘List’valueless

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

‘Old Basis File’int

Default $=0$

If $i>0$, the basis maps information will be obtained from this file. The file will usually have been output previously as a ‘New Basis File’ or ‘Backup Basis File’. A full description of information recorded in ‘New Basis File’ and ‘Backup Basis File’ is given in Gill et al. (2005a).

The file will not be acceptable if the number of rows or columns in the problem has been altered.

‘Partial Price’int

Default $=1$

This argument is recommended for large problems that have significantly more variables than constraints. It reduces the work required for each ‘pricing’ operation (where a nonbasic variable is selected to become superbasic). When $i=1$, all columns of the constraint matrix $\left(\begin{array}{cc}A& -I\end{array}\right)$ are searched. Otherwise, $A$ and $I$ are partitioned to give $i$ roughly equal segments $A_{\text{j}}$ and $I_{\text{j}}$, for $\text{j}=1,2,\dots ,i$. If the previous pricing search was successful on $A_{j-1}$ and $I_{j-1}$, the next search begins on the segments $A_j$ and $I_j$. (All subscripts here are modulo $i$.) If a reduced gradient is found that is larger than some dynamic tolerance, the variable with the largest such reduced gradient (of appropriate sign) is selected to become superbasic. If nothing is found, the search continues on the next segments $A_{j+1}$ and $I_{j+1}$, and so on.

For time-stage models having $t$ time periods, ‘Partial Price’ $t$ (or $t/2$ or $t/3$) may be appropriate.

‘Pivot Tolerance’float

Default $={ϵ}^{\frac{2}{3}}$

During the solution of QP subproblems, the pivot tolerance is used to prevent columns entering the basis if they would cause the basis to become almost singular.

When $x$ changes to $x+\alpha p$ for some search direction $p$, a ‘ratio test’ determines which component of $x$ reaches an upper or lower bound first. The corresponding element of $p$ is called the pivot element. Elements of $p$ are ignored (and, therefore, cannot be pivot elements) if they are smaller than the pivot tolerance $r$.

It is common for two or more variables to reach a bound at essentially the same time. In such cases, the ‘Minor Feasibility Tolerance’ (say, $t$) provides some freedom to maximize the pivot element and thereby improve numerical stability. Excessively small values of $t$ should, therefore, not be specified. To a lesser extent, the ‘Expand Frequency’ (say, $f$) also provides some freedom to maximize the pivot element. Excessively large values of $f$ should, therefore, not be specified.

‘Print File’int

Default $=0$

If $i>0$, the following information is output to a file associated with unit $i$ during the solution of each problem:

• a listing of the options;

• some statistics about the problem;

• the amount of storage available for the $LU$ factorization of the basis matrix;

• notes about the initial basis resulting from a Crash procedure or a Basis file;

• the iteration log;

• basis factorization statistics;

• the exit $\text{errno}$ condition and some statistics about the solution obtained;

• the printed solution, if requested.

These items are described in Further Comments and Monitoring Information. Further brief output may be directed to the ‘Summary File’.

‘Print Frequency’int

Default $=100$

If $i>0$, one line of the iteration log will be printed every $i$th iteration. A value such as $i=10$ is suggested for those interested only in the final solution. If $i\le 0$, the value of $i=99999999$ is used and effectively no checks are made.

‘Proximal Point Method’int

Default $=1$

$i=1\text{or}2$ specifies minimization of ${\parallel x-x_0\parallel }_1$ or $\frac{1}{2}{\parallel x-x_0\parallel }_2^2$ when the starting point $x_0$ is changed to satisfy the linear constraints (where $x_0$ refers to nonlinear variables).

‘Punch File’int

Default $=0$

The ‘Punch File’ from a previous run may be used as an ‘Insert File’ for a later run on the same problem. A full description of information recorded in ‘Insert File’ and ‘Punch File’ is given in Gill et al. (2005a).

If $i_1>0$, the final solution obtained will be output to the file. For linear programs, this format is compatible with various commercial systems.

If $i_2>0$ the ‘Insert File’ containing basis information will be read from unit $i_2$. The file will usually have been output previously as a ‘Punch File’. The file will not be accessed if ‘Old Basis File’ is specified.

‘Insert File’int

Default $=0$

The ‘Punch File’ from a previous run may be used as an ‘Insert File’ for a later run on the same problem. A full description of information recorded in ‘Insert File’ and ‘Punch File’ is given in Gill et al. (2005a).

If $i_1>0$, the final solution obtained will be output to the file. For linear programs, this format is compatible with various commercial systems.

If $i_2>0$ the ‘Insert File’ containing basis information will be read from unit $i_2$. The file will usually have been output previously as a ‘Punch File’. The file will not be accessed if ‘Old Basis File’ is specified.

‘Scale Option’int

Default $=0$

Three scale options are available as follows:

$i$	Meaning
0	No scaling. This is recommended if it is known that $x$ and the constraint matrix never have very large elements (say, larger than $100$).
1	The constraints and variables are scaled by an iterative procedure that attempts to make the matrix coefficients as close as possible to $1.0$ (see Fourer (1982)). This will sometimes improve the performance of the solution procedures.
2	The constraints and variables are scaled by the iterative procedure. Also, a certain additional scaling is performed that may be helpful if the right-hand side $b$ or the solution $x$ is large. This takes into account columns of $\left(\begin{array}{cc}A& -I\end{array}\right)$ that are fixed or have positive lower bounds or negative upper bounds.

Option ‘Scale Tolerance’ affects how many passes might be needed through the constraint matrix. On each pass, the scaling procedure computes the ratio of the largest and smallest nonzero coefficients in each column:

$${\rho }_j={max}_j\left|a_{ij}\right|/{min}_i\left|a_{ij}\right|(a_{ij}\ne 0)\text{.}$$

If ${max}_j{\rho }_j$ is less than $r$ times its previous value, another scaling pass is performed to adjust the row and column scales. Raising $r$ from $0.9$ to $0.99$ (say) usually increases the number of scaling passes through $A$. At most $10$ passes are made. The value of $r$ should lie in the range $0<r<1$.

‘Scale Print’ causes the row scales $r\left(i\right)$ and column scales $c\left(j\right)$ to be printed to ‘Print File’, if ‘System Information Yes’ has been specified. The scaled matrix coefficients are ${\overline{a}}_{ij}=a_{ij}c\left(j\right)/r\left(i\right)$, and the scaled bounds on the variables and slacks are ${\overline{l}}_j=l_j/c\left(j\right)$, ${\overline{u}}_j=u_j/c\left(j\right)$, where $c\left(j\right)=r(j-n)$ if $j>n$.

‘Scale Tolerance’float

Default $=0.9$

Three scale options are available as follows:

$i$	Meaning
0	No scaling. This is recommended if it is known that $x$ and the constraint matrix never have very large elements (say, larger than $100$).
1	The constraints and variables are scaled by an iterative procedure that attempts to make the matrix coefficients as close as possible to $1.0$ (see Fourer (1982)). This will sometimes improve the performance of the solution procedures.
2	The constraints and variables are scaled by the iterative procedure. Also, a certain additional scaling is performed that may be helpful if the right-hand side $b$ or the solution $x$ is large. This takes into account columns of $\left(\begin{array}{cc}A& -I\end{array}\right)$ that are fixed or have positive lower bounds or negative upper bounds.

Option ‘Scale Tolerance’ affects how many passes might be needed through the constraint matrix. On each pass, the scaling procedure computes the ratio of the largest and smallest nonzero coefficients in each column:

$${\rho }_j={max}_j\left|a_{ij}\right|/{min}_i\left|a_{ij}\right|(a_{ij}\ne 0)\text{.}$$

If ${max}_j{\rho }_j$ is less than $r$ times its previous value, another scaling pass is performed to adjust the row and column scales. Raising $r$ from $0.9$ to $0.99$ (say) usually increases the number of scaling passes through $A$. At most $10$ passes are made. The value of $r$ should lie in the range $0<r<1$.

‘Scale Print’ causes the row scales $r\left(i\right)$ and column scales $c\left(j\right)$ to be printed to ‘Print File’, if ‘System Information Yes’ has been specified. The scaled matrix coefficients are ${\overline{a}}_{ij}=a_{ij}c\left(j\right)/r\left(i\right)$, and the scaled bounds on the variables and slacks are ${\overline{l}}_j=l_j/c\left(j\right)$, ${\overline{u}}_j=u_j/c\left(j\right)$, where $c\left(j\right)=r(j-n)$ if $j>n$.

‘Scale Print’valueless

Three scale options are available as follows:

$i$	Meaning
0	No scaling. This is recommended if it is known that $x$ and the constraint matrix never have very large elements (say, larger than $100$).
1	The constraints and variables are scaled by an iterative procedure that attempts to make the matrix coefficients as close as possible to $1.0$ (see Fourer (1982)). This will sometimes improve the performance of the solution procedures.
2	The constraints and variables are scaled by the iterative procedure. Also, a certain additional scaling is performed that may be helpful if the right-hand side $b$ or the solution $x$ is large. This takes into account columns of $\left(\begin{array}{cc}A& -I\end{array}\right)$ that are fixed or have positive lower bounds or negative upper bounds.

Option ‘Scale Tolerance’ affects how many passes might be needed through the constraint matrix. On each pass, the scaling procedure computes the ratio of the largest and smallest nonzero coefficients in each column:

$${\rho }_j={max}_j\left|a_{ij}\right|/{min}_i\left|a_{ij}\right|(a_{ij}\ne 0)\text{.}$$

If ${max}_j{\rho }_j$ is less than $r$ times its previous value, another scaling pass is performed to adjust the row and column scales. Raising $r$ from $0.9$ to $0.99$ (say) usually increases the number of scaling passes through $A$. At most $10$ passes are made. The value of $r$ should lie in the range $0<r<1$.

‘Scale Print’ causes the row scales $r\left(i\right)$ and column scales $c\left(j\right)$ to be printed to ‘Print File’, if ‘System Information Yes’ has been specified. The scaled matrix coefficients are ${\overline{a}}_{ij}=a_{ij}c\left(j\right)/r\left(i\right)$, and the scaled bounds on the variables and slacks are ${\overline{l}}_j=l_j/c\left(j\right)$, ${\overline{u}}_j=u_j/c\left(j\right)$, where $c\left(j\right)=r(j-n)$ if $j>n$.

‘Solution File’int

Default $=0$

If $i>0$, the final solution will be output to file $i$ (whether optimal or not). All numbers are printed in 1pe16.6 format.

To see more significant digits in the printed solution, it will sometimes be useful to make $i$ refer to ‘Print File’.

‘Summary File’int

Default $=0$

If $i_1>0$, a brief log will be output to the file associated with unit $i_1$, including one line of information every $i_2$th iteration. In an interactive environment, it is useful to direct this output to the terminal, to allow a run to be monitored online. (If something looks wrong, the run can be manually terminated.) Further details are given in The Summary File.

‘Summary Frequency’int

Default $=100$

If $i_1>0$, a brief log will be output to the file associated with unit $i_1$, including one line of information every $i_2$th iteration. In an interactive environment, it is useful to direct this output to the terminal, to allow a run to be monitored online. (If something looks wrong, the run can be manually terminated.) Further details are given in The Summary File.

‘Superbasics Limit’int

Default $=n_1$

This option places a limit on the storage allocated for superbasic variables. Ideally, $i$ should be set slightly larger than the ‘number of degrees of freedom’ expected at an optimal solution.

For nonlinear problems, the number of degrees of freedom is often called the ‘number of independent variables’. Normally, $i$ need not be greater than $n+1$, where $n_1$ is the number of nonlinear variables. For many problems, $i$ may be considerably smaller than $n$. This will save storage if $n$ is very large.

‘Suppress Parameters’valueless

Normally nlp2_sparse_solve prints the options file as it is being read, and then prints a complete list of the available keywords and their final values. The option ‘Suppress Parameters’ tells nlp2_sparse_solve not to print the full list.

‘System Information No’valueless

Default

This option prints additional information on the progress of major and minor iterations, and Crash statistics. See Monitoring Information.

‘System Information Yes’valueless

This option prints additional information on the progress of major and minor iterations, and Crash statistics. See Monitoring Information.

‘Timing Level’int

Default $=0$

If $i>0$, some timing information will be output to the Print file, if $\text{‘Print File'}>0$.

‘Unbounded Objective’float

Default $=1.0e+15$

These arguments are intended to detect unboundedness in nonlinear problems. During a linesearch, $F_{obj}$ is evaluated at points of the form $x+\alpha p$, where $x$ and $p$ are fixed and $\alpha$ varies. If $\left|F_{obj}\right|$ exceeds $r_1$ or $\alpha$ exceeds $r_2$, iterations are terminated with the exit message $errno$ = 5.

If singularities are present, unboundedness in $F_{obj}\left(x\right)$ may be manifested by a floating-point overflow (during the evaluation of $F_{obj}(x+\alpha p)$), before the test against $r_1$ can be made.

Unboundedness in $x$ is best avoided by placing finite upper and lower bounds on the variables.

‘Unbounded Step Size’float

Default $=\text{infbnd}$

These arguments are intended to detect unboundedness in nonlinear problems. During a linesearch, $F_{obj}$ is evaluated at points of the form $x+\alpha p$, where $x$ and $p$ are fixed and $\alpha$ varies. If $\left|F_{obj}\right|$ exceeds $r_1$ or $\alpha$ exceeds $r_2$, iterations are terminated with the exit message $errno$ = 5.

If singularities are present, unboundedness in $F_{obj}\left(x\right)$ may be manifested by a floating-point overflow (during the evaluation of $F_{obj}(x+\alpha p)$), before the test against $r_1$ can be made.

Unboundedness in $x$ is best avoided by placing finite upper and lower bounds on the variables.

‘Verify Level’int

Default $=0$

This option refers to finite difference checks on the derivatives computed by the user-supplied functions. Derivatives are checked at the first point that satisfies all bounds and linear constraints.

$i$	Meaning
$0$	Only a ‘cheap’ test will be performed, requiring two calls to $usrfun$.
$1$	Individual gradients will be checked (with a more reliable test). A key of the form OK or Bad? indicates whether or not each component appears to be correct.
$2$	Individual columns of the problem Jacobian will be checked.
$3$	Options 2 and 1 will both occur (in that order).
$-1$	Derivative checking is disabled.

$\text{‘Verify Level'}=3$ should be specified whenever a new $usrfun$ is being developed.

‘Violation Limit’float

Default $=1.0e+6$

This keyword defines an absolute limit on the magnitude of the maximum constraint violation, $r$, after the linesearch. On completion of the linesearch, the new iterate $x_{k+1}$ satisfies the condition

$$v_i\left(x_{k+1}\right)\le rmax(1,v_i\left(x_0\right))\text{,}$$

where $x_0$ is the point at which the nonlinear constraints are first evaluated and $v_i\left(x\right)$ is the $i$th nonlinear constraint violation $v_i\left(x\right)=max(0,l_i-f_i\left(x\right),f_i\left(x\right)-u_i)$.

The effect of this violation limit is to restrict the iterates to lie in an expanded feasible region whose size depends on the magnitude of $r$. This makes it possible to keep the iterates within a region where the objective is expected to be well-defined and bounded below. If the objective is bounded below for all values of the variables, $r$ may be any large positive value.

Raises

NagValueError

(errno $1$)

The initialization function nlp2_sparse_init() has not been called.

(errno $2$)

On entry, one but not both of $\text{nxname}$ and $\text{nfname}$ is equal to $1$. $\text{nxname}=⟨value⟩$ and $\text{nfname}=⟨value⟩$.

(errno $2$)

On entry, $objrow=⟨value⟩$ and $\text{nf}=⟨value⟩$.

Constraint: $0\le objrow\le \text{nf}$.

(errno $2$)

On entry, bounds $xlow$ and $xupp$ for $⟨value⟩$ are equal and infinite. $xlow=xupp=⟨value⟩$ and $\text{infbnd}=⟨value⟩$.

(errno $2$)

On entry, bounds for $⟨value⟩$ are inconsistent. $xlow=⟨value⟩$ and $xupp=⟨value⟩$.

(errno $2$)

On entry, bounds $xlow$ and $xupp$ for variable $⟨value⟩$ are equal and infinite. $xlow=xupp=⟨value⟩$ and $\text{infbnd}=⟨value⟩$.

(errno $2$)

On entry, bounds for variable $⟨value⟩$ are inconsistent. $xlow=⟨value⟩$ and $xupp=⟨value⟩$.

(errno $2$)

On entry, bounds $flow$ and $fupp$ for $⟨value⟩$ are equal and infinite. $flow=fupp=⟨value⟩$ and $\text{infbnd}=⟨value⟩$.

(errno $2$)

On entry, bounds for $⟨value⟩$ are inconsistent. $flow=⟨value⟩$ and $fupp=⟨value⟩$.

(errno $2$)

On entry, bounds $flow$ and $fupp$ for variable $⟨value⟩$ are equal and infinite. $flow=fupp=⟨value⟩$ and $\text{infbnd}=⟨value⟩$.

(errno $2$)

On entry, bounds for variable $⟨value⟩$ are inconsistent. $flow=⟨value⟩$ and $fupp=⟨value⟩$.

(errno $2$)

Basis file dimensions do not match this problem.

(errno $2$)

Array element $igfun[⟨value⟩]=⟨value⟩$ is out of range $1$ to $\text{nf}=⟨value⟩$, or array element $jgvar[⟨value⟩]=⟨value⟩$ is out of range $1$ to $n=⟨value⟩$.

(errno $2$)

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

Constraint: $\text{lena}\ge 1$.

(errno $2$)

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

Constraint: $\text{leng}\ge 1$.

(errno $2$)

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

Constraint: $\text{nf}\ge 1$.

(errno $2$)

On entry, $n=⟨value⟩$.

Constraint: $n\ge 1$.

(errno $2$)

On entry, $start=⟨value⟩$.

Constraint: $start=0$, $1$ or $2$.

(errno $2$)

On entry, $nea=⟨value⟩$.

Constraint: $nea\ge 0$.

(errno $2$)

On entry, $neg=⟨value⟩$.

Constraint: $neg\ge 0$.

(errno $2$)

On entry, $nea=⟨value⟩$, $n=⟨value⟩$ and $\text{nf}=⟨value⟩$.

Constraint: $nea\le n\times \text{nf}$.

(errno $2$)

On entry, $neg=⟨value⟩$, $n=⟨value⟩$ and $\text{nf}=⟨value⟩$.

Constraint: $neg\le n\times \text{nf}$.

(errno $2$)

On entry, $\text{nxname}=⟨value⟩$ and $n=⟨value⟩$.

Constraint: $\text{nxname}=1$ or $\text{n}$.

(errno $2$)

On entry, $\text{nfname}=⟨value⟩$ and $\text{nf}=⟨value⟩$.

Constraint: $\text{nfname}=1$ or $\text{nf}$.

(errno $8$)

User-supplied function computes incorrect objective derivatives.

(errno $8$)

User-supplied function computes incorrect constraint derivatives.

(errno $11$)

Internal error: memory allocation failed when attempting to allocate workspace sizes $⟨value⟩$, $⟨value⟩$ and $⟨value⟩$. Please contact NAG.

(errno $12$)

Internal memory allocation was insufficient. Please contact NAG.

(errno $13$)

An error has occurred in the basis package. Check that arrays $iafun$, $javar$, $igfun$ and $jgvar$ contain values in the appropriate ranges and do not define duplicate elements of $a$ or $g$. Set the option ‘Print File’ and examine the output carefully for further information.

(errno $14$)

An unexpected error has occurred. Set the option ‘Print File’ and examine the output carefully for further information.

Warns

NagAlgorithmicWarning

(errno $3$)

The requested accuracy could not be achieved.

(errno $4$)

The linear constraints appear to be infeasible.

(errno $4$)

The problem appears to be infeasible. The linear equality constraints could not be satisfied.

(errno $4$)

The problem appears to be infeasible. Nonlinear infeasibilites have been minimized.

(errno $4$)

The problem appears to be infeasible. Infeasibilites have been minimized.

(errno $5$)

The problem appears to be unbounded. The objective function is unbounded.

(errno $5$)

The problem appears to be unbounded. The constraint violation limit has been reached.

(errno $10$)

User-supplied function requested termination.

NagAlgorithmicMajorWarning

(errno $6$)

Iteration limit reached.

(errno $6$)

Major iteration limit reached.

(errno $6$)

The value of the option ‘Superbasics Limit’ is too small.

(errno $7$)

Numerical difficulties have been encountered and no further progress can be made.

(errno $9$)

User-supplied function is undefined at the first feasible point.

(errno $9$)

User-supplied function is undefined at the initial point.

(errno $9$)

Unable to proceed into undefined region of user-supplied function.

Notes

nlp2_sparse_solve is designed to minimize a linear or nonlinear function subject to bounds on the variables and sparse linear or nonlinear constraints. It is suitable for large-scale linear and quadratic programming and for linearly constrained optimization, as well as for general nonlinear programs of the form

$$\begin{array}{c}{minimize}_x{f}_0\left(x\right)\text{subject to}l\le \left(\begin{array}{c}x\\ f\left(x\right)\\ A_{L}x\end{array}\right)\le u\text{,}\end{array}$$

where $x$ is an $n$-vector of variables, $l$ and $u$ are constant lower and upper bounds, $f_0\left(x\right)$ is a smooth scalar objective function, $A_L$ is a sparse matrix, and $f\left(x\right)$ is a vector of smooth nonlinear constraint functions $\left\{f_i\left(x\right)\right\}$. The option ‘Maximize’ specifies that $f_0\left(x\right)$ should be maximized instead of minimized.

Ideally, the first derivatives (gradients) of $f_0\left(x\right)$ and $f_i\left(x\right)$ should be known and coded by you. If only some of the gradients are known, nlp2_sparse_solve estimates the missing ones by finite differences.

If $f_0\left(x\right)$ is linear and $f\left(x\right)$ is absent, (1) is a linear program (LP) and nlp2_sparse_solve applies the primal simplex method (see Dantzig (1963)). Sparse basis factors are maintained by LUSOL (see Gill et al. (1987)) as in MINOS (see Murtagh and Saunders (1995)).

If only the objective is nonlinear, the problem is linearly constrained (LC) and tends to solve more easily than the general case with nonlinear constraints (NC). For both nonlinear cases, nlp2_sparse_solve applies a sparse Sequential Quadratic Programming (SQP) method (see Gill et al. (2002)), using limited-memory quasi-Newton approximations to the Hessian of the Lagrangian. The merit function for step-length control is an augmented Lagrangian, as in the dense SQP solver nlp2_solve() (see Gill et al. (1986) and Gill et al. (1992)).

nlp2_sparse_solve is suitable for nonlinear problems with thousands of constraints and variables, and is most efficient if only some of the variables enter nonlinearly, or there are relatively few degrees of freedom at a solution (i.e., many constraints are active). However, there is no limit on the number of degrees of freedom.

nlp2_sparse_solve allows linear and nonlinear constraints and variables to be entered in an arbitrary order, and uses one function to define all the nonlinear functions.

The optimization problem is assumed to be in the form

$${minimize}_x{F}_{obj}\left(x\right)\text{subject to}{l}_x\le x\le u_x\text{,}{l}_F\le F\left(x\right)\le u_F\text{,}$$

where the upper and lower bounds are constant, $F\left(x\right)$ is a vector of smooth linear and nonlinear constraint functions $\left\{F_i\left(x\right)\right\}$, and $F_{obj}\left(x\right)$ is one of the components of $F$ to be minimized, as specified by the input argument $objrow$. nlp2_sparse_solve reorders the variables and constraints so that the problem is in the form (1).

Upper and lower bounds are specified for all variables and functions. The $j$th constraint may be defined as an equality by setting $l_j=u_j$. If certain bounds are not present, the associated elements of $l$ or $u$ should be set to special values that are treated as $-\infty$ or $+\infty$. Free variables and free constraints (‘free rows’) have both bounds infinite.

In general, the components of $F$ are structured in the sense that they are formed from sums of linear and nonlinear functions of just some of the variables. This structure can be exploited by nlp2_sparse_solve.

In many cases, the vector $F\left(x\right)$ is a sum of linear and nonlinear functions. nlp2_sparse_solve allows these terms to be specified separately, so that the linear part is defined just once by the input arguments $iafun$, $javar$ and $a$. Only the nonlinear part is recomputed at each $x$.

Suppose that each component of $F\left(x\right)$ is of the form

$$F_i\left(x\right)=f_i\left(x\right)+\sum _{j=1}^n{A}_{ij}{x}_j\text{,}$$

where $f_i\left(x\right)$ is a nonlinear function (possibly zero) and the elements $A_{ij}$ are constant. The $nf\times n$ Jacobian of $F\left(x\right)$ is the sum of two sparse matrices of the same size: $F^{\prime }\left(x\right)=G\left(x\right)+A$, where $G\left(x\right)=f^{\prime }\left(x\right)$ and $A$ is the matrix with elements $\left\{A_{ij}\right\}$. The two matrices must be non-overlapping in the sense that each element of the Jacobian $F^{\prime }\left(x\right)=G\left(x\right)+A$ comes from $G\left(x\right)$ or $A$, but not both. The element cannot be split between $G\left(x\right)$ and $A$.

For example, the function

$$\begin{array}{c}F\left(x\right)=\left(\begin{array}{c}3x_1+e^{x_2}{x}_4+x_2^2+4x_4-x_3+x_5\\ x_2+x_3^2+\sin \left(x_4\right)-3x_5\\ x_1-x_3\end{array}\right)\end{array}$$

can be written as

$$\begin{array}{c}F\left(x\right)=f\left(x\right)+Ax=\left(\begin{array}{c}{e}^{x_2}{x}_4+x_2^2+4x_4\\ x_3^2+\sin \left(x_4\right)\\ 0\end{array}\right)+\left(\begin{array}{c}3x_1-x_3+x_5\\ x_2-3x_5\\ x_1-x_3\end{array}\right)\text{,}\end{array}$$

in which case

$$\begin{array}{c}{F}^{\prime }\left(x\right)=\left(\begin{array}{ccccc}3& e^{x_2}{x}_4+2x_2& -1& e^{x_2}+4& 1\\ 0& 1& 2x_3& \cos \left(x_4\right)& -3\\ 1& 0& -1& 0& 0\end{array}\right)\end{array}$$

can be written as $F^{\prime }\left(x\right)=f^{\prime }\left(x\right)+A=G\left(x\right)+A$, where

$$\begin{array}{c}G\left(x\right)=\left(\begin{array}{ccccc}0& e^{x_2}{x}_4+2x_2& 0& e^{x_2}+4& 0\\ 0& 0& 2x_3& \cos \left(x_4\right)& 0\\ 0& 0& 0& 0& 0\end{array}\right)\text{,}A=\left(\begin{array}{ccccc}3& 0& -1& 0& 1\\ 0& 1& 0& 0& -3\\ 1& 0& -1& 0& 0\end{array}\right)\text{.}\end{array}$$

Note: the element $e^{x_2}+4$ of $F^{\prime }\left(x\right)$ appears in $G\left(x\right)$ and is not split between $G\left(x\right)$ and $A$ although it contains a linear term.

The nonzero elements of $A$ and $G$ are provided to nlp2_sparse_solve in coordinate form. The elements of $A$ are entered as triples $(i,j,A_{ij})$ in the arrays $iafun$, $javar$ and $a$. The sparsity pattern $G$ is entered as pairs $(i,j)$ in the arrays $igfun$ and $jgvar$. The corresponding entries $G_{ij}$ (any that are known) are assigned to appropriate array elements $g\left(k\right)$ in $usrfun$.

The elements of $A$ and $G$ may be stored in any order. Duplicate entries are ignored. $igfun$ and $jgvar$ may be defined automatically by function nlp2_sparse_jacobian() when $\text{‘Derivative Option'}=0$ is specified and $usrfun$ does not provide any gradients.

Throughout this document the symbol $ϵ$ is used to represent the machine precision (see machine.precision).

nlp2_sparse_solve is based on SNOPTA, which is part of the SNOPT package described in Gill et al. (2005b).

References

Dantzig, G B, 1963, Linear Programming and Extensions, Princeton University Press

Eldersveld, S K, 1991, Large-scale sequential quadratic programming algorithms, PhD Thesis, Department of Operations Research, Stanford University, Stanford

Fourer, R, 1982, Solving staircase linear programs by the simplex method, Math. Programming (23), 274–313

Gill, P E, Murray, W and Saunders, M A, 2002, SNOPT: An SQP Algorithm for Large-scale Constrained Optimization (12), 979–1006, SIAM J. Optim.

Gill, P E, Murray, W and Saunders, M A, 2005, Users’ guide for SQOPT 7: a Fortran package for large-scale linear and quadratic programming, Report NA 05-1, Department of Mathematics, University of California, San Diego, https://www.ccom.ucsd.edu/~peg/papers/sqdoc7.pdf

Gill, P E, Murray, W and Saunders, M A, 2005, Users’ guide for SNOPT 7.1: a Fortran package for large-scale linear nonlinear programming, Report NA 05-2, Department of Mathematics, University of California, San Diego, https://www.ccom.ucsd.edu/~peg/papers/sndoc7.pdf

Gill, P E, Murray, W, Saunders, M A and Wright, M H, 1986, Users’ guide for NPSOL (Version 4.0): a Fortran package for nonlinear programming, Report SOL 86-2, Department of Operations Research, Stanford University

Gill, P E, Murray, W, Saunders, M A and Wright, M H, 1987, Maintaining $LU$ factors of a general sparse matrix, Linear Algebra and its Applics. (88/89), 239–270

Gill, P E, Murray, W, Saunders, M A and Wright, M H, 1992, Some theoretical properties of an augmented Lagrangian merit function, Advances in Optimization and Parallel Computing, (ed P M Pardalos), 101–128, North Holland

Hock, W and Schittkowski, K, 1981, Test Examples for Nonlinear Programming Codes. Lecture Notes in Economics and Mathematical Systems (187), Springer–Verlag

Murtagh, B A and Saunders, M A, 1978, Large-scale linearly constrained optimization, 14, 41–72, Math. Programming

Murtagh, B A and Saunders, M A, 1982, A projected Lagrangian algorithm and its implementation for sparse nonlinear constraints, Math. Program. Stud. (16), 84–118

Murtagh, B A and Saunders, M A, 1995, MINOS 5.4 users’ guide, Report SOL 83-20R, Department of Operations Research, Stanford University

See also

naginterfaces.library.examples.opt.nlp2_sparse_solve_ex.main()