NAG Toolbox: nag_opt_qp_dense_solve (e04nf)

Purpose

Syntax

Description

References

Parameters

Compulsory Input Parameters

1:     $\mathrm{a}\left(\mathit{lda},:\right)$ – double array

The first dimension of the array a must be at least $\max \left(1,\mathbf{nclin}\right)$.

The second dimension of the array a must be at least $\mathbf{n}$ if $\mathbf{nclin}>0$ and at least $1$ if $\mathbf{nclin}=0$.

The $\mathit{i}$th row of a must contain the coefficients of the $\mathit{i}$th general linear constraint, for $\mathit{i}=1,2,\dots ,m_L$.

If $\mathbf{nclin}=0$, a is not referenced.

2:     $\mathrm{bl}\left(\mathbf{n}+\mathbf{nclin}\right)$ – double array

3:     $\mathrm{bu}\left(\mathbf{n}+\mathbf{nclin}\right)$ – double array

bl must contain the lower bounds and bu the upper bounds, for all the constraints in the following order. The first $n$ elements of each array must contain the bounds on the variables, and the next $m_L$ elements the bounds for the general linear constraints (if any). To specify a nonexistent lower bound (i.e., $l_j=-\infty$), set $\mathbf{bl}\left(j\right)\le -\mathit{bigbnd}$, and to specify a nonexistent upper bound (i.e., $u_j=+\infty$), set $\mathbf{bu}\left(j\right)\ge \mathit{bigbnd}$; the default value of $\mathit{bigbnd}$ is ${10}^{20}$, but this may be changed by the optional parameter Infinite Bound Size. To specify the $j$th constraint as an equality, set $\mathbf{bl}\left(j\right)=\mathbf{bu}\left(j\right)=\beta$, say, where $\left|\beta \right|<\mathit{bigbnd}$.

Constraints:

• $\mathbf{bl}\left(\mathit{j}\right)\le \mathbf{bu}\left(\mathit{j}\right)$, for $\mathit{j}=1,2,\dots ,\mathbf{n}+\mathbf{nclin}$;
• if $\mathbf{bl}\left(j\right)=\mathbf{bu}\left(j\right)=\beta$, $\left|\beta \right|<\mathit{bigbnd}$.

4:     $\mathrm{cvec}\left(:\right)$ – double array

The dimension of the array cvec must be at least $\mathbf{n}$ if the problem is of type LP, QP2 (the default) or QP4, and at least $1$ otherwise

The coefficients of the explicit linear term of the objective function when the problem is of type LP, QP2 (the default) and QP4.

If the problem is of type FP, QP1, or QP3, cvec is not referenced.

5:     $\mathrm{h}\left(\mathit{ldh},\mathit{sdh}\right)$ – double array

ldh, the first dimension of the array, must satisfy the constraint

• if the problem is of type QP1, QP2 (the default), QP3 or QP4, $\mathit{ldh}\ge \mathbf{n}$ or at least the value of the optional parameter Hessian Rows;
• if the problem is of type FP or LP, $\mathit{ldh}\ge 1$.

May be used to store the quadratic term $H$ of the QP objective function if desired. In some cases, you need not use h to store $H$ explicitly (see the specification of function qphess). The elements of h are referenced only by function qphess. The number of rows of $H$ is denoted by $m$, whose default value is $n$. (The optional parameter Hessian Rows may be used to specify a value of $m<n$.)

If the default version of qphess is used and the problem is of type QP1 or QP2 (the default), the first $m$ rows and columns of h must contain the leading $m$ by $m$ rows and columns of the symmetric Hessian matrix $H$. Only the diagonal and upper triangular elements of the leading $m$ rows and columns of h are referenced. The remaining elements need not be assigned.

If the default version of qphess is used and the problem is of type QP3 or QP4, the first $m$ rows of h must contain an $m$ by $n$ upper trapezoidal factor of the symmetric Hessian matrix $H^{\mathrm{T}}H$. The factor need not be of full rank, i.e., some of the diagonal elements may be zero. However, as a general rule, the larger the dimension of the leading nonsingular sub-matrix of h, the fewer iterations will be required. Elements outside the upper trapezoidal part of the first $m$ rows of h need not be assigned.

