How to use the package¶

Overview¶

Calling sequences¶

Access to the package requires a USE statement

use ral_nlls_double

The user can then call one of the procedures:

nlls_solve() solves the non-linear least squares problem.

nlls_iterate() performs one iteration of the non-linear least squares solver.

The calling sequences of these subroutines are outlined in Argument lists and calling sequences.

The derived data types¶

For each problem, the user must employ the derived types defined by the module to declare scalars of the types nlls_inform and nlls_options. If nlls_iterate() is to be used, then a scalar of the type nlls_workspace must also be defined. The following pseudocode illustrates this.

use nlls_module
!...
type (NLLS_inform) :: inform
type (NLLS_options) :: options
type (NLLS_workspace) :: work ! needed if nlls_iterate to be called
!...

The components of nlls_options and nlls_inform are explained below in Data types.

Argument lists and calling sequences¶

We use square brackets to indicate arguments. In each call, optional arguments follow the argument inform. Since we reserve the right to add additional optional arguments in future releases of the code, we strongly recommend that all optional arguments be called by keyword, not by position.

The term package type is used to mean double precision.

To solve the non-linear least squares problem¶

subroutine nlls_solve(n, m, X, eval_r, eval_J, eval_Hf, params, options, inform[, weights, eval_HP])¶

Solves the non-linear least squares problem.

Parameters:

n [integer,in] :: holds the number $n$ of variables to be fitted; i.e., $n$ is the length of the unknown vector $\bm x$ . Restriction: n $>$ 0.
m [integer,in] :: holds the number $m$ of data points available; i.e., $m$ is the number of residuals $r_i$ . Restriction: m $\geq$ 0.
X (n) [real,inout] :: on entry, it must hold the initial guess for $\bm x$ , and on successful exit it holds the solution to the non-linear least squares problem.
eval_r [procedure] :: given a point ${\bm x} _{k}^{}$ , returns the vector ${\bm r} ({\bm x} _{k}^{})$ . Further details of the format required are given in eval_r() in User-supplied function evaluation routines.
eval_J [procedure] :: given a point ${\bm x} _{k}^{}$ , returns the $m \times n$ Jacobian matrix, ${\bm J} _{k}^{}$ , of ${\bm r}$ evaluated at ${\bm x} _{k}^{}$ . Further details of the format required are given in eval_J() in User-supplied function evaluation routines.
eval_Hf [procedure] :: given vectors ${\bm x} \in \mathbb{R}^n$ and ${\bm r} \in \mathbb{R}^m$ , returns the quantity $\sum_{i=1}^m ( {\bm r} )_i \nabla^2 {\bm r} _i ( {\bm x} )$ . Further details of the format required are given in eval_Hf() in User-supplied function evaluation routines. If exact_second_derivative = .false. in nlls_options, then this is not referenced.
params [params_base_type,in] :: holds parameters to be passed to the user-defined routines eval_r(), eval_J(), and eval_Hf(). Further details of its use are given in User-supplied function evaluation routines.
options [nlls_options,in] :: controls execution of algorithm.
inform [nlls_inform,out] :: components provide information about the execution of the subroutine.

Options:

weights (n) [real] :: If present, this holds the square-roots of the diagonal entries of the weighting matrix, ${\bm W}$ . If absent, then the norm in the least squares problem is taken to be the 2-norm, that is, ${\bm W} = I$ .
eval_HP [procedure] :: If present, this is a routine that, given vectors ${\bm x}, {\bm y} \in \mathbb{R}^m$ , returns the matrix $P({\bm x},{\bm y}) := ( H_1({\bm x}){\bm y} \dots H_m({\bm x}){\bm y})$ . Further details of the format required are given in eval_HP() in User-supplied function evaluation routines. This is only referenced if model = 4 in nlls_options.

To iterate once¶

subroutine nlls_iterate(n, m, X, eval_r, eval_J, eval_Hf, params, options, inform[, weights])¶

A call of this form allows the user to step through the solution process one iteration at a time.

n, m, eval_F, eval_J, eval_HF, params, info, options and weights are as in the desciption of nlls_solve().

