# The need for speed

Thanks for all the comments about the attributes of Python and Julia. It seems to me that the most prudent choice is to learn Python and Julia. However, what I would really like to know is just how fast these languages really are and here is the test. What I want to do is to fit networks of coupled ODEs (and  PDEs) to data using MCMC (see here). This means I need a language that loops fast. An example in pseudo-Matlab code would be

for n = 1:N

for i = 1:T

y(i+1) = M\y(i)

end

Compare to data and set new parameters

end

where h is a parameter and M is some matrix (say 1000 dimensional), which is sometimes a Toeplitz matrix but not always. Hence, in each time step I need to invert a matrix, which can depend on time so I can’t always precompute, and do a matrix multiplication. Then in each parameter setting step I need to sum an objective function like the mean square error over all the data points. The code to do this in C or Fortran can be pretty complicated because you have to keep track of all the indices and call linear algebra libraries. I thus want something that has the simple syntax of Matlab but is as fast as C. Python seems to be too slow for our needs but maybe we haven’t optimized the code. Julia seems like the perfect fit but let me know if I am just deluded.

## 8 thoughts on “The need for speed”

1. Rick G says:

I’d like to see your Python implementation, and then benchmark it against the way I would write it. I’m still skeptical that Julia will make an important difference (although, like MATLAB, that small block of code might be more readable).

Also, have you looked into STAN (http://mc-stan.org/)? I’m not sure if the likelihood can be expressed using the output of a system of differential equations, but in every other way it is supposed to dominate other MC implementations.

Like

2. @rick. Anything named after Ulam is definitely worth looking at. The manual says they are working on an ODE solver. It seems that I need to get out more often.

Like

3. robclewley says:

The fast loops might be easy to do with Cython, which has extensive access back to the high level Python objects that you might need to use in your calculation.

There was a recent “bake off” and investigation of MCMC in python on this blog: http://healthyalgorithms.com/2014/06/27/mcmc-in-python-another-thing-that-turns-out-to-be-hard/. The best part of the situation now is that it’s so easy to blend different tools together. For instance, the ODE solving tools in Matlab and Julia are extremely bare bones, but you can easily interface to Python for anything fancy. I can’t help you with PDEs though :O

Like

4. If you already know “R” and are doing mostly linear algebra, you can get orders of magnitude improvement in performance from a proper compilation with BLAS / ATLAS. (See Rick G’s comment on Python.)

Like

5. Josh Porter says:

It seems like you’re leaning away from Matlab, but just in case: a faster way to solve ODEs in Matlab is using the SUNDIALS software (http://computation.llnl.gov/casc/sundials/main.html), which has a Matlab interface. The speed is a huge improvement over Matlab’s built-in solvers.

Like