通过最新发布的Neanderthal，可以在CPU和GPU上直接生成随机向量和矩阵，基本可以在眨眼间生成数十亿高质量随机数。本文，作者对该功能的性能及实现动机进行了深入说明

如果想在机器上运行本文代码，需要创建一个包含Clojure依赖的项目。如果只是想试一下，可以直接复制Neanderthal项目中的Hello World例程。只需要稍微改动，就可以在Java项目中运行这些代码。

随机代码的Hello World

与往常一样，本文示例由动态代码自动生成。因此，为了能在本地环境正常运行，我需要展示一些必要的导入模块。我希望你能够在自己的开发环境中验证这些代码。如果只是想单纯地读一下，可以忽略下面的代码。

(require '[uncomplicate.commons.core :refer [with-release let-release]]
         '[uncomplicate.fluokitten.core :refer [fmap!]]
         '[uncomplicate.clojurecuda.core :refer [with-default synchronize!]]
         '[uncomplicate.clojurecl.core :refer [finish!] :as ocl]
         '[uncomplicate.neanderthal
           [native :refer [fv fge native-float]]
           [core :refer [submatrix native]]
           [math :refer [sqrt log sin pi sqr]]
           [cuda :refer [cuv cuge] :as cuda]
           [opencl :refer [clv clge] :as opencl]
           [random :refer [rand-normal! rand-uniform! rng-state]]]
         '[criterium.core :refer [quick-bench]])

让我们先从容易的代码开始入手，在内存中声明一个向量x。

(def x (fv 5))

我希望在这个向量中保存一些随机数。这样做的原因暂且不论，可能我只是想保存一些测试或演示数字，也可能正在创建一个需要随机数的仿真，不管怎么样，我需要一个随机向量。

通过调用函数rand-uniform可以很容易获得：

(rand-uniform! x)
nil#RealBlockVector[float, n:5, offset: 0, stride:1]
[   0.88    0.62    0.22    0.72    0.16 ]

那这个函数可以用到更复杂的结构上么，例如GE矩阵？答案是肯定的。

(def a (fge 3 2))

(rand-normal! 33 2.5 a)

nil#RealGEMatrix[float, mxn:4x3, layout:column, offset:0]
   ▥       ↓       ↓       ↓       ┓
   →      33.90   36.38   34.29
   →      31.09   32.56   33.44
   →      36.96   37.39   32.83
   →      33.46   36.02   34.23
   ┗                               ┛

那rand-normal又是什么呢？很多时候你可能需要均匀分布的数字，但有时候你需要一些正态的数字。rand-uniform生成均匀分布的随机数，默认是U(0,1)，而rand-normal则生成正态分布的随机数，默认是N（0,1）。

生成十亿随机数需要多久？

我在这里陷入了一个进退两难的情况，是应该先介绍实现这个功能的动机，还是先炫耀下该功能的性能？我决定先展示一下其性能，毕竟在你点击本文链接的时候就已经表明你对随机数生成很感兴趣，相信你也已经知道这有多麻烦。

我前面已经保证在一眨眼间可以做很多事情，那一眨眼究竟是多长呢？根据维基百科，眨一下眼需要平均约100～150毫秒，也有其它研究表明为100～140毫秒。

其实选择何种结果并没有很大不同，因为我们可以假设这个时间更短。我们可以把老电影中帧切换作为一次眨眼，这要比上面的任意一种结果都要快。

在CPU上：一颗五年前的Intel i7-4790

首先，让我们在内存里创建一个1000维的向量，然后测量生成这些随机均匀分布项的速度。

