4 ms·

Here's a gist that optimizes the author's Julia code to give a 4x speedup over the unvectorized code, and handily beat the vectorized codes too: https://gist.g

Here's a gist that optimizes the author's Julia code to give a 4x speedup over the unvectorized code, and handily beat the vectorized codes too: https://gist.github.com/jwmerrill/9715447 https://gist.github.com/jwmerrill/9715447 In general, carefully written devectorized Julia will be faster than equivalent vectorized code because it is possible to avoid allocating containers at intermediate stages of the computation. This is surprising for people coming from e.g. Matlab, R, or Numpy, because vectorized code is often faster in those environments. In those languages, direct scalar code is slow, and vectorization is about expressing the computation in a form that will call out to fast C operations, but the C still pays a price to allocate intermediates. In Julia, the JIT is able to emit fast machine code for scalar operations, so you don't need to express your code in a form that calls out to C. There are JIT compilers for Matlab and Python too, but Julia's JIT is very high quality and everyone uses it all the time. The main changes I made here were hoisting computations of constants out of the loop, and preallocating a 2D Float64 array to hold the results. These were worth roughly 2x each. I've chosen the storage format so that it is accessed in column order. Both optimizations have their own sections in the Julia Manual's performance tips section: http://julia.readthedocs.org/en/latest/manual/performance-tips/ http://julia.readthedocs.org/en/latest/manual/performance-ti...

shoo 13y ago
your point at the end is good - if we're only interested in the final price, then we should be able to do this in constant memory, independent of the number of timesteps and samples. edit: on my machine the running time seems dominated by the generation of random numbers, rather than anything else. perhaps another optimisation is to elimate the `exp` from the inner loop, and just do addition. if you're only interested in the final price, you can do one `exp` at the end. i.e. change paths[1, i] = st = s0 for j in 2:M st *= exp(a + b*randn()) paths[j, i] = st end to this paths[1, i] = log_st = log(s0) for j in 2:M log_st += (a + b*randn()) paths[j, i] = log_st end
- Retric 13y ago
  One solution to this is use a circular data structure filled with random numbers. Then dedicate a thread to updating that ring with new random numbers. Meanwhile you pull random numbers as fast as you want. It's not good for encryption, but for simulation it's fine. http://en.wikipedia.org/wiki/Circular_buffer http://en.wikipedia.org/wiki/Circular_buffer Another approach is to keep reseeding a fast RNG with a slow but higher quality RNG.
  - StefanKarpinski 13y ago
    You don't have to do anything quite that complicated. You can just pre-generate an array of Gaussian random samples and do table lookup in the inner loop: paths = zeros(Float64, M, I) rands = zeros(Float64, M-1, 1) for i in 1:I randn!(rands) paths[1, i] = st = s0 for j in 2:M st *= exp(a + b*rands[j-1]) paths[j, i] = st end end That shaves off another 35% percent of runtime on my machine.
    - jwmerrill 13y ago
      Nice! In case it isn't obvious, randn! takes an existing vector and fills it up with new random normal variates. This allows computing a bunch of random variates at once, which is a little more efficient, without having to allocate a new vector every time through the loop. Another option would be to generate all M-1xI random numbers ahead of time, but that would be less memory efficient.
      - StefanKarpinski 13y ago
        I tried that too, but it wasn't any faster – you've already gotten most of the gains to be had by lifting the random number generation out of the innermost loop. So I went with the option that uses less memory.
      - shoo 13y ago
        i tried the approach of generating them ahead of time and then loading them from disk. it only uses 80mb of memory for this example. the difference i see between generating numbers on-the-fly using boost and reading them from disk is 6.41s vs 0.32s, ie a 20x speedup, for 10 repetitions of 100,000 paths over 100 timesteps. edit: here's the code - https://gist.github.com/fcostin/3ed3512ad167adf9c785 https://gist.github.com/fcostin/3ed3512ad167adf9c785

shoo 13y ago

your point at the end is good - if we're only interested in the final price, then we should be able to do this in constant memory, independent of the number of timesteps and samples. edit: on my machine the running time seems dominated by the generation of random numbers, rather than anything else. perhaps another optimisation is to elimate the `exp` from the inner loop, and just do addition. if you're only interested in the final price, you can do one `exp` at the end. i.e. change paths[1, i] = st = s0 for j in 2:M st *= exp(a + b*randn()) paths[j, i] = st end to this paths[1, i] = log_st = log(s0) for j in 2:M log_st += (a + b*randn()) paths[j, i] = log_st end

Retric 13y ago
One solution to this is use a circular data structure filled with random numbers. Then dedicate a thread to updating that ring with new random numbers. Meanwhile you pull random numbers as fast as you want. It's not good for encryption, but for simulation it's fine. http://en.wikipedia.org/wiki/Circular_buffer http://en.wikipedia.org/wiki/Circular_buffer Another approach is to keep reseeding a fast RNG with a slow but higher quality RNG.
- StefanKarpinski 13y ago
  You don't have to do anything quite that complicated. You can just pre-generate an array of Gaussian random samples and do table lookup in the inner loop: paths = zeros(Float64, M, I) rands = zeros(Float64, M-1, 1) for i in 1:I randn!(rands) paths[1, i] = st = s0 for j in 2:M st *= exp(a + b*rands[j-1]) paths[j, i] = st end end That shaves off another 35% percent of runtime on my machine.
  - jwmerrill 13y ago
    Nice! In case it isn't obvious, randn! takes an existing vector and fills it up with new random normal variates. This allows computing a bunch of random variates at once, which is a little more efficient, without having to allocate a new vector every time through the loop. Another option would be to generate all M-1xI random numbers ahead of time, but that would be less memory efficient.
    - StefanKarpinski 13y ago
      I tried that too, but it wasn't any faster – you've already gotten most of the gains to be had by lifting the random number generation out of the innermost loop. So I went with the option that uses less memory.
    - shoo 13y ago
      i tried the approach of generating them ahead of time and then loading them from disk. it only uses 80mb of memory for this example. the difference i see between generating numbers on-the-fly using boost and reading them from disk is 6.41s vs 0.32s, ie a 20x speedup, for 10 repetitions of 100,000 paths over 100 timesteps. edit: here's the code - https://gist.github.com/fcostin/3ed3512ad167adf9c785 https://gist.github.com/fcostin/3ed3512ad167adf9c785