Optimization in R

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