NAG FL Interface
e04ybf (lsq_​check_​hessian)

1 Purpose

e04ybf checks that a user-supplied subroutine for evaluating the second derivative term of the Hessian matrix of a sum of squares is consistent with a user-supplied subroutine for calculating the corresponding first derivatives.

2 Specification

Fortran Interface
Subroutine e04ybf ( m, n, lsqfun, lsqhes, x, fvec, fjac, ldfjac, b, lb, iw, liw, w, lw, ifail)
Integer, Intent (In) :: m, n, ldfjac, lb, liw, lw
Integer, Intent (Inout) :: iw(liw), ifail
Real (Kind=nag_wp), Intent (In) :: x(n)
Real (Kind=nag_wp), Intent (Inout) :: fjac(ldfjac,n), w(lw)
Real (Kind=nag_wp), Intent (Out) :: fvec(m), b(lb)
External :: lsqfun, lsqhes
C Header Interface
#include <nag.h>
void  e04ybf_ (const Integer *m, const Integer *n,
void (NAG_CALL *lsqfun)(Integer *iflag, const Integer *m, const Integer *n, const double xc[], double fvec[], double fjac[], const Integer *ldfjac, Integer iw[], const Integer *liw, double w[], const Integer *lw),
void (NAG_CALL *lsqhes)(Integer *iflag, const Integer *m, const Integer *n, const double fvec[], const double xc[], double b[], const Integer *lb, Integer iw[], const Integer *liw, double w[], const Integer *lw),
const double x[], double fvec[], double fjac[], const Integer *ldfjac, double b[], const Integer *lb, Integer iw[], const Integer *liw, double w[], const Integer *lw, Integer *ifail)
The routine may be called by the names e04ybf or nagf_opt_lsq_check_hessian.

3 Description

Routines for minimizing a sum of squares of m nonlinear functions (or ‘residuals’), fix1,x2,,xn, for i=1,2,,m and mn, may require you to supply a subroutine to evaluate the quantities
bjk=i=1mfi 2fi xjxk  
for j=1,2,,n and k=1,2,,j. e04ybf is designed to check the bjk calculated by such user-supplied subroutines. As well as the routine to be checked (lsqhes), you must supply a subroutine (lsqfun) to evaluate the fi and their first derivatives, and a point x=x1,x2,,xnT at which the checks will be made. Note that e04ybf checks routines of the form required by e04hef. e04ybf is essentially identical to CHKLSH in the NPL Algorithms Library.
e04ybf first calls user-supplied subroutines lsqfun and lsqhes to evaluate the first derivatives and the bjk at x. Let J denote the m by n matrix of first derivatives of the residuals. The Hessian matrix of the sum of squares,
G=JTJ+B,  
is calculated and projected onto two orthogonal vectors y and z to give the scalars yTGy and zTGz respectively. The same projections of the Hessian matrix are also estimated by finite differences, giving
p=yTgx+hy-yTgx/h  and q=zTgx+hz-zTgx/h  
respectively, where g denotes the gradient vector of the sum of squares at the point in brackets and h is a small positive scalar. If the relative difference between p and yTGy or between q and zTGz is judged too large, an error indicator is set.

4 References

None.

5 Arguments

1: m Integer Input
2: n Integer Input
On entry: the number m of residuals, fix, and the number n of variables, xj.
Constraint: 1nm.
3: lsqfun Subroutine, supplied by the user. External Procedure
lsqfun must calculate the vector of values fix and their first derivatives fi xj at any point x. (e04hef gives you the option of resetting arguments of lsqfun to cause the minimization process to terminate immediately. e04ybf will also terminate immediately, without finishing the checking process, if the argument in question is reset.)
The specification of lsqfun is:
Fortran Interface
Subroutine lsqfun ( iflag, m, n, xc, fvec, fjac, ldfjac, iw, liw, w, lw)
Integer, Intent (In) :: m, n, ldfjac, liw, lw
Integer, Intent (Inout) :: iflag, iw(liw)
Real (Kind=nag_wp), Intent (In) :: xc(n)
Real (Kind=nag_wp), Intent (Inout) :: fjac(ldfjac,n), w(lw)
Real (Kind=nag_wp), Intent (Out) :: fvec(m)
C Header Interface
void  lsqfun_ (Integer *iflag, const Integer *m, const Integer *n, const double xc[], double fvec[], double fjac[], const Integer *ldfjac, Integer iw[], const Integer *liw, double w[], const Integer *lw)
1: iflag Integer Input/Output
On entry: to lsqfun, iflag will be set to 2.
On exit: if you reset iflag to some negative number in lsqfun and return control to e04ybf, the routine will terminate immediately with ifail set to your setting of iflag.
2: m Integer Input
On entry: the numbers m of residuals.
3: n Integer Input
On entry: the numbers n of variables.
4: xcn Real (Kind=nag_wp) array Input
On entry: the point x at which the values of the fi and the fi xj are required.
5: fvecm Real (Kind=nag_wp) array Output
On exit: unless iflag is reset to a negative number, fveci must contain the value of fi at the point x, for i=1,2,,m.
6: fjacldfjacn Real (Kind=nag_wp) array Output
On exit: unless iflag is reset to a negative number, fjacij must contain the value of fi xj at the point x, for i=1,2,,m and j=1,2,,n.
7: ldfjac Integer Input
On entry: the first dimension of the array fjac as declared in the (sub)program from which e04ybf is called.
8: iwliw Integer array Workspace
9: liw Integer Input
10: wlw Real (Kind=nag_wp) array Workspace
11: lw Integer Input
These arguments are present so that lsqfun will be of the form required by e04hef. lsqfun is called with e04ybf's arguments iw, liw, w, lw as these arguments. If the recommendation in e04hef is followed, you will have no reason to examine or change the elements of iw or w. In any case, lsqfun must not change the first 5×n+m+m×n+n×n-1/2 (or 6+2×m if n=1) elements of w.
lsqfun must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which e04ybf is called. Arguments denoted as Input must not be changed by this procedure.
Note: lsqfun should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by e04ybf. If your code inadvertently does return any NaNs or infinities, e04ybf is likely to produce unexpected results.
Note:  e04yaf should be used to check the first derivatives calculated by lsqfun before e04ybf is used to check the bjk since e04ybf assumes that the first derivatives are correct.
4: lsqhes Subroutine, supplied by the user. External Procedure
lsqhes must calculate the elements of the symmetric matrix
Bx=i=1mfixGix,  
at any point x, where Gix is the Hessian matrix of fix. (As with lsqfun, an argument can be set to cause immediate termination.)
The specification of lsqhes is:
Fortran Interface
Subroutine lsqhes ( iflag, m, n, fvec, xc, b, lb, iw, liw, w, lw)
Integer, Intent (In) :: m, n, lb, liw, lw
Integer, Intent (Inout) :: iflag, iw(liw)
Real (Kind=nag_wp), Intent (In) :: fvec(m), xc(n)
Real (Kind=nag_wp), Intent (Inout) :: w(lw)
Real (Kind=nag_wp), Intent (Out) :: b(lb)
C Header Interface
void  lsqhes_ (Integer *iflag, const Integer *m, const Integer *n, const double fvec[], const double xc[], double b[], const Integer *lb, Integer iw[], const Integer *liw, double w[], const Integer *lw)
1: iflag Integer Input/Output
On entry: is set to a non-negative number.
On exit: if lsqhes resets iflag to some negative number, e04ybf will terminate immediately, with ifail set to your setting of iflag.
2: m Integer Input
On entry: the numbers m of residuals.
3: n Integer Input
On entry: the numbers n of variables.
4: fvecm Real (Kind=nag_wp) array Input
On entry: the value of the residual fi at the point x, for i=1,2,,m, so that the values of the fi can be used in the calculation of the elements of b.
5: xcn Real (Kind=nag_wp) array Input
On entry: the point x at which the elements of b are to be evaluated.
6: blb Real (Kind=nag_wp) array Output
On exit: unless iflag is reset to a negative number b must contain the lower triangle of the matrix Bx, evaluated at the point in xc, stored by rows. (The upper triangle is not needed because the matrix is symmetric.) More precisely, bjj-1/2+k must contain i=1mfi 2fi xjxk evaluated at the point x, for j=1,2,,n and k=1,2,,j.
7: lb Integer Input
On entry: gives the length of the array b.
8: iwliw Integer array Workspace
9: liw Integer Input
10: wlw Real (Kind=nag_wp) array Workspace
11: lw Integer Input
As in lsqfun, these arguments correspond to the arguments iw, liw, w, lw of e04ybf. lsqhes must not change the first 5×n+m×n+n×n-1/2 (or 6+2×m if n=1) elements of w.
lsqhes must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which e04ybf is called. Arguments denoted as Input must not be changed by this procedure.
Note: lsqhes should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by e04ybf. If your code inadvertently does return any NaNs or infinities, e04ybf is likely to produce unexpected results.
5: xn Real (Kind=nag_wp) array Input
On entry: xj, for j=1,2,,n, must be set to the coordinates of a suitable point at which to check the bjk calculated by lsqhes. ‘Obvious’ settings, such as 0 or 1, should not be used since, at such particular points, incorrect terms may take correct values (particularly zero), so that errors could go undetected. For a similar reason, it is preferable that no two elements of x should have the same value.
6: fvecm Real (Kind=nag_wp) array Output
On exit: unless you set iflag negative in the first call of lsqfun, fveci contains the value of fi at the point supplied by you in x, for i=1,2,,m.
7: fjacldfjacn Real (Kind=nag_wp) array Output
On exit: unless you set iflag negative in the first call of lsqfun, fjacij contains the value of the first derivative fi xj at the point given in x, as calculated by lsqfun, for i=1,2,,m and j=1,2,,n.
8: ldfjac Integer Input
On entry: the first dimension of the array fjac as declared in the (sub)program from which e04ybf is called.
Constraint: ldfjacm.
9: blb Real (Kind=nag_wp) array Output
On exit: unless you set iflag negative in lsqhes, bj×j-1/2+k contains the value of bjk at the point given in x as calculated by lsqhes, for j=1,2,,n and k=1,2,,j.
10: lb Integer Input
On entry: the dimension of the array b as declared in the (sub)program from which e04ybf is called.
Constraint: lbn+1×n/2.
11: iwliw Integer array Workspace
12: liw Integer Input
This array appears in the argument list purely so that, if e04ybf is called by another library routine, the library routine can pass quantities to user-supplied subroutines lsqfun and lsqhes via iw. iw is not examined or changed by e04ybf. In general you must provide an array iw, but are advised not to use it.
On entry: the first dimension of the array iw as declared in the (sub)program from which e04ybf is called. The actual length of iw as declared in the subroutine from which e04ybf is called.
Constraint: liw1.
13: wlw Real (Kind=nag_wp) array Workspace
14: lw Integer Input
On entry: the first dimension of the array w as declared in the (sub)program from which e04ybf is called. The actual length of w as declared in the subroutine from which e04ybf is called.
Constraints:
  • if n>1, lw5×n+m+m×n+n×n-1/2;
  • if n=1, lw6+2×m.
