在 GPU 上生成随机数

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

随机数是许多仿真或估计算法的关键部分。您可以在 MATLAB ® 中使用 rand、randi 和 randn 函数生成随机数。这些函数可以使用几种不同的数字生成算法之一。

检查 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 上使用。这些生成器名称相同，在给定相同初始状态的情况下，会产生完全相同的结果。当您想在 GPU 和 CPU 上生成相同的随机数集时，这很有用。有关详细信息，请参阅GPU 上的随机数流。

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

更改默认随机数生成器

检查是否可用 GPU。

gpu = gpuDevice;
disp(gpu.Name + " GPU detected and available.")

NVIDIA RTX A5000 GPU detected and available.

gpurng 函数可以存储和重置 GPU 的生成器状态。您还可以使用 gpurng 在不同的生成器之间切换。在更换生成器之前，请保存当前状态，以便在测试结束时能够恢复该状态。

oldState = gpurng;
gpurng(0,"Philox");
disp(gpurng)

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

生成均匀分布的随机数

rand 和 randi 函数在 GPU 上生成均匀分布的随机数。这两个函数的行为非常相似，此处仅对 rand 的性能进行了测试。请使用 gputimeit 函数来测量执行时间，以确保计时结果的准确性，因为 gputimeit 会多次调用该函数以抵消初始化开销，并正确处理同步问题。

为了比较不同生成器的性能，使用 rand 使用每个生成器在 GPU 上生成大量随机数。对于每个生成器，使用 $1 0^{7}$ 生成 100 次随机数，并使用 gputimeit 测量每次函数调用的耗时。

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

绘制结果。

figure
hold on
histogram(gputimesU(:,1),BinWidth=1e-5)
histogram(gputimesU(:,2),BinWidth=1e-5)
histogram(gputimesU(:,3),BinWidth=1e-5)

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

Figure contains an axes object. The axes object with title Generating Samples in U(0,1) Using NVIDIA RTX A5000, xlabel Time to generate 10 toThePowerOf 7 baseline blank random blank numbers blank (sec), ylabel Frequency contains 3 objects of type histogram. These objects represent Philox, Threefry, CombRecursive.

新型发电机 ThreeFry 和 Philox 性能相似。两者都比 CombRecursive 更快。

生成正态分布的随机数

许多仿真依赖于从正态分布中采样的扰动。与上面的测试类似，使用 randn 来比较这三个生成器在生成正态分布随机数时的性能。

reset(gpu)
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

绘制结果。

figure
hold on
histogram(gputimesN(:,1),BinWidth=5e-5)
histogram(gputimesN(:,2),BinWidth=5e-5)
histogram(gputimesN(:,3),BinWidth=5e-5)

legend(generators)
xlabel("Time to generate 10^7 random numbers (sec)")
ylabel("Frequency")
title("Generating Samples in N(0,1) Using " + gpu.Name)
hold off

Figure contains an axes object. The axes object with title Generating Samples in N(0,1) Using NVIDIA RTX A5000, xlabel Time to generate 10 toThePowerOf 7 baseline blank random blank numbers blank (sec), ylabel Frequency contains 3 objects of type histogram. These objects represent Philox, Threefry, CombRecursive.

结果再次表明，Threefry 和 Philox 生成器的性能相当，且都比 CombRecursive 明显更快。生成正态分布值所需的额外计算，降低了每个生成器生成值的速率。

恢复发电机至初始状态。

gpurng(oldState);

结论

虽然具体结果会因您的显卡而异，但每个生成器都有其独特的特点。

Threefry

快
最长周期（双精度下的 $2^{514}$ ）
基于广为人知且经过充分验证的 Threefish 算法

Philox

快
中间阶段（双精度中的 $2^{192}$ ）

CombRecursive

最慢的
最短周期（双精度下的 $2^{191}$ ）
在实际应用中拥有悠久的历史

参考资料

[1] L'Ecuyer, P., and R. Simard."TestU01:A C library for empirical testing of random number generators."《ACM 数学软件汇刊》Vol. 33, No. 4, 2007, article 22.

另请参阅

gpurng | parallel.gpu.RandStream

主题

GPU 上的随机数流