如何为最大周期播种并行线性同余伪随机数生成器?

计算科学 并行计算 随机数生成
2021-12-10 09:22:48

通常,当我在 C 中播种顺序随机数生成器时,我使用调用

srand(time(NULL))

然后使用

rand() mod N

获得一个介于 0 和 N-1 之间的随机整数。但是,当我并行执行此操作时,对 time(NULL) 的调用彼此非常接近,以至于它们最终成为完全相同的数字。

我尝试使用线性同余随机数生成器:

xn+1:=axn+c(modm)用于某些整数a,c,m

我知道为一些大整数k选择m=2^k会产生快速的结果,因为模数运算符可以通过截断数字来计算。然而,我发现很难建立产生具有大周期并行随机序列的种子。我知道一个周期长度是最大的,如果 m=2kk

  1. c 和 m 互为素数
  2. a-1 能被 m 的所有素因子整除
  3. 如果 m 是 4 的倍数,则 a-1 也必须是 4 的倍数。

(来源:维基百科

但是如何确保所有随机数流都具有这个最大属性呢?就 MPI 而言,如何使用线性同余方法合并ranksize以产生最大周期?使用滞后 Fibonacci 或 Mersenne Twister 生成更长的并行随机流会更容易吗?

4个回答

诀窍是交错每个进程的 LCG 流:对于个进程,我们修改 LCGp

xn+1:=axn+c(modm),

成为

xn+p:=Apxn+Cp(modm),

其中有效地向前推进步。我们可以通过扩展原始 LCG 步骤来快速推导出它们:ApCpp

xn+2=a(axn+c)+c(modm)=a2xn+(a+1)c(modm),

并且模式是是从数字开始,依次乘以并加次,然后乘以,所有的结果。Apapmodm,Cp0a1 pcmodm

最后一步是确保每个进程的步长 LCG 不重叠:只需使用的进程,并且一个具有独立周期的并行 LCG已准备就绪,其中是原始周期,是进程数。如果每个进程的修改后的 LCG 均等使用,则整个周期将并行恢复。prxrN/pNp

我大约六个月前实现了这个(也许是天真的),你可以在这里找到代码。

Katzgrabber 有一篇非常不错的教程概述论文,科学计算中的随机数:简介,我指出那些想成为 PRNG 用户进行科学计算的人。线性同余生成器很快,但这就是它们的全部功能。它们的周期很短,很容易出错;即使您满足 a、c 和 m 之间的通常要求,a、c 和 m 看起来完全合理的组合最终也会产生可怕的相关输出。

更糟糕的是,在 m 是 2 的幂的一种常见情况下(因此 mod 操作很快),低位的周期比整个序列的周期短得多,所以如果你正在做 rand() % N,你的时间比你预期的还要短。

作为一般经验法则,滞后斐波那契、MT 和 WELL 生成器具有更好的属性,并且它们仍然相当快。

在并行播种方面,Jack Poulson 的方法很好,因为它给出了一个定义良好的数字序列,在处理器之间平均分配。如果这无关紧要,您可以做任何合理的事情来播种不同的 PRNG;同一篇论文提出了很多人独立提出的一些建议,将 PID 或 MPI 任务编号与时间进行散列。建议的特定公式是

long seedgen(void)  {
    long s, seed, pid;

    pid = getpid();
    s = time ( &seconds );

    seed = abs(((s*181)*((pid-83)*359))%104729); 
    return seed;
}

我对那个具体的实现没有什么特别的看法,但一般的做法肯定是合理的。

在相当数量的线程上传播典型的顺序 RNG 的一个简单想法是让单个线程尽可能快地推进种子,并且仅将每千分之一左右的种子发送到内存中。然后让您的每个其他线程拾取这些间隔的参考种子中的一个并处理该块中的 1000 个值,即再次重新生成块中的 1000 个种子,生成它们的伪随机绘制,然后执行您的任务涉及的任何其他处理。

这很有效,因为对于计算量不大的 RNG(LCG 肯定是其中之一,但许多其他应该属于这一类),真正的瓶颈是将种子发送到内存(也可能是后续处理)。如果您运行 LCG 时没有将任何内容发送到内存,那么整个事情应该留在 CPU 寄存器中并且速度非常快。即使对于更复杂的 RNG,您也应该留在 L1 缓存中并且非常快。

我已经在 LCG 中使用了这种非常简单的方法,出于遗留原因,我们必须保留它。在典型的多核工作站中,我们基本上可以线性加速到 4-8 个线程。但是现在我将尝试 Jack Poulson 的答案中的方法,并希望变得更快:)。

OTOH,我相信这个简单的技巧应该适用于其他固有的顺序 RNG。

这个答案来晚了,但你应该看看SPRNG它专为并行可扩展性而设计,并支持几种类型的 PRNG。