NLS examples

## multistart attempt gsl_nls
gsl_nls(
  fn = y ~ b1 / (1 + b2 * exp(-b3 * x)),
  data = hobbs_weed,
  start = list(b1 = c(0, 1000), b2 = c(0, 1000), b3 = c(0, 10))
)
#> Nonlinear regression model
#>   model: y ~ b1/(1 + b2 * exp(-b3 * x))
#>    data: hobbs_weed
#>       b1       b2       b3 
#> 196.1863  49.0916   0.3136 
#>  residual sum-of-squares: 2.587
#> 
#> Algorithm: multifit/levenberg-marquardt, (scaling: more, solver: qr)
#> 
#> Number of iterations to convergence: 3 
#> Achieved convergence tolerance: 2.842e-14

## multistart attempt 2 gsl_nls
gsl_nls(
  fn = y ~ b1 / (1 + b2 * exp(-b3 * x)),
  data = hobbs_weed,
  start = list(b1 = NA, b2 = NA, b3 = NA)
)
#> Nonlinear regression model
#>   model: y ~ b1/(1 + b2 * exp(-b3 * x))
#>    data: hobbs_weed
#>       b1       b2       b3 
#> 196.1863  49.0916   0.3136 
#>  residual sum-of-squares: 2.587
#> 
#> Algorithm: multifit/levenberg-marquardt, (scaling: more, solver: qr)
#> 
#> Number of iterations to convergence: 4 
#> Achieved convergence tolerance: 5.329e-15

gauss1 <- nls_test_problem("Gauss1")

## data
str(gauss1$data)
#> 'data.frame':	250 obs. of  2 variables:
#>  $ y: num  97.6 97.8 96.6 92.6 91.2 ...
#>  $ x: num  1 2 3 4 5 6 7 8 9 10 ...

## model + target parameters
gauss1[c("fn", "target")]
#> $fn
#> y ~ b1 * exp(-b2 * x) + b3 * exp(-(x - b4)^2/b5^2) + b6 * exp(-(x - 
#>     b7)^2/b8^2)
#> <environment: 0x562f10e8ceb8>
#> 
#> $target
#>           b1           b2           b3           b4           b5           b6 
#>  98.77821087   0.01049728 100.48990633  67.48111128  23.12977336  71.99450300 
#>           b7           b8 
#> 178.99805021  18.38938902

## initial attempt nls
nls(
  formula = gauss1$fn,
  data = gauss1$data,
  start = c(b1 = 100, b2 = 0, b3 = 100, b4 = 75, b5 = 25, b6 = 75, b7 = 125, b8 = 25),
  lower = 0,
  algorithm = "port"
)
#> Error in nls(formula = gauss1$fn, data = gauss1$data, start = c(b1 = 100, : Convergence failure: singular convergence (7)

## multistart attempt gsl_nls
gsl_nls(
  fn = gauss1$fn,
  data = gauss1$data,
  start = list(b1 = 100, b2 = c(0, 1), b3 = NA, b4 = NA, b5 = NA, b6 = NA, b7 = NA, b8 = NA),
  lower = 0
)
#> Nonlinear regression model
#>   model: y ~ b1 * exp(-b2 * x) + b3 * exp(-(x - b4)^2/b5^2) + b6 * exp(-(x -     b7)^2/b8^2)
#>    data: gauss1$data
#>       b1       b2       b3       b4       b5       b6       b7       b8 
#>  98.7782   0.0105 100.4899  67.4811  23.1298  71.9945 178.9981  18.3894 
#>  residual sum-of-squares: 1316
#> 
#> Algorithm: multifit/levenberg-marquardt, (scaling: more, solver: qr)
#> 
#> Number of iterations to convergence: 4 
#> Achieved convergence tolerance: 2.274e-13

lubricant <- nls_test_problem("Lubricant")

## data
str(lubricant$data)
#> 'data.frame':	53 obs. of  3 variables:
#>  $ y : num  5.11 6.39 7.39 5.79 5.11 ...
#>  $ x1: num  0 0 0 0 0 0 0 0 0 0 ...
#>  $ x2: num  0.001 0.741 1.407 0.363 0.001 ...

## model + target parameters
lubricant[c("fn", "target")]
#> $fn
#> y ~ b1/(b2 + x1) + b3 * x2 + b4 * x2^2 + b5 * x2^3 + (b6 * x2 + 
#>     b7 * x2^3) * exp(-x1/(b8 + b9 * x2^2))
#> <environment: 0x562f0ab01c30>
#> 
#> $target
#>            b1            b2            b3            b4            b5 
#> 1054.54053186  206.54577890    1.46031479   -0.25965483    0.02257371 
#>            b6            b7            b8            b9 
#>    0.40138428    0.03528405   57.40463143   -0.47672110

## initial attempt nls
nls(
  formula = lubricant$fn,
  data = lubricant$data,
  start = c(b1 = 1000, b2 = 200, b3 = 1, b4 = -0.01, b5 = 0.01, b6 = 0.01, b7 = 0.01, b8 = 50, b9 = -0.01),
  algorithm = "port"
)
#> Error in nls(formula = lubricant$fn, data = lubricant$data, start = c(b1 = 1000, : Convergence failure: false convergence (8)

## multistart attempt gsl_nls
gsl_nls(
  fn = lubricant$fn,
  data = lubricant$data,
  start = c(b1 = NA, b2 = NA, b3 = NA, b4 = NA, b5 = NA, b6 = NA, b7 = NA, b8 = NA, b9 = NA)
)
#> Nonlinear regression model
#>   model: y ~ b1/(b2 + x1) + b3 * x2 + b4 * x2^2 + b5 * x2^3 + (b6 * x2 +     b7 * x2^3) * exp(-x1/(b8 + b9 * x2^2))
#>    data: lubricant$data
#>         b1         b2         b3         b4         b5         b6         b7 
#> 1054.54080  206.54584    1.46031   -0.25965    0.02257    0.40138    0.03528 
#>         b8         b9 
#>   57.40467   -0.47672 
#>  residual sum-of-squares: 0.08702
#> 
#> Algorithm: multifit/levenberg-marquardt, (scaling: more, solver: qr)
#> 
#> Number of iterations to convergence: 6 
#> Achieved convergence tolerance: 5.69e-16

Multistart algorithm details

Before evaluating more NLS test problems, this section provides a more comprehensive overview of the implemented multistart algorithm. The implementation is primarily based on the global optimization algorithm in (Hickernell and Yuan 1997), but slightly modified to the context of trust-region nonlinear least squares. The multistart algorithm consists of multiple major iterations. At the start of major iteration \(M\), a set of pseudo-random starting points is sampled inside the grid of starting ranges. Starting points with an almost singular (approximate) Hessian matrix are discarded immediately, as these are unlikely to converge to a local optimum. For the remaining points, a few iterations of the NLS optimizer are applied in order to distinguish promising from unpromising starting points. At each major iteration, only the \(q\) most promising starting points are kept, and if a starting point is not discarded from this set for at least \(s\) major iterations, a full NLS optimization routine is executed using this starting point. If a new optimal solution is found, the number of optimal stationary points NOSP is incremented and the number of found worse stationary points NWSP is reset to zero. If the obtained solution has already been observed before, or if it converges to a local minimum that is worse than the current optimal solution, the number of worse stationary points NSWP is incremented. If the number of found worse stationary points (NWSP) becomes much larger than the number of optimal starting points (NOSP), then it is likely that we have exhausted the search space and are unable to further improve the current optimal solution. The following pseudo-code provides a high-level description of the multistart procedure. For simplicity, we first consider the scenario in which all \(p\) parameter starting ranges \([l_1, u_1] \times \ldots \times [l_{p}, u_{p}]\) are pre-defined and fixed.

1. INITIALIZE
   • Choose parameters \(N \geq q \geq 1\), \(L \geq \ell \geq 1\), \(s \geq 1\), \(r > 1\), maxstart \(\geq 1\), minsp \(\geq 1\). Set \(M = 0\), NOSP = NWSP = 0, and initialize a scaling vector \(\boldsymbol{D} = (0.75, \dots, 0.75) \in [0, 1]^p\).
2. SAMPLE
   • Sample \(N\) initial \(p\)-dimensional points \(\boldsymbol{\theta}_1,\ldots, \boldsymbol{\theta}_N\) inside the grid of starting ranges \([l_1, u_1] \times \ldots \times [l_{p}, u_{p}]\) from a pseudo-random Sobol sequence.
   • For \(i = 1,\ldots,N\) and \(k = 1,\ldots,p\), rescale the sampled points inside the grid of starting ranges favoring points near zero according to³: \(\theta^*_{i,k} = \frac{\text{sgn}(\theta_{i,k})}{D_k} \cdot ((|\theta_{i,k}| - l_k^+ - u_k^- + 1)^{D_k} - 1) + l_k^+ - u_k^-\).
3. REDUCE 1 AND CONCENTRATE
   • Discard all points for which \(\det\left(J(\boldsymbol{\theta}^*_i)^TJ(\boldsymbol{\theta}^*_i)\right) < \epsilon\), with \(\epsilon > 0\) small and \(J(\boldsymbol{\theta}^*_i)\) the Jacobian evaluated at \(\boldsymbol{\theta}^*_i\).
   • For each remaining point apply (a small number) \(\ell\) iterations of the NLS optimizer, obtaining concentrated points \(\boldsymbol{\theta}_1^{**},\ldots,\boldsymbol{\theta}_{N_1}^{**}\)
4. REDUCE 2
   • Identify the first \(q\) order statistics of \(F(\boldsymbol{\theta}_1^{**}),\ldots,F(\boldsymbol{\theta}_{N_1}^{**})\), with \(F(\boldsymbol{\theta}_i^{**})\) the NLS objective evaluated at \(\boldsymbol{\theta}_i^{**}\). Denote \(I_q\) for the set of indices (of size \(q\)).
   • Update counter \(S(i) = S(i) + 1\) if \(i \in I_q\), otherwise set \(S(i) = 0\).
5. OPTIMIZE
   • For each \(i\) with \(S(i) \geq s\), if \(F(\boldsymbol{\theta}_i^{**}) < (1 + \delta) \cdot F(\boldsymbol{\theta}^{opt})\), apply \(L\) iterations of the NLS optimizer, obtaining \(\boldsymbol{\theta}_i^{***}\), and set \(S(i) = 0\).
   • If \(F(\boldsymbol{\theta}_i^{***}) < F(\boldsymbol{\theta}^{opt})\), set \(\boldsymbol{\theta}^{opt} = \boldsymbol{\theta}_i^{***}\), and update the scaling vector \(\boldsymbol{D}\) by \(D_k = \left(\frac{\min_\kappa \Delta_\kappa}{\Delta_k}\right)^\alpha\), with exponent \(\alpha \in [0, 1]\) and \(\boldsymbol{\Delta} = \text{diag}(J(\boldsymbol{\theta}^{***}_i)^TJ(\boldsymbol{\theta}^{***}_i))\), i.e. the diagonal damping matrix as used by Marquardt’s scaling method. Update the number of found optimal stationary points NOSP = NOSP + 1, and reset the number of found worse stationary points NWSP = 0.
   • If \(F(\boldsymbol{\theta}_i^{***}) \geq F(\boldsymbol{\theta}^{opt})\), set NWSP = NWSP + 1.
     
6. REPEAT
   • If NWSP \(\geq\) \(r \cdot\) NSP and NSP \(\geq\) minsp (minimum required number of stationary points) then stop with success. Otherwise, if \(M <\) maxstart, increment \(M = M + 1\) and return to step 2. sampling new pseudo-random points for each \(i\) with \(S(i) = 0\), reusing \(\boldsymbol{\theta}_i = \boldsymbol{\theta}_i^{**}\) if \(S(i) > 0\).

It is important to point out that there is no guarantee that the NLS objective \(F(\boldsymbol{\theta}^{opt})\) at the solution returned by the multistart procedure indeed evaluates to the global minimum inside the grid of starting ranges. This may for instance be due to the rescaling function in step 2., which is a type of inverse logistic function scaled to the starting range \([l_k, u_k]\). The scaling exponent \(D_k\) is calculated from the damping matrix in Marquardt’s scaling method, thereby rescaling each parameter differently based on an approximate measure of its order of magnitude. If the damping matrix that is used to calculate the \(D_k\)’s is highly –but incorrectly– confident about the order of magnitudes of certain parameters, we may not explore the parameter starting ranges \([l_k, u_k]\) sufficiently broadly in the subsequent major iterations. Increasing the number of sampled points \(N\) at the start of each major iteration may help to overcome such issues. If we suspect that the returned optimal solution \(\boldsymbol{\theta}^{opt}\) is only a local optimizing solution, we can always force the multistart procedure to continue searching until a better optimum is found by increasing minsp (minimum required number of stationary points).

Example regression problem (Misra1a)

## nls model + data
(misra1a <- nls_test_problem("Misra1a"))
#> $data
#>        y     x
#> 1  10.07  77.6
#> 2  14.73 114.9
#> 3  17.94 141.1
#> 4  23.93 190.8
#> 5  29.61 239.9
#> 6  35.18 289.0
#> 7  40.02 332.8
#> 8  44.82 378.4
#> 9  50.76 434.8
#> 10 55.05 477.3
#> 11 61.01 536.8
#> 12 66.40 593.1
#> 13 75.47 689.1
#> 14 81.78 760.0
#> 
#> $fn
#> y ~ b1 * (1 - exp(-b2 * x))
#> <environment: 0x562f11445798>
#> 
#> $start
#>    b1    b2 
#> 5e+02 1e-04 
#> 
#> $target
#>           b1           b2 
#> 2.389421e+02 5.501564e-04 
#> 
#> attr(,"class")
#> [1] "nls_test_formula"

## nls model fit
with(misra1a, 
     gsl_nls(
       fn = fn,
       data = data,
       start = start
     )
)
#> Nonlinear regression model
#>   model: y ~ b1 * (1 - exp(-b2 * x))
#>    data: data
#>        b1        b2 
#> 2.389e+02 5.502e-04 
#>  residual sum-of-squares: 0.1246
#> 
#> Algorithm: multifit/levenberg-marquardt, (scaling: more, solver: qr)
#> 
#> Number of iterations to convergence: 16 
#> Achieved convergence tolerance: 6.745e-15

Example optimization problem (Rosenbrock)

## nls model fit
with(nls_test_problem("Rosenbrock"), 
     gsl_nls(
       fn = fn,
       y = y,
       start = start,
       jac = jac
     )
)
#> Nonlinear regression model
#>   model: y ~ fn(x)
#> x1 x2 
#>  1  1 
#>  residual sum-of-squares: 9.762e-19
#> 
#> Algorithm: multifit/levenberg-marquardt, (scaling: more, solver: qr)
#> 
#> Number of iterations to convergence: 16 
#> Achieved convergence tolerance: 7.696e-14

Benchmark NLS fits

To conclude, we benchmark the performance of the multistart algorithm by computing NLS model fits for each of the 59 test problems using the multistart algorithm with no starting values provided, i.e. all starting values are set to NA. As trust region method we choose respectively: the default Levenberg-Marquardt algorithm (algorithm = "lm"); the double dogleg algorithm (algorithm = "ddogleg"); and the Levenberg-Marquardt algorithm with geodesic acceleration (algorithm = "lmaccel"). The maximum number of allowed iterations maxiter is set to \(10\ 000\), all other tuning parameters in the control argument are kept at their default values according to gsl_nls_control(). For comparison, we also compute single-start NLS model fits using the default Levenberg Marquardt algorithm (algorithm = "lm"), with as naive choice of starting values a vector of all ones \((1, \ldots, 1)\), similar to nls() when argument start is missing.

The table below displays the NLS model fit results for each individual test problem using the following status colors:

• success ; the NLS routine converged successfully and the residual sum-of-squares (SSR) under the fitted parameters coincides (up to a small numeric tolerance) with the SSR under the target parameter values⁴. The total runtime is displayed in seconds, timed on a modern laptop computer (Intel i7-8550U CPU, 1.80GHz) using a single core.
  
• false convergence ; the NLS routine converged successfully, but the residual SSR under the fitted parameters is larger than the SSR under the target parameter values.
• max. iterations ; the NLS routine failed to converge within the maximum number of allowed iterations.
• failed/non-zero exit ; the NLS routine failed to converge and returns an error or an NLS object with a non-zero exit code.

We observe that the naive single-start model fits manage to correctly fit about half of the test problems (27 out of 59), suggesting that these test problems are straightforward to optimize and do not require well-informed starting values. The multistart model fits using the double dogleg method improve upon the naive single-start model fits achieving correct convergence for 51 out of 59 test problems. The multistart Levenberg-Marquardt model fits correctly converge for a few more test problems (56 out of 59). Finally, the most robust results are obtained with the multistart model fits using the Levenberg-Marquardt algorithm with geodesic acceleration, which correctly fit all 59 test problems without initializing proper starting values or starting ranges! 🎉