Haskell: Efficiency iterative map with added noise

Question

I am wondering how to improve the time performance of the white noise addition to the logistic map? The noise is only allowed to be added after calculating the values (as it is an iterative map.)

 module Generating

import System.Random (Random,randomR,random,mkStdGen,StdGen)
import Data.Random.Normal (mkNormals,normal)
import qualified Data.Vector.Unboxed as U 
import Control.Monad.State

genR :: (Random a, Fractional a) => Int -> (a, StdGen)
genR x = randomR (0,1.0) (mkStdGen x)


new ::Double-> Double ->Int -> (Double,Int) -> U.Vector (Double,Int)
new skal r n = U.unfoldrN n go
  where  
   go (x0,g0)  = let  !eins= (1.0-x0)
                      !x=x0 `seq` eins `seq` r*x0*eins
                      !g=g0+1
                      !noise= skal*(fst $ genR g)
             in Just ((x0+noise,g0),(x,g))

fs :: (t, t1) -> t
fs (x,y)=x

first :: U.Vector (Double,Int)->U.Vector Double
first  =U.map (\(x,y)->x)  

As you can see, I actually only want the first value of the tuple, but the generator needs to be updated.

Any suggestions? Maybe State Monads?

Answer 1

tl;dr: Don't try to optimize a Haskell program without using profiling and benchmarking. Adding random exclamation marks and seqs is almost never going to work. The big problem here, in fact, is that StdGen is an incredibly slow random number generator and it's completely dominating the execution time of your program. You need to replace it to make any significant progress.

Here's the longer answer: A good first step is to install a benchmarking library, like criterion, and write a test case:

import Criterion.Main

...your program above...

vect1 :: (Double, Int) -> U.Vector Double
vect1 = first . new 0.5 1 10000

main = defaultMain [
  bench "vect1" $ nf vect1 (0,1)
  ]

In my case, the results look like:

benchmarking vect1
time                 8.097 ms   (8.071 ms .. 8.125 ms)
                     1.000 R²   (1.000 R² .. 1.000 R²)
mean                 8.140 ms   (8.124 ms .. 8.162 ms)
std dev              52.90 μs   (36.32 μs .. 91.72 μs)

so we have about 8 milliseconds per run to generate a 10000-element vector.

Now, let's get rid of all the bangs, seqs and intermediate calculations you added to try to speed things up:

new :: Double-> Double -> Int -> (Double, Int) -> U.Vector (Double,Int)
new skal r n = U.unfoldrN n go
  where  
   go (x0,g0)  = let  x = r * x0 * (1-x0)
                      g = g0 + 1
                      noise = skal * (fst $ genR g)
                 in Just ((x0+noise, g0), (x,g))

Rerunning, here are the results:

time                 8.195 ms   (8.168 ms .. 8.235 ms)

Ah, so they had no effect at all. Glad we got rid of them.

Now, it's worth noting that unfoldrN is carrying along the accumulator for you which holds on to your g. You don't also need to include g in the result if you're going to throw it away anyway, so we can simplify new to:

new :: Double-> Double -> Int -> (Double, Int) -> U.Vector Double
new skal r n = U.unfoldrN n go
  where  
   go (x0,g0)  = let  x = r * x0 * (1-x0)
                      g = g0 + 1
                      noise = skal * (fst $ genR g)
                 in Just (x0+noise, (x,g))

and drop the first call from the definition of vect1:

vect1 :: (Double, Int) -> U.Vector Double
vect1 = new 0.5 1 10000

This gives:

time                 8.289 ms   (8.238 ms .. 8.373 ms)

so it didn't really make a difference. Undoubtedly, the compiler was able to optimize away the useless extra Doubles anyway, so changing the code had no effect.

A more serious problem with the algorithm is that it uses generators in a really weird way. A StdGen is meant to be seeded and then reused to generate multiple random numbers, not to be generated fresh from a seed based on a counter. We really ought to rewrite new to use the generator properly:

new :: Double-> Double -> Int -> (Double, Int) -> U.Vector Double
new skal r n (x0, seed) = U.unfoldrN n go (x0, g0)
  where
   g0 = mkStdGen seed  -- create initial generator from seed
   go (x0,g0)  = let  (eps, g) = randomR (0, 1.0) g0 -- use generator properly
                      x = r * x0 * (1-x0)
                      noise = skal * eps
                 in Just (x0 + noise, (x, g))

though again, this makes almost no difference to our benchmarking times. I'll admit that this one surprised me. I thought it would have a significant effect. Good thing I was benchmarking these changes so I had actual objective evidence of the effect (or lack of effect) of this change!

Now, it seems like it's probably time to profile our program and see what it's spending its time doing.

$ stack ghc -- -prof -fprof-auto -O2 Generating.hs
$ ./Generating -n 100 +RTS -p   # run 100 iterations

If you look at the Generating.prof file that's output, you'll see that the majority of time is spent in System.Random, like so:

COST CENTRE               MODULE                            SRC                                                    %time %alloc

randomR                   System.Random                     System/Random.hs:409:3-27                               21.7   24.0
stdNext                   System.Random                     System/Random.hs:(518,1)-(528,64)                       15.4   16.6
randomIvalInteger         System.Random                     System/Random.hs:(468,1)-(489,76)                       12.2   12.0
randomIvalInteger.f       System.Random                     System/Random.hs:(486,8)-(489,76)                       11.0    4.8
randomIvalInteger.f.v'    System.Random                     System/Random.hs:489:25-76                               7.0    8.6

It turns out that Haskell's standard random number generator is appallingly slow, and we'll need to replace it with something faster to make any more progress.

The mersenne-random-pure64 package provides a fast Mersenne Twister implementation that produces high quality random numbers, and we can rewrite new to use it. Note that randomDouble returns a uniform random number in the interval [0,1):

import System.Random.Mersenne.Pure64
new :: Double-> Double -> Int -> (Double, Int) -> U.Vector Double
new skal r n (x0, seed) = U.unfoldrN n go (x0, g0)
  where
   g0 = pureMT (fromIntegral seed)
   go (x0,g0)  = let (eps, g) = randomDouble g0
                     x = r * x0 * (1-x0)
                     noise = skal * eps
                 in Just (x0 + noise, (x, g))

Re-benchmarking (recompiled without profiling) gives:

benchmarking vect1
time                 106.7 μs   (106.4 μs .. 107.0 μs)
                     1.000 R²   (1.000 R² .. 1.000 R²)
mean                 107.1 μs   (106.7 μs .. 107.7 μs)
std dev              1.415 μs   (842.3 ns .. 2.377 μs)

Note that that's 107 microseconds, so it's about 75 times faster.

That's where I'd stop, but if you decide to continue optimizing, make sure you refer frequently to the profiling and benchmarking results to make sure your changes are having an effect.

I highly recommend Googling for "profiling haskell programs" and for the "criterion" library and taking some time to learn how to use these tools.

For reference, the final program is:

import Criterion.Main
import qualified Data.Vector.Unboxed as U 
import System.Random.Mersenne.Pure64

new :: Double-> Double -> Int -> (Double, Int) -> U.Vector Double
new skal r n (x0, seed) = U.unfoldrN n go (x0, g0)
  where
   g0 = pureMT (fromIntegral seed)
   go (x0,g0)  = let (eps, g) = randomDouble g0
                     x = r * x0 * (1-x0)
                     noise = skal * eps
                 in Just (x0 + noise, (x, g))

vect1 :: (Double, Int) -> U.Vector Double
vect1 = new 0.5 1 10000

main = defaultMain [
  bench "vect1" $ nf vect1 (0,1)
  ]