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.
nlls_setup_bounds()
sets up a bound-constrained non-linear least squares problem.
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 of variables to be fitted; i.e., is the length of the unknown vector . Restriction: n 0.
- m [integer,in] :: holds the number of data points available; i.e., is the number of residuals . Restriction: m 0.
- X (n) [real,inout] :: on entry, it must hold the initial guess for , and on successful exit it holds the solution to the non-linear least squares problem.
- eval_r [procedure] :: given a point , returns the vector . Further details of the format required are given in
eval_r()
in User-supplied function evaluation routines. - eval_J [procedure] :: given a point , returns the Jacobian matrix, , of evaluated at . Further details of the format required are given in
eval_J()
in User-supplied function evaluation routines. - eval_Hf [procedure] :: given vectors and , returns the quantity . Further details of the format required are given in
eval_Hf()
in User-supplied function evaluation routines. Ifexact_second_derivative = .false.
innlls_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()
, andeval_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, . If absent, then the norm in the least squares problem is taken to be the 2-norm, that is, .
- eval_HP [procedure] :: If present, this is a routine that, given vectors , returns the matrix . Further details of the format required are given in
eval_HP()
in User-supplied function evaluation routines. This is only referenced ifmodel = 4
innlls_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, inform, 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 . On return it holds the value of 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.
- X (n) [real,inout] :: on the first call it must hold the initial guess for . On return it holds the value of at the current iterate, and must be passed unaltered to any subsequent call to
The user may use the components convergence_normf
and
convergence_normg
and converge_norms
in nlls_inform
to determine whether the iteration has
converged.
To set up box constraints¶
-
subroutine
nlls_setup_bounds
(params, n, blx, bux, options, info)¶ A call of this form passes in information about box contraints that the solution should satisfy.
n, params, options, info are as in the description of
nlls_solve()
Parameters: - blx (n) [real,in] :: is used to pass in lower bounds on the solution.
- bux (n) [real,in] :: is used to pass in upper bounds on the solution.
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 ¶
A subroutine must be supplied to calculate for a given vector . 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, thennlls_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 at which to evaluate .
- r (m) [real,out] :: must be set by the routine to hold the residual function evaluated at the current point , .
- params [params_base_type,in] :: is passed unchanged as provided in the call to
nlls_solve()
/nlls_iterate()
.
- status [integer,inout] :: is initialised to
For evaluating the function ¶
A subroutine must be supplied to calculate for a given vector . 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, thennlls_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 at which to evaluate .
- J (m*n) [real,out] :: must be set by the routine to hold the Jacobian of the residual function evaluated at the current point , .
J(i*m+j)
must be set to hold . - params [params_base_type,in] :: is passed unchanged as provided in the call to
nlls_solve()
/nlls_iterate()
.
- status [integer,inout] :: is initialised to
For evaluating the function ¶
A subroutine must be supplied to calculate for given vectors and ; here denotes the th component of the vector . 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, thennlls_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 at which to evaluate .
- r (m) [real,in] :: holds , 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 , held by columns as a vector, where denotes the th component of , the vector passed to the routine.
- params [params_base_type,in] :: is passed unchanged as provided in the call to
nlls_solve()
/nlls_iterate()
.
- status [integer,inout] :: is initialised to
For evaluating the function ¶
A subroutine may be supplied to calculate for given vectors . 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, thennlls_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 at which to evaluate the Hessians .
- y (n) [real,in] :: holds , the vector which multiplies each Hessian.
- HP (n*m) [real,out] :: must be set by the routine to hold the matrix , held by columns as a vector.
- params [params_base_type,in] :: is passed unchanged as provided in the call to
nlls_solve()
/nlls_iterate()
.
- status [integer,inout] :: is initialised to
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 everyprint_header
iterations whenprint_level > 1
.
Choice of Algorithm
Type fields: %
model
[integer,default=3] ::specifies the model, , 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
’sDTRS
method iftype_of_method=1
, orGalahad
’sDRQS
method iftype_of_method=2
.See Subproblem solves for details.
%
exact_second_derivatives
[logical,default=false] :: iftrue
, signifies that the exact second derivatives are available (and, iffalse
, 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 ifnlls_iterate()
is used.%
stop_g_absolute
[real,default=1e-5] :: specifies the absolute tolerance used in the convergence test on .%
stop_g_relative
[real,default=1e-8] :: specifies the relative tolerance used in the convergence test on .%
stop_f_absolute
[real,default=1e-5] :: specifies the absolute tolerance used in the convergence test on .%
stop_f_relative
[real,default=1e-8] :: specifies the relative tolerance used in the convergence test on .%
stop_s
[real,default=eps] :: specifies the tolerance used in the convergence test on .
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 ifrelative_tr_radius = 1
.%
initial_radius
[real,default=100.0] :: specifies the initial trust-region radius, .%
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 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 after which the trust-region radius is increased.%
eta_too_successful
[real,default=2.0] :: specifies that value of 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 . Possible values are:
1 use the usual step function. 2 use a the continuous method. %
reg_order
[real,default=0.0] :: iftype_of_method = 2
, the order of the regularization used ( in ((7))). Ifreg_order = 0.0
, then the algorithm chooses an appropriate value of .
Scaling options
Type fields: %
scale
[integer,default=1] :: specifies how, if at all, we scale the Jacobian. We calculate a diagonal scaling matrix, , as follows:%
scale_trim_max
[logical,default=true] :: specifies whether or not to trim large values of the scaling matrix, . Iftrue
, .%
scale_max
[real,default=1e11] :: specifies the maximum value allowed ifscale_trim_max = true
.%
scale_trim_min
[logical,default=true] :: specifies whether or not to trim small values of the scaling matrix, . Iftrue
, .%
scale_min
[real,default=1e-11] :: specifies the minimum value allowed ifscale_trim_max = true
.%
scale_require_increase
[logical,default=false] :: specifies whether or not to require 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, ifmodel=3
, at which second derivatives are used.%
hybrid_tol
[real,default=2.0] :: ifmodel=3
, specifies the value such that if the method switches to a quasi-Newton method.%
hybrid_switch_its
[integer,default=1] :: ifmodel=3
, sets how many iterates in a row must the condition in the definition ofhybrid_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] :: ifnlls_method = 3
, specifies the maximum number of iterations allowed in the More-Sorensen method.%
more_sorensen_maxits
:: ifnlls_method = 3
, specifies the maximum number of iterations allowed in the More-Sorensen method.%
more_sorensen_shift
[real,default=1e-13] :: ifnlls_method = 3
, specifies the shift to be used in the More-Sorensen method.%
more_sorensen_tiny
[real,default=10.0*eps] :: ifnlls_method = 3
, specifies the value below which numbers are considered to be essentially zero.%
more_sorensen_tol
[real,default=1e-3] :: ifnlls_method = 3
, specifies the tolerance to be used in the More-Sorensen method.
Box Bound Options These options are used if
nlls_setup_bounds
is called to set up box constraints.Type fields: %
box_nFref_max
[integer,default=4] :: Memory size for the non-monotone linesearch.%
box_gamma
[real,default=0.9995] :: Kanzow sufficient decrease ratio (see equtaion 25 in Kanzow [2004])%
box_decmin
[real,default=2.0e-16] :: FIXME! macheps%
box_bigbnd
[real,default=1.0e20] :: Magic number to consider box bound (+/-) infinity%
box_wolfe_descent
[real,default=1.0e-4] :: Wolfe descent condition%
box_wolfe_curvature
[real,default=0.9] :: Wolfe curvature condition%
box_kanzow_power
[real,default=2.1] :: Tolerance to consider projected direction a descent direction. See LS STEP, Section 4, p392, Kanzow [2014]%
box_kanzow_descent
[real,default=1.0e-8] :: FIXME! sqrt(mcheps)%
box_quad_model_descent
[real,default=1.0e-8] :: FIXME! sqrt(mcheps)%
box_tr_test_step
[logical,default=true] :: Take a projected trust region step when trust region test is OK? If false, force a linesearch or projected gradient step.%
box_wolfe_test_step
[logical,default=true] :: Take a project trust region step when Wolfe test is OK? If false, force a linesearch or projected gradient step.%
box_tau_min
[real,default=0.25] :: Threshold to determine if the projection of the trust region direction is too severe (:math:0 < tau_min 1)%
box_tau_descent
[real,default=1.0e-4] ::tau_descent
in order to test for descent%
box_max_ntrfail
[integer,default=2] :: Max times trust region iterations can fail without passing the various descent tests. Ignored when .%
box_quad_match
[integer,default=1] :: Number of consecutive times quadratic model matches required before setting initial alpha step for projected gradient step equal toscale_alpha*alpha_k-1
(FIXME)%
box_alpha_scale
[real,default=1.0] :: Initial step scale (if )%
box_Delta_scale
[real,default=2.0] :: Scaling factor to use when updating Delta from linesearch/project gradient step%
box_tau_wolfe
[real,default=0.3] :: FIXME%
box_tau_tr_step
[real,default=0.3] :: FIXME%
box_ls_step_maxit
[integer,default=20] :: FIXME%
box_lineseach_type
[integer,default=1] ::Linesearch type – available options are:
1 Dennis-Schnabel linesearch. 2 Hager-Zhang linesearch. For more information, see Bound constraints
Other options
Type fields: %
output_progress_vectors
[logical,default=false] :: if true, outputs the progress vectorsnlls_inform%resvec
andnlls_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 , and an unchanged problem is solved 1 a non-linear least-squares problem of size is solved implicitly. Can only be used if . 2 a non-linear least-squares problem of size is solved implicitly. See Incorporating the regularization term for details.
%
regularization_term
[real,default=0.0] :: specifies the regularization weight, , used when implicitly solving the least-squares problem.%
regularization_power
[real,default=0.0] :: specifies the regularization index, , 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 usingeval_hf
.%
hp_eval
[intgeer] :: gives the total number of evaluations of the Hessian of the objective function usingeval_hp
.%
convergence_normf
[integer] :: tells us if the test on the size of 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] :: ifoutput_progress_vectors=true
innlls_options
, holds the vector of residuals.%
resvec
:: ifoutput_progress_vectors=true
innlls_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. |
-16 | Initial iterate x0 is not usable as starting point |
-17 | Unsupported value of linesearch type (box_linesearch_type ) |
-18 | Bad bound constraints (blx <= bux) |
-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 . |
-501 | Linesearch in projected gradient direction failed |
-900 | Illegal value of print_level in options array. Valid range is 0 to 5. |
-999 | Unexpected error occured |