e04dgf/e04dga : NAG Library, Mark 26

Note: this routine 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.

e04dgf/e04dga minimizes an unconstrained nonlinear function of several variables using a pre-conditioned, limited memory quasi-Newton conjugate gradient method. First derivatives (or an ‘acceptable’ finite difference approximation to them) are required. It is intended for use on large scale problems.

e04dga is a version of e04dgf that has additional arguments in order to make it safe for use in multithreaded applications (see Section 5). The initialization routine e04wbf must have been called before calling e04dga.

Fortran Interface

Subroutine e04dgf (

n, objfun, iter, objf, objgrd, x, iwork, work, iuser, ruser, ifail)

Integer, Intent (In)	::	n
Integer, Intent (Inout)	::	iuser(*), ifail
Integer, Intent (Out)	::	iter, iwork(n+1)
Real (Kind=nag_wp), Intent (Inout)	::	x(n), ruser(*)
Real (Kind=nag_wp), Intent (Out)	::	objf, objgrd(n), work(13*n)
External	::	objfun

C Header Interface

#include nagmk26.h

void

e04dgf_ ( const Integer *n,
void (NAG_CALL *objfun)( Integer *mode, const Integer *n, const double x[], double *objf, double objgrd[], const Integer *nstate, Integer iuser[], double ruser[]),
Integer *iter, double *objf, double objgrd[], double x[], Integer iwork[], double work[], Integer iuser[], double ruser[], Integer *ifail)

Fortran Interface

Subroutine e04dga (

n, objfun, iter, objf, objgrd, x, iwork, work, iuser, ruser, lwsav, iwsav, rwsav, ifail)

Integer, Intent (In)	::	n
Integer, Intent (Inout)	::	iuser(*), iwsav(610), ifail
Integer, Intent (Out)	::	iter, iwork(n+1)
Real (Kind=nag_wp), Intent (Inout)	::	x(n), ruser(*), rwsav(475)
Real (Kind=nag_wp), Intent (Out)	::	objf, objgrd(n), work(13*n)
Logical, Intent (Inout)	::	lwsav(120)
External	::	objfun

C Header Interface

#include nagmk26.h

void

e04dga_ ( const Integer *n,
void (NAG_CALL *objfun)( Integer *mode, const Integer *n, const double x[], double *objf, double objgrd[], const Integer *nstate, Integer iuser[], double ruser[]),
Integer *iter, double *objf, double objgrd[], double x[], Integer iwork[], double work[], Integer iuser[], double ruser[], logical lwsav[], Integer iwsav[], double rwsav[], Integer *ifail)

Before calling e04dga, or either of the option setting routines e04dja or e04dka, routine e04wbf must be called. The specification for e04wbf is:

Fortran Interface

Subroutine e04wbf (

rname, cwsav, lcwsav, lwsav, llwsav, iwsav, liwsav, rwsav, lrwsav, ifail)

Integer, Intent (In)	::	lcwsav, llwsav, liwsav, lrwsav
Integer, Intent (Inout)	::	ifail
Integer, Intent (Out)	::	iwsav(liwsav)
Real (Kind=nag_wp), Intent (Out)	::	rwsav(lrwsav)
Logical, Intent (Out)	::	lwsav(llwsav)
Character (*), Intent (In)	::	rname
Character (80), Intent (Out)	::	cwsav(lcwsav)

C Header Interface

#include nagmk26.h

void	e04wbf ( const char rname, char cwsav[], const Integer lcwsav, logical lwsav[], const Integer llwsav, Integer iwsav[], const Integer liwsav, double rwsav[], const Integer lrwsav, Integer ifail, const Charlen length_rname, const Charlen length_cwsav)

e04wbf should be called with

rname ='e04dga'

. lcwsav, llwsav, liwsav and lrwsav, the declared lengths of cwsav, lwsav, iwsav and rwsav respectively, must satisfy:

$lcwsav \geq 1$
$llwsav \geq 120$
$liwsav \geq 610$
$lrwsav \geq 475$

The contents of the arrays cwsav, lwsav, iwsav and rwsav must not be altered between calling routines e04dga, e04dja, e04dka and e04wbf.

e04dgf/e04dga is designed to solve unconstrained minimization problems of the form

\underset{x \in R^{n}}{minimize} F (x) subject to - \infty \leq x \leq \infty,

where

x

is an

n

-element vector.

You must supply an initial estimate of the solution.

For maximum reliability, it is preferable to provide all first partial derivatives. If all of the derivatives cannot be provided, you are recommended to obtain approximate values (using finite differences) by calling e04xaf/e04xaa from within objfun. This is illustrated in Section 10 in e04djf/e04dja.

The method used by e04dgf/e04dga is described in Section 11.

Gill P E and Murray W (1979) Conjugate-gradient methods for large-scale nonlinear optimization Technical Report SOL 79-15 Department of Operations Research, Stanford University

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

If on entry

ifail = 0

or

- 1

, explanatory error messages are output on the current error message unit (as defined by x04aaf).

Note: e04dgf/e04dga may return useful information for one or more of the following detected errors or warnings.

Errors or warnings detected by the routine:

On successful exit (

ifail = 0

) the accuracy of the solution will be as defined by the optional parameter Optimality Tolerance (

default value = ε^{0.8}

).

e04dgf/e04dga makes calls to BLAS and/or LAPACK routines, which may be threaded within the vendor library used by this implementation. Consult the documentation for the vendor library for further information.

Please consult the X06 Chapter Introduction for information on how to control and interrogate the OpenMP environment used within this routine. Please also consult the Users' Note for your implementation for any additional implementation-specific information.

To evaluate an ‘acceptable’ set of finite difference intervals using e04xaf/e04xaa requires

2

function evaluations per variable for a well-scaled problem and up to

6

function evaluations per variable for a badly scaled problem.

This example finds a minimum of the function

F = e^{x_{1}} (4 x_{1}^{2} + 2 x_{2}^{2} + 4 x_{1} x_{2} + 2 x_{2} + 1) .

The initial point is

x_{0} = {(- 1.0, 1.0)}^{T},

and

F (x_{0}) = 1.8394

(to five figures).

The optimal solution is

x^{*} = {(0.5, - 1.0)}^{T},

and

F (x^{*}) = 0

.

The document for e04djf/e04dja includes an example program to solve the same problem using some of the optional parameters described in Section 12.

Note: the following programs illustrate the use of e04dgf and e04dga.

Note: the remainder of this document is intended for more advanced users. Section 11 contains a detailed description of the algorithm which may be needed in order to understand Section 12. Section 12 describes the optional parameters which may be set by calls to e04djf/e04dja and/or e04dkf/e04dka.

This section contains a description of the method used by e04dgf/e04dga.

e04dgf/e04dga uses a pre-conditioned conjugate gradient method and is based upon algorithm PLMA as described in Section 4.8.3 of Gill and Murray (1979) and Gill et al. (1981).

The algorithm proceeds as follows:

Let

x_{0}

be a given starting point and let

k

denote the current iteration, starting with

k = 0

. The iteration requires

g_{k}

, the gradient vector evaluated at

x_{k}

, the

k

th estimate of the minimum. At each iteration a vector

p_{k}

(known as the direction of search) is computed and the new estimate

x_{k + 1}

is given by

x_{k} + α_{k} p_{k}

where

α_{k}

(the step length) minimizes the function

F (x_{k} + α_{k} p_{k})

with respect to the scalar

α_{k}

. A choice of initial step

α_{0}

is taken as

α_{0} = \min \{1, 2 \times |F_{k} - F_{est}| / g_{k}^{T} g_{k}\}

where

F_{est}

is a user-supplied estimate of the function value at the solution. If

F_{est}

is not specified, the software always chooses the unit step length for

α_{0}

. Subsequent step length estimates are computed using cubic interpolation with safeguards.

A quasi-Newton method can be used to compute the search direction

p_{k}

by updating the inverse of the approximate Hessian

(H_{k})

and computing

p_{k + 1} = - H_{k + 1} g_{k + 1} .

(1)

The updating formula for the approximate inverse is given by

H_{k + 1} = H_{k} - \frac{1}{y_{k}^{T} s_{k}} (H_{k} y_{k} s_{k}^{T} + s_{k} y_{k}^{T} H_{k}) + \frac{1}{y_{k}^{T} s_{k}} (1 + \frac{y_{k}^{T} H_{k} y_{k}}{y_{k}^{T} s_{k}}) s_{k} s_{k}^{T},

(2)

where

y_{k} = g_{k - 1} - g_{k}

and

s_{k} = x_{k + 1} - x_{k} = α_{k} p_{k}

.

The method used to obtain the search direction is based upon computing

p_{k + 1}

as

- H_{k + 1} g_{k + 1}

where

H_{k + 1}

is a matrix obtained by updating the identity matrix with a limited number of quasi-Newton corrections. The storage of an

n

by

n

matrix is avoided by storing only the vectors that define the rank two corrections – hence the term ‘limited-memory’ quasi-Newton method. The precise method depends upon the number of updating vectors stored. For example, the direction obtained with the ‘one-step’ limited memory update is given by (1) using (2) with

H_{k}

equal to the identity matrix, viz.

p_{k + 1} = - g_{k + 1} + \frac{1}{y_{k}^{T} s_{k}} (s_{k}^{T} g_{k + 1} y_{k} + y_{k}^{T} g_{k + 1} s_{k}) - \frac{s_{k}^{T} g_{k + 1}}{y_{k}^{T} s_{k}} (1 + \frac{y_{k}^{T} y_{k}}{y_{k}^{T} s_{k}}) s_{k} .

Using a limited-memory quasi-Newton formula, such as the one above, guarantees

p_{k + 1}

to be a descent direction if all the inner products

y_{k}^{T} s_{k}

are positive for all vectors

y_{k}

and

s_{k}

used in the updating formula.

Several optional parameters in e04dgf/e04dga define choices in the problem specification or the algorithm logic. In order to reduce the number of formal arguments of e04dgf/e04dga 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 following is a list of the optional parameters available. A full description of each optional parameter is provided in Section 12.1.

Optional parameters may be specified by calling one, or both, of the routines e04djf/e04dja and e04dkf/e04dka before a call to e04dgf/e04dga.

e04djf/e04dja reads options from an external options file, with Begin and End as the first and last lines respectively and each intermediate line defining a single optional parameter. For example,

Begin
    Print Level = 1
End

The call

Call e04djf/e04dja (ioptns, inform)

can then be used to read the file on unit ioptns. inform will be zero on successful exit. e04djf/e04dja should be consulted for a full description of this method of supplying optional parameters.

e04dkf/e04dka can be called to supply options directly, one call being necessary for each optional parameter. For example,

Call e04dkf/e04dka ('Print Level = 1')

e04dkf/e04dka should be consulted for a full description of this method of supplying optional parameters.

All optional parameters not specified by you are set to their default values. Optional parameters specified by you are unaltered by e04dgf/e04dga (unless they define invalid values) and so remain in effect for subsequent calls unless altered by you.

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

The summary line contains:

Keywords and character values are case and white space insensitive.

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

This value of

r

specifies the user-supplied guess of the optimum objective function value

F_{est}

. This value is used to calculate an initial step length

α_{0}

(see Section 11). If the value of

r

is not specified (the default), then this has the effect of setting

α_{0}

to unity. It should be noted that for badly scaled functions a unit step along the steepest descent direction will often compute the objective function at very large values of

x

.

The parameter defines

ε_{r}

, which is intended to be a measure of the accuracy with which the problem function

