nag_opt_lsq_deriv (e04gbc) : NAG Library, Mark 26

nag_opt_lsq_deriv (e04gbc) is a comprehensive algorithm for finding an unconstrained minimum of a sum of squares of

m

nonlinear functions in

n

variables

(m \geq n)

. First derivatives are required.

nag_opt_lsq_deriv (e04gbc) is intended for objective functions which have continuous first and second derivatives (although it will usually work even if the derivatives have occasional discontinuities).

nag_opt_lsq_deriv (e04gbc) is applicable to problems of the form:

Minimize ​ F (x) = \sum_{i = 1}^{m} {[f_{i} (x)]}^{2}

where

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

and

m \geq n

. (The functions

f_{i} (x)

are often referred to as ‘residuals’.) You must supply a function to calculate the values of the

f_{i} (x)

and their first derivatives

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

at any point

x

.

From a starting point

x^{(1)}

nag_opt_lsq_deriv (e04gbc) generates a sequence of points

x^{(2)}, x^{(3)}, \dots,

which is intended to converge to a local minimum of

F (x)

. The sequence of points is given by

x^{(k + 1)} = x^{(k)} + α^{(k)} p^{(k)}

where the vector

p^{(k)}

is a direction of search, and

α^{(k)}

is chosen such that

F (x^{(k)} + α^{(k)} p^{(k)})

is approximately a minimum with respect to

α^{(k)}

.

The vector

p^{(k)}

used depends upon the reduction in the sum of squares obtained during the last iteration. If the sum of squares was sufficiently reduced, then

p^{(k)}

is the Gauss–Newton direction; otherwise the second derivatives of the

f_{i} (x)

are taken into account using a quasi-Newton updating scheme.

The method is designed to ensure that steady progress is made whatever the starting point, and to have the rapid ultimate convergence of Newton's method.

Gill P E and Murray W (1978) Algorithms for the solution of the nonlinear least squares problem SIAM J. Numer. Anal. 15 977–992

Intermediate and final results are printed out by default. The level of printed output can be controlled with the option

options . print_level

(see Section 11.2). The default,

options . print_level = Nag_Soln_Iter

, provides a single line of output at each iteration and the final result. The line of results printed at each iteration gives:

Itn	the current iteration number $k$ .
Nfun	the cumulative number of calls to lsqfun.
Objective	the current value of the objective function, $F (x^{(k)})$ .
Norm g	the Euclidean norm of the gradient of $F (x^{(k)})$ .
Norm x	the Euclidean norm of $x^{(k)}$ .
Norm(x(k-1)-x(k))	the Euclidean norm of $x^{(k - 1)} - x^{(k)}$ .
Step	the step $α^{(k)}$ taken along the computed search direction $p^{(k)}$ .

The printout of the final result consists of:

x	the final point $x^{*}$ .
g	the gradient of $F$ at the final point.
Residuals	the values of the residuals $f_{i}$ at the final point.
Sum of squares	the value of $F (x^{*})$ , the sum of squares of the residuals at the final point.

If the problem is reasonably well scaled and a successful exit is made, then, for a computer with a mantissa of

t

decimals, one would expect to get about

t / 2 - 1

decimals accuracy in the components of

x

and between

t - 1

(if

F (x)

is of order 1 at the minimum) and

2 t - 2

(if

F (x)

is close to zero at the minimum) decimals accuracy in

F (x)

.

A successful exit (

fail . code = NE_NOERROR

) is made from nag_opt_lsq_deriv (e04gbc) when (B1, B2 and B3) or B4 or B5 hold, where

\begin{array}{l} B1 \equiv α^{(k)} \times ‖p^{(k)}‖ < (options . optim_tol + ε) \times (1.0 + ‖x^{(k)}‖) \\ B2 \equiv |F^{(k)} - F^{(k - 1)}| < {(options . optim_tol + ε)}^{2} \times (1.0 + F^{(k)}) \\ B3 \equiv ‖g^{(k)}‖ < ε^{1 / 3} \times (1.0 + F^{(k)}) \\ B4 \equiv F^{(k)} < ε^{2} \\ B5 \equiv ‖g^{(k)}‖ < {(ε \times \sqrt{F^{(k)}})}^{1 / 2} \end{array}