If a non-default version of qphess is supplied, then in some cases it may be desirable to use a one-dimensional array to transmit data to qphess. h is then declared as an ldh by $1$ array, where $\mathit{ldh}\ge \mathbf{n}\times \left(\mathbf{n}+1\right)/2$.

In other situations, it may be desirable to compute $Hx$ or $H^{\mathrm{T}}Hx$ without accessing h – for example, if $H$ or $H^{\mathrm{T}}H$ is sparse or has special structure. The arguments h and ldh may then refer to any convenient array.

If the problem is of type FP or LP, h is not referenced.

6:     $\mathrm{qphess}$ – function handle or string containing name of m-file

In general, you need not provide a version of qphess, because a ‘default’ function with name nag_opt_qp_dense_sample_qphess (e54nfu) is included in the Library. However, the algorithm of nag_opt_qp_dense_solve (e04nf) requires only the product of $H$ or $H^{\mathrm{T}}H$ and a vector $x$; and in some cases you may obtain increased efficiency by providing a version of qphess that avoids the need to define the elements of the matrices $H$ or $H^{\mathrm{T}}H$ explicitly.

qphess is not referenced if the problem is of type FP or LP, in which case qphess may be the function nag_opt_qp_dense_sample_qphess (e54nfu).

[hx, user, iwsav] = qphess(n, jthcol, h, ldh, x, user, iwsav)

Input Parameters

1:     $\mathrm{n}$ – int64int32nag_int scalar

This is the same argument as supplied to this function. See the description for the top level argument n.

2:     $\mathrm{jthcol}$ – int64int32nag_int scalar

Specifies whether or not the vector $x$ is a column of the identity matrix.

$\mathbf{jthcol}=j>0$

The vector $x$ is the $j$th column of the identity matrix, and hence $Hx$ or $H^{\mathrm{T}}Hx$ is the $j$th column of $H$ or $H^{\mathrm{T}}H$, respectively. This may in some cases require very little computation and qphess may be coded to take advantage of this. However special code is not necessary because $x$ is always stored explicitly in the array x.

$\mathbf{jthcol}=0$

$x$ has no special form.

3:     $\mathrm{h}\left(\mathit{ldh},\mathit{sdh}\right)$ – double array

This is the same argument as supplied to this function. See the description for the top level argument h.

4:     $\mathrm{ldh}$ – int64int32nag_int scalar

This is the same argument as supplied to this function. See the description for the top level argument ldh.

5:     $\mathrm{x}\left(\mathbf{n}\right)$ – double array

The vector $x$.

6:     $\mathrm{user}$ – Any MATLAB object

qphess is called from nag_opt_qp_dense_solve (e04nf) with the object supplied to nag_opt_qp_dense_solve (e04nf).

7:     $\mathrm{iwsav}\left(610\right)$ – int64int32nag_int array

iwsav contains information that is required by the default function 'e54nfu'.

Output Parameters

1:     $\mathrm{hx}\left(\mathbf{n}\right)$ – double array

The product $Hx$ if the problem is of type QP1 or QP2 (the default), or the product $H^{\mathrm{T}}Hx$ if the problem is of type QP3 or QP4.

2:     $\mathrm{user}$ – Any MATLAB object

3:     $\mathrm{iwsav}\left(610\right)$ – int64int32nag_int array

7:     $\mathrm{istate}\left(\mathbf{n}+\mathbf{nclin}\right)$ – int64int32nag_int array

Need not be set if the (default) optional parameter Cold Start is used.

If the optional parameter Warm Start has been chosen, istate specifies the desired status of the constraints at the start of the feasibility phase. More precisely, the first $n$ elements of istate refer to the upper and lower bounds on the variables, and the next $m_L$ elements refer to the general linear constraints (if any). Possible values for $\mathbf{istate}\left(j\right)$ are as follows:

$\mathbf{istate}\left(j\right)$	Meaning
0	The corresponding constraint should not be in the initial working set.
1	The constraint should be in the initial working set at its lower bound.
2	The constraint should be in the initial working set at its upper bound.
3	The constraint should be in the initial working set as an equality. This value must not be specified unless $\mathbf{bl}\left(j\right)=\mathbf{bu}\left(j\right)$.

The values $-2$, $-1$ and $4$ are also acceptable but will be reset to zero by the function. If nag_opt_qp_dense_solve (e04nf) has been called previously with the same values of n and nclin, istate already contains satisfactory information. (See also the description of the optional parameter Warm Start.) The function also adjusts (if necessary) the values supplied in x to be consistent with istate.

Constraint: $-2\le \mathbf{istate}\left(\mathit{j}\right)\le 4$, for $\mathit{j}=1,2,\dots ,\mathbf{n}+\mathbf{nclin}$.

8:     $\mathrm{x}\left(\mathbf{n}\right)$ – double array

An initial estimate of the solution.

9:     $\mathrm{lwsav}\left(120\right)$ – logical array

10:   $\mathrm{iwsav}\left(610\right)$ – int64int32nag_int array

11:   $\mathrm{rwsav}\left(475\right)$ – double array

The arrays lwsav, iwsav and rwsav must not be altered between calls to any of the functions nag_opt_qp_dense_solve (e04nf), nag_opt_qp_dense_option_string (e04nh).

Optional Input Parameters

1:     $\mathrm{n}$ – int64int32nag_int scalar

Default: the dimension of the array x.

$n$, the number of variables.

Constraint: $\mathbf{n}>0$.

2:     $\mathrm{nclin}$ – int64int32nag_int scalar

Default: the first dimension of the array a.

$m_L$, the number of general linear constraints.

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

3:     $\mathrm{user}$ – Any MATLAB object

user is not used by nag_opt_qp_dense_solve (e04nf), but is passed to qphess. Note that for large objects it may be more efficient to use a global variable which is accessible from the m-files than to use user.

Output Parameters

1:     $\mathrm{istate}\left(\mathbf{n}+\mathbf{nclin}\right)$ – int64int32nag_int array

The status of the constraints in the working set at the point returned in x. The significance of each possible value of $\mathbf{istate}\left(j\right)$ is as follows:

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

2:     $\mathrm{x}\left(\mathbf{n}\right)$ – double array

The point at which nag_opt_qp_dense_solve (e04nf) terminated. If $\mathbf{ifail}=\mathbf{0}$, $\mathbf{1}$ or $\mathbf{4}$, x contains an estimate of the solution.

3:     $\mathrm{iter}$ – int64int32nag_int scalar

The total number of iterations performed.

4:     $\mathrm{obj}$ – double scalar

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

5:     $\mathrm{ax}\left(\max \left(1,\mathbf{nclin}\right)\right)$ – double array

The final values of the linear constraints $Ax$.

If $\mathbf{nclin}=0$, ax is not referenced.

6:     $\mathrm{clamda}\left(\mathbf{n}+\mathbf{nclin}\right)$ – double array

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

7:     $\mathrm{user}$ – Any MATLAB object

8:     $\mathrm{lwsav}\left(120\right)$ – logical array

9:     $\mathrm{iwsav}\left(610\right)$ – int64int32nag_int array

10:   $\mathrm{rwsav}\left(475\right)$ – double array

11:   $\mathrm{ifail}$ – int64int32nag_int scalar

$\mathbf{ifail}=\mathbf{0}$ unless the function detects an error (see Error Indicators and Warnings).

nag_opt_qp_dense_solve (e04nf) returns with $\mathbf{ifail}=\mathbf{0}$ if $x$ is a strong local minimizer, i.e., the reduced gradient (Norm Gz; see Printed output) is negligible, the Lagrange multipliers (Lagr Mult; see Printed output) are optimal and $H_R$ (the reduced Hessian of $f\left(x\right)$; see Definition of Search Direction) is positive semidefinite.

Error Indicators and Warnings

Cases prefixed with W are classified as warnings and do not generate an error of type NAG:error_n. See nag_issue_warnings.

W  $\mathbf{ifail}=1$

