Optimization in R
On 04/08/2007 10:30 AM, Andrew Clausen wrote:
Hi Duncan, On Sat, Aug 04, 2007 at 09:25:39AM -0400, Duncan Murdoch wrote:
This is interesting work; thanks for doing it. Could I make a suggestion? Why not put together a package containing those test optimization problems, and offer to include other interesting ones as they arise?
Good idea.
You could also include your wrapper for the gsl function and your own improvements to optim.
I have submitted my gsl multimin() wrapper for inclusion into the Rgsl package, which seems to be the "right" place for it. (Its maintainer is currently enjoying a vacation :) It would be nice if all of these methods could be accessed with the existing optim() interface, so that users could easily try different algorithms. That is, users could specify method="BFGS" or "BFGS-R" or "BFGS-GSL", etc. Is there a convenient mechanism for packages registering new methods?
No there isn't, but I don't really see it as necessary. The R optim() function is reasonably short, mainly setting up defaults specific to each of the supported optimization methods. The C do_optim() function has a bit more code common to the methods, but still only about 50 lines. You could copy this common part into your own function following the same argument and return value conventions and a user would just need to change the function name in addition to specifying your method, not really that much harder. (You could call your function optim and use it as a wrapper for stats::optim if you wanted to avoid even this, but I wouldn't recommend it, since then behaviour would depend on the search list ordering.) There's a small risk that the argument list to optim() will change incompatibly sometime in the future, but I think that's unlikely.
One incompatibility with my BFGS implementation is that it returns the *inverse* Hessian, which is a natural by-product of the BFGS algorithm. Indeed, R's existing BFGS implementation also calculates the inverse Hessian, and discards it! (Disclaimer: as far as I know, there are no theorems that say that BFGS's inverse Hessians are any good. In practice, they seem to be.)
If you return more than optim() returns it shouldn't have any serious ill effects.
The inverse Hessian is more useful than the Hessian for statistics because it gives the variance-covariance matrix for maximum likelihood estimators. When the Hessian is close to being singular (aka "computationally singular"), solve() sometimes fails to invert it. I think this means we should change the optim() interface. For example, an extra logical parameter, "inv.hessian" could control whether an inv.hessian matrix is returned.
I'd suggest always returning it if it's useful and the calculation is reliable and cheap. Adding "inv.hessian" as a general parameter would be troublesome with some of the other methods, where the inverse Hessian isn't already calculated, because of the inversion problem you mention above. Duncan Murdoch
On your first point: I agree that a prototype implementation in R makes sense, but I suspect there exist problems where the overhead would not be negligible (e.g. ones where the user has written the objective function in C for speed). So I think you should keep in mind the possibility of moving the core of your improvements to C once you are happy with them.
Fair enough. Cheers, Andrew
______________________________________________ R-devel at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel