Question: I am trying to use R to estimate a multinomial logit model with a manual specification. I have found a few packages that allow you to estimate MNL models [here](http://hosho.ees.hokudai.ac.jp/~kubo/Rdoc/library/VGAM/html/multinomial.html) or [here](http://cran.r-project.org/web/packages/mlogit/index.html). 

I've found some other writings on "rolling" your own MLE function [here](http://www.mayin.org/ajayshah/KB/R/documents/mle/mle.html). However, from my digging around - all of these functions and packages rely on the internal `optim` function. 

In my benchmark tests, `optim` is the bottleneck. Using a simulated dataset with ~16000 observations and 7 parameters, R takes around 90 seconds on my machine. The equivalent model in [Biogeme](http://transp-or.epfl.ch/page63023.html) takes ~10 seconds. A colleague who writes his own code in [Ox](http://www.doornik.com/) reports around 4 seconds for this same model.

Does anyone have experience with writing their own MLE function or can point me in the direction of something that is optimized beyond the default `optim` function (no pun intended)? 

If anyone wants the R code to recreate the model, let me know - I'll glady provide it. I haven't provided it since it isn't directly relevant to the problem of optimizing the `optim` function and to preserve space...

*EDIT: Thanks to everyone for your thoughts. Based on a myriad of comments below, we were able to get R in the same ballpark as Biogeme for more complicated models, and R was actually faster for several smaller / simpler models that we ran. I think the long term solution to this problem is going to involve writing a separate maximization function that relies on a fortran or C library, but am certainly open to other approaches.*

Answer: 
Tried with the nlm() function already? Don't know if it's much faster, but it does improve speed. Also check the options. optim uses a slow algorithm as the default. You can gain a > 5-fold speedup by using the Quasi-Newton algorithm (method="BFGS") instead of the default. If you're not concerned too much about the last digits, you can also set the tolerance levels higher of nlm() to gain extra speed. 

```
f <- function(x) sum((x-1:length(x))^2)

a <- 1:5

system.time(replicate(500,
     optim(a,f)
))
   user  system elapsed 
   0.78    0.00    0.79 

system.time(replicate(500,
     optim(a,f,method="BFGS")
))
   user  system elapsed 
   0.11    0.00    0.11 

system.time(replicate(500,
     nlm(f,a)
))
   user  system elapsed 
   0.10    0.00    0.09 

system.time(replicate(500,
      nlm(f,a,steptol=1e-4,gradtol=1e-4)
))
   user  system elapsed 
   0.03    0.00    0.03 

```
Did you consider the material on the [CRAN Task View for Optimization](http://cran.r-project.org/web/views/Optimization.html) ?