Variable results depending on optimisation algorithms

[Dominique Fauteux]

Hello everyone,

 

I have met a peculiar situation with our lemming dataset on Bylot Island where we are fitting full-session models where one session is one trapping session (6 secondary periods; PP1=June_2004, PP2=June_2005, PP3=June_2006, etc.). The reason why we are fitting full-session models is because we subsequently add habitat covariates to test their effects directly on density. With these full-session models, we isolate the natural yearly variation in density and obtain coefficient, which is much easier to interpret than interpreting outputs of 10 single-session models.

 

Because we sometimes ended up with NaNs in the outputs, we tried different starting values without much success. Then, we decided to try different optimisation algorithms and we ended up with this AICc comparison table:

 

### Test of different optimization algorithms for full-session models

> AIC(model.set.full)

 

Newton_Raphson       D~session g0~session sigma~session halfnormal   21 -1426.652 2895.305 2900.438   0.000 0.1335

BFGS                 D~session g0~session sigma~session halfnormal   21 -1426.652 2895.305 2900.438   0.000 0.1335

BFGS_Autoini5        D~session g0~session sigma~session halfnormal   21 -1426.652 2895.305 2900.438   0.000 0.1335

BFGS_autoini2        D~session g0~session sigma~session halfnormal   21 -1426.652 2895.305 2900.438   0.000 0.1335

L_BFGS_B             D~session g0~session sigma~session halfnormal   21 -1426.653 2895.306 2900.439   0.001 0.1334

L_BFGS_B_autoini2    D~session g0~session sigma~session halfnormal   21 -1426.656 2895.312 2900.445   0.007 0.1330

L_BFGS_B_autoini5    D~session g0~session sigma~session halfnormal   21 -1426.663 2895.325 2900.458   0.020 0.1322

CG                   D~session g0~session sigma~session halfnormal   21 -1427.337 2896.673 2901.807   1.369 0.0673

Nelder_Mead_Autoini5 D~session g0~session sigma~session halfnormal   21 -1436.004 2914.008 2919.141  18.703 0.0000

Nelder_Mead          D~session g0~session sigma~session halfnormal   21 -1447.865 2937.730 2942.863  42.425 0.0000

Nelder_Mead_Autoini2 D~session g0~session sigma~session halfnormal   21 -1454.306 2950.612 2955.746  55.308 0.0000

SANN                 D~session g0~session sigma~session halfnormal   21 -1709.528 3461.056 3466.190 565.752 0.0000

 

Based on these results, we discovered that the Nelder-Mead algorithm that we initially used was less robust (falling more easily in local minima) to variable starting values than the BFGS and L-BFGS-B algorithms that seem to have all reached the same minimum. In contrast, all the Nelder-Mead models we tested here gave very different log-likelihood values. Here and in other tests with our dataset, the Newton-Raphson kept giving NaNs in the outputs (being our first evidence of future difficulties...). Seeing the stability of the results with the BFGS compared to Nelder-Mead, we switched to BFGS for our analyses and L-BFGS-B as a "rescue" algorithm in case of problems with BFGS due to the highly similar outputs.

 

Is this robustness of BFGS to variable starting values been observed by other users of the forum? Anyone may have a clue why the Nelder-Mead may be less robust in some circumstances?

 

For the sake of exploration, I’ve done the same comparison with a simpler model with the Huggins paraterization (CL=T):

 

### Test of different optimization algorithms for Huggins parameterization (CL=TRUE) models

> AIC(model.set.huggins)

                                            model   detectfn npar    logLik      AIC     AICc  dAICc AICcwt

Newton.Raphson        g0~session sigma~session halfnormal   14 -1409.076 2846.151 2848.397  0.000 0.1320

BFGS                  g0~session sigma~session halfnormal   14 -1409.076 2846.151 2848.397  0.000 0.1320

L.BFGS.B              g0~session sigma~session halfnormal   14 -1409.076 2846.151 2848.397  0.000 0.1320

BFGS.autoini5         g0~session sigma~session halfnormal   14 -1409.076 2846.151 2848.397  0.000 0.1320