and where

‖.‖

,

ε

and the optional parameter

options . optim_tol

are as defined in Section 11.2, while

F^{(k)}

and

g^{(k)}

are the values of

F (x)

and its vector of first derivatives at

x^{(k)}

.

If

fail . code = NE_NOERROR

then the vector in x on exit,

x_{sol}

, is almost certainly an estimate of

x_{true}

, the position of the minimum to the accuracy specified by

options . optim_tol

.

If

fail . code = NW_COND_MIN

, then

x_{sol}

may still be a good estimate of

x_{true}

, but to verify this you should make the following checks. If

(a)	the sequence $\{F (x^{(k)})\}$ converges to $F (x_{sol})$ at a superlinear or a fast linear rate, and
(b)	$g {(x_{sol})}^{T} g (x_{sol}) < 10 ε$ ,

where

T

denotes transpose, then it is almost certain that

x_{sol}

is a close approximation to the minimum. When (b) is true, then usually

F (x_{sol})

is a close approximation to

F (x_{true})

.

Further suggestions about confirmation of a computed solution are given in the e04 Chapter Introduction.

nag_opt_lsq_deriv (e04gbc) is not threaded in any implementation.

The number of iterations required depends on the number of variables, the number of residuals, the behaviour of

F (x)

, the accuracy demanded and the distance of the starting point from the solution. The number of multiplications performed per iteration of nag_opt_lsq_deriv (e04gbc) varies, but for

m > > n

is approximately

n \times m^{2} + O (n^{3})

. In addition, each iteration makes at least one call of lsqfun. So, unless the residuals can be evaluated very quickly, the run time will be dominated by the time spent in lsqfun.

Ideally, the problem should be scaled so that, at the solution,

F (x)

and the corresponding values of the

x_{j}

are each in the range

(- 1, + 1)

, and so that at points one unit away from the solution,

F (x)

differs from its value at the solution by approximately one unit. This will usually imply that the Hessian matrix of

F (x)

at the solution is well-conditioned. It is unlikely that you will be able to follow these recommendations very closely, but it is worth trying (by guesswork), as sensible scaling will reduce the difficulty of the minimization problem, so that nag_opt_lsq_deriv (e04gbc) will take less computer time.

When the sum of squares represents the goodness-of-fit of a nonlinear model to observed data, elements of the variance-covariance matrix of the estimated regression coefficients can be computed by a subsequent call to nag_opt_lsq_covariance (e04ycc), using information returned in the arrays

options . s

and

options . v

. See nag_opt_lsq_covariance (e04ycc) for further details.

This example finds the 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.

$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

The program uses (0.5, 1.0, 1.5) as the initial guess at the position of the minimum.

The program shows the use of certain optional parameters, with some option values being assigned directly within the program text and by reading values from a data file. The options structure is declared and initialized by nag_opt_init (e04xxc). A value is then assigned directly to options

options . outfile

and three further options are read from the data file by use of nag_opt_read (e04xyc). The memory freeing function nag_opt_free (e04xzc) is used to free the memory assigned to the pointers in the option structure. You must not use the standard C function free() for this purpose.

A number of optional input and output arguments to nag_opt_lsq_deriv (e04gbc) are available through the structure argument options, type Nag_E04_Opt. An argument may be selected by assigning an appropriate value to the relevant structure member; those arguments not selected will be assigned default values. If no use is to be made of any of the optional parameters you should use the NAG defined null pointer, E04_DEFAULT, in place of options when calling nag_opt_lsq_deriv (e04gbc); the default settings will then be used for all arguments.

Before assigning values to options directly the structure must be initialized by a call to the function nag_opt_init (e04xxc). Values may then be assigned to the structure members in the normal C manner.

