From ggrothendieck at gmail.com  Sun Aug  5 00:56:47 2007
From: ggrothendieck at gmail.com (Gabor Grothendieck)
Date: Sat, 4 Aug 2007 18:56:47 -0400
Subject: [Rd] Optimization in R
In-Reply-To: <46B4CFCD.2000202@stats.uwo.ca>
References: <20070804051231.GB3016@econ.upenn.edu>
	<46B47E53.2040801@stats.uwo.ca>
	<20070804143032.GC16431@econ.upenn.edu> <46B4A200.10109@stats.uwo.ca>
	<971536df0708041005h55c2488by76709cf4fb989690@mail.gmail.com>
	<46B4C2B2.3030109@stats.uwo.ca>
	<971536df0708041123n5d59b1ebw18618d6340b76300@mail.gmail.com>
	<46B4C56D.8020408@stats.uwo.ca>
	<971536df0708041153o174d73f5vd2c28f483a6bf6b6@mail.gmail.com>
	<46B4CFCD.2000202@stats.uwo.ca>
Message-ID: <971536df0708041556q7faf64ber38bf0da574927adc@mail.gmail.com>

I don't have an example of that but that does not make it less
desirable.  If one wants to use method 1, 2 or 3 then one can
use optim with a method= but if one wants to use methods 4
or 5 then one must use an entirely different function.  Surely
it would be better to be consistent from the user's viewpoint
and allow all of them to work consistently through the same
interface.

On 8/4/07, Duncan Murdoch <murdoch at stats.uwo.ca> wrote:
> 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
>
>



[Rd] Optimization in R
Gabor Grothendieck