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(n)
      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(n)
      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: