[Rd] Optimization in R

[Duncan Murdoch]

On 04/08/2007 2:53 PM, Gabor Grothendieck wrote:
> The example of generic functions.

Show me an example where we have a list of ways to do a calculation 
passed as an argument (analogous to the method argument of optim), where 
the user is allowed to add his own function to the list.

Duncan Murdoch
> 
> On 8/4/07, Duncan Murdoch <murdoch at stats.uwo.ca> wrote:
>> 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
>>
> 
> ______________________________________________
> R-devel at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel