在 GPU 上生成随机数

此示例展示如何在 GPU 支持的不同随机数生成器之间切换。

随机数是许多仿真或估计算法的关键部分。通常，这些数字是使用函数 rand、randi 和 randn 生成的。Parallel Computing Toolbox™ 提供了三个相应的函数用于直接在 GPU 上生成随机数：rand、randi 和 randn。这些函数可以使用几种不同的数字生成算法之一。

d = gpuDevice;
fprintf("This example is run on a " + d.Name + " GPU.")

This example is run on a GeForce GTX 1080 GPU.

探索 GPU 随机数生成器

函数parallel.gpu.RandStream.list提供了可用生成器的简短描述。

parallel.gpu.RandStream.list

 
The following random number generator algorithms are available:
 
MRG32K3A:         Combined multiple recursive generator (supports parallel streams)
Philox4x32_10:    Philox 4x32 generator with 10 rounds (supports parallel streams)
Threefry4x64_20:  Threefry 4x64 generator with 20 rounds (supports parallel streams)

每个生成器在设计时都考虑了并行使用，可提供多个独立的随机数流。但是，它们各自都有一些优点和缺点：

CombRecursive（也称为 MRG32k3a）：该生成器于 1999 年推出，并已得到广泛的测试和使用。
Philox（也称为 Philox4x32_10）：2011 年推出的新生成器，专为在 GPU 等高度并行系统中实现高性能而设计。
Threefry（也称为 Threefry4x64_20）：2011 年推出的基于现有加密 ThreeFish 算法的新生成器，该算法经过了广泛的测试和使用。该生成器旨在在 GPU 等高度并行系统中提供良好的性能。这是 GPU 计算的默认生成器。

GPU 上可用的三个生成器也可用于 MATLAB® 中的 CPU。MATLAB 生成器具有相同的名称，并且在相同的初始状态下产生相同的结果。当您想在 GPU 和 CPU 上生成相同的随机数集时，这很有用。有关详细信息，请参阅GPU 上的随机数流。

所有这些生成器都通过了标准 TestU01 测试套件 [1]。

更改默认随机数生成器

函数gpurng可以存储和重置 GPU 的生成器状态。您还可以使用 gpurng 在提供的不同生成器之间切换。在更换生成器之前，存储现有状态，以便在这些测试结束时恢复。

oldState = gpurng;

gpurng(0, "Philox4x32-10");
disp(gpurng)

     Type: 'philox'
     Seed: 0
    State: [7×1 uint32]

生成均匀分布的随机数

使用 rand 或 randi 在 GPU 上生成均匀分布的随机数。在性能方面，这两个函数的行为非常相似，这里只测量 rand。gputimeit用于测量性能，以确保准确的计时结果，自动多次调用该函数并正确处理同步和其他计时问题。

为了比较不同生成器的性能，使用 rand 使用每个生成器在 GPU 上生成大量随机数。在下面的代码中，rand 生成 $1 0^{7}$ 随机数，每个生成器被调用 100 次。每次运行都使用 gputimeit 来计时。生成大量随机数样本可能需要几分钟。结果表明了 GPU 上可用的三个随机数生成器之间的性能比较。

generators = ["Philox","Threefry","CombRecursive"];
gputimesU = nan(100,3);
for g=1:numel(generators)
    % Set the generator
    gpurng(0, generators{g});
    % Perform calculation 100 times, timing the generator 
    for rep=1:100
        gputimesU(rep,g) = gputimeit(@() rand(10000,1000,"gpuArray"));
    end
end

% Plot the results
figure
hold on
histogram(gputimesU(:,1),"BinWidth",1e-4);
histogram(gputimesU(:,2),"BinWidth",1e-4);
histogram(gputimesU(:,3),"BinWidth",1e-4)

legend(generators)
xlabel("Time to generate 10^7 random numbers (sec)")
ylabel("Frequency")
title("Generating samples in U(0,1) using " + d.Name)
hold off

较新的生成器 Threefry 和 Philox 具有类似的性能。两者都比 CombRecursive 更快。

生成正态分布的随机数

许多仿真依赖于从正态分布中采样的扰动。与均匀测试类似，使用randn比较三个生成器在生成正态分布的随机数时的性能。生成大量随机数样本可能需要几分钟。

generators = ["Philox","Threefry","CombRecursive"];
gputimesN = nan(100,3);
for g=1:numel(generators)
    % Set the generator
    gpurng(0, generators{g});
    % Perform calculation 100 times, timing the generator 
    for rep=1:100
        gputimesN(rep,g) = gputimeit(@() randn(10000,1000,"gpuArray"));
    end
end

% Plot the results
figure
hold on
histogram(gputimesN(:,1),"BinWidth",1e-4);
histogram(gputimesN(:,2),"BinWidth",1e-4)
histogram(gputimesN(:,3),"BinWidth",1e-4)
legend(generators)
xlabel("Time to generate 10^7 random numbers (sec)")
ylabel("Frequency")
title("Generating samples in N(0,1) using " + d.Name)
hold off

再次，结果显示 Threefry 和 Philox 生成器的性能相似，并且都明显比 CombRecursive 快。产生正态分布值所需的额外工作会降低每个生成器产生值的速率。

在完成之前，恢复原始生成器状态。

gpurng(oldState);

结论

在此示例中，比较了三个 GPU 随机数生成器。具体结果取决于您的 GPU 和计算平台。每个生成器都有一些优点（+）并且也有一些注意事项（-）。

Threefry

(+) 快
(+) 基于著名且经过充分测试的 Threefish 算法
（-）在实际使用中相对较新

Philox

(+) 快
（-）在实际使用中相对较新

CombRecursive

(+) 在实际使用中有着长期的记录
(-) 最慢

参考

[1] L'Ecuyer, P., and R. Simard."TestU01:A C library for empirical testing of random number generators."ACM Transactions on Mathematical Software.Vol. 33, No. 4, 2007, article 22.

另请参阅

gpurng | parallel.gpu.RandStream