[Rd] Optimization in R

[Duncan Murdoch]

On 04/08/2007 1:12 AM, Andrew Clausen wrote:
> Hi all,
> 
> I've been working on improving R's optim() command, which does general purpose
> unconstrained optimization.  Obviously, this is important for many statistics
> computations, such as maximum likelihood, method of moments, etc.  I have
> focused my efforts of the BFGS method, mainly because it best matches my
> current projects.
> 
> Here's a quick summary of what I've done:
>  * implemented my own version of BFGS in R,
> 	http://www.econ.upenn.edu/~clausen/computing/bfgs.zip
>  * written a wrapper for the GNU Scientific Library's optimization function,
> multimin(),
> 	http://www.econ.upenn.edu/~clausen/computing/multimin.zip
>  * written some tricky functions to compare implementations,
> 	http://www.econ.upenn.edu/~clausen/computing/tests.zip
> 
> My own implementation has several advantages over optim()'s implementation
> (which you can see in the vmmin() function in
> 	https://svn.r-project.org/R/trunk/src/main/optim.c)
>  * the linesearch algorithm (More-Thuente) quickly finds a region of interest
> to zoom into.  Moreover, it strikes a much better balance between finding
> a point that adequately improves upon the old point, but doesn't waste too
> much time finding a much better point.  (Unlike optim(), it uses the standard
> Wolfe conditions with weak parameters.)
>  * the linesearch algorithm uses interpolation, so it finds an acceptable
> point more quickly.
>  * implements "box" constraints.
>  * easier to understand and modify the code, partly because it's written in R.
> 
> Of course, this comes at the (slight?) overhead cost of being written in R.
> 
> The test suite above takes the first few functions from the paper
> 
>         Moré, Garbow, and Hillstrom, "Testing Unconstrained
>         Optimization Software", ACM Trans Math Softw 7:1 (March 1981)
> 
> The test results appear below, where "*" means "computed the right solution",
> and "!" means "got stuck".
> 
> test                    optim           clausen         gsl
> --------------------------------------------------------------
> bard                                                    !
> beale
> brown-scaled
> freudenstein-roth
> gaussian                                                *
> helical-valley          *               *
> jennrich-sampson                                        *
> meyer                                                   *
> powell-scaled                           *
> rosenbrock                              *
> 
> The table indiciates that all three implementations of BFGS failed to compute
> the right answer in most cases.  I suppose this means they are all quite
> deficient.  Of course, this doesn't imply that they perform badly on real
> statistics problems -- but in my limited experience with my crude econometric
> models, they do perform badly.   Indeed, that's why I started investigating in
> the first place.
> 
> For what it's worth, I think:
>  * the optimization algorithms should be written in R -- the overhead is
> small compared to the cost of evaluating likelihood functions anyway, and is
> easily made up by the better algorithms that are possible.
>  * it would be useful to keep a repository of interesting optimization
> problems relevant to R users.  Then R developers can evaluate "improvements".

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?  You could also include your wrapper for the gsl function 
and your own improvements to optim.

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.

Duncan Murdoch

P.S. I dropped the help-gsl mailing list from the distribution, since my 
post is mainly about R.