[Rd] Optimization in R

[Gabor Grothendieck]

The example of generic functions.

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
> >>
>
>