Namespace: nagcpp

Namespace: opt

Function: handle_solve_bounds_foas

NAG CPP Interface
nagcpp::opt::handle_solve_bounds_foas (e04kf)

Note: this function uses optional parameters to define choices in the problem specification and in the details of the algorithm. If you wish to use default settings for all of the optional parameters, you need only read Sections 1 to 10 of this document. If, however, you wish to reset some or all of the settings please refer to Section 11 for a detailed description of the algorithm and to Section 12 for a detailed description of the specification of the optional parameters.

Keyword Search:

NAG Library Manual, Mark 29.2

Interfaces: FL CL CPP AD PY MB

NAG CPP Interface Introduction

E04 (Opt) Chapter Contents

e04kf: FL CL CPP AD PY MB

Namespace: nagcpp

Namespace: opt

Function: handle_solve_bounds_foas

▸▿ Contents

Settings help

CPP Name Style:
NAG (a00aa)
Short (impl_details)
Long (info::impl_details)
Full (nagcpp::info::impl_details)

1 Purpose

handle_solve_bounds_foas is a solver from the NAG optimization modelling suite for bound-constrained large-scale Nonlinear Programming (NLP) problems. It is a first-order active-set method (FOAS) that has low memory requirements and thus is suitable for very large-scale problems.

2 Specification

#include "e04/nagcpp_e04kf.hpp"
#include "e04/nagcpp_class_CommE04RA.hpp"

template <typename COMM, typename OBJFUN, typename OBJGRD, typename MONIT, typename X, typename RINFO, typename STATS>

void function handle_solve_bounds_foas(COMM &comm, OBJFUN &&objfun, OBJGRD &&objgrd, MONIT &&monit, X &&x, RINFO &&rinfo, STATS &&stats, OptionalE04KF opt)

template <typename COMM, typename OBJFUN, typename OBJGRD, typename MONIT, typename X, typename RINFO, typename STATS>

void function handle_solve_bounds_foas(COMM &comm, OBJFUN &&objfun, OBJGRD &&objgrd, MONIT &&monit, X &&x, RINFO &&rinfo, STATS &&stats)

3 Description

handle_solve_bounds_foas solves large-scale bound-constrained nonlinear optimization problems of the form

\underset{x \in R^{n}}{minimize} f (x) subject to ℓ_{x} \leq x \leq u_{x},

(1)

where

n

is the number of decision variables,

ℓ_{x}

u_{x}

, and

x

are

n

-dimensional vectors, and the nonlinear objective function

f (x)

is assumed to be sufficiently smooth.

The solver is a first-order method (i.e., uses only first derivatives) that has very low memory requirements and, therefore, is suitable for very large bound-constrained problems. It is based on an active-set method coupled to a nonmonotone projected gradient algorithm (NPG), nonlinear conjugate gradient method (CG) and its limited-memory variant (LCG). The active-set method is based on alternating between both solvers, the NPG step handles the box constraints and identifies a suitable search space while the CG step explores it for a solution.

For a detailed description of the algorithm see Section 11. Under standard assumptions on the problem (smoothness of the first derivative of the objective) the algorithm converges to a local solution or to a critical point.

handle_solve_bounds_foas serves as a solver for problems stored as a handle. The handle points to an internal data structure which defines the problem and serves as a means of communication for functions in the NAG optimization modelling suite. First, the problem handle is initialized by calling handle_init. Then some of the functions handle_set_linobj, handle_set_quadobj, handle_set_nlnobj and handle_set_simplebounds may be called to formulate the objective and to define (simple) box constraints for the problem. Once the problem is fully described, the handle may be passed to the solver handle_solve_bounds_foas. When the handle is no longer needed, handle_free should be called to destroy it and deallocate the memory held within. See Section 3.1 in the E04 Chapter Introduction for more details about the NAG optimization modelling suite.

The algorithm behaviour can be modified by various optional parameters (see Section 12) which can be set by handle_opt_set and e04zpf (no CPP interface) anytime between the initialization of the handle by handle_init and a call to the solver. Once the solver has finished, options may be modified for the next solve. The solver may be called repeatedly with various starting points and/or optional parameters. Option getter handle_opt_get can be called to retrieve the current value of any option.

The optional parameter Task may be used to switch the problem to maximization, while FOAS Estimate Derivatives can be used to complete missing elements from the gradient. Optional parameter Verify Derivatives may help verify the correctness of the gradient vector before starting to solve a problem.

Several options may have significant impact on the performance of the solver. Even if the defaults were chosen to suit the majority of problems, it is recommended that you experiment in order to find the most suitable set of options for a particular problem, see Sections 11 and 12 for further details.

4 References

Dai Y-H and Kou C-X (2013) A Nonlinear Conjugate Gradient Algorithm with an Optimal Property and an Improved Wolfe Line Search SIAM J. Optim. 23(1) 296–320

Gill P E and Leonard M W (2003) Limited-Memory Reduced-Hessian Methods for Large-Scale Unconstrained Optimization SIAM J. Optim. 14(2) 380–401

Hager W W and Zhang H (2005) A New Conjugate Gradient Method with Guaranteed Descent and an Efficient Line Search SIAM J. Optim. 16(1) 170–192

Hager W W and Zhang H (2006a) Algorithm 851: CG DESCENT, a Conjugate Gradient Method with Guaranteed Descent ACM Trans. Math. Software 32(1) 113–137

Hager W W and Zhang H (2006b) A New Active Set Algorithm for Box Constrained Optimization SIAM J. Optim. 17(2) 525–557

Hager W W and Zhang H (2013) The Limited Memory Conjugate Gradient Method SIAM J. Optim. 23(4) 2150–2168

Nocedal J and Wright S J (2006) Numerical Optimization (2nd Edition) Springer Series in Operations Research, Springer, New York

5 Arguments

1: $comm$ – CommE04RA Input/Output

Communication structure. An object of either the derived class CommE04RA or its base class NoneCopyableComm can be supplied. It is recommended that the derived class is used. If the base class is supplied it must first be initialized via a call to opt::handle_init (e04ra).

2: $objfun$ – void function Function

objfun must calculate the value of the nonlinear objective function

f (x)

at a specified point

x

. If there is no nonlinear objective (e.g., handle_set_linobj and handle_set_quadobj was called to define a linear or quadratic objective function), objfun will never be called by handle_solve_bounds_foas and may be the dummy function bounds_dummy_objfun (included in the NAG Library.)

void function objfun(const utility::array1D<double,data_handling::ArgIntent::IntentIN> &x, double &fx, types::f77_integer &inform)

1: $x (nvar)$ – double array Input: On entry: the vector $x$ of variable values at which the objective function is to be evaluated.
2: $fx$ – double Output: On exit: the value of the objective function at $x$ .
3: $inform$ – types::f77_integer Input/Output: On entry: a non-negative value.
In some cases, it is known beforehand that the evaluations of the objective function and its gradient are required at the same point $x$ , in such cases, $inform = 1$ . This may help to optimize your code in order to avoid recalculations of common quantities when evaluating both the objective function and gradient; the objective function is always evaluated before the objective gradient. This notification parameter may be safely ignored if such optimization is not required.

On exit: may be used to indicate that the function cannot be evaluated at the requested point $x$ by setting $inform < 0$ . Returning NaN or $\pm \infty$ in fx has the same effect. The algorithm will try to recover if the objective cannot be evaluated; if recovery is not possible it will stop with $errorid = 25$ .
4: $nvar$ – types::f77_integer Input: On entry: $n$ , the current number of decision variables $x$ in the model.

3: $objgrd$ – void function Function

objgrd must calculate the values of the nonlinear objective function gradient

\frac{\partial f}{\partial x}

at a specified point

x

. Every call to objgrd is preceded by a call to objfun at the same point, if this is known in advance, both functions will be notified via

inform = 1

. If there is no nonlinear objective (e.g., handle_set_linobj and handle_set_quadobj was called to define a linear or quadratic objective function), objgrd will never be called by handle_solve_bounds_foas and objgrd may be the dummy function bounds_dummy_objgrd (included in the NAG Library.)

If the optional parameter

FOAS Estimate Derivatives = YES

, then after returning from objgrd the gradient vector is checked for missing entries which you have not supplied. Missing entries are estimated using the finite difference method, see optional parameter FOAS Estimate Derivatives description for more details.

void function objgrd(const utility::array1D<double,data_handling::ArgIntent::IntentIN> &x, utility::array1D<double,data_handling::ArgIntent::IntentINOUT> &fdx, types::f77_integer &inform)

1: $x (nvar)$ – double array Input: On entry: the vector $x$ of variable values at which the objective function gradient is to be evaluated.
2: $fdx (nnzfd)$ – double array Input/Output: On entry: the elements should only be assigned and not referenced.

On exit: the values of the nonzero elements in the sparse gradient vector of the objective function, in the order specified by idxfd in a previous call to handle_set_nlnobj. $fdx (i - 1)$ will store the gradient element $\frac{\partial f}{\partial x_{j}}$ , where $j = idxfd (i - 1)$ .
3: $inform$ – types::f77_integer Input/Output: On entry: a non-negative value.
If $inform = 1$ , then the previous call to objfun was also made at the same point $x$ with $inform = 1$ . This may help to optimize your code in order to avoid recalculations of common quantities when evaluating both the objective function and gradient. This notification parameter may be safely ignored if such optimization is not required.

On exit: may be used to inform that the gradient cannot be evaluated at the requested point $x$ by setting $inform < 0$ . Returning NaN or $\pm \infty$ in any element of fdx has the same effect. The algorithm will try to recover if the gradient cannot be evaluated; if recovery is not possible it will stop with $errorid = 25$ .
4: $nvar$ – types::f77_integer Input: On entry: $n$ , the current number of decision variables $x$ in the model.
5: $nnzfd$ – types::f77_integer Input: On entry: the number of nonzero elements in the sparse gradient vector of the objective function, as was set in a previous call to handle_set_nlnobj.

4: $monit$ – void function Function

monit is provided to enable you to monitor the progress of the optimization and optionally to terminate the solver early if necessary. It is invoked at the end of every

i

th iteration where

i

is given by the FOAS Monitor Frequency (the default is

0

, monit is not called).

monit may be the dummy function bounds_dummy_monit (included in the NAG Library).

void function monit(const utility::array1D<double,data_handling::ArgIntent::IntentIN> &x, const utility::array1D<double,data_handling::ArgIntent::IntentIN> &rinfo, const utility::array1D<double,data_handling::ArgIntent::IntentIN> &stats)

1: $x (nvar)$ – double array Input: On entry: the vector $x$ of decision variables at the current iteration.
2: $rinfo (100)$ – double array Input: On entry: error measures and various indicators at the end of the current iteration as described in rinfo.
3: $stats (100)$ – double array Input: On entry: solver statistics at the end of the current iteration as described in stats.
4: $nvar$ – types::f77_integer Input: On entry: $n$ , the current number of decision variables $x$ in the model.

5: $x (nvar)$ – double array Input/Output

On entry:

x_{0}

, the initial estimates of the variables,

x

On exit: the final values of the variables,

x

6: $rinfo (100)$ – double array Output

On exit: error measures and various indicators at the end of the final iteration as given in the table below:

$0$	Objective function value $f (x)$ .
$1$	Norm of inactive gradient, the objective gradient over the current search space. If the problem is unconstrained, then elements $1$ – $3$ coincide. See Section 11.4 for details.
$2$	Norm of projected direction, used in (3), see Section 11.2.
$3$	Norm of objective gradient.
$4$	Last step size ( $α_{k}$ ) used in (2) and (5), see Section 11.2 and Section 11.3.
$5$	Progress score, a positive value that measures the progress of the solver. A low score close to zero indicates poor progress, see (8) in Section 11.4.
$6$ – $99$	Reserved for future use.

7: $stats (100)$ – double array Output

On exit: solver statistics at the end of the final iteration as given in the table below:

$0$	Number of function evaluations performed by NPG, see Section 11.2.
$1$	Number of gradient evaluations performed by NPG.
$2$	Number of function evaluations performed by CG, see Section 11.3.
$3$	Number of gradient evaluations performed by CG.
$4$	Number of function evaluations performed by LCG.
$5$	Number of gradient evaluations performed by LCG.
$6$	Number of function evaluations used by the finite difference method to estimate missing components of the gradient.
$7$	Number of iterations.
$8$	Total time spent in the solver (including user-supplied function calls).
$9$	Total time spent in user-supplied objective function.
$10$	Total time spent in user-supplied objective gradient.
$11$ – $99$	Reserved for future use.

8: $opt$ – OptionalE04KF Input/Output

Optional parameter container, derived from Optional.

5.1Additional Quantities

1: $nvar$: $n$ , the current number of decision variables $x$ in the model.

6 Exceptions and Warnings

Errors or warnings detected by the function:

Note: in some cases handle_solve_bounds_foas may return useful information.

All errors and warnings have an associated numeric error code field, errorid, stored either as a member of the thrown exception object (see errorid), or as a member of opt.ifail, depending on how errors and warnings are being handled (see Error Handling for more details).

Raises: ErrorException

$errorid = 1$: comm::handle has not been initialized.

$errorid = 1$: comm::handle does not belong to the NAG optimization modelling suite,
has not been initialized properly or is corrupted.

$errorid = 1$: comm::handle has not been initialized properly or is corrupted.

$errorid = 2$: This solver does not support the model defined in the handle.

$errorid = 2$: The problem is already being solved.

$errorid = 4$: On entry, $nvar = ⟨ value ⟩$ ,
expected $value = ⟨ value ⟩$ .
Constraint: $nvar$ must match the current number of variables
of the model in the comm::handle.

$errorid = 7$: Please provide a proper objfun function.

$errorid = 7$: Please provide a proper objgrd function.

$errorid = 21$: The current starting point is unusable.

$errorid = 26$: User-provided gradient is likely to be incorrect.

$errorid = 10601$: On entry, argument $⟨ value ⟩$ must be a vector of size $⟨ value ⟩$ array.
Supplied argument has $⟨ value ⟩$ dimensions.

$errorid = 10601$: On entry, argument $⟨ value ⟩$ must be a vector of size $⟨ value ⟩$ array.
Supplied argument was a vector of size $⟨ value ⟩$ .

$errorid = 10601$: On entry, argument $⟨ value ⟩$ must be a vector of size $⟨ value ⟩$ array.
The size for the supplied array could not be ascertained.

$errorid = 10602$: On entry, the raw data component of $⟨ value ⟩$ is null.

$errorid = 10603$: On entry, unable to ascertain a value for $⟨ value ⟩$ .

$errorid = 10605$: On entry, the communication class $⟨ value ⟩$ has not been initialized correctly.

$errorid = 10703$: An exception was thrown during IO (writing).

$errorid = −99$: An unexpected error has been triggered by this routine.

$errorid = −399$: Your licence key may have expired or may not have been installed correctly.

$errorid = −999$: Dynamic memory allocation failed.

Raises: CallbackEarlyTermination

$errorid = 20$: User requested termination during a monitoring step.

Raises: WarningException

$errorid = 22$: Maximum number of iterations reached.

$errorid = 23$: The solver terminated after the maximum time allowed was exhausted.

$errorid = 24$: The solver was terminated because no further progress could be achieved.

$errorid = 25$: Invalid number detected in user-supplied function and recovery failed.

$errorid = 25$: Invalid number detected in user-supplied gradient and recovery failed.

$errorid = 50$: Problem was solved to an acceptable level; full accuracy was not achieved.

$errorid = 54$: The problem seems to be unbounded and the algorithm was stopped.

Raises: CallbackException

$errorid = 10701$: An exception was thrown in a callback.

$errorid = 10702$: The memory address for an array in a callback has changed.

7 Accuracy

The accuracy of the solution is determined by optional parameters FOAS Stop Tolerance and FOAS Rel Stop Tolerance. If

errorid = 0

on exit, the returned point satisfies (6) or (7) to the defined accuracy.

Please refer to Section 11.4 and the description of the particular options in Section 12.

8 Parallelism and Performance

Please see the description for the underlying computational routine in this section of the FL Interface documentation.

9 Further Comments

9.1 Description of the Printed Output

The solver can print information to give an overview of the problem and of the progress of the computation. The output may be sent to two independent streams (files) which are set by optional parameters Print File and Monitoring File. Optional parameters Print Level, Monitoring Level, Print Solution and Print Options determine the exposed level of detail. This allows, for example, a detailed log file to be generated while the condensed information is displayed on the screen.

By default (

Print File = 6

Print Level = 2

), the following sections are printed to the standard output:

Header
Optional parameters list (if $Print Options = YES$ )
Problem statistics
Verification of derivatives (if $Verify Derivatives = YES$ )
Iteration log
Summary
Solution (if $Print Solution = YES$ )

Header

The header is a message indicating the start of the solver. It should look like:

-------------------------------------------------------------------------------
 E04KF, First-order method for bound-constrained problems
-------------------------------------------------------------------------------

Optional parameters list

The list shows all options of the solver, each displayed on one line. The output contains the option name, its current value and an indicator for how it was set. The options unchanged from the default setting are noted by ‘d’, options you set are noted by ‘U’, and options reset by the solver are noted by ‘S’. Note that the output format is compatible with the file format expected by e04zpf (no CPP interface). The output might look as follows:

Begin of Options
    Print File                    =                   6     * d
    Print Level                   =                   2     * U
    Print Options                 =                 Yes     * d
    Print Solution                =                 All     * U
    Monitoring File               =                   9     * U
    Monitoring Level              =                   3     * U
    Foas Monitor Frequency        =                   0     * d
    Foas Print Frequency          =                   5     * U
End of Options

Problem statistics

Print Level \geq 2

, statistics on the problem is printed, for example:

 Problem Statistics
   No of variables                  8 (+1 disabled, +2 fixed)
     free (unconstrained)           3
     bounded                        2
   Objective function       Nonlinear

Verification of derivatives

Verify Derivatives = YES

, then the user-supplied nonlinear objective gradient is verified. If

Print Level < 5

, then it will only check those entries related to the actual search space and compare to a finite-difference approximation, it might look as follows:

 Derivative checker
     idx test        fdx     approx  rel error
       3 Fail   3.56E+02   3.57E+02   2.80E-03
       4  Ok    1.80E+02   1.80E+02   6.05E-09
       5  Ok    0.00E+00   0.00E+00   0.00E+00

On the other hand, if

Print Level = 5

, it will print one line for each gradient entry and the output should be similar to:

 Derivative checker
     idx test        fdx     approx  rel error
       3 Fail   3.56E+02   3.57E+02   2.80E-03
       4  Ok    1.80E+02   1.80E+02   6.05E-09
       5  Ok    0.00E+00   0.00E+00   0.00E+00
       7 Skip         FD         --         --
       2 Skip    BLX=BUX         --         --
       6 Skip    BLX=BUX         --         --
       1 Skip   disabled         --         --

Iteration log

Print Level \geq 2

, the solver prints the status of every

k

th iteration, specified using the optional parameter

FOAS Print Frequency = k

, with default value of

k = 1

. It will also print status information regarding the switch between NPG, CG and LCG. It will also print status information regarding the switch between NPG, CG and LCG.

Print Level = 2

, the output shows the iteration number (

0

represents the starting point), the current objective value, and optimality measures (norm of inactive gradient and the norm of the projected direction). Note that all these values for the last iteration are also available in rinfo. The output might look as follows:

-------------------------------------------------------------------------------
  iters |  objective |  optim  |   dir
-------------------------------------------------------------------------------
      20  7.71339E-01  1.22E-01  1.70E+00
      25  7.06709E-01  1.75E+00  1.54E+00
      30  7.06709E-01  1.75E+00  1.54E+00
      35  6.82989E-01  1.35E+00  1.35E+00
      30  6.55834E-01  3.26E+00  1.67E+00

Print Level = 3

, the solver also prints for each iteration the progress score, iteration type, size or type of performed step, accumulated number of objective function and gradient calls. The output takes the following form:

-------------------------------------------------------------------------------
  iters |  objective |  optim  |   dir   | progrss | it|   step  |    nf|    ng
-------------------------------------------------------------------------------
       0  1.10354E+02  0.00E+00  6.81E+01  1.00E+00       Start        1      1
       1  3.29016E+01  1.45E+01  1.45E+01  5.36E+00 NPG  1.39E-02      8      2
       2  3.29016E+01  1.45E+01  1.45E+01  5.36E+00     Switch CG      8      2
       3  1.84090E+01  2.48E+01  2.48E+01  2.85E+00  CG  2.42E-02     11      4
       4  1.84090E+01  2.48E+01  2.48E+01  2.85E+00      Restart      11      4
       5  9.48087E+00  1.25E+01  1.25E+01  1.45E+00  CG  1.49E-02     13      5
       6  4.20527E+00  8.28E+00  8.22E-01  2.51E+00  CG  3.61E-02     15      6
       7  7.20453E-01  3.08E+00  1.62E+00  9.70E+00 LCG  7.21E-01     17      7

Print Level > 3

, each iteration produces more information that expands over several lines. This additional information can contain:

The sizes of gradient in full space and subspace;
The quadratic model prediction error and counter;
The amount of search direction restarts;
The extreme eigenvalues of the basis matrix for the subspace;
Information and size of the active-set;
Whether the orthogonal property holds;
The action to take for the next iteration (switch solver);
The search space size;
The amount of memory vectors used.

The output might look as follows:

-----  Iteration details (CG)  ---------
Size of gradient                2.64E+01
Size of subspace gradient       2.64E+01
Eigenvalues R range   6.80E-01  6.80E-01
Orthogonal property                 lost
Action                     Switch to LCG
----------------------------------------
-----  LCG solver details  -------------
Search space has            2 element(s)
Memory vectors                   2 (  2)
----------------------------------------

Summary

Once the solver finishes, a detailed summary is produced:

