Billions of Random Numbers in a Blink of an Eye

I promised to show you a lot of work done in a blink of an eye, but how long is that? According to Wikipedia, it does take an average of 100-150 milliseconds, or 100-400 milliseconds, depending on the source of the research.

I'll show you that it is not important which one of these we choose as the reference, since we can go even faster and assume that the blink of an eye is whatever it takes for a frame in an old movie to switch, which was 1/24 of a second in the good old days of analog theaters. That is 41.97 milliseconds.

First, I'll create a 1000-dimension vector in the main memory and measure the speed of generating random uniformly distributed entries.

(def x1K (fv 1000))
(quick-bench (rand-uniform! x1K))

Evaluation count : 1420254 in 6 samples of 236709 calls.
             Execution time mean : 422.105972 ns
    Execution time std-deviation : 1.347492 ns
   Execution time lower quantile : 420.687916 ns ( 2.5%)
   Execution time upper quantile : 423.978301 ns (97.5%)
                   Overhead used : 1.369142 ns

It takes 422 /nano/seconds per 1000 entries, which is 0.42 /nano/seconds per entry.

Let's create a larger vector, one of million entries, and try that.

(def x1M (fv 1000000))
(quick-bench (rand-uniform! x1M))

Evaluation count : 1890 in 6 samples of 315 calls.
             Execution time mean : 319.874178 µs
    Execution time std-deviation : 4.399502 µs
   Execution time lower quantile : 317.444041 µs ( 2.5%)
   Execution time upper quantile : 327.527853 µs (97.5%)
                   Overhead used : 1.369142 ns

Even faster per entry, 0.32 nanoseconds! This looks promising.

The next step is a billion entries. I'll create two vectors of 500 million entries each, which together takes a billion. The reason I'm doing this instead of one vector, is that each entry takes 4 bytes, so a billion entries requires 4GBs of bytes. Java (and Intel MKL) buffers are indexed with integers, and the largest integer is 2147483647.

So, generating 1 billion random numbers on the CPU takes…

(def x500M (fv 500000000))
(def y500M (fv 500000000))
(quick-bench
 (do (rand-uniform! x500M)
     (rand-uniform! y500M)))

Evaluation count : 6 in 6 samples of 1 calls.
             Execution time mean : 347.586774 ms
    Execution time std-deviation : 5.861722 ms
   Execution time lower quantile : 341.168665 ms ( 2.5%)
   Execution time upper quantile : 356.432065 ms (97.5%)
                   Overhead used : 1.369142 ns

347 milliseconds. Even that is a billion random numbers in a blink of an eye, if we choose the slowest eyes according to Wikipedia. But, this is Clojure and Neanderthal, so why not try GPU? Of course!

But before that, let's also mention that we have only demonstrated uniformly distributed random numbers. We often need normally distributed numbers, which are a bit more difficult to generate, and rand-normal expected to be slower. Let's see…

(quick-bench
 (do (rand-normal! x500M)
     (rand-normal! y500M)))

Evaluation count : 6 in 6 samples of 1 calls.
             Execution time mean : 1.870732 sec
    Execution time std-deviation : 4.685253 ms
   Execution time lower quantile : 1.865571 sec ( 2.5%)
   Execution time upper quantile : 1.876130 sec (97.5%)
                   Overhead used : 1.369142 ns

Almost 2 seconds. Although this is extremely fast (see the next section about motivation) this is slower than a blink of an eye!

I hope that I don't need to convince you that using GPU computing with Neanderthal is almost effortless (read about this in many dedicated articles on this blog).

Other than creating a GPU context, everything else is polymorphic and transparent, and uses the same functions as the CPU examples, so we're going to the benchmarks right away.

(with-default
  (cuda/with-default-engine
    (with-release [x500M (cuv 500000000)
                   y500M (cuv 500000000)]
      (time (do (rand-uniform! x500M)
                (rand-uniform! y500M)
                (synchronize!))))))

"Elapsed time: 39.319381 msecs"

That's right! The most restrictive goal has been met right away! Neanderthal generated 1 billion random numbers in 39 milliseconds, which is even faster than motion picture frames! Certainly faster than a blink of an eye.

Fantastic! But, what about normally distributed numebers? They are more demanding…

(with-default
  (cuda/with-default-engine
    (with-release [x500M (cuv 500000000)
                   y500M (cuv 500000000)]
      (time (do (rand-normal! x500M)
                (rand-normal! y500M)
                (synchronize!))))))

"Elapsed time: 34.108762 msecs"

Even better! 34 milliseconds for 1 billion random numbers from Normal distribution. That's something!

Now, I like supporting open source, and CUDA is a proprietary platform, etc. etc. Of course, Neanderthal's code is open-source, but the CUDA platform itself is owned and controlled by Nvidia.

Here's my AMD GPU with open-source ROCm drivers that AMD provides.

(ocl/with-default
  (opencl/with-default-engine
    (with-release [x500M (clv 500000000)
                   y500M (clv 500000000)]
      (time (do (rand-uniform! x500M)
                (rand-uniform! y500M)
                (finish!))))))

"Elapsed time: 37.707126 msecs"