The iterations were terminated at a dead-point. The necessary conditions for optimality are satisfied but the sufficient conditions are not. (The reduced gradient is negligible, the Lagrange multipliers are optimal, but $H_R$ is singular or there are some very small multipliers.) If ${\nabla }^{2}f\left(x\right)$ is not positive definite, $x$ is not necessarily a local solution of the problem and verification of optimality requires further information. If ${\nabla }^{2}f\left(x\right)$ is positive semidefinite or the problem is of type LP, $x$ gives the global minimum value of the objective function, but the final $x$ is not unique.

W  $\mathbf{ifail}=2$

The solution appears to be unbounded, i.e., the objective function is not bounded below in the feasible region. This value of ifail occurs if a step larger than Infinite Step Size ($\text{default value}={10}^{20}$) would have to be taken in order to continue the algorithm, or the next step would result in an element of $x$ having magnitude larger than Infinite Bound Size ($\text{default value}={10}^{20}$).

W  $\mathbf{ifail}=3$

No feasible point was found, i.e., it was not possible to satisfy all the constraints to within the feasibility tolerance. In this case, the constraint violations at the final $x$ will reveal a value of the tolerance for which a feasible point will exist – for example, when the feasibility tolerance for each violated constraint exceeds its Slack (see Printed output) at the final point. The modified problem (with an altered feasibility tolerance) may then be solved using a Warm Start. You should check that there are no constraint redundancies. If the data for the constraints are accurate only to the absolute precision $\sigma$, you should ensure that the value of the optional parameter Feasibility Tolerance ($\text{default value}=\sqrt{\epsilon }$, where $\epsilon$ is the machine precision) is greater than $\sigma$. For example, if all elements of $A$ are of order unity and are accurate only to three decimal places, the Feasibility Tolerance should be at least ${10}^{-3}$.

   $\mathbf{ifail}=4$

The limiting number of iterations was reached before normal termination occurred.

The values of the optional parameters Feasibility Phase Iteration Limit ($\text{default value}=\max \left(50,5\left(n+m_L\right)\right)$) and Optimality Phase Iteration Limit ($\text{default value}=\max \left(50,5\left(n+m_L\right)\right)$) may be too small. If the method appears to be making progress (e.g., the objective function is being satisfactorily reduced), either increase the iterations limit and rerun nag_opt_qp_dense_solve (e04nf) or, alternatively, rerun nag_opt_qp_dense_solve (e04nf) using the Warm Start facility to specify the initial working set.

   $\mathbf{ifail}=5$

The reduced Hessian exceeds its assigned dimension. The algorithm needed to expand the reduced Hessian when it was already at its maximum dimension, as specified by the optional parameter Maximum Degrees of Freedom ($\text{default value}=n$).

The value of the optional parameter Maximum Degrees of Freedom is too small. Rerun nag_opt_qp_dense_solve (e04nf) with a larger value (possibly using the Warm Start facility to specify the initial working set).

   $\mathbf{ifail}=6$

An input argument is invalid.

   $\mathbf{ifail}=7$

The designated problem type was not FP, LP, QP1, QP2, QP3 or QP4. Rerun nag_opt_qp_dense_solve (e04nf) with the optional parameter Problem Type set to one of these values.

   $\mathbf{\text{Overflow}}$

If the printed output before the overflow error contains a warning about serious ill-conditioning in the working set when adding the $j$th constraint, it may be possible to avoid the difficulty by increasing the magnitude of the Feasibility Tolerance ($\text{default value}=\sqrt{\epsilon }$, where $\epsilon$ is the machine precision) and rerunning the program. If the message recurs even after this change, the offending linearly dependent constraint (with index ‘$j$’) must be removed from the problem.

   $\mathbf{ifail}=-99$

An unexpected error has been triggered by this routine. Please contact NAG.

   $\mathbf{ifail}=-399$

Your licence key may have expired or may not have been installed correctly.

   $\mathbf{ifail}=-999$

Dynamic memory allocation failed.

Accuracy

Further Comments

Scaling

Description of the Printed Output

Example

Open in the MATLAB editor: e04nf_example

function e04nf_example


fprintf('e04nf example results\n\n');

% Problem specification
n     = 7;
m     = 7;
nclin = 7;
a     = [ 1,     1,       1,       1,       1,      1,       1;
          0.15,  0.04,    0.02,    0.04,    0.02,   0.01,    0.03;
          0.03,  0.05,    0.08,    0.02,    0.06,   0.01,    0;
          0.02,  0.04,    0.01,    0.02,    0.02,   0.00,    0;
          0.02,  0.03,    0.00,    0.00,    0.01,   0.00,    0;
          0.70,  0.75,    0.80,    0.75,    0.80,   0.97,    0;
          0.02,  0.06,    0.08,    0.12,    0.02,   0.01,    0.97];
bl    = [-0.01; -0.10;   -0.01;   -0.04;   -0.10;  -0.01;   -0.01;...
         -0.13; -1e25;   -1e25;   -1e25;   -1e25;  -0.0992; -0.003];
bu    = [ 0.01;  0.15;    0.03;    0.02;    0.05;   1e25;    1e25;...
         -0.13; -0.0049; -0.0064; -0.0037; -0.0012; 1e25;    0.002];
cvec  = [-0.02; -0.2;    -0.2;    -0.2;    -0.2;    0.04;    0.04];
x     = [-0.01; -0.03;    0;      -0.01;   -0.1;    0.02;    0.01];

h     = zeros(m,n);
h(1,  1)   = 2;
h(2,  2)   = 2;
h(3:4,3:4) = 2;
h(5,  5)   = 2;
h(6:7,6:7) = -2;

% Initialize and set options
istate     = zeros(14, 1, 'int64');
[cwsav,lwsav,iwsav,rwsav,ifail] = e04wb('e04nf');
[lwsav, iwsav, rwsav, inform] = e04nh(...
                                      'NoList', lwsav, iwsav, rwsav);
[lwsav, iwsav, rwsav, inform] = e04nh(...
                                      'Print Level = 0', lwsav, iwsav, rwsav);

% Solve
[istate, x, iter, obj, ax, clamda, user, lwsav, iwsav, rwsav, ifail] = ...
    e04nf(...
          a, bl, bu, cvec, h, @qphess, istate, x, lwsav, iwsav, rwsav);
fprintf('Minimum value     :  %9.4f\n\n',obj);
fprintf('Found after %3d iterations at x:\n   ',iter);
fprintf(' %9.4f',x);
fprintf('\nLinear contrained values Ax:\n   ');
fprintf(' %9.4f',ax);
fprintf('\n');



function [hx, user, iwsav] = qphess(n, jthcol, h, ldh, x, user, iwsav)

  if (iwsav(365) == 3 || iwsav(365) == 4)
    % Problem type QP1 and QP2
    hx = h*x;
  elseif(iwsav(365) == 5 || iwsav(365) == 6)
    hx = transpose(h)*h*x;
  else
    hx = zeros(n,1);
  end

e04nf example results

Minimum value     :     0.0370

Found after   7 iterations at x:
      -0.0100   -0.0699    0.0183   -0.0243   -0.0620    0.0138    0.0041
Linear contrained values Ax:
      -0.1300   -0.0059   -0.0064   -0.0045   -0.0029   -0.0992   -0.0030

Algorithmic Details

Overview

Definition of Search Direction

Main Iteration

Choosing the Initial Working Set

Optional Parameters

• Check Frequency
• Cold Start
• Crash Tolerance
• Defaults
• Expand Frequency
• Feasibility Phase Iteration Limit
• Feasibility Tolerance
• Hessian Rows
• Infinite Bound Size
• Infinite Step Size
• Iteration Limit
• Iters
• Itns
• List
• Maximum Degrees of Freedom
• Minimum Sum of Infeasibilities
• Monitoring File
• Nolist
• Optimality Phase Iteration Limit
• Optimality Tolerance
• Print Level
• Problem Type
• Rank Tolerance
• Warm Start

[lwsav, iwsav, rwsav, inform] = e04nh('Print Level = 5', lwsav, iwsav, rwsav);

Description of the Optional Parameters

• the keywords, where the minimum abbreviation of each keyword is underlined (if no characters of an optional qualifier are underlined, the qualifier may be omitted);
• a parameter value, where the letters $a$, $i$ and $r$ denote options that take character, integer and real values respectively;
• the default value, where the symbol $\epsilon$ is a generic notation for machine precision (see nag_machine_precision (x02aj)).

Description of Monitoring Information