User Documentation for RALFit (Fortran version)¶

RALFit computes a solution $\vx$ to the non-linear least-squares problem

$\min_\vx \ F(\vx) := \frac{1}{2}\| \vr(\vx) \|_{\vW}^2 + \frac{\sigma}{p}\| \vx\|_2^p,$

where $\vW\in\mathbb{R}^{m\times m}$ is a diagonal, non-negative, weighting matrix, and $\vr(\vx) =(\comp[1]{r}(\vx), \comp[2]{r}(\vx),...,\comp[m]{r}(\vx))^T$ is a non-linear function.

A typical use may be to fit a function $f(\vx)$ to the data $y_i, \ t_i$ , weighted by the uncertainty of the data, $\sigma_i$ , so that

$r_i(\vx) := y_i - f(\vx;t_i),$

and $\vW$ is the diagonal matrix such that $\vW_{ii} = (1/\sqrt{\sigma_i}).$ For this reason we refer to the function $\vr$ as the residual function.

Various algorithms for solving this problem are implemented – see Description of the method used.

Contents¶

Installation¶

Obtaining the code¶

The latest version of the code can be downloaded from GitHub by issuing the command

git clone https://github.com/ralna/RALFit.git

Building the library¶

From the RALFit/libRALFit/ directory, issue the commands:

mkdir build
cd build
cmake ..
make

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:

% error [integer,default=6] :: the Fortran unit number for error messages. If it is negative, these messages will be suppressed.
% 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:

< 1	No informational output will occur.
1	Gives a one-line summary for each iteration.
2	As 1, plus gives a summary of the inner iteration for each iteration.
3	As 2, plus gives more verbose (debugging) output.

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

Solving a regularized problem

Type fields:

% regularization [integer,default=0] ::

specifies if a regularized non-linear least squares problem needs to be solved, and if so, what method is used to solve it. Possible values are:

0	$\sigma = 0$ , and an unregularized 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 in the least-squares problem.
% regularization_power [real,default=0.0] :: specifies the regularization index, $p$ , used in the least-squares problem.

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.

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:

% reg_order [real,default=0.0] :: if nlls_method = 4, the order of the regularization used ( $p$ in TODO (eq:: reg_subproblem)). If reg_order = 0.0, then the algorithm chooses an appropriate value of $p$ .