-------------------------------------------------------------------------------
Status: converged, an optimal solution was found
-------------------------------------------------------------------------------
Value of the objective             3.55353E-14
Norm of inactive gradient          2.53812E-07
Norm of projected direction        1.68481E-07
Iterations                                  32
Function evaluations                        80
FD func. evaluations                         5
Gradient evaluations                        48
  NPG function calls                         0
  NPG gradient calls                         0
  CG function calls                         11
  CG gradient calls                          9
  LCG function calls                        69
  LCG gradient calls                        39
-------------------------------------------------------------------------------

It starts with a status line of the overall result, followed by the final objective value as well as the gradient (unconstrained problem) or inactive gradient and projected direction (constrained problem) norms. If

Print Level > 1

Monitoring Level > 1

, it will additionally report iteration count, objective function and gradient call information.

Optionally, if

Stats Time = YES

, the timings for the user-supplied objective are displayed. It might look as follows:

Timing
  Total time spent                       12.54 sec
  Total time in obj function              3.01 sec ( 24.0%)
  Total time in obj gradient              5.28 sec ( 42.1%)

Solution

Print Solution = YES

, the values of the primal variables are printed, furthermore if the problem is constrained, also the dual variables (see Lagrangian Multipliers) and their bounds are reported. It might look as follows:

 Primal variables:
   idx   Lower bound       Value       Upper bound
     1     Disabled         NaN          Disabled
     2   5.00000E+00    5.00000E+00    5.00000E+00
     3  -1.00000E+00    8.00000E-01    8.00000E-01
     4  -2.00000E+00    6.40000E-01    2.00000E+00
     5   2.00000E+00    2.00000E+00    9.00000E+00
     6  -1.00000E+00   -1.00000E+00   -1.00000E+00
     7   0.00000E+00    0.00000E+00    1.00000E+00
 
 Box bounds dual variables:
   idx   Lower bound       Value       Upper bound       Value
     1     Disabled         NaN          Disabled         NaN
     2   5.00000E+00    0.00000E+00    5.00000E+00    0.00000E+00
     3  -1.00000E+00    0.00000E+00    8.00000E-01    4.00000E-01
     4  -2.00000E+00    0.00000E+00    2.00000E+00    0.00000E+00
     5   2.00000E+00    0.00000E+00    9.00000E+00    0.00000E+00
     6  -1.00000E+00    0.00000E+00   -1.00000E+00    0.00000E+00
     7   0.00000E+00    0.00000E+00    1.00000E+00    0.00000E+00

10 Example

In this example, we minimize in

ℝ^{2}

the Rosenbrock function over simple bounds on the variables. The problem to solve is

\begin{array}{rcl} \underset{x \in ℝ^{2}}{minimize} f (x) = {(1 - x_{1})}^{2} + 100 {(x_{2} - x_{1}^{2})}^{2} & subject to & −1 \leq x_{1} \leq 0.8, \\ −2 \leq x_{2} \leq 2 . \end{array}

The initial guess is

x_{0} = {(- 1.5, 1.9)}^{T}

, and the expected solution point is

x^{*} = {(0.8, 0.64)}^{T}

with objective value

f (x^{*}) = 0.04

10.1 Example Program

Source File	Data	Results
ex_e04kf.cpp	None	ex_e04kf.r

11 Algorithmic Details

This section contains the description of the underlying algorithms used in handle_solve_bounds_foas, a first-order active-set (FOAS) method with very low memory requirements suitable for large-scale bound-constrained nonlinear optimization problems. For further details, see Hager and Zhang (2006b) and references therein.

11.1 Active-Set Method and Algorithm Outline

Active-set is a useful method to tackle problems with bound constraints. It derives its name from the partitioning of the search space: elements of

x

that are fixed at a bound (

x_{i} = ℓ_{i}

x_{i} = u_{i}

), are said to be active, while the rest of the elements of

x

are said to be inactive. The goal of this method is to estimate which variables will be active at the solution point while optimizing over the inactive components.

handle_solve_bounds_foas consists of a constrained solver, an unconstrained solver and a set of rules to switch between them. The constrained solver (nonmonotone projected gradient, NPG) has two purposes: while solving (1) it tries to guess which components of

x

are active at a solution. Once a reasonable guess is available, control is transferred to an unconstrained solver (conjugate gradient method, CG, and its limited-memory variant, LCG) that operates only over the inactive elements of

x

but convergence is much faster that NPG.

The following is an outline of the implemented algorithm.

FOAS Algorithm

(i)Make initial point $x_{0}$ feasible.
(ii)Loop until stopping criteria are satisfied
1. (a)Constrained problem (NPG): start solving constrained problem (1) and try to identify a suitable inactive search space; if found, then go to (ii)(b).
2. (b)Unconstrained subproblem (CG): solve unconstrained problem (4) over the elements of $x$ marked inactive. If elements of $x$ become newly active or it is deemed that some active elements should be explored, then go back to (ii)(a).

11.2 Constrained Subproblem (NPG)

By setting

Ω = {x \in ℝ^{n} : ℓ_{x} \leq x \leq u_{x}}

, the constrained problem (1) can be stated as

\min_{x \in Ω} f (x)

and can be solved using NPG which is an iterative method of the form

x_{k + 1} = x_{k} + α_{k} d_{Ω} (x_{k}),

(2)

with step size

α_{k}

and projection direction

d_{Ω} (x)

The projected direction,

d_{Ω} (x)

, is obtained by a projection of the objective gradient

g (x)

onto

Ω

d_{Ω} (x) = P_{Ω} (x - g (x)) - x,

(3)

with

P_{Ω}

the Euclidean projection operator over the box

Ω

. Note that if

Ω = ℝ^{n}

, then

d_{Ω} (x_{k}) = - g (x_{k})

and the method reduces to nonmonotone steepest descent.

The step size

α_{k} > 0

is chosen to guarantee global convergence, i.e., guarantee sufficient progress at each iteration using the nonmonotone Armijo condition

f (x_{k + 1}) \leq f_{k}^{R} + c α_{k} {g (x_{k})}^{T} d_{Ω} (x_{k}),   with ​ 0 < c < 1,

where

f_{k}^{R}

is a reference objective function value for iteration

k

, possibly

f_{k}^{R} = f (x_{k})

. The step size is estimated using a nonmonotone backtracking technique. As soon as a suitable inactive set is identified (see Hager and Zhang (2006b)), the control is transferred to CG (LCG) to accelerate the convergence.

11.3 Unconstrained Subproblem (CG)

The exploration of the inactive search space is performed using the CG method which is an iterative method to solve

\min_{x \in ℝ^{n}} f (x)

(4)

and each iterate is updated using the expression

x_{k + 1} = x_{k} + α_{k} d_{k},

(5)

where

d_{k}

is a direction of descent that combines the current gradient

g (x)

and a scaled version of the previous direction,

β_{k} d_{k - 1}

, with

β_{k} \in ℝ

d_{k + 1} = - g (x_{k + 1}) + β_{k} d_{k}, d_{0} = - g (x_{0}) .

The step size

α_{k} > 0

is chosen to ensure global convergence and is achieved by satisfying the weak Wolfe conditions

f (x_{k + 1}) - f (x_{k}) \leq δ α_{k} {g (x_{k})}^{T} d_{k}, {g (x_{k + 1})}^{T} d_{k} \geq σ {g (x_{k})}^{T} d_{k},

with

0 < δ < σ < 1

or the approximate Wolfe conditions (see Hager and Zhang (2005))

ε ε

\begin{matrix} (2 δ - 1) {g (x_{k})}^{T} d_{k} \geq {g (x_{k + 1})}^{T} d_{k}, & {g (x_{k + 1})}^{T} d_{k} \geq σ {g (x_{k})}^{T} d_{k},   and \\ {g (x_{k + 1})}^{T} d_{k} \leq {g (x_{k})}^{T} d_{k} + ε_{k}, \end{matrix}

with

δ < \min (0.5, σ)

and

ε_{k} > 0

. The line search performed by handle_solve_bounds_foas will always satisfy either/or both stated conditions. The line search builds a series of nested intervals that contains a point satisfying the above sufficient decrease requirements.

It is important to note that when

Ω \neq ℝ^{n}

, then the used gradient is actually the inactive gradient, it is the gradient,

g (x_{k})

, but with zeros in the same positions where the elements of

x_{k}

are marked as active.

The performance of the CG method can degrade when orthogonality is lost between consecutive search directions. Therefore, handle_solve_bounds_foas uses a limited number of previous search directions to detect and restore orthogonality. When the current search direction is no longer orthogonal, it is discarded and a quasi-Newton variant known as limited-memory CG (LCG) is used to build a new search direction orthogonal to the explored subspace.

The available memory is used to build an approximation of the Hessian and the new search direction is estimated using

d_{k} = {- H}_{k} g (x_{k}) .

The quasi-Newton Hessian approximation

H_{k}

, is updated at each iteration using the BFGS update formula, see Nocedal and Wright (2006).

Experiments in literature indicate that an adequate range for the quasi-Newton (amount of vectors used to build

H_{k}

) lies between

7

and

20

. The default value for the maximum amount of vectors used in handle_solve_bounds_foas is

11

and can be changed via the FOAS Memory.

11.4 Stopping Criteria

A point is considered a solution when there are no feasible descent directions to use. Under this circumstance the function will stop, declaring to have found a solution.

