On Thu, Mar 14, 2013 at 12:08 PM, Ben Bolker <bbolker@gmail.com> wrote:

> Ross Boylan <ross@...> writes:
>
> >
> > lme4 appears to use Nelder Mead for the final stage of its
> > optimization.  The archives show various concerns about stability, but
> > it looks as if even the prior optimizer was in the simplex class.
>
>   Depends how far you go back.  nlme and older versions of lme4 used
> nlminb.  This is derivative-based, but you'd have to look at the guts
> of nlminb and at the PORT documentation for a detailed description ...
>
>  http://netlib.bell-labs.com/cm/cs/cstr/153.pdf
>
> I think it's quasi-Newton ... digging around more in the R code
> guts shows it's NL2SOL -- you could dig some more on Netlib to
> find the details
>

Not really.  NL2SOL is part of the Port collection of optimizers but nlminb
uses another of the Port Fortran routines.

Originally the R version of the nlme package used the optimizer that Ross
wrote based on Dennis and Schnabel's code and the S-PLUS version of nlme
used nlminb based on the Port code from Bell Labs.  There were some
embarassing instances where the S-PLUS version would converge when the R
version either didn't converge or converged to a false optimum.

When the Port code became publicly available I wrote the nlminb function
using it.

The optimization in lme4 requires a constrained optimizer (some of the
parameters must be positive) whereas nlme can use an unconstrained
optimizer.  The non-negativity constraints were transformed away in nlme by
taking logarithms but it turns out that is not a good thing to do.
 Initially I tried the L-BFGS algorithm in optim but it has problems.  It
may have been subsequently fixed but at the time I was using it there were
cases where it converged to a minimum for which the function value was
greater than the value at the starting values.  So we used nlminb for the
constrained optimization.  Then that had problems and we switched to bobyqa
from the minqa package.  Then that had problems so I coded up an
Nelder-Mead based on the C code in the nlopt package from MIT.  It has
generally proven useful although none of these optimizers are perfect.

>
> > This surprised me, since I have found Nelder Mead to be relatively slow
> > and imprecise, and in a more complex mixed model with sampling I've been
> > using optim with L-BFGS-B, which is quasi-newton.  One of the parameter
> > estimates kept creeping up over iterations; the bounds were necessary to
> > cut it off.
>
>   Not quite sure what you mean here.  I have personally found L-BFGS-B
> to be somewhat unstable in my experience ...  The optimplus package
> on R-forge gives a lot more optimizer choices ...
>
> > If anyone can give me more insight into why Nelder Mead is in use, and
> > perhaps whether I should be using it myself, I'd appreciate it.
>
>   Nelder-Mead *is* fairly slow, but it tends to be slightly more robust
> over a broad class of problems.  The current development version of lme4
> lets you pick your own optimizer, so you can make the speed vs
> robustness tradeoff however you like ...
>
> > The current behavior, in which about 10% of the simulated datasets have
> > convergence problems, with one of the parameters heading toward an
> > implausible value (a correlation of 1.0, where atanh(rho) is the actual
> > parameter being estimated) is certainly not ideal.  I tried simulated
> > annealing to see if the algorithm had wandered mistakenly toward a local
> > rather than global optimum; there was no evidence that it had (the SA
> > ended up in the same neighborhood as BFGS, and when that point was the
> > starting value for BFGS it ended on essentially the same point as
> before).
>
>    Can you send me a simulated example that experiences this problem?
>
>   I can't tell at the moment whether you are using lme4 or your
> own hand-rolled solution to the problem.
>
>   There are two possibilities here: (1) the best-fit solution really
> is in the limit as rho -> Inf / correlation -> 1.0; this is similar to
> the case where the MLE variance of a random effect is zero, even when
> the true model (from which the data were simulated) has a non-zero
> variance; (2) the best-fit solution is very close to +/- 1.0.
>
>   It's going to be tough either way, but in case #1 you would
> definitely be better off with a parameter space in which the
> singular cases corresponded to _simple_ box constraints on the
> parameters.   For the Cholesky parameterization this is true
> for the variances (for a 2-by-2 matrix, theta_1=0 and theta_2=theta_3=0
> are the conditions for the variances to be 0), and (mostly) for
> the correlations -- the correlations go to +/- 1 when theta_3=0
> *or* when theta_2 -> infinity, but in the latter case the variance
> will blow up to infinity as well (so probably not a case worth
> worrying about).
>
> _______________________________________________
> R-sig-mixed-models@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>

	[[alternative HTML version deleted]]