% 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:	`%` `calculate_svd_J` [logical,default=false] :: specifies whether or not to calculate the singular value decomposition of ${\tt J}$ at each iteration. `%` `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.

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

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

Description of the method used¶

RALFit computes a solution $\vx$ to the non-linear least-squares problem

(1)¶ $\min_\vx \ F(\vx) := \frac{1}{2}\| \vr(\vx) \|_{\vW}^2 + \frac{\sigma}{p}\| \vx\|_2^p,$

Here we describe the method used to solve (1). RALFit implements an iterative method that, at each iteration, calculates and returns a step $\vs$ that reduces the model by an acceptable amount by solving (or approximating a solution to) a subproblem, as detailed in Subproblem solves.

The algorithm is iterative. At each point, $\iter{\vx}$ , the algorithm builds a model of the function at the next step, $F({\iter{\vx}+\iter{\vs}})$ , which we refer to as $m_k(\cdot)$ . We allow either a Gauss-Newton model, a (quasi-)Newton model, or a Newton-tensor model; see The models for more details.

Once the model has been formed we find a candidate for the next step by solving a subitable subproblem. The quantity

(2)¶ $\rho = \frac{F(\iter{\vx}) - F(\iter{\vx} + \iter{\vs})}{\iter{m}(\iter{\vx}) - \iter{m}(\iter{\vx} + \iter{\vs})}$

is then calculated. If this is sufficiently large we accept the step, and $\iter[k+1]{\vx}$ is set to $\iter{\vx} + \iter{\vs}$ ; if not, the parameter $\Delta_k$ is reduced and the resulting new trust-region sub-problem is solved. If the step is very successful – in that $\rho$ is close to one – $\Delta_k$ is increased. Details are explained in Accepting the step and updating the parameter.

This process continues until either the residual, $\|\vr(\iter{\vx})\|_\vW$ , or a measure of the gradient, $\|{\iter{\vJ}}^T\vW\vr(\iter{\vx})\|_2 / \|\vr(\iter{\vx})\|_\vW$ , is sufficiently small.

The models¶

A vital component of the algorithm is the choice of model employed. There are four choices available, controlled by the parameter model of nlls_options.

model = 1

this implements the Gauss-Newton model. Here we replace $\vr(\iter[k]{\vx} + \vs)$ by its first-order Taylor approximation, $\vr(\iter{\vx}) + \iter{\vJ}\vs$ . The model is therefore given by

(3)¶ $m_k^{GN}(\vs) = \frac{1}{2} \|\vr(\iter{\vx}) + \iter{\vJ}\vs\|_\vW^2.$

model = 2

this implements the Newton model. Here, instead of approximating the residual, $\vr(\cdot)$ , we take as our model the second-order Taylor approximation of the function, $F(\iter[k+1]{\vx}).$ Namely, we use

(4)¶

where $\iter{\vg} = {\iter{\vJ}}^T\vW \vr(\iter{\vx})$ and $\iter{\vH} = \sum_{i=1}^m\iter[i]{r}(\iter{\vx}) \vW \nabla^2 \iter[i]{r}(\iter{\vx}).$ Note that $m_k^{N}(\vs) = m_k^{GN}(\vs) + \frac{1}{2}\vs^T\iter{\vH} \vs$ .

If the second derivatives of $\vr(\cdot)$ are not available (i.e., the option exact_second_derivatives is set to false, then the method approximates the matrix $\iter{\vH}$ ; see Approximating the Hessian.

model = 3

This implements a hybrid model. In practice the Gauss-Newton model tends to work well far away from the solution, whereas Newton performs better once we are near to the minimum (particularly if the residual is large at the solution). This option will try to switch between these two models, picking the model that is most appropriate for the step. In particular, we start using $m_k^{GN}(\cdot)$ , and switch to $m_k^{N}(\cdot)$ if $\|{\iter{\vg}}\|_2 \leq \mathtt{hybrid\_tol} \frac{1}{2}\|\vr(\iter{\vx})\|^2_\vW$ for more than hybrid_switch_its iterations in a row. If, in subsequent iterations, we fail to get a decrease in the function value, then the algorithm interprets this as being not sufficiently close to the solution, and thus switches back to using the Gauss-Newton model.

The exact method used is described below:

$& \mathbf{if } \texttt{ use\_second\_derivatives} \qquad \textit{ // previous step used Newton model} \\ & \qquad \mathbf{if } \|\iter[k+1]{\tg}\| > \|\iter[k]{\tg} \| \\ & \qquad \qquad \texttt{use\_second\_derivatives = false} \qquad \textit{ // Switch back to Gauss-Newton} \\ & \qquad \qquad {\iter[temp]{\thess}} = \iter[k]{\thess}, \ \iter[k]{\thess} = 0 \qquad \textit{ // Copy Hessian back to temp array} \\ & \qquad \mathbf{end if } \\ & \mathbf{else} \\ & \qquad \mathbf{if } \|\iter[k+1]{\tg}\| / \texttt{normF}_{k+1} < \texttt{hybrid\_tol} \\ & \qquad \qquad \texttt{hybrid\_count = hybrid\_count + 1} \qquad \textit{ // Update the no of successive failures} \\ & \qquad \qquad \textbf{if } \texttt{hybrid\_count = hybrid\_count\_switch\_its} \\ & \qquad \qquad \qquad \texttt{use\_second\_derivatives = true} \\ & \qquad \qquad \qquad \texttt{hybrid\_count = 0} \\ & \qquad \qquad \qquad \iter[temp]{\thess} = {\iter[k]{\thess}} \textit{// Copy approximate Hessian back} \\ & \qquad \qquad \textbf{end if} \\ & \qquad \textbf{end if} \\ & \textbf{end if} \\$

model = 4

this implements a Newton-tensor model. This uses a second order Taylor approximation to the residual, namely

$r_i(\iter{\vx} + \vs) \approx (\iter{\vt}(\vs))_i := r_i(\iter{\vx}) + (\iter{\vJ})_i\vs + \frac{1}{2}\vs^T B_{ik}\vs,$

where $(\iter{\vJ})_i$ is the ith row of $\iter{\vJ}$ , and $B_{ik}$ is $\nabla^2 r_i(\iter{\vx})$ . We use this to define our model

(5)¶ $m_k^{NT}(\vs) = \frac{1}{2}\|\vt_k(\vs)\|_\vW^2.$

Approximating the Hessian¶

If the exact Hessian is not available, we approximate it using the method of Dennis, Gay, and Welsch [4]. The method used is given as follows:

$& \textbf{function} \ \iter[k+1]{\thess} = \mathtt{rank\_one\_update}(\td ,\iter[k]{\tg},\iter[k+1]{\tg}, \iter[k+1]{\tr},\iter[k]{\tJ},\iter[k]{\thess}) \\ & \ty = \iter[k]{\tg} - \iter[k+1]{\tg} \\ & \widehat{\ty} = {\iter[k]{\tJ}}^T \iter[k+1]{\tr} - \iter[k+1]{\tg} \\ & \widehat{\iter[k]{\thess}} = \min\left( 1, \frac{|\td^T\widehat{\ty}|}{|\td^T\iter[k]{\thess}\td|}\right) \iter[k]{\thess} \\ & \iter[k+1]{\thess} = \widehat{\iter[k]{\thess}} + \left(({\iter[k+1]{\widehat{\ty}}} - \iter[k]{\thess}\td )^T\td\right)/\ty^T\td$

It is sometimes the case that this approximation becomes corrupted, and the algorithm may not recover from this. To guard against this, if model = 3 in nlls_options and we are using this approximation to the Hessian in our (quasi-Newton) model, we test against the Gauss-Newton model if the first step is unsuccessful. If the Gauss-Newton step would have been successful, we discard the approximate Hessian information, and recompute the step using Gauss-Newton.

In the case where model=3, the approximation to the Hessian is updated at each step whether or not it is needed for the current calcuation.

Subproblem solves¶

The main algorithm calls a number of subroutines. The most vital is the subroutine calculate_step, which finds a step that minimizes the model chosen, subject to a globalization strategy. The algorithm supports the use of two such strategies: using a trust-region, and regularization. If Gauss-Newton, (quasi-)Newton, or a hybrid method is used (model = 1,2,3 in nlls_options), then the model function is quadratic, and the methods available to solve the subproblem are described in The trust region method and Regularization. If the Newton-Tensor model is selected (model = 4 in nlls_options), then this model is not quadratic, and the methods available are described in Newton-Tensor subproblem.

Note that, when calculating the step, if the initial regularization parameter $\sigma$ in (1) is non-zero, then we must modify ${\iter[k]{\tJ}}^T\iter[k]{\tJ}$ to take into account the Jacobian of the modified least squares problem being solved. Practically, this amounts to making the change

${\iter[k]{\tJ}}^T\iter[k]{\tJ} = {\iter[k]{\tJ}}^T\iter[k]{\tJ} + \begin{cases} \sigma I & \text{if }p = 2\\ \frac{\sigma p}{2} \|\iter[k]{\vx}\|^{p-4}\iter[k]{\vx}{\iter[k]{\vx}}^T & \text{otherwise} \end{cases}.$

The trust region method¶

If model = 1, 2, or 3, and type_of_method=1, then we solve the subproblem

(6)¶ $\iter{\vs} = \arg \min_{\vs} \ \iter{m} (\vs) \quad \mathrm{s.t.} \quad \|\vs\|_B \leq \Delta_k,$

and we take as our next step the minimum of the model within some radius of the current point. The method used to solve this is dependent on the control parameter optionsnlls_method. The algorithms called for each of the options are listed below:

nlls_method = 1

approximates the solution to (6) by using Powell’s dogleg method. This takes as the step a linear combination of the Gauss-Newton step and the steepest descent step, and the method used is described here:

$& \textbf{function} \ \texttt{dogleg}(\tt \tJ, {\tr}, \thess, \tg,\Delta)\\ & \alpha = \|\tg\|^2 / \|\tJ * \tg\|^2 \\ & \td_{\rm sd} = \alpha \,\tg \\ & \text{solve } \td_{\rm gn} = \arg \min_{\tx}\|\tJ \tx- \tr\|_2 \\ & \textbf{if } \|\td_{\rm gn}\| \leq \Delta \textbf{then} \\ & \qquad \td = \td_{\rm gn} \\ & \textbf{else if } \|\alpha \, \td_{\rm sd}\| \geq \Delta \\ & \qquad \td = (\Delta / \|\td_{\rm sd}\|) \td_{\rm sd} \\ & \textbf{else} \\ & \qquad \td = \alpha \, \td_{\rm sd} + \beta\, (\td_{\rm gn} - \alpha \td_{\rm sd}), \ \text{where } \beta \text{ is chosen such that } \|\td\| = \Delta \\ & \textbf{end if}$

nlls_method = 2

solves the trust region subproblem using the trust region solver of Adachi, Iwata, Nakatsukasa, and Takeda. This reformulates the problem (6) as a generalized eigenvalue problem, and solves that. See [1] for more details.

nlls_method = 3

this solves (6) using a variant of the More-Sorensen method. In particular, we implement Algorithm 7.3.6 in Trust Region Methods by Conn, Gould and Toint [2].

nlls_method = 4

this solves (6) by first converting the problem into the form

$\min_\vp \vw^T \vp + \frac{1}{2} \vp^T \vD \vp \quad {\rm s.t.} \quad \|\vp\| \leq \Delta,$

where $\vD$ is a diagonal matrix. We do this by performing an eigen-decomposition of the Hessian in the model. Then, we call the Galahad routine DTRS; see the Galahad [3] documentation for further details.

Regularization¶

If model = 1, 2, or 3, and type_of_method=2, then the next step is taken to be the minimum of the model with a regularization term added:

(7)¶ $\iter{\vs} = \arg \min_{\vs} \ \iter{m} (\vs) + \frac{1}{\Delta_k}\cdot \frac{1}{p} \|\vs\|_B^p,$

At present, only one method of solving this subproblem is supported:

nlls_method = 4:

this solves (7) by first converting the problem into the form

$\min_\vp \vw^T \vp + \frac{1}{2} \vp^T \vD \vp + \frac{1}{p}\|\vp\|_2^p,$

where $\vD$ is a diagonal matrix. We do this by performing an eigen-decomposition of the Hessian in the model. Then, we call the Galahad routine DRQS; see the Galahad [3] documentation for further details.

Newton-Tensor subproblem¶

If model=4, then the non-quadratic Newton-Tensor model is used. As such, none of the established subproblem solvers described in The trust region method or Regularization can be used.

If we use regularization (with $p=2$ ), then the subproblem we need to solve is of the form

(8)¶ $\min_\vs \frac{1}{2}\sum_{i=1}^mW_{ii}{(\vt_k(\vs))_i}^2 + \frac{1}{2\Delta_k}\|\vs\|_2^2$

Note that (8) is a sum-of-squares, and as such can be solved by calling nlls_solve() recursively. We support two options:

inner_method = 1

if this option is selected, then nlls_solve() is called to solve (5) directly. The current regularization parameter of the ‘outer’ method is used as a base regularization in the ‘inner’ method, so that the (quadratic) subproblem being solved in the ‘inner’ call is of the form

$\min_\vs \, m_k(\vs) + \frac{1}{2}\left(\frac{1}{\Delta_k} + \frac{1}{\delta_k}\right)\|\vs\|_B^2,$

where $m_k(\vs)$ is a quadratic model of (5), $\Delta_k$ is the (fixed) regularization parameter of the outer iteration, and $\delta_k$ the regularization parameter of the inner iteration, which is free to be updated as required by the method.

inner_method = 2

in this case we use nlls_solve() to solve the regularized model (8) directly. The number of parameters for this subproblem is $n+m$ . Specifically, we have a problem of the form

$\min_\vs \frac{1}{2} \|\widehat{\vr}(\vs)\|_\vW^2, \quad \text{where } (\widehat{\vr}(\vs))_i = \begin{cases} (\vt_k(\vs))_i & 1 \leq i \leq m \\ \frac{1}{\sqrt{\Delta_k}}s_i& m+1 \leq i \leq n+m \end{cases}.$

This subproblem can then be solved using any of the methods described in The trust region method or Regularization.

inner_method = 3

In this case, nlls_solve() is called recursively with the inbuilt feature of solving a regularized problem, as described in Incorporating the regularization term

Accepting the step and updating the parameter¶

Once a step has been suggested, we must decide whether or not to accept the step, and whether the trust region radius or regularization parameter, as appropriate, should grow, shrink, or remain the same.

These decisions are made with reference to the parameter, $\rho$ (2), which measures the ratio of the actual reduction in the model to the predicted reduction in the model. If this is larger than eta_successful in nlls_options, then the step is accepted.

The value of $\Delta_k$ then needs to be updated, if appropriate. The package supports two options:

tr_update_strategy = 1: a step-function is used to decide whether or not to increase or decrease $\Delta_k$ , as described here:

$& \textbf{if } \rho \leq \texttt{eta\_success\_but\_reduce} \textbf{ then} \\ & \qquad \Delta = \texttt{radius\_reduce} * \Delta \qquad \textit{// \ reduce }\mathit{ \Delta} \\ & \textbf{else if } \rho \leq \texttt{eta\_very\_successful} \\ & \qquad \Delta = \Delta \qquad \textit{// } \mathit{\Delta} \textit{ stays unchanged} \\ & \textbf{else if } \rho \leq \texttt{eta\_too\_successful} \\ & \qquad \Delta = \texttt{radius\_increase} * \Delta \qquad \textit{// increase }\mathit{\Delta} \\ & \textbf{else if } \rho > \texttt{eta\_too\_successful}\\ & \qquad \Delta = \Delta \qquad \textit{// too successful: accept step, but don't change }\mathit{\Delta}\\ & \textbf{end if }$
tr_update_strategy = 2: a continuous function is used to make the decision [5], as described below. On the first call, the parameter $\nu$ is set to $2.0$ .

$& \textbf{if } \rho \geq \texttt{eta\_too\_successful} \\ & \qquad \Delta = \Delta \qquad \textit{// }\mathit{\Delta}\textit{ stays unchanged} \\ & \textbf{else if } \rho > \texttt{eta\_successful} \\ & \qquad \Delta = \Delta * \min\left(\texttt{radius\_increase}, 1 - \left( (\texttt{radius\_increase} -1)*((1 - 2*\rho)^3) \right)\right) \\ & \qquad \nu = \texttt{radius\_reduce} \\ & \textbf{else if } \rho \leq \texttt{eta\_successful} \\ & \qquad \Delta = \nu * \Delta \\ & \qquad \nu = 0.5 * \nu \\ & \textbf{end if }$

Incorporating the regularization term¶

If a non-zero regularization term is required in (1), then this is handled by transforming the problem internally into a new non-linear least-squares problem. The formulation used will depend on the value of regularization in nlls_options.

regularization = 1

This is only supported if $\bf p = 2$ . We solve a least squares problem with $n$ additional degrees of freedom. The new function, $\widehat{\vr} : \mathbb{R}^{n}\rightarrow\mathbb{R}^{m+n}$ , is defined as

$\widehat{\vr}_i(\vx) = \begin{cases} \vr_i(\vx) & \text{for } i = 1,\dots, m \\ \sqrt{\sigma}[\vx]_j & \text{for } i = m+j, \ j = 1,\dots,n \end{cases}$

where $[\vx]_j$ denotes the $j$ th component of $\vx$ .

This problem is now in the format of a standard non-linear least-squares problem. In addition to the function values, the we also need a Jacobian and some more information about the Hessian. For our modified function, the Jacobian is

$\widehat{\vJ}(\vx) = \begin{bmatrix} \vJ(\vx) \\ \sqrt{\sigma} I \end{bmatrix},$

and the other function that needs to be supplied is given by

$\widehat{\vH}_k(\vx) = \sum_{i=1}^{n+m} \widehat{\vr}_i(\vx) \nabla^2 \widehat{\vr}_i(\vx) = \sum_{i=1}^{n} {\vr}_i(\vx) \nabla^2 {\vr}_i(\vx) = {\vH}_k(\vx).$

We solve these problems implicitly by modifing the code so that the user does not need do any additional work. We can simply note that

$\|\widehat{\vr}(\vx)\|^2 = \|\vr(\vx) \|^2 + \sigma \|\vx\|^2,$

$\widehat{\vJ}^T\widehat{\vr} = \vJ^T\vr + \sigma \vx,$

and that

$\widehat{\vJ}^T\widehat{\vJ} = \vJ^T\vJ + \sigma I.$

We also need to update the value of the model. Since the Hessian vanishes, we only need to be concerned with the Gauss-Newton model. We have that

$\widehat{m}_k^{GN}(\vs)& = \frac{1}{2} \|\widehat{\vr}(\vx_k) + \widehat{\vJ}_k\vs\|^2\\ &=\frac{1}{2} \left(\widehat{\vr}^T \widehat{\vr} + 2 \vs^T \widehat{\vJ}_k^T\widehat{\vr} + \vs^T\widehat{\vJ}_k^T \widehat{\vJ}_k \vs \right) \\ &= \frac{1}{2} \left(\vr^T\vr + \sigma \vx^T\vx + 2(\vs^T{\vJ_k}^T\vr + \sigma \vs^T\vx) + \vs^T {\vJ_k}^T{\vJ_k}\vs + \sigma \vs^T\vs \right)\\ &= m_k^{GN}(\vs) + \frac{1}{2}\sigma(\vx^T\vx + 2\vs^T\vx + \vs^T\vs) \\ &= m_k^{GN}(\vs) + \frac{\sigma}{2} \| \vx + \vs \|^2$

regularization=2

We solve a non-linear least-squares problem with one additional degree of freedom.

Since the term $\frac{\sigma}{p}\|\vx\|_2^p$ is non-negative, we can write

$F_\sigma(\vx) = \frac{1}{2} \left( \|\vr(\vx)\|^2 + \left(\left(\frac{2 \sigma}{p}\right)^{1/2} \|\vx\|^{p/2}\right)^2 \right),$

thereby defining a new non-linear least squares problem involving the function $\vr:\mathbb{R}^{n} \rightarrow \mathbb{R}^{m+1}$ such that

$\bar{r}_i(\vx) = \begin{cases} r_i(\vx) & 1 \leq i \leq m \\ \frac{2\sigma}{p} \|\vx\|^{p/2}& i = m+1 \end{cases}.$

The Jacobian for this new function is given by

$\bar{\vJ}(\vx) = \begin{bmatrix} \vJ(\vx) \\ \left(\frac{\sigma p}{2}\right)^{1/2} \|\vx\|^{(p-4)/2}\vx^T \end{bmatrix},$

and we get that

$\nabla^2 \bar{r}_{m+1} = \left(\frac{\sigma p}{2}\right)^{1/2} \|\vx\|^{(p-4)/2}\left(I + \frac{\vx\vx^T}{\|\vx\|^2}\right).$

As for the case where regularization=1, we simply need to update quantities in our non-linear least squares code to solve this problem, and the changes needed in this case are

$\|\bar{\vr}(\vx)\|^2 = \|\vr(\vx)\|^2 + \frac{2\sigma}{p} \|\vx\|^p,$

$\bar{\vJ}^T\bar{\vr} = \vJ^T \vr + \sigma \|\vx\|^{p-2}\vx,$

$\bar{\vJ}^T\bar{\vJ} = \vJ^T\vJ + \frac{\sigma p }{2}\|\vx\|^{p-4}\vx\vx^T,$

$\sum_{i=1}^{m+1} \bar{r}_i(\vx)\bar{\vH}_i(\vx) = \sigma \|\vx\|^{p-4} \left(\|\vx\|^2 I + \vx\vx^T\right) + \sum_{i=1}^{m} {r}_i(\vx)\vH_i(\vx)$

We also need to update the model. Here we must consider the Gauss-Newton and Newton models separately.

$\bar{m}_k^{GN}(\vs)& = \frac{1}{2} \|\bar{\vr}(\vx_k) + \bar{\vJ}_k\vs\|^2\\ &=\frac{1}{2} \left(\bar{\vr}^T \bar{\vr} + 2 \vs^T \bar{\vJ}_k^T\bar{\vr} + \vs^T\bar{\vJ}_k^T \bar{\vJ}_k \vs \right) \\ &= \frac{1}{2} \left(\vr^T\vr + \frac{2\sigma}{p} \|\vx\|^p + 2(\vs^T{\vJ_k}^T\vr + \sigma\|\vx\|^{p-2} \vs^T\vx) + \vs^T {\vJ_k}^T{\vJ_k}\vs + \frac{\sigma p}{2}\|\vx\|^{p-4} (\vs^T\vx)^2 \right)\\ &= m_k^{GN}(\vs) + \sigma\left( \frac{1}{p}\|\vx\|^p + \|\vx\|^{p-2}\vs^T\vx + \frac{p}{4} \|\vx\|^{p-4}(\vs^T\vx)^2 \right).$

If we use a Newton model then

$\bar{m}_k^N(\vs) &= \bar{m}_k^{GN}(\vs) + \frac{1}{2} \vs^T \bar{\vH_k}\vs\\ & = \bar{m}_k^{GN}(\vs) + \frac{1}{2} \vs^T \left( \vH_k + \sigma\|\vx\|^{p-2}\left(I + \frac{\vx\vx^T}{\|\vx\|^2}\right)\right)\vs\\ & = \bar{m}_k^{GN}(\vs) + \frac{1}{2} \vs^T \vH_k \vs + \frac{\sigma}{2}\|\vx\|^{p-4}\vs^T\left(\vx^T\vx I + \vx\vx^T \right)\vs \\ & = \bar{m}_k^{GN}(\vs) + \frac{1}{2} \vs^T \vH_k \vs + \frac{\sigma}{2}\|\vx\|^{p-4}\left((\vx^T\vx)(\vs^T\vs) + (\vx^T\vs)^2\right)$

[1]	Adachi, Satoru and Iwata, Satoru and Nakatsukasa, Yuji and Takeda, Akiko (2015). Solving the trust region subproblem by a generalized eigenvalue problem. Technical report, Mathematical Engineering, The University of Tokyo.

[2]	Conn, A. R., Gould, N. I., & Toint, P. L. (2000). Trust region methods. SIAM.

[3]	(1, 2) Gould, N. I., Orban, D., & Toint, P. L. (2003). GALAHAD, a library of thread-safe Fortran 90 packages for large-scale nonlinear optimization. ACM Transactions on Mathematical Software (TOMS), 29(4), 353-372.

[4]	Nocedal, J., & Wright, S. (2006). Numerical optimization. Springer Science & Business Media.

[5]	Nielsen, Hans Bruun (1999). Damping parameter in Marquadt’s MEthod. Technical report TR IMM-REP-1999-05, Department of Mathematical Modelling, Technical University of Denmark (http://www2.imm.dtu.dk/documents/ftp/tr99/tr05_99.pdf)

Examples¶

Consider fitting the function $y(t) = x_1e^{x_2 t}$ to data $(\bm{t}, \bm{y})$ using a non-linear least squares fit.

The residual function is given by

$r_i(\vx) = x_1 e^{x_2 t_i} - y_i.$

We can calculate the Jacobian and Hessian of the residual as

$\nabla r_i(\vx) = \left(\begin{array}{cc} e^{x_2 t_i} & t_i x_1 e^{x_2 t_i} \end{array}\right),$

$\nabla^2 r_i(\vx) = \left(\begin{array}{cc} 0 & t_i e^{x_2 t_i} \\ t_i e^{x_2 t_i} & x_1 t_i^2 e^{x_2 t_i} \end{array}\right).$

For some data points, $y_i$ , $t_i$ , $(i = 1,\dots,m)$ the user must return

$\vr(\vx) = \begin{bmatrix} r_1(\vx)\\ \vdots \\ r_m(\vx) \end{bmatrix}, \quad \vJ(\vx) = \begin{bmatrix} \nabla r_1(\vx) \\ \vdots \\ \nabla r_m(\vx) \\ \end{bmatrix}, \quad Hf(\vx) = \sum_{i=1}^m (\vr)_i \nabla^2 r_i(\vx),$

where, in the case of the Hessian, $(\vr)_i$ is the $i\text{th}$ component of a residual vector passed to the user.

$i$	1	2	3	4	5
$t_i$	1	2	4	5	8
$y_i$	3	4	6	11	20

and initial guess $\vx = (2.5, 0.25)$ , the following code performs the fit (with no weightings, i.e., $\vW = \vI$ ).

! examples/Fortran/nlls_example.f90
! 
! Attempts to fit the model y_i = x_1 e^(x_2 t_i)
! For parameters x_1 and x_2, and input data (t_i, y_i)
module fndef_example
   use ral_nlls_double, only : params_base_type
   implicit none

   integer, parameter :: wp = kind(0d0)

   type, extends(params_base_type) :: params_type
      real(wp), dimension(:), allocatable :: t ! The m data points t_i
      real(wp), dimension(:), allocatable :: y ! The m data points y_i
   end type

contains
   ! Calculate r_i(x; t_i, y_i) = x_1 e^(x_2 * t_i) - y_i
   subroutine eval_r(status, n, m, x, r, params)
      integer, intent(out) :: status
      integer, intent(in) :: n
      integer, intent(in) :: m
      real(wp), dimension(*), intent(in) :: x
      real(wp), dimension(*), intent(out) :: r
      class(params_base_type), intent(inout) :: params

      real(wp) :: x1, x2

      x1 = x(1)
      x2 = x(2)
      select type(params)
      type is(params_type)
         r(1:m) = x1 * exp(x2*params%t(:)) - params%y(:)
      end select

      status = 0 ! Success
   end subroutine eval_r
   ! Calculate:
   ! J_i1 = e^(x_2 * t_i)
   ! J_i2 = t_i x_1 e^(x_2 * t_i)
   subroutine eval_J(status, n, m, x, J, params)
      integer, intent(out) :: status
      integer, intent(in) :: n
      integer, intent(in) :: m
      real(wp), dimension(*), intent(in) :: x
      real(wp), dimension(*), intent(out) :: J
      class(params_base_type), intent(inout) :: params

      real(wp) :: x1, x2

      x1 = x(1)
      x2 = x(2)
      select type(params)
      type is(params_type)
         J(  1:  m) = exp(x2*params%t(1:m))                     ! J_i1
         J(m+1:2*m) = params%t(1:m) * x1 * exp(x2*params%t(1:m))! J_i2
      end select

      status = 0 ! Success
   end subroutine eval_J
   ! Calculate
   ! HF = sum_i r_i H_i
   ! Where H_i = [ 0                t_i e^(x_2 t_i)    ]
   !             [ t_i e^(x_2 t_i)  t_i^2 x_1 e^(x_2 t_i)  ]
   subroutine eval_HF(status, n, m, x, r, HF, params)
      integer, intent(out) :: status
      integer, intent(in) :: n
      integer, intent(in) :: m
      real(wp), dimension(*), intent(in) :: x
      real(wp), dimension(*), intent(in) :: r
      real(wp), dimension(*), intent(out) :: HF
      class(params_base_type), intent(inout) :: params

      real(wp) :: x1, x2

      x1 = x(1)
      x2 = x(2)
      select type(params)
      type is(params_type)
         HF(    1) = sum(r(1:m) * 0)                                      ! H_11
         HF(    2) = sum(r(1:m) * params%t(1:m) * exp(x2*params%t(1:m)))  ! H_21
         HF(1*n+1) = HF(2)                                                ! H_12
         HF(1*n+2) = sum(r(1:m) * (params%t(1:m)**2) * x1 * exp(x2*params%t(1:m)))! H_22
      end select

      status = 0 ! Success
   end subroutine eval_HF
end module fndef_example

program nlls_example
   use ral_nlls_double
   use fndef_example
   implicit none

   type(nlls_options) :: options
   type(nlls_inform) :: inform

   integer :: m,n
   real(wp), allocatable :: x(:)
   type(params_type) :: params

   ! Data to be fitted
   m = 5
   allocate(params%t(m), params%y(m))
   params%t(:) = (/ 1.0, 2.0, 4.0,  5.0,  8.0 /)
   params%y(:) = (/ 3.0, 4.0, 6.0, 11.0, 20.0 /)
   
   ! Call fitting routine
   n = 2
   allocate(x(n))
   x = (/ 2.5, 0.25 /) ! Initial guess
   call nlls_solve(n, m, x, eval_r, eval_J, eval_HF, params, options, inform)
   if(inform%status.ne.0) then
      print *, "ral_nlls() returned with error flag ", inform%status
      stop
   endif

   ! Print result
   print *, "Found a local optimum at x = ", x
   print *, "Took ", inform%iter, " iterations"
   print *, "     ", inform%f_eval, " function evaluations"
   print *, "     ", inform%g_eval, " gradient evaluations"
   print *, "     ", inform%h_eval, " hessian evaluations"
end program nlls_example

This returns the following output:

User Documentation for RALFit (Fortran version)¶

Contents¶

Installation¶

Obtaining the code¶

Building the library¶

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¶

Description of the method used¶

The models¶

Approximating the Hessian¶

Subproblem solves¶

The trust region method¶

Regularization¶

Newton-Tensor subproblem¶

Accepting the step and updating the parameter¶

Incorporating the regularization term¶

Examples¶

Indices and tables¶

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})$ ¶