Optional parameter settings may also be read from a text file using the function nag_opt_read (e04xyc) in which case initialization of the options structure will be performed automatically if not already done. Any subsequent direct assignment to the options structure must not be preceded by initialization.

If assignment of functions and memory to pointers in the options structure is required, this must be done directly in the calling program. They cannot be assigned using nag_opt_read (e04xyc).

For easy reference, the following list shows the members of options which are valid for nag_opt_lsq_deriv (e04gbc) together with their default values where relevant. The number

ε

is a generic notation for machine precision (see nag_machine_precision (X02AJC)).

Boolean list	Nag_TRUE
Nag_PrintType print_level	$Nag_Soln_Iter$
char outfile[80]	stdout
void (*print_fun)()	NULL
Boolean deriv_check	Nag_TRUE
Integer max_iter	$\max (50, 5 n)$
double optim_tol	$\sqrt{ε}$
Nag_LinFun minlin	$Nag_Lin_Deriv$
double linesearch_tol	0.9 (0.0 if $n = 1$ )
double step_max	100000.0
double *s	size n
double *v	size $n \times n$
Integer tdv	n
Integer grade
Integer iter
Integer nf

On entry: if

options . list = Nag_TRUE

the argument settings in the call to nag_opt_lsq_deriv (e04gbc) will be printed.

On entry: the level of results printout produced by nag_opt_lsq_deriv (e04gbc). The following values are available:

$Nag_NoPrint$	No output.
$Nag_Soln$	The final solution.
$Nag_Iter$	One line of output for each iteration.
$Nag_Soln_Iter$	The final solution and one line of output for each iteration.
$Nag_Soln_Iter_Full$	The final solution and detailed printout at each iteration.

Details of each level of results printout are described in Section 11.3.

Constraint:

options . print_level = Nag_NoPrint

,

Nag_Soln

,

Nag_Iter

,

Nag_Soln_Iter

or

Nag_Soln_Iter_Full

.

On entry: the name of the file to which results should be printed. If

options . outfile [0] ='
  \0
'

then the stdout stream is used.

On entry: printing function defined by you; the prototype of

options . print_fun

is

void (*print_fun)(const Nag_Search_State *st, Nag_Comm *comm);

See Section 11.3.1 for further details.

On entry: if

options . deriv_check = Nag_TRUE

a check of the derivatives defined by lsqfun will be made at the starting point x. The derivative check is carried out by a call to nag_opt_lsq_check_deriv (e04yac). A starting point of

x = 0

or

x = 1

should be avoided if this test is to be meaningful, but if either of these starting points is necessary then nag_opt_lsq_check_deriv (e04yac) should be used to check lsqfun at a different point prior to calling nag_opt_lsq_deriv (e04gbc).

On entry: the limit on the number of iterations allowed before termination.

Constraint:

options . max_iter \geq 0

.

On entry: the accuracy in

x

to which the solution is required. If

x_{true}

is the true value of

x

at the minimum, then

x_{sol}

, the estimated position prior to a normal exit, is such that

‖x_{sol} - x_{true}‖ < options . optim_tol \times (1.0 + ‖x_{true}‖),

where

‖y‖ = \sqrt{\sum_{j = 1}^{n} y_{j}^{2}}

. For example, if the elements of

x_{sol}

are not much larger than 1.0 in modulus and if

options . optim_tol = 1.0 \times 10^{- 5}

, then

x_{sol}

is usually accurate to about five decimal places. (For further details see Section 7.) If

F (x)

and the variables are scaled roughly as described in Section 9 and

ε

is the machine precision, then a setting of order

options . optim_tol = \sqrt{ε}

will usually be appropriate.

Constraint:

10 ε \leq options . optim_tol < 1.0

.

On entry:

options . minlin

specifies whether the linear minimizations (i.e., minimizations of

F (x^{(k)} + α^{(k)} p^{(k)})

with respect to

α^{(k)}

) are to be performed by a function which just requires the evaluation of the

f_{i} (x)

,

Nag_Lin_NoDeriv

, or by a function which also requires the first derivatives of the

f_{i} (x)

,

Nag_Lin_Deriv

.

It will often be possible to evaluate the first derivatives of the residuals in about the same amount of computer time that is required for the evaluation of the residuals themselves – if this is so then nag_opt_lsq_deriv (e04gbc) should be called with

options . minlin

set to

Nag_Lin_Deriv

. However, if the evaluation of the derivatives takes more than about four times as long as the evaluation of the residuals, then a setting of

Nag_Lin_NoDeriv

will usually be preferable. If in doubt, use the default setting

Nag_Lin_Deriv

as it is slightly more robust.

Constraint:

options . minlin = Nag_Lin_Deriv

or

Nag_Lin_NoDeriv

.

If

options . minlin = Nag_Lin_NoDeriv

then the default value of

options . linesearch_tol

will be changed from 0.9 to 0.5 if

n > 1

.

On entry:

options . linesearch_tol

specifies how accurately the linear minimizations are to be performed.

Every iteration of nag_opt_lsq_deriv (e04gbc) involves a linear minimization, i.e., minimization of

F (x^{(k)} + α^{(k)} p^{(k)})

with respect to

α^{(k)}

. The minimum with respect to

α^{(k)}

will be located more accurately for small values of

options . linesearch_tol

(say 0.01) than for large values (say 0.9). Although accurate linear minimizations will generally reduce the number of iterations performed by nag_opt_lsq_deriv (e04gbc), they will increase the number of calls of lsqfun made each iteration. On balance it is usually more efficient to perform a low accuracy minimization.

Constraint:

0.0 \leq options . linesearch_tol < 1.0

.

On entry: an estimate of the Euclidean distance between the solution and the starting point supplied. (For maximum efficiency, a slight overestimate is preferable.) nag_opt_lsq_deriv (e04gbc) will ensure that, for each iteration,

\sum_{j = 1}^{n} {(x_{j}^{(k)} - x_{j}^{(k - 1)})}^{2} \leq {(options . step_max)}^{2}

where

k

is the iteration number. Thus, if the problem has more than one solution, nag_opt_lsq_deriv (e04gbc) is most likely to find the one nearest to the starting point. On difficult problems, a realistic choice can prevent the sequence

x^{(k)}

entering a region where the problem is ill-behaved and can help avoid overflow in the evaluation of

F (x)

. However, an underestimate of

options . step_max

can lead to inefficiency.

Constraint:

options . step_max \geq options . optim_tol

.

On entry: n values of memory will be automatically allocated by nag_opt_lsq_deriv (e04gbc) and this is the recommended method of use of

options . s

. However, you may supply memory from the calling program.

On exit: the singular values of the Jacobian matrix at the final point. Thus

options . s

may be useful as information about the structure of your problem.

On entry:

n \times n

values of memory will be automatically allocated by nag_opt_lsq_deriv (e04gbc) and this is the recommended method of use of

options . v

. However, you may supply memory from the calling program.

On exit: the matrix

V

associated with the singular value decomposition

J = {U S V}^{T}

of the Jacobian matrix at the final point, stored by rows. This matrix may be useful for statistical purposes, since it is the matrix of orthonormalized eigenvectors of

J^{T} J

.

On entry: if memory is supplied then

options . tdv

must contain the last dimension of the array assigned to

options . tdv

as declared in the function from which nag_opt_lsq_deriv (e04gbc) is called.

On exit: the trailing dimension used by

options . v

. If the NAG default memory allocation has been used this value will be n.

Constraint:

options . tdv \geq n

.

On exit: the grade of the Jacobian at the final point. nag_opt_lsq_deriv (e04gbc) estimates the dimension of the subspace for which the Jacobian matrix can be used as a valid approximation to the curvature (see Gill and Murray (1978)); this estimate is called the grade.

On exit: the number of iterations which have been performed in nag_opt_lsq_deriv (e04gbc).

On exit: the number of times the residuals have been evaluated (i.e., the number of calls of lsqfun).

The level of printed output can be controlled with the structure members

options . list

and

options . print_level

(see Section 11.2). If

options . list = Nag_TRUE

then the argument values to nag_opt_lsq_deriv (e04gbc) are listed, whereas the printout of results is governed by the value of

options . print_level

. The default of

options . print_level = Nag_Soln_Iter

provides a single line of output at each iteration and the final result. This section describes all of the possible levels of results printout available from nag_opt_lsq_deriv (e04gbc).

When

options . print_level = Nag_Iter

or

Nag_Soln_Iter

a single line of output is produced on completion of each iteration, this gives the following values:

Itn	the current iteration number $k$ .
Nfun	the cumulative number of calls to lsqfun.
Objective	the value of the objective function, $F (x^{(k)})$ .
Norm g	the Euclidean norm of the gradient of $F (x^{(k)})$ .
Norm x	the Euclidean norm of $x^{(k)}$ .
Norm(x(k-1)-x(k))	the Euclidean norm of $x^{(k - 1)} - x^{(k)}$ .
Step	the step $α^{(k)}$ taken along the computed search direction $p^{(k)}$ .

When

options . print_level = Nag_Soln_Iter_Full

more detailed results are given at each iteration. Additional values output are:

Grade	the grade of the Jacobian matrix. (See description of $options . grade$ , Section 9.)
x	the current point $x^{(k)}$ .
g	the current gradient of $F (x^{(k)})$ .
Singular values	the singular values of the current approximation to the Jacobian matrix.

If

options . print_level = Nag_Soln

,

Nag_Soln_Iter

or

Nag_Soln_Iter_Full

the final result consists of:

x	the final point $x^{*}$ .
g	the gradient of $F$ at the final point.
Residuals	the values of the residuals $f_{i}$ at the final point.
Sum of squares	the value of $F (x^{*})$ , the sum of squares of the residuals at the final point.

If

options . print_level = Nag_NoPrint

then printout will be suppressed; you can print the final solution when nag_opt_lsq_deriv (e04gbc) returns to the calling program.

You may also specify your own print function for output of iteration results and the final solution by use of the

options . print_fun

function pointer, which has prototype

void (*print_fun)(const Nag_Search State *st, Nag_Comm *comm);

The rest of this section can be skipped if the default printing facilities provide the required functionality.

When a user-defined function is assigned to

options . print_fun

this will be called in preference to the internal print function of nag_opt_lsq_deriv (e04gbc). Calls to the user-defined function are again controlled by means of the

options . print_level

member. Information is provided through st and comm, the two structure arguments to

options . print_fun

. The structure member

comm \to it_prt

is relevant in this context. If

comm \to it_prt = Nag_TRUE

then the results from the last iteration of nag_opt_lsq_deriv (e04gbc) are in the following members of st:

The relevant members of the structure comm are:

NAG Library Function Document

nag_opt_lsq_deriv (e04gbc)

▸▿ Contents

1 Purpose

2 Specification

3 Description

4 References

5 Arguments

5.1 Description of Printed Output

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 Optional Parameters

11.1 Optional Parameter Checklist and Default Values

11.2 Description of the Optional Parameters

11.3 Description of Printed Output

11.3.1 Output of results via a user-defined printing function

NAG Library Function Documentnag_opt_lsq_deriv (e04gbc)

▸▿ Contents

1 Purpose

2 Specification

3 Description

4 References

5 Arguments

5.1 Description of Printed Output

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 Optional Parameters

11.1 Optional Parameter Checklist and Default Values

11.2 Description of the Optional Parameters

11.3 Description of Printed Output

11.3.1 Output of results via a user-defined printing function

NAG Library Function Document

nag_opt_lsq_deriv (e04gbc)