R: solving a small system of equations fast
2011-02-20 · Original WordPress address
My google-fu failed me on this one, so for the next person googling "stupid fast library for solving a system of equations in R such that the trade off between speed and 'robustness' is way on the speed side", you're welcome.
I'm on the final leg of my dissertation and I have some code that loops through a discreet state-space with about 300k points. Each of these iterations requires a two-equation non-linear system to be solved. Because all of this is contained in a loop that is (supposed to be) converging on a pair of value functions, its not super critical that the solutions to this system are "correct" every time. I want speed, baby.
Here's the R libraries I tried:
- nleqslv with the functions written to have a high penalty: each iteration took 0.020 seconds with a standard deviation of 0.007 seconds
- optim using the "make the system a single quadratic" trick and having a high penalty: mean 0.253 seconds with a standard deviation of 0.130 seconds
- optim using the "make the system a single quadratic" trick and using the "L-BFGS-B" algorithm to deal with boundedness: mean 0.014 seconds with a standard deviation of 0.032 seconds
- pso (partical swarm optimization) with bounds: mean 20 seconds!
I was kind of excited about particle swarm optimization but it was waaaaaay slow (at least with the default settings). The best was optim (available with the R core installation) using the quadratic trick. The quadratic trick is just defining all your functions in the system with the right hand side equal to zero and then having the algorithm minimize the sum of squares of those equations. Interestingly the nleqslv library, which is built for solving "medium-sized" systems, did almost as well on average but had had a much better variance. I'm not sure if this is because nleqslv gives up earlier or if it does a better job staying out of pathological parts of the parameter space relative to optim.
Anyway, there you go. If you find your self solving a small system of equations lots and lots of times over and over again, optim should do the trick for you in terms of speed.
PS - R seems to be the wrong environment for this task. Fortran code that does a very similar thing runs about 10 times faster.