[Rd] Optimization in R

[Duncan Murdoch]

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