Re: [Rd] Optimization in R

From: Duncan Murdoch <>
Date: Sat, 04 Aug 2007 15:13:17 -0400

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 <> 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 <> 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 <> 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
>>>>>>> ______________________________________________
>>>>>>> mailing list
>>>>>> ______________________________________________
>>>>>> mailing list
>>>>> ______________________________________________
>>>>> mailing list
> ______________________________________________
> mailing list

______________________________________________ mailing list Received on Sat 04 Aug 2007 - 19:16:43 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Tue 07 Aug 2007 - 14:37:38 GMT.

Mailing list information is available at Please read the posting guide before posting to the list.