NAG Library Routine Document
E04XAF/E04XAA computes an approximation to the gradient vector and/or the Hessian matrix for use in conjunction with, or following the use of an optimization routine (such as E04UFF/E04UFA
E04XAA is a version of E04XAF that has additional parameters in order to make it safe for use in multithreaded applications (see Section 5
). The initialization routine E04WBF must
have been called before calling E04XAA.
2.1 Specification for E04XAF
|SUBROUTINE E04XAF (
||MSGLVL, N, EPSRF, X, MODE, OBJFUN, LDH, HFORW, OBJF, OBJGRD, HCNTRL, H, IWARN, WORK, IUSER, RUSER, INFO, IFAIL)
||MSGLVL, N, MODE, LDH, IWARN, IUSER(*), INFO(N), IFAIL
||EPSRF, X(N), HFORW(N), OBJF, OBJGRD(N), HCNTRL(N), H(LDH,*), WORK(*), RUSER(*)
2.2 Specification for E04XAA
|SUBROUTINE E04XAA (
||MSGLVL, N, EPSRF, X, MODE, OBJFUN, LDH, HFORW, OBJF, OBJGRD, HCNTRL, H, IWARN, WORK, IUSER, RUSER, INFO, LWSAV, IWSAV, RWSAV, IFAIL)
||MSGLVL, N, MODE, LDH, IWARN, IUSER(*), INFO(N), IWSAV(1), IFAIL
||EPSRF, X(N), HFORW(N), OBJF, OBJGRD(N), HCNTRL(N), H(LDH,*), WORK(*), RUSER(*), RWSAV(1)
E04XAF/E04XAA is similar to routine FDCALC described in Gill et al. (1983a)
. It should be noted that this routine aims to compute sufficiently accurate estimates of the derivatives for use with an optimization algorithm. If you require more accurate estimates you should refer to Chapter D04
E04XAF/E04XAA computes finite difference approximations to the gradient vector and the Hessian matrix for a given function. The simplest approximation involves the forward-difference formula, in which the derivative
of a univariate function
is approximated by the quantity
for some interval
, where the subscript 'F' denotes ‘forward-difference’ (see Gill et al. (1983b)
To summarise the procedure used by E04XAF/E04XAA (for the case when the objective function is available and you require estimates of gradient values and Hessian matrix diagonal values, i.e.,
) consider a univariate function
at the point
. (In order to obtain the gradient of a multivariate function
-vector, the procedure is applied to each component of
, keeping the other components fixed.) Roughly speaking, the method is based on the fact that the bound on the relative truncation error in the forward-difference approximation tends to be an increasing function of
, while the relative condition error bound is generally a decreasing function of
, hence changes in
will tend to have opposite effects on these errors (see Gill et al. (1983b)
The ‘best’ interval
is given by
is an estimate of
is an estimate of the relative error associated with computing the function (see Chapter 8 of Gill et al. (1981)
). Given an interval
is defined by the second-order approximation
The decision as to whether a given value of
is acceptable involves
, the following bound on the relative condition error in
is taken as an arbitrary large number.)
The procedure selects the interval
(to be used in computing
) from a sequence of trial intervals
. The initial trial interval is taken as
unless you specify the initial value to be used.
The value of
for a trial value
is defined as ‘acceptable’ if it lies in the interval
. In this case
is taken as
, and the current value of
is used to compute
is unacceptable, the next trial interval is chosen so that the relative condition error bound will either decrease or increase, as required. If the bound on the relative condition error is too large, a larger interval is used as the next trial value in an attempt to reduce the condition error bound. On the other hand, if the relative condition error bound is too small,
The procedure will fail to produce an acceptable value of in two situations. Firstly, if is extremely small, then may never become small, even for a very large value of the interval. Alternatively, may never exceed , even for a very small value of the interval. This usually implies that is extremely large, and occurs most often near a singularity.
As a check on the validity of the estimated first derivative, the procedure provides a comparison of the forward-difference approximation computed with
(as above) and the central-difference approximation computed with
. Using the central-difference formula the first derivative can be approximated by
. If the values
do not display some agreement, neither can be considered reliable.
When both function and gradients are available and you require the Hessian matrix (i.e.,
) E04XAF/E04XAA follows a similar procedure to the case above with the exception that the gradient function
is substituted for the objective function and so the forward-difference interval for the first derivative of
with respect to variable
is computed. The
th column of the approximate Hessian matrix is then defined as in Chapter 2 of Gill et al. (1981)
is the best forward-difference interval associated with the
th component of
is the vector with unity in the
th position and zeros elsewhere.
When only the objective function is available and you require the gradients and Hessian matrix (i.e.,
) E04XAF/E04XAA again follows the same procedure as the case for
except that this time the value of
for a trial value
is defined as acceptable if it lies in the interval
and the initial trial interval is taken as
The approximate Hessian matrix
is then defined as in Chapter 2 of Gill et al. (1981)
Gill P E, Murray W, Saunders M A and Wright M H (1983a) Documentation for FDCALC and FDCORE Technical Report SOL 83–6 Stanford University
Gill P E, Murray W, Saunders M A and Wright M H (1983b) Computing forward-difference intervals for numerical optimization SIAM J. Sci. Statist. Comput. 4 310–321
Gill P E, Murray W and Wright M H (1981) Practical Optimization Academic Press
- 1: MSGLVL – INTEGERInput
: must indicate the amount of intermediate output desired (see Section 8.1
for a description of the printed output). All output is written on the current advisory message unit (see X04ABF
|Value ||Definition |
|0 ||No printout |
|1 ||A summary is printed out for each variable plus any warning messages. |
|Other ||Values other than and should normally be used only at the direction of NAG. |
- 2: N – INTEGERInput
On entry: the number of independent variables.
- 3: EPSRF – REAL (KIND=nag_wp)Input
: must define
, which is intended to be a measure of the accuracy with which the problem function
can be computed. The value of
should reflect the relative precision of
, i.e., acts as a relative precision when
is large, and as an absolute precision when
is small. For example, if
is typically of order
and the first six significant digits are known to be correct, an appropriate value for
A discussion of EPSRF
is given in Chapter 8 of Gill et al. (1981)
. If EPSRF
is either too small or too large on entry a warning will be printed if
, the parameter IWARN
set to the appropriate value on exit and E04XAF/E04XAA will use a default value of
is the machine precision
If on entry then E04XAF/E04XAA will use the default value internally. The default value will be appropriate for most simple functions that are computed with full accuracy.
- 4: X(N) – REAL (KIND=nag_wp) arrayInput
On entry: the point at which the derivatives are to be computed.
- 5: MODE – INTEGERInput/Output
: indicates which derivatives are required.
- The gradient and Hessian diagonal values having supplied the objective function via OBJFUN.
- The Hessian matrix having supplied both the objective function and gradients via OBJFUN.
- The gradient values and Hessian matrix having supplied the objective function via OBJFUN.
: is changed only
if you set MODE
negative in OBJFUN
, i.e., you have requested termination of E04XAF/E04XAA.
- 6: OBJFUN – SUBROUTINE, supplied by the user.External Procedure
must calculate the objective function; otherwise if
must calculate the objective function and the gradients.
The specification of OBJFUN
||MODE, N, NSTATE, IUSER(*)
||X(N), OBJF, OBJGRD(N), RUSER(*)
- 1: MODE – INTEGERInput/Output
indicates which parameter values within OBJFUN
need to be set.
: to OBJFUN
is always set to the value that you set it to before the call to E04XAF/E04XAA.
: its value must not be altered unless you wish to indicate a failure within OBJFUN
, in which case it should be set to a negative value. If MODE
is negative on exit from OBJFUN
, the execution of E04XAF/E04XAA is terminated with IFAIL
set to MODE
- 2: N – INTEGERInput
On entry: the number of variables as input to E04XAF/E04XAA.
- 3: X(N) – REAL (KIND=nag_wp) arrayInput
On entry: the point at which the objective function (and gradients if ) is to be evaluated.
- 4: OBJF – REAL (KIND=nag_wp)Output
On exit: must be set to the value of the objective function.
- 5: OBJGRD(N) – REAL (KIND=nag_wp) arrayOutput
must contain the value of the first derivative with respect to
need not be set.
- 6: NSTATE – INTEGERInput
: will be set to
on the first call of OBJFUN
by E04XAF/E04XAA, and is
for all subsequent calls. Thus, if you wish, NSTATE
may be tested within OBJFUN
in order to perform certain calculations once only. For example you may read data.
- 7: IUSER() – INTEGER arrayUser Workspace
- 8: RUSER() – REAL (KIND=nag_wp) arrayUser Workspace
is called with the parameters IUSER
as supplied to E04XAF/E04XAA. You are free to use the arrays IUSER
to supply information to OBJFUN
as an alternative to using COMMON global variables.
must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which E04XAF/E04XAA is called. Parameters denoted as Input
be changed by this procedure.
- 7: LDH – INTEGERInput
: the first dimension of the array H
as declared in the (sub)program from which E04XAF/E04XAA is called.
- 8: HFORW(N) – REAL (KIND=nag_wp) arrayInput/Output
: the initial trial interval for computing the appropriate partial derivative to the
, then the initial trial interval is computed by E04XAF/E04XAA (see Section 3
On exit: is the best interval found for computing a forward-difference approximation to the appropriate partial derivative for the th variable.
- 9: OBJF – REAL (KIND=nag_wp)Output
: the value of the objective function evaluated at the input vector in X
- 10: OBJGRD(N) – REAL (KIND=nag_wp) arrayOutput
contains the best estimate of the first partial derivative for the
contains the first partial derivative for the
th variable evaluated at the input vector in X
- 11: HCNTRL(N) – REAL (KIND=nag_wp) arrayOutput
On exit: is the best interval found for computing a central-difference approximation to the appropriate partial derivative for the th variable.
- 12: H(LDH,) – REAL (KIND=nag_wp) arrayOutput
the second dimension of the array H
must be at least
and at least
, the estimated Hessian diagonal elements are contained in the first column of this array.
If or , the estimated Hessian matrix is contained in the leading by part of this array.
- 13: IWARN – INTEGEROutput
on successful exit.
If the value of EPSRF
on entry is too small or too large then IWARN
is set to
respectively on exit and the default value for EPSRF
is used within E04XAF/E04XAA.
then warnings will be printed if EPSRF
is too small or too large.
- 14: WORK() – REAL (KIND=nag_wp) arrayWorkspace
the dimension of the array WORK
must be at least
and at least
- 15: IUSER() – INTEGER arrayUser Workspace
- 16: RUSER() – REAL (KIND=nag_wp) arrayUser Workspace
are not used by E04XAF/E04XAA, but are passed directly to OBJFUN
and may be used to pass information to this routine as an alternative to using COMMON global variables.
- 17: INFO(N) – INTEGER arrayOutput
represents diagnostic information on variable
. (See Section 6
for more details.)
- 18: IFAIL – INTEGERInput/Output
Note: for E04XAA, IFAIL does not occur in this position in the parameter list. See the additional parameters described below
must be set to
. If you are unfamiliar with this parameter you should refer to Section 3.3
in the Essential Introduction for details.
For environments where it might be inappropriate to halt program execution when an error is detected, the value
is recommended. If the output of error messages is undesirable, then the value
is recommended. Otherwise, if you are not familiar with this parameter, the recommended value is
. When the value is used it is essential to test the value of IFAIL on exit.
unless the routine detects an error or a warning has been flagged (see Section 6
- Note: the following are additional parameters for specific use with E04XAA. Users of E04XAF therefore need not read the remainder of this description.
- 18: LWSAV() – LOGICAL arrayCommunication Array
- 19: IWSAV() – INTEGER arrayCommunication Array
- 20: RWSAV() – REAL (KIND=nag_wp) arrayCommunication Array
These parameters are no longer required by E04XAF/E04XAA.
- 21: IFAIL – INTEGERInput/Output
see the parameter description for IFAIL
6 Error Indicators and Warnings
On exit from E04XAF/E04XAA both diagnostic parameters INFO
should be tested. IFAIL
represents an overall diagnostic indicator, whereas the integer array INFO
represents diagnostic information on each variable.
If on entry
, explanatory error messages are output on the current error message unit (as defined by X04AAF
Errors or warnings detected by the routine:
A negative value of IFAIL
indicates an exit from E04XAF/E04XAA because you set MODE
negative in OBJFUN
. The value of IFAIL
will be the same as your setting of MODE
On entry, one or more of the following conditions are satisfied: , is invalid.
One or more variables have a nonzero INFO
value. This may not necessarily represent an unsuccessful exit – see diagnostic information on INFO
Diagnostic information returned via INFO
is as follows:
The appropriate function appears to be constant.
is set to the initial trial interval value (see Section 3
) corresponding to a well-scaled problem and Error est
. in the printed output is set to zero. This value occurs when the estimated relative condition error in the first derivative approximation is unacceptably large for every value of the finite difference interval. If this happens when the function is not constant the initial interval may be too small; in this case, it may be worthwhile to rerun E04XAF/E04XAA with larger initial trial interval values supplied in HFORW
(see Section 3
). This error may also occur if the function evaluation includes an inordinately large constant term or if EPSRF
is too large.
The appropriate function appears to be linear or odd.
is set to the smallest interval with acceptable bounds on the relative condition error in the forward- and backward-difference estimates. In this case, the estimated relative condition error in the second derivative approximation remained large for every trial interval, but the estimated error in the first derivative approximation was acceptable for at least one interval. If the function is not linear or odd the relative condition error in the second derivative may be decreasing very slowly, it may be worthwhile to rerun E04XAF/E04XAA with larger initial trial interval values supplied in HFORW
(see Section 3
The second derivative of the appropriate function appears to be so large that it cannot be reliably estimated (i.e., near a singularity). is set to the smallest trial interval.
This value occurs when the relative condition error estimate in the second derivative remained very small for every trial interval.
If the second derivative is not large the relative condition error in the second derivative may be increasing very slowly. It may be worthwhile to rerun E04XAF/E04XAA with smaller initial trial interval values supplied in HFORW
(see Section 3
). This error may also occur when the given value of EPSRF
is not a good estimate of a bound on the absolute error in the appropriate function (i.e., EPSRF
is too small).
The algorithm terminated with an apparently acceptable estimate of the second derivative. However the forward-difference estimates of the appropriate first derivatives (computed with the final estimate of the ‘optimal’ forward-difference interval) and the central difference estimates (computed with the interval used to compute the final estimate of the second derivative) do not agree to half a decimal place. The usual reason that the forward- and central-difference estimates fail to agree is that the first derivative is small.
If the first derivative is not small, it may be helpful to execute the procedure at a different point.
If on exit the algorithm terminated successfully, i.e., the forward-difference estimates of the appropriate first derivatives (computed with the final estimate of the ‘optimal’ forward-difference interval ) and the central-difference estimates (computed with the interval used to compute the final estimate of the second derivative) agree to at least half a decimal place.
In short word length implementations when computing the full Hessian matrix given function values only (i.e., ) the elements of the computed Hessian will have at best to figures of accuracy.
To evaluate an acceptable set of finite difference intervals for a well-scaled problem, the routine will require around two function evaluations per variable; in a badly scaled problem however, as many as six function evaluations per variable may be needed.
If you request the full Hessian matrix supplying both function and gradients (i.e.,
) or function only (i.e.,
) then a further N
function evaluations respectively are required.
8.1 Description of the Printed Output
The following is a description of the printed output from E04XAF/E04XAA as controlled by the parameter MSGLVL
is as follows:
||number of variable for which the difference interval has been computed.
||th variable of as set by you.
|F. dif. int.
||the best interval found for computing a forward-difference approximation to the appropriate partial derivative with respect to the th variable.
|C. dif. int.
||the best interval found for computing a central-difference approximation to the appropriate partial derivative with respect to the th variable.
||a bound on the estimated error in the final forward-difference approximation. When , Error est. is set to zero.
||best estimate of the first partial derivative with respect to the th variable.
|Hess diag est.
||best estimate of the second partial derivative with respect to the th variable.
||the number of function evaluations used to compute the final difference intervals for the th variable.
||the value of INFO for the th variable.
This example computes the gradient vector and the Hessian matrix of the following function:
at the point
9.1 Program Text
Note: the following programs illustrate the use of E04XAF and E04XAA.
Program Text (e04xafe.f90)
Program Text (e04xaae.f90)
9.2 Program Data
9.3 Program Results
Program Results (e04xafe.r)
Program Results (e04xaae.r)