F (x)

can be computed. If

r < ε

or

r \geq 1

, the default value is used.

The value of

ε_{r}

should reflect the relative precision of

1 + |F (x)|

; i.e.,

ε_{r}

acts as a relative precision when

|F|

is large, and as an absolute precision when

|F|

is small. For example, if

F (x)

is typically of order

1000

and the first six significant digits are known to be correct, an appropriate value for

ε_{r}

would be

10^{- 6}

. In contrast, if

F (x)

is typically of order

10^{- 4}

and the first six significant digits are known to be correct, an appropriate value for

ε_{r}

would be

10^{- 10}

. The choice of

ε_{r}

can be quite complicated for badly scaled problems; see Chapter 8 of Gill et al. (1981) for a discussion of scaling techniques. The default value is appropriate for most simple functions that are computed with full accuracy. However when the accuracy of the computed function values is known to be significantly worse than full precision, the value of

ε_{r}

should be large enough so that no attempt will be made to distinguish between function values that differ by less than the error inherent in the calculation.

The value of

i

specifies the maximum number of iterations allowed before termination. If

i < 0

, the default value is used.

Problems whose Hessian matrices at the solution contain sets of clustered eigenvalues are likely to be minimized in significantly fewer than

n

iterations. Problems without this property may require anything between

n

and

5 n

iterations, with approximately

2 n

iterations being a common figure for moderately difficult problems.

The value

r

controls the accuracy with which the step

α

taken during each iteration approximates a minimum of the function along the search direction (the smaller the value of

r

, the more accurate the linesearch). The default value

r = 0.9

requests an inaccurate search, and is appropriate for most problems. A more accurate search may be appropriate when it is desirable to reduce the number of iterations – for example, if the objective function is cheap to evaluate. If

r < 0

or

r \geq 1

, the default value is used.

Normally each optional parameter specification is printed as it is supplied. Optional parameter Nolist may be used to suppress the printing and optional parameter List may be used to restore printing.

If

r > 0

, the maximum allowable step length for the linesearch is taken as

\min (\frac{1}{x02amf ()}, \frac{r}{‖p_{k}‖})

. If

r \leq 0

, the default value is used.

The parameter

r

specifies the accuracy to which you wish the final iterate to approximate a solution of the problem. Broadly speaking,

r

indicates the number of correct figures desired in the objective function at the solution. For example, if

r

is

10^{- 6}

and termination occurs with

ifail = 0

(see Section 5), then the final point satisfies the termination criteria, where

τ_{F}

represents Optimality Tolerance. If

r < ε_{r}

or

r \geq 1

, the default value is used. If Optimality Tolerance is chosen below a certain threshold, it will automatically be reset to another value.

The value

i

controls the amount of printout produced by e04dgf/e04dga, as indicated below. A detailed description of the printout is given in Section 12.2 (summary output at each iteration and the final solution).

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

These keywords take effect only if

Verify Level > 0

. They may be used to control the verification of gradient elements computed by objfun. For example, if the first

30

elements of the objective gradient appeared to be correct in an earlier run, so that only element

31

remains questionable, it is reasonable to specify

Start Objective Check at Variable = 31

. If the first

30

variables appear linearly in the objective, so that the corresponding gradient elements are constant, the above choice would also be appropriate.

If

i_{1} \leq 0

or

i_{1} > \max (1, \min (n, i_{2}))

, the default value is used. If

i_{2} \leq 0

or

i_{2} > n

, the default value is used.

These keywords refer to finite difference checks on the gradient elements computed by objfun. Gradients are verified at the user-supplied initial estimate of the solution. The possible choices for

i

are as follows:

$i$	Meaning
$- 1$	No checks are performed.
$0$	Only a ‘cheap’ test will be performed, requiring one call to objfun.
$1$	In addition to the ‘cheap’ test, individual gradient elements will also be checked using a reliable (but more expensive) test.

For example, the objective gradient will be verified if Verify,

