[Rd] Optimization in R

[Duncan Murdoch]

On 04/08/2007 2:23 PM, Gabor Grothendieck wrote:
> For the same reason that generic functions exist.  They don't have
> a lot of common code but it makes easier to use.  Perhaps the argument
> is not as strong here since the class tends to be implicit whereas the
> method is explicit but it would still be a convenience.

Can you give other examples where we do this?  The ones I can think of 
(graphics drivers and finalizers) involve a large amount of common 
machinery that it's difficult or impossible for the user to duplicate. 
That's not the case here.

Duncan Murdoch

> On 8/4/07, Duncan Murdoch <murdoch at stats.uwo.ca> wrote:
>> On 04/08/2007 1:05 PM, Gabor Grothendieck wrote:
>>> I think it would be desirable for optim to have a dispatching mechanism
>>> that allows users to add their own optimization techniques to those
>>> provided without having to modify optim and without having to come
>>> up with a new visible function.  For example, if we call optim as
>>> optim(...whatever..., method = "X") then that might dispatch
>>> optim.X consistent with S3 except here X is explicitly given rather
>>> than being a class.
>> Why?  This would make sense if optim() had a lot of code common to all
>> methods, but it doesn't.
>>
>> Duncan Murdoch
>>
>>> On 8/4/07, Duncan Murdoch <murdoch at stats.uwo.ca> wrote:
>>>> 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
>>>> ______________________________________________
>>>> R-devel at r-project.org mailing list
>>>> https://stat.ethz.ch/mailman/listinfo/r-devel
>>>>
>>> ______________________________________________
>>> R-devel at r-project.org mailing list
>>> https://stat.ethz.ch/mailman/listinfo/r-devel
>>