E04YBF : NAG Library, Mark 24

3 Description

Routines for minimizing a sum of squares of

m

nonlinear functions (or ‘residuals’),

f_{i} (x_{1}, x_{2}, \dots, x_{n})

, for

i = 1, 2, \dots, m

and

m \geq n

, may require you to supply a subroutine to evaluate the quantities

b_{j k} = \sum_{i = 1}^{m} f_{i} \frac{\partial^{2} f_{i}}{\partial x_{j} \partial x_{k}}

for

j = 1, 2, \dots, n

and

k = 1, 2, \dots, j

. E04YBF is designed to check the

b_{j k}

calculated by such user-supplied subroutines. As well as the routine to be checked (LSQHES), you must supply a subroutine (LSQFUN) to evaluate the

f_{i}

and their first derivatives, and a point

x = {(x_{1}, x_{2}, \dots, x_{n})}^{T}

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

b_{j k}

x

. Let

J

denote the

m

n

matrix of first derivatives of the residuals. The Hessian matrix of the sum of squares,

G = J^{T} J + B,

is calculated and projected onto two orthogonal vectors

y

and

z

to give the scalars

y^{T} G y

and

z^{T} G z

respectively. The same projections of the Hessian matrix are also estimated by finite differences, giving

\begin{matrix} p = (y^{T} g (x + h y) - y^{T} g (x)) / h and \\ q = (z^{T} g (x + h z) - z^{T} g (x)) / h \end{matrix}

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

y^{T} G y

or between

q

and

z^{T} G z

is judged too large, an error indicator is set.

5 Parameters

1: M – INTEGERInput

2: N – INTEGERInput

On entry: the number

m

of residuals,

f_{i} (x)

, and the number

n

of variables,

x_{j}

Constraint:

1 \leq N \leq M

3: LSQFUN – SUBROUTINE, supplied by the user.External Procedure

LSQFUN must calculate the vector of values

f_{i} (x)

and their first derivatives

\frac{\partial f_{i}}{\partial x_{j}}

at any point

x

. (E04HEF gives you the option of resetting parameters of LSQFUN to cause the minimization process to terminate immediately. E04YBF will also terminate immediately, without finishing the checking process, if the parameter in question is reset.)

The specification of LSQFUN is:

SUBROUTINE LSQFUN (

IFLAG, M, N, XC, FVEC, FJAC, LDFJAC, IW, LIW, W, LW)

INTEGER	IFLAG, M, N, LDFJAC, IW(LIW), LIW, LW
REAL (KIND=nag_wp)	XC(N), FVEC(M), FJAC(LDFJAC,N), W(LW)

1: IFLAG – INTEGERInput/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 – INTEGERInput: On entry: the numbers $m$ of residuals.
3: N – INTEGERInput: On entry: the numbers $n$ of variables.
4: XC(N) – REAL (KIND=nag_wp) arrayInput: On entry: the point $x$ at which the values of the $f_{i}$ and the $\frac{\partial f_{i}}{\partial x_{j}}$ are required.
5: FVEC(M) – REAL (KIND=nag_wp) arrayOutput: On exit: unless IFLAG is reset to a negative number, $FVEC (i)$ must contain the value of $f_{i}$ at the point $x$ , for $i = 1, 2, \dots, m$ .
6: FJAC(LDFJAC,N) – REAL (KIND=nag_wp) arrayOutput: On exit: unless IFLAG is reset to a negative number, $FJAC (i, j)$ must contain the value of $\frac{\partial f_{i}}{\partial x_{j}}$ at the point $x$ , for $i = 1, 2, \dots, m$ and $j = 1, 2, \dots, n$ .
7: LDFJAC – INTEGERInput: On entry: the first dimension of the array FJAC as declared in the (sub)program from which E04YBF is called.
8: IW(LIW) – INTEGER arrayWorkspace
9: LIW – INTEGERInput
10: W(LW) – REAL (KIND=nag_wp) arrayWorkspace
11: LW – INTEGERInput: These parameters are present so that LSQFUN will be of the form required by E04HEF. LSQFUN is called with E04YBF's parameters IW, LIW, W, LW as these parameters. 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 \times N + M + M \times N + N \times (N - 1) / 2$ (or $6 + 2 \times 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. Parameters denoted as Input must not be changed by this procedure.

Note: E04YAF should be used to check the first derivatives calculated by LSQFUN before E04YBF is used to check the

b_{j k}

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

B (x) = \sum_{i = 1}^{m} f_{i} (x) G_{i} (x),

at any point

x

, where

G_{i} (x)

is the Hessian matrix of

f_{i} (x)

. (As with LSQFUN, a parameter can be set to cause immediate termination.)

The specification of LSQHES is:

SUBROUTINE LSQHES (

IFLAG, M, N, FVEC, XC, B, LB, IW, LIW, W, LW)

INTEGER	IFLAG, M, N, LB, IW(LIW), LIW, LW
REAL (KIND=nag_wp)	FVEC(M), XC(N), B(LB), W(LW)

1: IFLAG – INTEGERInput/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 – INTEGERInput: On entry: the numbers $m$ of residuals.
3: N – INTEGERInput: On entry: the numbers $n$ of variables.
4: FVEC(M) – REAL (KIND=nag_wp) arrayInput: On entry: the value of the residual $f_{i}$ at the point $x$ , for $i = 1, 2, \dots, m$ , so that the values of the $f_{i}$ can be used in the calculation of the elements of B.
5: XC(N) – REAL (KIND=nag_wp) arrayInput: On entry: the point $x$ at which the elements of B are to be evaluated.
6: B(LB) – REAL (KIND=nag_wp) arrayOutput: On exit: unless IFLAG is reset to a negative number B must contain the lower triangle of the matrix $B (x)$ , evaluated at the point in XC, stored by rows. (The upper triangle is not needed because the matrix is symmetric.) More precisely, $B (j (j - 1) / 2 + k)$ must contain $\sum_{i = 1}^{m} f_{i} \frac{\partial^{2} f_{i}}{\partial x_{j} \partial x_{k}}$ evaluated at the point $x$ , for $j = 1, 2, \dots, n$ and $k = 1, 2, \dots, j$ .
7: LB – INTEGERInput: On entry: gives the length of the array B.
8: IW(LIW) – INTEGER arrayWorkspace
9: LIW – INTEGERInput
10: W(LW) – REAL (KIND=nag_wp) arrayWorkspace
11: LW – INTEGERInput: As in LSQFUN, these parameters correspond to the parameters IW, LIW, W, LW of E04YBF. LSQHES must not change the first $5 \times N + M \times N + N \times (N - 1) / 2$ (or $6 + 2 \times 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. Parameters denoted as Input must not be changed by this procedure.

5: X(N) – REAL (KIND=nag_wp) arrayInput

On entry:

X (j)

, for

j = 1, 2, \dots, n

, must be set to the coordinates of a suitable point at which to check the

b_{j k}

calculated by LSQHES. ‘Obvious’ settings, such as

0

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: FVEC(M) – REAL (KIND=nag_wp) arrayOutput

On exit: unless you set IFLAG negative in the first call of LSQFUN,

FVEC (i)

contains the value of

f_{i}

at the point supplied by you in X, for

i = 1, 2, \dots, m

7: FJAC(LDFJAC,N) – REAL (KIND=nag_wp) arrayOutput

On exit: unless you set IFLAG negative in the first call of LSQFUN,

FJAC (i, j)

contains the value of the first derivative

\frac{\partial f_{i}}{\partial x_{j}}

at the point given in X, as calculated by LSQFUN, for

i = 1, 2, \dots, m

and

j = 1, 2, \dots, n

8: LDFJAC – INTEGERInput

On entry: the first dimension of the array FJAC as declared in the (sub)program from which E04YBF is called.

Constraint:

LDFJAC \geq M

9: B(LB) – REAL (KIND=nag_wp) arrayOutput

On exit: unless you set IFLAG negative in LSQHES,

B (j \times (j - 1) / 2 + k)

contains the value of

b_{j k}

at the point given in X as calculated by LSQHES, for

j = 1, 2, \dots, n

and

k = 1, 2, \dots, j

10: LB – INTEGERInput

On entry: the dimension of the array B as declared in the (sub)program from which E04YBF is called.

Constraint:

LB \geq (N + 1) \times N / 2

11: IW(LIW) – INTEGER arrayWorkspace

This array appears in the parameter 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.

12: LIW – INTEGERInput

On entry: the actual length of IW as declared in the subroutine from which E04YBF is called.

Constraint:

LIW \geq 1

13: W(LW) – REAL (KIND=nag_wp) arrayWorkspace

14: LW – INTEGERInput

On entry: the actual length of W as declared in the subroutine from which E04YBF is called.

Constraints:

if $N > 1$ , $LW \geq 5 \times N + M + M \times N + N \times (N - 1) / 2$ ;
if $N = 1$ , $LW \geq 6 + 2 \times M$ .

15: IFAIL – INTEGERInput/Output

On entry: IFAIL must be set to

0

- 1 ​ or ​ 1

. 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

- 1 ​ or ​ 1

is recommended. If the output of error messages is undesirable, then the value

1

is recommended. Otherwise, because for this routine the values of the output parameters may be useful even if

IFAIL \neq 0

on exit, the recommended value is

- 1

. 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

- 1

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

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

Errors or warnings detected by the routine:

IFAIL < 0

A negative value of IFAIL indicates an exit from E04YBF because you have set IFLAG negative in user-supplied subroutines LSQFUN or LSQHES. The setting of IFAIL will be the same as your setting of IFLAG. The check on LSQHES will not have been completed.

IFAIL = 1

On entry,	$M < N$ ,
or	$N < 1$ ,
or	$LDFJAC < M$ ,
or	$LB < (N + 1) \times N / 2$ ,
or	$LIW < 1$ ,
or	$LW < 5 \times N + M + M \times N + N \times (N - 1) / 2$ , if $N > 1$ ,
or	$LW < 6 + 2 \times M$ , if $N = 1$ .

IFAIL = 2

You should check carefully the derivation and programming of expressions for the

b_{j k}

, because it is very unlikely that LSQHES is calculating them correctly.

9 Example

Suppose that it is intended to use E04HEF to find least squares estimates of

x_{1}, x_{2}

and

x_{3}

in the model

y = x_{1} + \frac{t_{1}}{x_{2} t_{2} + x_{3} t_{3}}

using the

15

sets of data given in the following table.

\begin{array}{r} y & t_{1} & t_{2} & t_{3} \\ 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 \end{array}

This example program could be used to check the

b_{j k}

calculated by LSQHES required. (The call of E04YBF is preceded by a call of E04YAF to check LSQFUN which calculates the first derivatives.)

9.1 Program Text

9.2 Program Data

9.3 Program Results

NAG Library Routine Document

E04YBF

+− Contents

1 Purpose

2 Specification

3 Description

4 References

5 Parameters

6 Error Indicators and Warnings

7 Accuracy

8 Further Comments

9 Example

9.1 Program Text

9.2 Program Data

9.3 Program Results

NAG Library Routine DocumentE04YBF

+− Contents

1 Purpose

2 Specification

3 Description

4 References

5 Parameters

6 Error Indicators and Warnings

7 Accuracy

8 Further Comments

9 Example

9.1 Program Text

9.2 Program Data

9.3 Program Results

NAG Library Routine Document

E04YBF