15: ifail Integer Input/Output
On entry: ifail must be set to 0, -1 or 1 to set behaviour on detection of an error; these values have no effect when no error is detected.
A value of 0 causes the printing of an error message and program execution will be halted; otherwise program execution continues. A value of -1 means that an error message is printed while a value of 1 means that it is not.
If halting is not appropriate, the value -1 or 1 is recommended. If message printing is undesirable, then the value 1 is recommended. Otherwise, the value -1 is recommended since useful values can be provided in some output arguments even when ifail0 on exit. When the value -1 or 1 is used it is essential to test the value of ifail on exit.
On exit: ifail=0 unless the routine detects an error or a warning has been flagged (see Section 6).

6 Error Indicators and Warnings

If on entry ifail=0 or -1, explanatory error messages are output on the current error message unit (as defined by x04aaf).
Errors or warnings detected by the routine:
Note: in some cases e04ybf may return useful information.
ifail=1
On entry, lb=value and n+1×n/2=value.
Constraint: lbn+1×n/2.
On entry, ldfjac=value and m=value.
Constraint: ldfjacm.
On entry, liw=value.
Constraint: liw1.
On entry, m=value and n=value.
Constraint: mn.
On entry, n=value.
Constraint: n1.
On entry, n=1 and lw=value.
Constraint: if n=1 then lw6+2×m; that is, value.
On entry, n>1 and lw=value.
Constraint: if n>1 then lw5×n+m+m×n+n-1×n/2; that is, value.
ifail=2
It is very likely that you have made an error in setting up b in lsqhes.
ifail<0
User requested termination by setting iflag negative in lsqfun or lsqhes.
ifail=-99
An unexpected error has been triggered by this routine. Please contact NAG.
See Section 7 in the Introduction to the NAG Library FL Interface for further information.
ifail=-399
Your licence key may have expired or may not have been installed correctly.
See Section 8 in the Introduction to the NAG Library FL Interface for further information.
ifail=-999
Dynamic memory allocation failed.
See Section 9 in the Introduction to the NAG Library FL Interface for further information.

7 Accuracy

ifail is set to 2 if
yTGy-phyTGy+1.0 or zTGz-qhzTGz+1.0  
where h is set equal to ε (ε being the machine precision as given by x02ajf) and other quantities are defined as in Section 3.

8 Parallelism and Performance

e04ybf 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.

9 Further Comments

e04ybf calls lsqhes once and lsqfun three times.

10 Example

Suppose that it is intended to use e04hef to find least squares estimates of x1,x2 and x3 in the model
y = x1 + t1 x2 t2+ x3 t3  
using the 15 sets of data given in the following table.
y t1 t2 t3 0.14 1.0 15.0 1.0 0.18 2.0 14.0 2.0 0.22 3.0 13.0 3.0 0.25 4.0 12.0 4.0 0.29 5.0 11.0 5.0 0.32 6.0 10.0 6.0 0.35 7.0 9.0 7.0 0.39 8.0 8.0 8.0 0.37 9.0 7.0 7.0 0.58 10.0 6.0 6.0 0.73 11.0 5.0 5.0 0.96 12.0 4.0 4.0 1.34 13.0 3.0 3.0 2.10 14.0 2.0 2.0 4.39 15.0 1.0 1.0  
This example program could be used to check the bjk calculated by lsqhes required. (The call of e04ybf is preceded by a call of e04yaf to check lsqfun which calculates the first derivatives.)

10.1 Program Text

Program Text (e04ybfe.f90)

10.2 Program Data

Program Data (e04ybfe.d)

10.3 Program Results

Program Results (e04ybfe.r)