Verify = YES

, Verify Gradients, Verify Objective Gradients or

Verify Level = 1

is specified.

This section describes the intermediate printout and final printout produced by e04dgf/e04dga. You can control the level of printed output (see the description of the optional parameter Print Level). Note that the intermediate printout and final printout are produced only if

Print Level \geq 10

(the default for e04dgf, by default no output is produced by e04dga).

The following line of summary output (

< 80

characters) is produced at every iteration. In all cases, the values of the quantities are those in effect on completion of the given iteration.

Itn	is the iteration count.
Step	is the step $α_{k}$ taken along the computed search direction. On reasonably well-behaved problems, the unit step (i.e., $α_{k} = 1$ ) will be taken as the solution is approached.
Nfun	is the cumulated number of evaluations of the objective function needed for the linesearch. Evaluations needed for the verification of the gradients by finite differences are not included. Nfun is printed as a guide to the amount of work required for the linesearch. e04dgf/e04dga will perform at most $11$ function evaluations per iteration.
Objective	is the value of the objective function at $x_{k}$ .
Norm G	is the Euclidean norm of the gradient of the objective function at $x_{k}$ .
Norm X	is the Euclidean norm of $x_{k}$ .
Norm (X(k-1)-X(k))	is the Euclidean norm of $x_{k - 1} - x_{k}$ .

The following describes the printout for each variable.

Variable	gives the name (Varbl) and index $j$ , for $j = 1, 2, \dots, n$ of the variable.
Value	is the value of the variable at the final iteration.
Gradient Value	is the value of the gradient of the objective function with respect to the $j$ th variable at the final iteration.

Numerical values are output with a fixed number of digits; they are not guaranteed to be accurate to this precision.

(i)	$F_{k - 1} - F_{k} < τ_{F} (1 + \|F_{k}\|)$
(ii)	$‖x_{k - 1} - x_{k}‖ < \sqrt{τ_{F}} (1 + ‖x_{k}‖)$
(iii)	$‖g_{k}‖ \leq \sqrt[3]{τ_{F}} (1 + \|F_{k}\|)$ or $‖g_{k}‖ < ε_{A}$

NAG Library Routine Document

e04dgf (uncon_conjgrd_comp_old)
e04dga (uncon_conjgrd_comp)

▸▿ Contents

1

Purpose

2

Specification

2.1

Specification for e04dgf

2.2

Specification for e04dga

3

Description

4

References

5

Arguments

6

Error Indicators and Warnings

7

Accuracy

8

Parallelism and Performance

9

Further Comments

10

Example

10.1

Program Text

10.2

Program Data

10.3

Program Results

11

Algorithmic Details

12

Optional Parameters

12.1

Description of the Optional Parameters

12.2

Description of Printed Output

NAG Library Routine Document

e04dgf (uncon_conjgrd_comp_old)e04dga (uncon_conjgrd_comp)

▸▿ Contents

1 Purpose

2 Specification

2.1 Specification for e04dgf

2.2 Specification for e04dga

3 Description

4 References

5 Arguments

6 Error Indicators and Warnings

7 Accuracy

8 Parallelism and Performance

9 Further Comments

10 Example

10.1 Program Text

10.2 Program Data

10.3 Program Results

11 Algorithmic Details

12 Optional Parameters

12.1 Description of the Optional Parameters

12.2 Description of Printed Output

e04dgf (uncon_conjgrd_comp_old)
e04dga (uncon_conjgrd_comp)

1

Purpose

2

Specification

2.1

Specification for e04dgf

2.2

Specification for e04dga

3

Description

4

References

5

Arguments

6

Error Indicators and Warnings

7

Accuracy

8

Parallelism and Performance

9

Further Comments

10

Example

10.1

Program Text

10.2

Program Data

10.3

Program Results

11

Algorithmic Details

12

Optional Parameters

12.1

Description of the Optional Parameters

12.2

Description of Printed Output