(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

这个过程显示每创建一千个单元使用422纳秒，也就是0.42 纳秒每单元。

接下来让我们创建一个百万级别的向量。

(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

生成速度更快，只需要32纳秒！

下一步是十亿个单元。我会分别创建两个五亿单元，之所以使用两个，主要原因是因为只一个的话会使内存溢出。每个单元为4个字节，十亿就是4GB。Java（以及Intel的MKL）缓冲是用整数类型来做索引的，而最大的整数是2147483647（4G）。

所以，生成10亿个随机数需要的时间是：

(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毫秒。如果我们选择维基百科中最慢的眨眼时间，我们勉强可以称在一眨眼间生成了十亿个随机数。但是，我们使用了Clojure和Neanderthal，那为什么不用一下GPU呢？

在开始以前，我需要先说明一下，我们这里只是展示了均匀分布的随机数。更多的时候我们需要正态分布的随机数，但这更困难，而且rand-normal也会更慢一点。

(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

差不多2秒钟，尽管这已经很快了，但这仍然比眨一下眼睛要慢。

在GPU上：Nvidia GTX 1080Ti

我相信我已经在很多地方强调过，在GPU计算中使用Neanderthal是非常容易的，所以也就不在这里重复了。

除了创建GPU上下文，Neanderthal中所有的东西都是多态且透明的，并且使用同上面CPU例子一样的函数，让我们马上来测试一下吧。

(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"

怎么样，这个十分苛刻的目标在这里被瞬间实现。Neanderthal用39毫秒生成了10亿个随机数，这要比动态图片帧的切换都要快，当然也比眨一下眼更快。

那正态分布随机数会怎么样呢？毕竟正态分布是更想要的。

(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"

结果只需要34毫秒。

在GPU上：AMD Vega 64

现在，我更喜欢支持开源的东西，Neanderthal的代码就是开源的，CUDA虽然开源但仍然被Nvidia控制，更像是一个私有的开发平台。

下面是AMD GPU以及他们开源的ROCm驱动。

(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"

在一个完全开源的平台上只需要37毫秒！

(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"

当然，现在你也希望rand-normal能够具有一样的性能。

为什么需要这个功能？

几乎每个大型开发框架或编程语言都提供了随机生成函数。这不就可以了嘛？

默认的随机数函数一般具有三个问题：

通过它们生成的随机数质量很差。
生成过程很慢。
只能生成均匀分布的随机数并且质量很差。

第三个问题比较容易解决，可以通过实现一个函数将rand的输出转变为正态分布，但其他两个就很麻烦了。

内置RNG（随机数生成器）

下面让我们來看一下如何使用内置rand函数。

(rand)

nil0.13648137992696296

如果我们需要一个随机数序列，生成过程仍然很简单。

(map rand (range 5))

nil(0.0 0.1955719758267589 1.0403752058141191 2.0357222475376053 1.7017948352691565)

这里有个bug。注意这些数是逐步增长的，这是因为rand的第一个参数决定了它的上限。如果需要整个序列具有相同固定的界限，或者默认值1，我们需要另外创建一个lambda函数。

(map (fn [_] (rand)) (range 5))
nil(0.027135352571733606 0.5705001188754536 0.3043893027311386 0.9307741452527688 0.42642378633952305)

如果需要用随机数填充一个更复杂的结构会怎么样呢？我们只能通过自己编程实现，或者希望手边的类库具有对映射的多态支持。

幸运的是，Neanderthal能够支持fmap!，这是Fluokitten提供的其中一种Clojure映射的多态替代方案。

(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
   ┗                               ┛

而遗憾的是，如果需要这些随机数正态分布，我们就需要特别创建一个函数。我在这里使用之前文章中所讨论的结果。

(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
   ┗                               ┛

尽管结果已经可以满足大多数情况，但一个非常重要的技术障碍就是rand只能在CPU上运行！如果我们需要在GPU上初始化一个随机向量，这时候rand就无能为力了。

之前提到过：通过CPU生成随机数然后转移到GPU所需要的时间要远比直接在GPU上生成这些随机数的时间要长，而且这种设计也很丑陋。

内置RNG质量很差

这个问题其实很难说清楚。如果不是事前已经了解的话，一般很难注意到它对程序产生的影响，毕竟需要的数量可能并不多，因此很难注意到这些数其实并没有很好的分布。

当然，这些产生的随机数也并不是真正的随机数，它们又被称作“伪随机数”，因为它们通过一个确定的算法生成。这就是为什么当你提供一些相同的随机种子，你会得到一些相同序列的随机数。

通过计算机生成真实的随机数是很困难的，很多时候需要特殊的硬件，或者是真实世界中的一些随机源。这对于密码学领域尤其重要，随机数生成很困难也很慢。当然本文讨论的并不是这种情况。

幸运的是，只要不是跟密码相关，伪随机数就已经可以了。例如，你只是需要生成42个随机数用于测试，rand函数就可以满足，但是如果你需要一些随机数用于探索MCMC（马尔可夫链蒙特卡洛），rand函数就不可以了。

Neanderthal使用了Philox和ARS5 RNG，这要远远好于那些内建的rand函数。

内置RNG很慢

这个就很容易解释了。

(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

对比CPU上的0.32纳秒以及GPU上的0.034纳秒，15纳秒看起来并不太慢。

正态分布情况下会怎么样呢？

(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

十亿个随机数的话应该需要15或者67秒，但是同时需要考虑保存到内存所需要的时间，让我们正式测试一下。

(time
 (do (fmap! rand-gaussian x500M)
     (fmap! rand-gaussian y500M)))
"Elapsed time: 69576.897458 msecs"

我承认rand函数对大多数情况来说已经足够快了，唯一不足的是它产生的随机数的质量并不太好，而且它也不能在GPU上运行。

事实上，Neanderthal在CPU上可以比rand函数快50倍，在GPU上可以很轻松地达到200倍。

其他免费功能

除了性能上的优势，Neanderthal函数还可以很容易地适应不同的数据结构。Neanderthal不仅能支持矩阵，也可以支持更棘手的亚矩阵结构。

例如下面这个rand-uniform函数，它可以让主矩阵体同亚矩阵的数据完全分开。

(def a (fge 4 3))
(def sub-a (submatrix a 0 1 4 2))
(rand-uniform! 0 2 sub-a)
nil#RealGEMatrix[float, mxn:4x2, layout:column, offset:4]
   ▥       ↓       ↓       ┓
   →       1.78    0.09
   →       0.32    0.32
   →       1.86    1.94
   →       1.43    0.95
   ┗                       ┛
a
nil#RealGEMatrix[float, mxn:4x3, layout:column, offset:0]
   ▥       ↓       ↓       ↓       ┓
   →      33.90    1.78    0.09
   →      31.09    0.32    0.32
   →      36.96    1.86    1.94
   →      33.46    1.43    0.95
   ┗                               ┛

当然，该函数也可以在GPU上运行。

(with-default
  (cuda/with-default-engine
    (with-release [a (cuge 5 4)
                   sub-a (submatrix a 1 1 3 2)]
      (rand-normal! 55 3 sub-a)
      (native a))))
: nil#RealGEMatrix[float, mxn:5x4, layout:column, offset:0]
:    ▥       ↓       ↓       ↓       ↓       ┓
:    →       0.00    0.00    0.00    0.00
:    →       0.00   57.59   63.07    0.00
:    →       0.00   53.33   55.46    0.00
:    →       0.00   52.48   58.21    0.00
:    →       0.00    0.00    0.00    0.00
:    ┗                                       ┛

如果你需要代码能够被复现，并且在每次运行测试时都能得到相同的随机数，那么你可以提供一个特殊的随机种子。

(rand-uniform! (rng-state native-float 11) 2 3 (fv 3))
nil#RealBlockVector[float, n:3, offset: 0, stride:1]
[   2.73    2.40    2.89 ]
(rand-uniform! (rng-state native-float 11) 2 3 (fv 3))
nil#RealBlockVector[float, n:3, offset: 0, stride:1]
[   2.73    2.40    2.89 ]

英文原文：https://dragan.rocks/articles/19/Billion-random-numbers-blink-eye-Clojure

创作场景

如何在眨眼间生成数十亿随机数?