Supercool! 37 milliseconds on a completely open GPU platform!

(ocl/with-default
  (opencl/with-default-engine
    (with-release [x500M (clv 500000000)
                   y500M (clv 500000000)]
      (time (do (rand-normal! x500M)
                (rand-normal! y500M)
                (finish!))))))

"Elapsed time: 37.836357 msecs"

Of course, by now you expect that rand-normal has the same performance…

I declare victory!

Here's how we use the built-in rand function.

(rand)

nil0.13648137992696296

Nice and simple! That part is good.

And, if we need a sequence of these, it's still nice.

(map rand (range 5))

nil(0.0 0.1955719758267589 1.0403752058141191 2.0357222475376053 1.7017948352691565)

There's a bug here. Note how the numbers progressively increase. This is because the first argument of rand controls the upper bound. If we want the whole sequence to have the same fixed bounds, or the default 1, we need to create a lambda function in-place.

(map (fn [_] (rand)) (range 5))

nil(0.027135352571733606 0.5705001188754536 0.3043893027311386 0.9307741452527688 0.42642378633952305)

What if we need to fill a more complex structure with random numebers? We need to program it ourselves, or to hope that the library at hand has a polymorphic support for mapping.

Fortunately, Neanderthal does support fmap!, one of Fluokitten's polymorphic alternatives to Clojure's map.

(fmap! (fn [_] (rand)) (fge 2 3))

nil#RealGEMatrix[float, mxn:2x3, layout:column, offset:0]
   ▥       ↓       ↓       ↓       ┓
   →       0.48    0.59    0.33
   →       0.13    0.49    0.32
   ┗                               ┛

Unfortunately, if we need these numbers to be normally distributed, we have to create a specialized function. I use the result discussed in the earlier article, Deep Learning from Scratch to GPU - 13 - Initializing Weights.

(defn rand-gaussian ^double [^double _]
  (double (* (sqrt (* -2.0 (log (double (rand)))))
             (sin (* 2.0 pi (double (rand)))))))

(fmap! rand-gaussian (fge 2 3))

nil#RealGEMatrix[float, mxn:2x3, layout:column, offset:0]
   ▥       ↓       ↓       ↓       ┓
   →       1.33   -0.66    1.04
   →       1.23   -0.97   -0.36
   ┗                               ┛

While this is good enough for lots of purposes, one important technical obstacle is that rand works only on the CPU! If we need to initialize random vectors on the GPU, rand is of no help.

As I demonstrated in Deep Learning from Scratch to GPU - 13 - Initializing Weights, generating random numbers on the CPU, and transferring them to the GPU takes much more time that computing the main algorithm that processes these random numbers on the GPU. That is ugly and not good enough.

This is actually and issue that I can't demonstrate easily. If you are not already aware of this, it probably didn't matter much in your applications, since you didn't need that many random numbers, so you couldn't discover that the numbers that built in random functions generate are not that well distributed.

Now, of course, the "random" numbers that the computer generates are not really random. They are "pseudorandom", meaning that there is a deterministic algorithm that generates them. This is why when you provide the same seed, you get the same sequence of "random" numbers.

Producing real random numbers on a computer is difficult, and often requires specialized hardware, or some other source of randomness from the real world. This is very important for cryptography, and is very difficult and slow. I am not talking about that level of randomness here.

Fortunately, for almost everything outside cryptography, various pseudorandom numbers are good enough. However, while the built-in RNG is good enough if you need to create a list of 42 random numbers for testing purposes, it is not good enough if you need to explore a distribution with MCMC.

I won't bother you again other than stating that Neanderthal uses Philox and/or ARS5 RNG which is much, much, better than the stuff you get from the built-in rand.

This is easy to demonstrate.

(quick-bench (rand))

Evaluation count : 35329860 in 6 samples of 5888310 calls.
             Execution time mean : 15.789613 ns
    Execution time std-deviation : 0.114434 ns
   Execution time lower quantile : 15.674532 ns ( 2.5%)
   Execution time upper quantile : 15.964610 ns (97.5%)
                   Overhead used : 1.364400 ns

While 15 ns does not seem too slow, compare it to 0.32 ns that we get on the CPU, and 0.034 ns that we got on the GPU.

What about normal distribution?

(quick-bench (rand-gaussian Double/NaN))

Evaluation count : 8797110 in 6 samples of 1466185 calls.
             Execution time mean : 67.180339 ns
    Execution time std-deviation : 0.511289 ns
   Execution time lower quantile : 66.824600 ns ( 2.5%)
   Execution time upper quantile : 67.927071 ns (97.5%)
                   Overhead used : 1.364400 ns

Billion numbers should take 15 or 67 seconds, but we should also consider that these numbers have to be written to the memory, so let's benchmark it.

(time
 (do (fmap! rand-gaussian x500M)
     (fmap! rand-gaussian y500M)))

"Elapsed time: 69576.897458 msecs"

Now, I admit that rand is fast enough for most stuff, and that the main reason why rand is not good enough is that the numbers it generates are poor, and that it does not work on the GPU.

The fact that Neanderthal can do this RNG 50 times faster on the CPU, and 200 times faster on the GPU is just an icing on the cake.