R devel - Aug 2007 - Optimization in R

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

Cheers,
Andrew

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.

Hi Pat,

On Sat, Aug 04, 2007 at 09:59:57AM +0100, Patrick Burns
wrote:
> Sounds like a good project.
Thanks :)
> How much extra overhead are you getting from the
> algorithm being in R?
On the Rosenbrock function (which is very quick to evaluate), here are the
system.time() results:
> system.time(bfgs(x0, f, g))
[1] 0.148 0.000 0.149 0.000 0.000
> system.time(optim(x0, f, g, method="BFGS"))
[1] 0.008 0.000 0.008 0.000 0.000

and the function evaluation counts:
> bfgs(x0, f, g)$counts
function gradient 
      95       58 
> optim(x0, f, g, method="BFGS")$counts
function gradient 
     318      100 

So the overhead is clearly much bigger, but still too small to matter for
most (?) applications.

Cheers,
Andrew

PS, my computer is a "Intel(R) Pentium(R) 4 CPU 2.80GHz" with a
1024 KB cache, according to /proc/cpuinfo.

Hi Andrew,
I have been working quite a bit with optim and friends on automated
nonlinear fitting, mainly for calibration. All of the optimizers
seem to have trouble in the tricky situation of fitting a log-logistic 
model when the upper asymptote is not well defined, and trying to
estimate a variance
parameter. Granted this is not necessarily the optimizer, but a
combination
of objective function, convergence criteria, scaling, ... 
But my point is that for automated fitting where many starting values,
iterations with
pseudo likelihood for the variance function, over multiple curves, the
overhead of
the function written in R vs C can become non-negligible for "simple"
univariate functions. 
So I have to agree with Duncan, that R is very good for prototyping but
C would be preferable.

I can contribute some "hard" problems if you start the package.

Nicholas
> Date: Sat, 4 Aug 2007 01:12:31 -0400
> From: Andrew Clausen <clausen at econ.upenn.edu>
> Subject: [Rd] Optimization in R
> To: r-devel at r-project.org
> Cc: help-gsl at gnu.org
> Message-ID: <20070804051231.GB3016 at econ.upenn.edu>
> Content-Type: text/plain; charset=unknown-8bit
> 
> 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".
> 
> Cheers,
> Andrew
> 
> 
> 
> ------------------------------
> 
> _______________________________________________
> R-devel at r-project.org mailing list  DIGESTED
> https://stat.ethz.ch/mailman/listinfo/r-devel
> 
> 
> End of R-devel Digest, Vol 54, Issue 4
> **************************************

Hi,
In my earlier post I eluded to a situation where that would be useful. 
In nlme, there is a choice of optimizers, minpack.lm has
Levenberg-Marquardt,
while nlminb has the port routines. For the same starting values,
different
optimizers will present different solutions, having a common interface
would make fitting with multiple optimizers very attractive. 

Also the inverse Hessian would be useful, for cases where the Hessian is
ill conditioned
a little regularization goes a long way. I believe the package accuracy 
has a nice solution.

Nicholas
 

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

Hi guys,

I ran into the same problem and solved it by changing one line in Andrew's
Makefile:

instead of:
gsl.so: vector.o multimin.o
    gcc -g -Wall -shared $^ -lgsl -o $@

use:
gsl.so: vector.o multimin.o
    gcc -g -Wall -shared $^  -o $@ -lgsl -lgslcblas -lm

Christophe


Thanks for the functions!
>
> I tried installing the multimin function. To get it to compile, I needed
> to change the Makefile to reflect my path and by adding the flags fPIC
> in response to the error: "/usr/bin/ld: vector.o: relocation
R_X86_64_32
> against `a local symbol' can not be used when making a shared object;
> recompile with -fPIC"
>
> However, I get the following running test.R:
>
> Error in dyn.load(x, as.logical(local), as.logical(now)) :
>        unable to load shared library '/home/mmorales/Desktop
> /multimin/gsl.so':
>  /usr/lib64/libgsl.so.0: undefined symbol: cblas_ctrmv
>
> I'm running R-2.5.1 compiled for x86_64 with a custom built ATLAS.
>

-- 
A Master Carpenter has many tools and is expert with most of them. If you
only know how to use a hammer, every problem starts to look like a nail.
Stay away from that trap.
Richard B Johnson.
--

Christophe Pouzat
Laboratoire de Physiologie Cerebrale
CNRS UMR 8118
UFR biomedicale de l'Universite Paris V
45, rue des Saints Peres
75006 PARIS
France

tel: +33 (0)1 42 86 38 28
fax: +33 (0)1 42 86 38 30

	[[alternative HTML version deleted]]

I am sorry for omitting a citation in my previous post.
The complete message is as follows (my text unchanged). PS

I would like to add a remark and a question.

Remark.

There is a part of R, which allows the user to select
among several methods for the same task and also to add
his own C code: random number generation. However, the interface
for optimization is more complex. In my opinion, looking
for a unified interface for this is desirable, but it is
a research problem, not a suggestion for an immediate
code modification.

Question.

Is there a way how to optimize a function written in C
using optim? This would be very useful, if the optimization
needs a lot of iterations. This may be done by defining 
an R function, which does nothing more than calling .C with
appropriate parameters,
but this looses efficiency. A more efficient solution
could be adding a specified entry point (or several, if derivatives
are also available), similar as in the user defined random number
generator. Then, a parameter of optim could control, whether
the function to be optimized is fn or the C entry point.

Petr Savicky.

On Sat, Aug 04, 2007 at 06:56:47PM -0400, Gabor Grothendieck
wrote:
> 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
> > >

Those are small parts of the calculation, not the whole thing.  The 
original point was that optim() is a very thin wrapper around the code 
to do the optimization.  I just don't see a need to make it more 
complicated so it can be used to wrap other methods.  Authors of new 
optimization methods can just create new functions, following the 
pattern set by optim(), and it will be easier for almost everyone.

Duncan Murdoch

On 06/08/2007 11:18 PM, Andrew Robinson wrote:
>  ... Variance and correlation model classes in nlme.
> 
> Cheers
> 
> Andrew
> 
> On Mon, Aug 06, 2007 at 09:55:38PM -0500, hadley wickham wrote:
>> 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.
>> Bin width selection in hist?  Family functions for glm?  Those come
>> quickly to my mind, but I'm sure there are others.
>>
>> Hadley
>>
>> ______________________________________________
>> R-devel at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-devel
>