通过最新发布的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。
我希望在这个向量中保存一些随机数。这样做的原因暂且不论，可能我只是想保存一些测试或演示数字，也可能正在创建一个需要随机数的仿真，不管怎么样，我需要一个随机向量。
通过调用函数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 ┗ ┛anil#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
评论