BFGS.autoini2         g0~session sigma~session halfnormal   14 -1409.076 2846.151 2848.397  0.000 0.1320

L.BFGS.B.autoini5     g0~session sigma~session halfnormal   14 -1409.076 2846.151 2848.397  0.000 0.1320

L.BFGS.B.autoini2     g0~session sigma~session halfnormal   14 -1409.076 2846.151 2848.397  0.000 0.1320

CG                    g0~session sigma~session halfnormal   14 -1409.631 2847.262 2849.508  1.111 0.0758

Nelder.Mead           g0~session sigma~session halfnormal   14 -1423.227 2874.454 2876.700 28.303 0.0000

Nelder.Mead.autoini2  g0~session sigma~session halfnormal   14 -1423.352 2874.703 2876.949 28.552 0.0000

Nelder.Mead.autoini5  g0~session sigma~session halfnormal   14 -1425.903 2879.806 2882.052 33.655 0.0000

 

I basically end up in a similar situation with the Newton-Raphson algorithm that leads to NaNs in the output, BFGS and L-BFGS-B with no NaNs in the outputs and being more robust to starting values than the Nelder-Mead algorithm. Here, I did not run the SANN based on the bad performance for the more complex models above and because it took much longer to run than the other models.

 

I do not have an answer to explain these differences, but the discrepancies may be specific to the dataset. WIth lemmings, we are stuck with few individuals per trapping sessions being between 7 to 51 individuals (21 to 64 captures respectively) in this particular case. 

Could the mask be a problem for some algorithms? Here are the characteristics of our custom mask:

traps.m <- make.mask(traps(lemming.ch), type = 'trapbuffer', buffer = 100, nx=128, ny=128)

traps.m.cov <- addCovariates(traps.m, traps(lemming.ch), columns = NULL, strict = FALSE, replace = FALSE)

summary(traps.m.cov)

Object class      mask

Mask type         trapbuffer

Number of points  15884

Spacing m         4.140625

Cell area ha      0.001714478

Total area ha     27.23276

x-range m         -97.92969 427.9297

y-range m         -97.92969 427.9297

Bounding box     

     x    y

1 -100 -100

2  430 -100

4  430  430

3 -100  430

 

If needed, I will provide more information about the dataset, no problem.

 

Thank you very much in advance for your time!

Dominique

[Murray Efford]

Hello Dominique

Thanks for laying out the problem. I can't add much. You do seem to have a dataset and model for which Nelder-Mead doesn't perform well. I give a counter example below just to reassure people that the different optimization methods can give the same answer (and as skeleton code for testing other datasets, anyone?).

library(secr)
CH <- subset(ovenCHp, sessions = 1:3) # a small example
msk20 <- make.mask(traps(CH[[1]]), buffer = 300, spacing = 20, type = 'trapbuffer')
methods <- c("Newton-Raphson", "BFGS", "Nelder-Mead")
fits <- list()
for (m in 1:length(methods)) {
    fits[[m]] <- secr.fit(CH, model = list(g0~session, sigma~session), CL = TRUE,
             method = methods[m], trace = FALSE, mask = msk20)
}
names(fits) <- methods
fits <- secrlist(fits)
AIC(fits, criterion = 'AIC')[,-c(2,6)]
#                                   model npar    logLik     AIC dAIC  AICwt
# Newton-Raphson g0~session sigma~session    6 -371.0973 754.195    0 0.3333
# BFGS           g0~session sigma~session    6 -371.0973 754.195    0 0.3333
# Nelder-Mead    g0~session sigma~session    6 -371.0973 754.195    0 0.3333

Some minor comments:

Given that the problem appears with CL  = TRUE it makes sense to troubleshoot that simpler model.

If the NaN values from Newton-Raphson are just in the variances then it may help to use a better algorithm (fdhess) for that final variance calculation. Do this with e.g.,

secr.fit(CH, model = list(g0~session, sigma~session), CL = TRUE,
         method = "Newton-Raphson", trace = FALSE, mask = msk20,
         details = list(hessian = "fdhess"))

The argument ‘ny’ is not used in make.mask() for type = 'trapbuffer'. If you do not specify spacing (recommended) then it is computed from ‘nx’ and applied also in the y-direction. Also, you are using a mask with a large number of points - is this really necessary? See secr-habitatmasks.pdf.

Murray

[Sebastien Comte]

Hi everyone,

Murray, I tried your code with one of my dataset on Tasmanian devils. 

I used a three session dataset with varying effort:

session1: 75obs  20traps 9occ

session2: 92obs  20traps 8occ

session3: 94obs  23traps 8occ

I have the same outputs with the three optimization methods:

names(fits) <- methods
fits <- secrlist(fits)
AIC(fits, criterion = 'AIC')[,-c(2,6)]
#                                   model npar    logLik     AIC dAIC  AICwt

# Newton-Raphson g0~session sigma~session    6 -799.9428 1611.886    0 0.3333
# BFGS           g0~session sigma~session    6 -799.9428 1611.886    0 0.3333
# Nelder-Mead    g0~session sigma~session    6 -799.9428 1611.886    0 0.3333

Cheers,
Sebastien

---

From: secr...@googlegroups.com <secr...@googlegroups.com> on behalf of Murray Efford <murray...@gmail.com>
Sent: Wednesday, June 3, 2020 8:27 AM
To: secr <secr...@googlegroups.com>
Subject: Re: Variable results depending on optimisation algorithms

 

--
You received this message because you are subscribed to the Google Groups "secr" group.
To unsubscribe from this group and stop receiving emails from it, send an email to secrgroup+...@googlegroups.com.
To view this discussion on the web, visit https://groups.google.com/d/msgid/secrgroup/a3cbd3f6-a15e-4fbd-b49e-aa8936ee2e33%40googlegroups.com.

University of Tasmania Electronic Communications Policy (December, 2014).
This email is confidential, and is for the intended recipient only. Access, disclosure, copying, distribution, or reliance on any of it by anyone outside the intended recipient organisation is prohibited and may be a criminal offence. Please delete if obtained in error and email confirmation to the sender. The views expressed in this email are not necessarily the views of the University of Tasmania, unless clearly intended otherwise.

[Dominique Fauteux]

Dear Dr. Efford,

I also tried your code with our lemming data with a tiny tweak as I used our mask but without any given nx value and I added a 4th model using L-BFGS-B for consistency with my previous message. Seems Nelder-Mead still gives a different logLikelihood compared to the other methods, which confirms that the problem is specific to the dataset:

                                       model npar    logLik      AIC  dAIC  AICwt

fits.Newton-Raphson g0~session sigma~session    6 -548.7399 1109.480 0.000 0.2568

fits.BFGS           g0~session sigma~session    6 -548.7399 1109.480 0.000 0.2568

fits.L-BFGS-B       g0~session sigma~session    6 -548.7399 1109.480 0.000 0.2568

fits.Nelder-Mead    g0~session sigma~session    6 -548.8513 1109.703 0.223 0.2297

As for the reason why we've given nx values (thanks for the correction about the ny!) is purely for graphical purposes to have smoother countours. I have tested our mask with and without given nx values on our full-session model with a habitat covariate and it did not change the results.

Also, thank you for the advice about changing the algorithm to calcule the hessian, it did help eliminate all NaNs in my few tests for now when using the Newton-Raphson optimisation.

Thank you very much for your help!

Dominique

Code for the results above:

CH <- subset(lemming.ch, sessions = 1:3) # a small example

msk20 <- make.mask(traps(lemming.ch[[1]]), buffer = 100, type = 'trapbuffer')

methods <- c("Newton-Raphson", "BFGS", "Nelder-Mead", "L-BFGS-B")

[Murray Efford]

Yes, it's good to extend the list of methods for comparison. There was nothing special about my code - your table was convincing enough!

Thanks to Sebastien for trying another dataset: we need to keep looking for data that trip up Nelder-Mead.

My reference to nx in make.mask() was ambiguous: I meant that it's good to specify 'spacing' directly (and report in papers), rather than rely on 'nx'.

I'm mystified by your reference to contours. Maybe that relates to plotting covariate effects later; usually it is possible to substitute a finer mask at the plotting stage.

Murray