Parameters:	X (n) [real,inout] :: on the first call it must hold the initial guess for $\bm x$ . On return it holds the value of $\bm x$ at the current iterate, and must be passed unaltered to any subsequent call to `nlls_iterate()`. w [nlls_workspace,inout] :: is used to store the current state of the iteration and should not be altered by the user.

The user may use the components convergence_normf and convergence_normg and converge_norms in nlls_inform to determine whether the iteration has converged.

User-supplied function evaluation routines¶

The user must supply routines to evaluate the residual, Jacobian and Hessian at a point. RALFit will call these routines internally.

In order to pass user-defined data into the evaluation calls, params_base_type is extended to a user_type, as follows:

type, extends( params_base_type ) :: user_type
   ! code declaring components of user_type
end type user_type

We recommend this type is wrapped in a module with the user-defined routines for evaluating the function, Jacobian, and Hessian.

The components of the extended type are accessed through a select type construct:

select type(params)
type is(user_type)
  ! code that accesses components of params that were defined within user_type
end select

For evaluating the function ${\bm r} ( {\bm x} )$ ¶

A subroutine must be supplied to calculate ${\bm r} ( {\bm x} )$ for a given vector ${\bm x}$ . It must implement the following interface:

abstract interface
   subroutine eval_r(status, n, m, x, r, params)
      integer, intent(inout) :: status
      integer, intent(in) :: n
      integer, intent(in) :: m
      double precision, dimension(n), intent(in) :: x
      double precision, dimension(m), intent(out) :: r
      class(params_base_type), intent(in) :: params
   end subroutine eval_r
end interface

subroutine eval_r(status, n, m, x, r, params)¶

Parameters:

status [integer,inout] :: is initialised to 0 before the routine is called. If it is set to a non-zero value by the routine, then nlls_solve() / nlls_iterate() will exit with error.
n [integer,in] :: is passed unchanged as provided in the call to nlls_solve()/nlls_iterate().
m [integer,in] :: is passed unchanged as provided in the call to nlls_solve()/nlls_iterate().
X (n) [real,in] :: holds the current point ${\bm x}_{k}^{}$ at which to evaluate ${\bm r} ( {\bm x} _{k}^{})$ .
r (m) [real,out] :: must be set by the routine to hold the residual function evaluated at the current point ${\bm x} _{k}^{}$ , ${\bm r} ({\bm x} _{k}^{})$ .
params [params_base_type,in] :: is passed unchanged as provided in the call to nlls_solve()/nlls_iterate().

For evaluating the function ${\bm J} = \nabla {\bm r} ( {\bm x} )$ ¶

A subroutine must be supplied to calculate ${\bm J} = \nabla {\bm r} ( {\bm x} )$ for a given vector ${\bm x}$ . It must implement the following interface:

abstract interface
   subroutine eval_J(status, n, m, x, J, params)
      integer, intent(inout) :: status
      integer, intent(in) :: n
      integer, intent(in) :: m
      double precision, dimension(n), intent(in)  :: x
      double precision, dimension(n*m), intent(out) :: J
      class(params_base_type), intent(in) :: params
  end subroutine eval_J
end interface

subroutine eval_J(status, n, m, x, J, params)¶

Parameters:

status [integer,inout] :: is initialised to 0 before the routine is called. If it is set to a non-zero value by the routine, then nlls_solve() / nlls_iterate() will exit with error.
n [integer,in] :: is passed unchanged as provided in the call to nlls_solve()/nlls_iterate().
m [integer,in] :: is passed unchanged as provided in the call to nlls_solve()/nlls_iterate().
X (n) [real,in] :: holds the current point ${\bm x}_{k}^{}$ at which to evaluate ${\bm J} ( {\bm x} _{k}^{})$ .
J (m*n) [real,out] :: must be set by the routine to hold the Jacobian of the residual function evaluated at the current point ${\bm x}_{k}^{}$ , ${\bm r} ( {\bm x} _{k}^{})$ . J(i*m+j) must be set to hold $\nabla_{x_j} r_i( {\bm x} _{k}^{})$ .
params [params_base_type,in] :: is passed unchanged as provided in the call to nlls_solve()/nlls_iterate().

For evaluating the function $Hf = \sum_{i=1}^m r_i( {\bm x} ) {\bm W} \nabla^2 r_i( {\bm x} )$ ¶

A subroutine must be supplied to calculate $Hf = \sum_{i=1}^m ( {\bm r} )_i \nabla^2 r_i( {\bm x} )$ for given vectors ${\bm x} \in \mathbb{R}^n$ and ${\bm r} \in \mathbb{R}^m$ ; here $( {\bm r} )_i$ denotes the $i$ th component of the vector ${\bm r}$ . The subroutine must implement the following interface:

abstract interface
   subroutine eval_Hf_type(status, n, m, x, r, Hf, params)
       integer, intent(inout) :: status
       integer, intent(in) :: n
       integer, intent(in) :: m
       double precision, dimension(n), intent(in)  :: x
       double precision, dimension(m), intent(in)  :: r
       double precision, dimension(n*n), intent(out) :: Hf
       class(params_base_type), intent(in) :: params
     end subroutine eval_Hf_type
end interface
:language: fortran

subroutine eval_Hf(status, n, m, x, r, Hf, params)¶

Parameters:

status [integer,inout] :: is initialised to 0 before the routine is called. If it is set to a non-zero value by the routine, then nlls_solve() / nlls_iterate() will exit with error.
n [integer,in] :: is passed unchanged as provided in the call to nlls_solve()/nlls_iterate().
m [integer,in] :: is passed unchanged as provided in the call to nlls_solve()/nlls_iterate().
X (n) [real,in] :: holds the current point ${\bm x}_{k}^{}$ at which to evaluate $\sum_{i=1}^m ( {\bm r} )_i \nabla^2 r_i( {\bm x} )$ .
r (m) [real,in] :: holds ${\bm W} {\bm r} ( {\bm x} )$ , the (weighted) residual, as computed from vector returned by the last call to eval_r().
Hf (n*n) [real,out] :: must be set by the routine to hold the matrix $\sum_{i = 1}^m ( {\bm r} )_{i}\nabla^2 r_{i}^{}( {\bm x} _{k}^{})$ , held by columns as a vector, where $( {\bm r} )_i$ denotes the $i$ th component of $\texttt{r}$ , the vector passed to the routine.
params [params_base_type,in] :: is passed unchanged as provided in the call to nlls_solve()/nlls_iterate().

For evaluating the function $P({\bm x},{\bm y}) := ( H_1({\bm x}){\bm y} \dots H_m({\bm x}){\bm y})$ ¶

A subroutine may be supplied to calculate $P({\bm x},{\bm y}) := ( H_1({\bm x}){\bm y} \dots H_m({\bm x}){\bm y})$ for given vectors ${\bm x}, {\bm y} \in \mathbb{R}^n$ . The subroutine must implement the following interface:

abstract interface
   subroutine eval_HP_type(status, n, m, x, y, HP, params)
       integer, intent(inout) :: status
       integer, intent(in) :: n
       integer, intent(in) :: m
       double precision, dimension(n), intent(in)  :: x
       double precision, dimension(n), intent(in)  :: y
       double precision, dimension(n*m), intent(out) :: HP
       class(params_base_type), intent(in) :: params
     end subroutine eval_HP_type
end interface
:language: fortran

subroutine eval_HP(status, n, m, x, y, HP, params)¶

Parameters:

status [integer,inout] :: is initialised to 0 before the routine is called. If it is set to a non-zero value by the routine, then nlls_solve() / nlls_iterate() will exit with error.
n [integer,in] :: is passed unchanged as provided in the call to nlls_solve()/nlls_iterate().
m [integer,in] :: is passed unchanged as provided in the call to nlls_solve()/nlls_iterate().
x (n) [real,in] :: holds the current point ${\bm x}_{k}^{}$ at which to evaluate the Hessians $\nabla^2 r_i( {\bm x_k} )$ .
y (n) [real,in] :: holds ${\bm y}$ , the vector which multiplies each Hessian.
HP (n*m) [real,out] :: must be set by the routine to hold the matrix $P({\bm x},{\bm y}) := ( H_1({\bm x}){\bm y} \dots H_m({\bm x}){\bm y})$ , held by columns as a vector.
params [params_base_type,in] :: is passed unchanged as provided in the call to nlls_solve()/nlls_iterate().

Data types¶

The derived data type for holding options¶

type nlls_options¶

This is used to hold controlling data. The components are automatically given default values in the definition of the type.

Printing Controls

Type fields:

% out [integer,default=6] :: the Fortran unit number for general messages. If it is negative, these messages will be suppressed.

% print_level [integer,default=0] ::

controls the level of output required. Options are:

0	No informational output will occur.
1	Prints a brief summary when finished.
2	Gives a one-line summary for each iteration.
3	As 2, but with more details.
4	As 3, plus gives a summary of the inner iteration for each iteration.
5	As 4, plus gives more verbose (debugging) output.

% print_options [logical,default=false] :: determines whether to print a list of all options and their values at the beggining of the solve.
% print_header [integer,default=30] :: prints the column header every print_header iterations when print_level > 1.

Choice of Algorithm

Type fields:

% model [integer,default=3] ::

specifies the model, $m_k(\cdot)$ , used. Possible values are:

1	Gauss-Newton (no Hessian).
2	(Quasi-)Newton (uses exact Hessian if `exact_second_derivatives` is true, otherwise builds an approximation to the Hessian).
3	Hybrid method (mixture of Gauss-Newton/(Quasi-)Newton, as determined by the package).
4	Newton-tensor method.

See The models for details.

% type_of_method [integer,default=1] ::
specifies the type of globalization method used. Possible values are:

1 Trust-region method.

2 Regularization.

See Subproblem solves for details.

% nlls_method [integer,default=4] ::

specifies the method used to solve (or approximate the solution to) the trust-region sub problem. Possible values are:

1	Powell’s dogleg method (approximates the solution).
2	The Adachi-Iwata-Nakatsukasa-Takeda (AINT) method.
3	The More-Sorensen method.
4	`Galahad`’s `DTRS` method if `type_of_method=1`, or `Galahad`’s `DRQS` method if `type_of_method=2`.

See Subproblem solves for details.

% exact_second_derivatives [logical,default=false] :: if true, signifies that the exact second derivatives are available (and, if false, approximates them using a secant method).

Stopping rules

Type fields:

% maxit [integer,default=100] :: gives the number of iterations the algorithm is allowed to take before being stopped. This is not accessed if nlls_iterate() is used.
% stop_g_absolute [real,default=1e-5] :: specifies the absolute tolerance used in the convergence test on $\|{\iter{\vJ}}^T\vr(\iter{\vx}))\|/\|\vr(\iter{\vx})\|$ .
% stop_g_relative [real,default=1e-8] :: specifies the relative tolerance used in the convergence test on $\|{\iter{\vJ}}^T\vr(\iter{\vx})\|/\|\vr(\iter{\vx})\|$ .
% stop_f_absolute [real,default=1e-5] :: specifies the absolute tolerance used in the convergence test on $\|\vr(\iter{\vx})\|$ .
% stop_f_relative [real,default=1e-8] :: specifies the relative tolerance used in the convergence test on $\|\vr(\iter{\vx})\|$ .
% stop_s [real,default=eps] :: specifies the tolerance used in the convergence test on $\|\iter{\vs}\|$ .

Trust region radius/regularization behaviour

Type fields:

% relative_tr_radius [integer,default=0] :: specifies whether the initial trust region radius should be scaled.
% initial_radius_scale [real,default=1.0] :: specifies the scaling parameter for the initial trust region radius, which is only used if relative_tr_radius = 1.
% initial_radius [real,default=100.0] :: specifies the initial trust-region radius, $\Delta$ .
% maximum_radius [real,default=1e8] :: specifies the maximum size permitted for the trust-region radius.
% eta_successful [real,default=1e-8] :: specifies the smallest value of $\rho$ such that the step is accepted. .. success_but_reduce is also available, but not documented
% eta_very_successful [real,default=0.9] :: specifies the value of $\rho$ after which the trust-region radius is increased.
% eta_too_successful [real,default=2.0] :: specifies that value of $\rho$ after which the step is accepted, but keep the trust-region radius unchanged.
% radius_increase [real,default=2.0] :: specifies the factor to increase the trust-region radius by.
% radius_reduce [real,default=0.5] :: specifies the factor to decrease the trust-region radius by.
% tr_update_strategy [integer,default=1] ::
specifies the strategy used to update $\Delta_k$ . Possible values are:

1 use the usual step function.

2 use a the continuous method.
% reg_order [real,default=0.0] :: if type_of_method = 2, the order of the regularization used ( $q$ in ((7))). If reg_order = 0.0, then the algorithm chooses an appropriate value of $q$ .

Scaling options

Type fields:

% scale [integer,default=1] :: specifies how, if at all, we scale the Jacobian. We calculate a diagonal scaling matrix, ${\tt D}$ , as follows:
% scale_trim_max [logical,default=true] :: specifies whether or not to trim large values of the scaling matrix, $D$ . If true, ${\tt D}_{i,i} \leftarrow min({\tt D}_{i,i}, {\tt scale\_max})$ .
% scale_max [real,default=1e11] :: specifies the maximum value allowed if scale_trim_max = true.
% scale_trim_min [logical,default=true] :: specifies whether or not to trim small values of the scaling matrix, ${\tt D}$ . If true, ${\tt D}_{i,i} \leftarrow max({\tt D}_{i,i}, {\tt scale_max})$ .
% scale_min [real,default=1e-11] :: specifies the minimum value allowed if scale_trim_max = true.
% scale_require_increase [logical,default=false] :: specifies whether or not to require ${\tt D}_{i,i}$ to increase before updating it.

Hybrid method options These options are used if model=3

Type fields:

% hybrid_switch [real,default=0.1] :: specifies the value, if model=3, at which second derivatives are used.
% hybrid_tol [real,default=2.0] :: if model=3, specifies the value such that if $\|{\iter{\vJ}}^T \vW \vr(\vx_k) \|_2 < \mathtt{hybrid\_tol} * 0.5 \|\vr(\vx_k)\|_\vW^2$ the method switches to a quasi-Newton method.
% hybrid_switch_its [integer,default=1] :: if model=3, sets how many iterates in a row must the condition in the definition of hybrid_tol hold before a switch.

Newton-Tensor options These options are used if model=4

Type fields:

% inner_method [integer,default=2] ::

if nlls_method = 4, specifies the method used to pass in the regularization parameter to the inner non-linear least squares solver. Possible values are:

1	The current regularization parameter is passed in as a base regularization parameter.
2	A larger non-linear least squares problem is explicitly formed to be solved as the subproblem.
3	The regularization is handled implicitly by calling RALFit recursively.

More-Sorensen options These options are used if nlls_method=3

Type fields:

% more_sorensen_maxits [integer,default=3] :: if nlls_method = 3, specifies the maximum number of iterations allowed in the More-Sorensen method.
% more_sorensen_maxits :: if nlls_method = 3, specifies the maximum number of iterations allowed in the More-Sorensen method.
% more_sorensen_shift [real,default=1e-13] :: if nlls_method = 3, specifies the shift to be used in the More-Sorensen method.
% more_sorensen_tiny [real,default=10.0*eps] :: if nlls_method = 3, specifies the value below which numbers are considered to be essentially zero.
% more_sorensen_tol [real,default=1e-3] :: if nlls_method = 3, specifies the tolerance to be used in the More-Sorensen method.

Other options

Type fields:	`%` `output_progress_vectors` [logical,default=false] :: if true, outputs the progress vectors `nlls_inform%resvec` and `nlls_inform%gradvec` at the end of the routine.

Internal options to help solving a regularized problem implicitly

Type fields:

% regularization [integer,default=0] ::

specifies the method by which a regularized non-linear least squares problem is solved implicitly. Is designed to be used when solving the nonlinear least-squares problem recursively. Possible values are:

0	$\sigma = 0$ , and an unchanged problem is solved
1	a non-linear least-squares problem of size $n+m$ is solved implicitly. Can only be used if $p=2$ .
2	a non-linear least-squares problem of size $n+1$ is solved implicitly.

See Incorporating the regularization term for details.

% regularization_term [real,default=0.0] :: specifies the regularization weight, $\sigma$ , used when implicitly solving the least-squares problem.
% regularization_power [real,default=0.0] :: specifies the regularization index, $p$ , used when implicitly solving the least-squares problem.

The derived data type for holding information¶

type nlls_inform¶

This is used to hold information about the progress of the algorithm.

Type fields:

% status [integer] :: gives the exit status of the subroutine. See Warning and error messages for details.
% error_message (80) [character] :: holds the error message corresponding to the exit status.
% alloc_status [integer] :: gives the status of the last attempted allocation/deallocation.
% bad_alloc (80) [character] :: holds the name of the array that was being allocated when an error was flagged.
% iter [integer] :: gives the total number of iterations performed.
% f_eval [integer] :: gives the total number of evaluations of the objective function.
% g_eval [integer] :: gives the total number of evaluations of the gradient of the objective function.
% h_eval [integer] :: gives the total number of evaluations of the Hessian of the objective function.
% convergence_normf [integer] :: tells us if the test on the size of $\vr$ is satisfied.
% convergence_normf :: that tells us if the test on the size of the gradient is satisfied.
% convergence_normf :: that tells us if the test on the step length is satisfied.
% resvec (iter+1) [real] :: if output_progress_vectors=true in nlls_options, holds the vector of residuals.
% resvec :: if output_progress_vectors=true in nlls_options, holds the vector of gradients.
% obj [real] :: holds the value of the objective function at the best estimate of the solution determined by the algorithm.
% norm_g [real] :: holds the gradient of the objective function at the best estimate of the solution determined by the package.
% scaled_g [real] :: holds a scaled version of the gradient of the objective function at the best estimate of the solution determined by the package.
% external_return [integer] :: gives the error code that was returned by a call to an external routine.
% external_name (80) [character] :: holds the name of the external code that flagged an error.
% step [real] :: holds the size of the last step taken.

The workspace derived data type¶

type nlls_workspace¶: This is used to save the state of the algorithm in between calls to nlls_iterate(), and must be used if that subroutine is required. It’s components are not designed to be accessed by the user.

Warning and error messages¶

A successful return from a subroutine in the package is indicated by status in nlls_inform having the value zero. A non-zero value is asscociated with an error message, which will be output on error in nlls_inform.

Possible values are:

-1	Maximum number of iterations reached without convergence.
-2	Error from evaluating a function/Jacobian/Hessian.
-3	Unsupported choice of model.
-4	Error return from an external routine.
-5	Unsupported choice of method.
-6	Allocation error.
-7	Maximum number of reductions of the trust radius reached.
-8	No progress being made in the solution.
-9	n $>$ m.
-10	Unsupported trust region update strategy.
-11	Unable to valid step when solving trust region subproblem.
-12	Unsupported scaling method.
-13	Error accessing pre-allocated workspace.
-14	Unsupported value in `type_of_method`.
-15	Unsupported value of `inner_method` passed in options.
-101	Unsupported model in dogleg (`nlls_method=1`).
-201	All eigenvalues are imaginary (`nlls_method=2`).
-202	Matrix with odd number of columns sent to `max_eig` subroutine (`nlls_method=2`).
-301	`more_sorensen_max_its` is exceeded in `more_sorensen` subroutine (`nlls_method=3`).
-302	Too many shifts taken in `more_sorensen` subroutine (`nlls_method=3`).
-303	No progress being made in `more_sorensen` subroutine (`nlls_method=3`).
-401	`model = 4` selected, but `exact_second_derivatives` is set to `false`.
-900	Illegal value of `print_level` in options array. Valid range is 0 to 5.
-950	Combination of method/regularization options not yet implemented
-999	Unexpected error occured

How to use the package¶

Overview¶

Calling sequences¶

The derived data types¶

Argument lists and calling sequences¶

To solve the non-linear least squares problem¶

To iterate once¶

User-supplied function evaluation routines¶

For evaluating the function ¶

For evaluating the function ¶

For evaluating the function ¶

For evaluating the function ¶

Data types¶

The derived data type for holding options¶

The derived data type for holding information¶

The workspace derived data type¶

Warning and error messages¶

For evaluating the function ${\bm r} ( {\bm x} )$ ¶

For evaluating the function ${\bm J} = \nabla {\bm r} ( {\bm x} )$ ¶

For evaluating the function $Hf = \sum_{i=1}^m r_i( {\bm x} ) {\bm W} \nabla^2 r_i( {\bm x} )$ ¶

For evaluating the function $P({\bm x},{\bm y}) := ( H_1({\bm x}){\bm y} \dots H_m({\bm x}){\bm y})$ ¶