If the problem is unconstrained, the function declares to have found a solution and stops when the first-order optimality condition is met within the defined absolute or relative tolerances (

ε_{tol} \geq 0

ε_{rel} \geq 0

{‖ g (x_{k}) ‖}_{p} \leq \max (ε_{tol}, ε_{rel} {‖ g (x_{0}) ‖}_{p}),

(6)

where

p

can designate either the Euclidean or Infinity (default) norms.

On the other hand, if the problem is constrained, the function characterizes to have found a solution when the first-order condition is satisfied for the projected direction,

d_{Ω}

, to the defined tolerances,

{‖ d_{Ω} (x_{k}) ‖}_{p} \leq \max (ε_{tol}, ε_{rel} {‖ d_{Ω} (x_{0}) ‖}_{p}) .

(7)

In the unconstrained case we have

Ω = ℝ^{n}

and both stopping criteria (6) and (7) coincide.

The stopping tolerances can be changed using the optional parameters FOAS Stop Tolerance and FOAS Rel Stop Tolerance, while the norm used can be set by using the optional parameter FOAS Tolerance Norm, see Section 12 for details. If these parameters are set too small in relation to the complexity and scaling of the problem, the function can terminate with

errorid = 22

24

50

Progress of the solver towards a solution is monitored using two criteria. The first one evaluates how poor the actual step has been and is estimated via

{- α}_{k} \frac{\nabla f {(x_{k})}^{T} d_{k}}{| f (x_{k}) |} < ε_{prog},

(8)

where

α_{k}

is given by the line search. If the above relation is consistently satisfied, then the solver stops with

errorid = 24

. The tolerance

ε_{prog}

is set using the optional parameter FOAS Progress Tolerance. The second criteria monitors the rate of convergence using

\frac{| f (x_{k - 1}) - f (x_{k}) |}{| f (x_{k}) |} < α_{k} ‖ \nabla f (x_{k}) ‖ ε_{slow} .

(9)

As with the first criteria, if the previous relation is deemed permanent, then the solver stops with

errorid = 50

. The tolerance

ε_{slow}

can be changed using the optional parameter FOAS Slow Tolerance.

11.5 A Note About Lagrangian Multipliers

It is often useful to have access to the Lagrangian multipliers (dual variables) associated with the constraints if there are any defined. In the case where only simple bounds are present, the multipliers directly relate to the values of the gradient at the solution. The multipliers of the active bounds are the absolute values of the associated elements of the gradient. The multipliers of the inactive bounds are always zero.

The multipliers based on the final gradient value (or its finite-difference approximation) can be retrieved by calling e04rxf (no CPP interface) with the command string

cmdstr =' Dual Variables'

. The format is the same as for other functions, see Section 3.1 in handle_solve_pennon (no CPP interface in the current release). Note that if the problem has not fully converged, the provided approximation might be quite crude.

12 Optional Parameters

Several optional parameters in handle_solve_bounds_foas define choices in the problem specification or the algorithm logic. In order to reduce the number of formal arguments of handle_solve_bounds_foas these optional parameters have associated default values that are appropriate for most problems. Therefore, you need only specify those optional parameters whose values are to be different from their default values.

The remainder of this section can be skipped if you wish to use the default values for all optional parameters.

The optional parameters can be changed by calling handle_opt_set anytime between the initialization of the handle and the call to the solver. Modification of the optional parameters during intermediate monitoring steps is not allowed. Once the solver finishes, the optional parameters can be altered again for the next solve.

The option values may be retrieved by handle_opt_get.

The following is a list of the optional parameters available. A full description of each optional parameter is provided in Section 12.1.

12.1 Description of the Optional Parameters

For each option, we give a summary line, a description of the optional parameter and details of constraints.

The summary line contains:

the keywords;
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 $ε$ is a generic notation for machine precision (see precision).

All options accept the value

DEFAULT

to return single options to their default states.

Keywords and character values are case and white space insensitive.

Defaults

This special keyword may be used to reset all optional parameters to their default values. Any value given with this keyword will be ignored.

FOAS Estimate Derivatives

a

Default

= NO

This option indicates whether to check for and estimate missing entries of the user-supplied gradient vector. Since the associated cost of estimating missing elements can be high, this option should only be used if strictly necessary. In general terms, if the gradient is not provided (in its entirety or for arbitrary points) potential degradation in the progress of the solver is to be expected. Depending on the complexity of the objective, the function may not achieve the desired optimality accuracy or even terminate with no possible further progress error

errorid = 24

, it is advisable to increase the values of FOAS Stop Tolerance and FOAS Rel Stop Tolerance when using this option.

Missing elements from the gradient vector are estimated by finite-differences using the perturbation interval specified by the optional parameter FOAS Finite Diff Interval.

FOAS Estimate Derivatives = NO

, the entries are not checked and all derivative elements need to be provided.

Constraint:

FOAS Estimate Derivatives = YES

NO

FOAS Finite Diff Interval

r

Default

= \sqrt{ε}

Specifies the relative perturbation size used to estimate a derivative using the forward (or backward) finite-difference method. Setting the value too small or too big may lead handle_solve_bounds_foas to terminated with

errorid = 24

25

Constraint:

10^{−12} \leq FOAS Finite Diff Interval \leq 10^{−1}

FOAS Iteration Limit

i

Default

= 10^{7}

This parameter sets the maximum number of iterations to be performed by handle_solve_bounds_foas. Setting the option too low might lead to

errorid = 22

Constraint:

FOAS Iteration Limit \geq 1

FOAS Memory

i

Default

= 11

This parameter specifies the maximum number of memory vectors to use in the LCG solver.

Constraint:

0 \leq FOAS Memory \leq 100

FOAS Monitor Frequency

i

Default

= 0

This parameter specifies the frequency on which to call the monitor function monit. If zero, the monitor function will not be called.

Constraint:

FOAS Monitor Frequency \geq 0

FOAS Print Frequency

i

Default

= 1

This parameter specifies the frequency with which to print information regarding each iteration to Print File and/or Monitoring File. By default, it will print information of every iteration.

Constraint:

FOAS Print Frequency \geq 1

FOAS Progress Tolerance

r

Default

= ε^{\frac{3}{4}}

Specifies the tolerance for

ε_{prog}

(see (8)) for which the function characterises a poor rate of progress given that it deems to be far from a solution. If this behaviour is persistent, then the function asserts that no substantial further progress can be achieved and the process is terminated with

errorid = 24

. Setting a high tolerance can lead to misinterpret reasonable progress for unsatisfactory progress or even issue a premature stop, see (8) in Section 11.4.

Constraint:

0 < FOAS Progress Tolerance < 1

FOAS Rel Stop Tolerance

r

Default

= ε^{\frac{3}{4}}

This parameter sets the value of

ε_{rel}

which specifies the relative tolerance for the convergence measures in the stopping criteria, see (6) and (7) in Section 11.4.

Constraint:

0 \leq ε_{rel} < 1

FOAS Restart Factor

r

Default

= 6.0

This factor specifies the frequency

nvar \times FOAS Restart Factor

with which the CG/LCG directions are replaced by the steepest descent direction (

d_{k} = - g_{k}

). Setting the value too small can potentially slow the convergence speed.

Constraint:

FOAS Restart Factor \geq 0

FOAS Slow Tolerance

r

Default

= ε^{\frac{1}{8}}

Specifies the tolerance for

ε_{slow}

(see (9)) for which the function characterises a slow rate of convergence. If this behaviour is deemed permanent, then the function asserts that no substantial improvement can be achieved and the process is terminated with

errorid = 50

. Setting a large tolerance can lead to incorrectly identifying a suboptimal solution, see (9) in Section 11.4.

Constraint:

FOAS Slow Tolerance > 0

FOAS Stop Tolerance

r

Default

= \max (10^{−6}, \sqrt{ε})

This parameter sets the value of

ε_{tol}

which specifies the tolerance for the convergence measures in the stopping criteria, see (6) and (7) in Section 11.4.

Constraint:

0 \leq ε_{tol} < 1

FOAS Tolerance Norm

a

Default

= INFINITY

This parameter specifies the norm used to measure some stopping metrics, such as optimality tolerances (see Section 11.4). It is possible to choose between 2-norm and ∞-norm. Solving problems using ∞-norm generally has lower computational costs than those based on 2-norm.

Constraint:

FOAS Tolerance Norm = INFINITY

TWO

Infinite Bound Size

r

Default

= 10^{20}

This defines the ‘infinite’ bound

bigbnd

in the definition of the problem constraints. Any upper bound greater than or equal to

bigbnd

will be regarded as

+ \infty

(and similarly any lower bound less than or equal to

- bigbnd

will be regarded as

- \infty

). Note that a modification of this optional parameter does not influence constraints which have already been defined; only the constraints formulated after the change will be affected.

Constraint:

Infinite Bound Size \geq 1000

Monitoring File

i

Default

= −1

i \geq 0

, the unit number for the secondary (monitoring) output. If set to

−1

, no secondary output is provided. The following information is output to the unit:

–a listing of the optional parameters;
–problem statistics, the iteration log and the final status as set by Monitoring Level;
–the solution if set by Print Solution.

Constraint:

Monitoring File \geq −1

Monitoring Level

i

Default

= 4

This parameter sets the amount of information detail that will be printed by the solver to the secondary output. The meaning of the levels is the same as Print Level.

Constraint:

0 \leq Monitoring Level \leq 5

Print File

i

Default

= advisory message unit number

i \geq 0

, the unit number for the primary output of the solver. If

Print File = −1

, the primary output is completely turned off independently of other settings. The default value is the advisory message unit number as defined by x04abf (no CPP interface) at the time of the optional parameter's initialization, e.g., at the initialization of the handle. The following information is output to the unit:

–a listing of optional parameters if set by Print Options;
–problem statistics, the iteration log and the final status from the solver as set by Print Level;
–the solution if set by Print Solution.

Constraint:

Print File \geq −1

Print Level

i

Default

= 2

This parameter defines how detailed information should be printed by the solver to the primary output.

$i$	Output
$0$	No output from the solver
$1$	Only the final status and the objective value
$2$	Problem statistics, one line per iteration showing the progress of the solution with respect to the convergence measures, final status and statistics
$3$	As level $2$ but each iteration line is longer and includes step length and progress measure
$4, 5$	As level $3$ but further details of each iteration are presented

Constraint:

0 \leq Print Level \leq 5

Print Options

a

Default

= YES

Print Options = YES

, a listing of optional parameters will be printed to the primary output.

Constraint:

Print Options = YES

NO

Print Solution

a

Default

= NO

Print Solution = YES

, the final values of the solution vector are printed on the primary and secondary outputs.

Constraint:

Print Solution = YES

NO

Stats Time

a

Default

= NO

This parameter allows you to turn on timings of various parts of the algorithm to give a better overview of where most of the time is spent. This might be helpful for a choice of different solving approaches. It is possible to choose between CPU and wall clock time. Choice

YES

is equivalent to

WALL CLOCK

Constraint:

Stats Time = YES

NO

CPU

WALL CLOCK

Task

a

Default

= MINIMIZE

This parameter specifies the required direction of the optimization. If

Task = FEASIBLE POINT

, the objective function (if set) is ignored and the algorithm stops as soon as a feasible point is found. If no objective function is set, Task reverts to

FEASIBLE POINT

automatically.

Constraint:

Task = MINIMIZE

MAXIMIZE

FEASIBLE POINT

Time Limit

r

Default

= 10^{6}

This parameter specifies a limit in seconds that the solver can use to solve one problem. If during the convergence check this limit is exceeded, the solver will terminate with

errorid = 23

error message.

Constraint:

Time Limit > 0

Verify Derivatives

a

Default

= NO

This parameter specifies whether the function should perform numerical checks on the consistency of the user-supplied gradient function. If any discrepancies are found,

errorid = 26

is returned. It is recommended that such checks are enabled when first developing the formulation of the problem, however, the verification process results in a significant increase of the number of the function evaluations and thus it shouldn't be used in production code.

Constraint:

Verify Derivatives = YES

NO

NAG Library Manual, Mark 29.2

Interfaces: FL CL CPP AD PY MB

NAG CPP Interface Introduction

E04 (Opt) Chapter Contents

e04kf: FL CL CPP AD PY MB

Namespace: nagcpp

Namespace: opt

Function: handle_solve_bounds_foas

NAG CPP Interfacenagcpp::opt::handle_solve_bounds_foas (e04kf)

▸▿ Contents

1 Purpose

2 Specification

3 Description

4 References

5 Arguments

5.1Additional Quantities

6 Exceptions and Warnings

7 Accuracy

8 Parallelism and Performance

9 Further Comments

9.1 Description of the Printed Output

10 Example

10.1 Example Program

11 Algorithmic Details

11.1 Active-Set Method and Algorithm Outline

11.2 Constrained Subproblem (NPG)

11.3 Unconstrained Subproblem (CG)

11.4 Stopping Criteria

11.5 A Note About Lagrangian Multipliers

12 Optional Parameters

12.1 Description of the Optional Parameters

NAG CPP Interface
nagcpp::opt::handle_solve_bounds_foas (e04kf)