用于蒙特卡罗模拟的mt19937_64种子生成的最佳方法是什么?

31

我正在编写一个运行Monte Carlo模拟的程序;具体来说,我使用Metropolis算法。该程序需要生成可能达数十亿的“随机”数。我知道Mersenne twister在Monte Carlo模拟中非常流行,但我想确保以最佳方式初始化发生器。

目前,我正在使用以下方法计算32位种子:

mt19937_64 prng; //pseudo random number generator
unsigned long seed; //store seed so that every run can follow the same sequence
unsigned char seed_count; //to help keep seeds from repeating because of temporal proximity

unsigned long genSeed() {
    return (  static_cast<unsigned long>(time(NULL))      << 16 )
         | ( (static_cast<unsigned long>(clock()) & 0xFF) << 8  )
         | ( (static_cast<unsigned long>(seed_count++) & 0xFF) );
}

//...

seed = genSeed();
prng.seed(seed);

我有一种感觉,有更好的方法来确保不重复的新种子,而且我相当确定mt19937_64可以用超过32位的数值进行初始化。是否有人有什么建议吗?


3
为什么这很重要?为什么你需要确保每次模拟运行时使用不同的种子?为什么你需要费心去做这件事?虽然这样做不会给你带来“更好”的随机数。 - jalf
3
可以,但使用时间戳等简单方式进行种子生成可以确保这一点。为什么需要像NASA级别的复杂性来绝对“保证”……我甚至不知道你试图保证什么。听起来过于矫枉过正了。 - jalf
1
然而你在问题中没有提到这一点。你问的是“最好的种子”,这是一个毫无意义的问题。你显然想要知道的是“我如何选择种子,以便不同的线程(或进程?),即使它们同时启动,也有最小的选择相同种子的机会”。这是一个合理的问题。但它与选择“最好的种子”没有任何关系。你应该更新你的问题,询问你想要回答的实际问题。 - jalf
我没有得到我希望的答案,但希望@Mathhead200得到了一些有用的答案。如果你有任何偏好,想让谁获得奖励,请留下评论。 - dyp
@dyp 没关系。我得到了一些非常有用的信息和对问题的洞察力。此外,人们仍在发布帖子,我很有希望会有比我更懂的人提供一个优美且看似平台无关的解决方案,这是我迄今为止所没有想到的。 (到目前为止,你的解决方案仍然是最优雅和看似平台无关的解决方案。) - Mathhead200
显示剩余4条评论
6个回答

28

使用std::random_device来生成种子。如果您的实现支持,它将提供不确定性随机数。否则,可以使用其他随机数引擎。

std::mt19937_64 prng;
seed = std::random_device{}();
prng.seed(seed);

std::random_deviceoperator() 返回一个无符号整型unsigned int,所以如果你的平台的int为32位,并且你需要64位的种子(seed),你需要调用它两次。

std::mt19937_64 prng;
std::random_device device;
seed = (static_cast<uint64_t>(device()) << 32) | device();
prng.seed(seed);

另一个可用的选项是使用std::seed_seq来为PRNG提供种子。这允许PRNG调用seed_seq::generate,产生一个非偏置序列,其范围为[0 ≤ i < 232),并且具有足够大的输出范围以填充其整个状态。

std::mt19937_64 prng;
std::random_device device;
std::seed_seq seq{device(), device(), device(), device()};
prng.seed(seq);

我调用random_device 4次来创建一个包含4个元素的序列,作为seed_seq的初始序列。然而,我不确定在初始序列的长度或元素来源方面,哪种做法是最佳实践。


1
据我所知,这仅使用了64位种子。我不确定这对每个目的都足够。我仍然想知道如何优雅地使用random_device提供PRNG可以利用的尽可能多的种子。 - dyp
1
@dyp 你可以使用std::seed_seq,将其与多次调用random_device一起使用,然后将其传递给mt19937_64::seed。根据cppreference的说明,seed_seq将生成分布在[0≤i<2^32)范围内的结果,但我不知道这样做是否比位移更好,或者您需要调用多少次random_device来构建输入范围以使其被认为是好的(也就是说,输入应该有多少个元素)。 - Praetorian
1
哎?我不太明白你的意思。我发布的示例代码以624个32位数量为种子(请参见http://coliru.stacked-crooked.com/a/5e754fc72c89ed59),这是19968位,因此可能填满了19937位的整个状态。 - dyp
@Mathhead200 就像我在结尾处所说的那样,我不知道什么是最佳方法,我认为使用seed_seq是一种的方法,但是我不知道要使用多少初始化程序,或者是否应该使用random_device。如果您将初始化程序保存到seed_seq中,则每次使用该序列进行初始化时都会产生相同的输出,因为它是确定性的。我只使用过std::mt19937一次,那次我用一个random_device的单个调用对其进行了初始化。 - Praetorian
@Mathhead200 你拒绝听我的意见 :) 我对这个东西的了解还不够,不能推荐使用单次调用 random_device 还是 seed_seq 来进行种子生成。如果生成的随机数质量对您的应用程序至关重要,您应该联系该领域的专家。您应该能够查看 random_device 的 MSVC 实现。如果我没记错的话,它会忽略可选参数并调用一些 Windows 加密API函数。至于它在云计算平台上的可用性,我一无所知。我想,如果这些平台有C++标准库实现,它就会可用。 - Praetorian
显示剩余7条评论

7
让我们回顾一下(包括注释),我们想要生成不同的种子以在以下每种情况下获得独立的随机数序列:
  1. 程序稍后在同一台机器上重新启动,
  2. 同时在同一台机器上启动两个线程,
  3. 程序同时在两台不同的机器上启动。
1使用自纪元以来的时间解决,2使用全局原子计数器解决,3使用平台相关的ID解决(请参见如何以跨平台方式获取(几乎)唯一的系统标识符?
现在问题是什么是将它们组合以获得uint_fast64_tstd::mt19937_64的种子类型)的最佳方法?我假设我们在此处不知道每个参数的范围或它们太大,因此我们不能通过位移轻松地获得唯一种子。 std::seed_seq是简单的方法,但是其返回类型uint_least32_t不是我们的最佳选择。
一个好的64位哈希函数是更好的选择。STL在functional头文件下提供了std :: hash,一种可能性是将上述三个数字连接成一个字符串,然后传递给哈希函数。返回类型是一个size_t,在64位机器上很可能满足我们的要求。
碰撞是不太可能但当然有可能的,如果你想确保不会建立包含重复序列的统计数据,你只能存储种子并丢弃重复的运行。
std :: random_device也可以用于生成种子(仍然可能发生碰撞,难以确定更频繁还是更少),但由于实现依赖于库并且可能降至伪随机生成器,因此必须检查设备的熵并避免使用零熵设备进行此目的,因为您可能会打破上述点(尤其是点3)。不幸的是,只有当您将程序带到特定的机器上并使用安装的库进行测试时才能发现熵。

谢谢,总结得非常好。我将查看关于获取唯一系统标识符的SO链接。我的主要担忧是我不知道未来这个程序将在哪些机器上运行。(random_device也有同样的问题。)-- 我认为,通常情况下,Mersenne Twister的种子比64位更多...? - Mathhead200
1
我认为seed_seq是为了提供任意长度的种子而设计的。请查看seed_seq::generategenerate生成多个32位值,但mt19937_64不需要使用基于64位数据类型的状态;即使如此,您仍然可以使用32位值来填充它(通过分布/适配器)。 - dyp
1
MT19937的状态有19937位。mt19937_64的一个构造函数接受一个64位种子,并从该种子计算初始状态。另一个构造函数接受一个SeedSequence和它的generate方法来填充它的状态。一些实现从种子序列中提取19968位;足以填满整个状态。 - dyp
@dyp 非常清楚,我的担忧只是关于seed_seq仅使用32位数字。糟糕的实现可能会导致高碰撞率(请参见上文)。 - DarioP
1
seed_seq :: generate算法在标准中已经完全规定。我不知道该算法的确切数学属性。它是专门用于从(有偏)任意长度的数组生成MT种子的。而std::mersenne_twister_engine使用的算法来从生成的值初始化其状态也已经充分指定(这是我刚刚学到的)。就我所知,不会出现高碰撞率:这些工具是专门为此目的设计的。 - dyp
显示剩余6条评论

4
据你的评论,我认为你关心的是确保如果一个进程同时启动多个模拟,则它们将获得不同的种子。您当前方法唯一的显着问题是竞态条件:如果要同时启动多个模拟,则必须从单独的线程中执行。如果从单独的线程中执行,则需要以线程安全的方式更新seed_count,否则多个模拟可能会使用相同的seed_count。您可以将其简单地设置为std :: atomic 来解决。除此之外,它似乎比必须复杂。您通过使用两个单独的计时器所获得的好处是什么?您只需执行以下简单操作即可:1.在程序启动时(使用高分辨率计时器)获取当前系统时间,并将其存储。2.为每个模拟分配唯一ID(这可以只是初始化为0的整数,应该生成而没有任何竞争条件,如上所述),每次启动模拟时增加此值,类似于seed_count。3.当播种模拟时,只需使用最初生成的时间戳+唯一标识符。如果您这样做,过程中的每个模拟都将获得唯一的种子。

这就是我试图使用上面的代码实现的内容。然而,该代码是从一个名为“MSD”的类中提取出来的片段,这些实际上是成员变量。(也许将seed_count设置为静态变量是有必要的。)在多线程情况下,每个线程都有自己的MSD实例,因此,在这种情况下,seed_count并没有起到帮助作用。(再次将seed_count设置为静态变量应该可以解决这个问题。)另一个问题是这些任务可能会分配给多台计算机,这种情况下,上述解决方案不适用。 - Mathhead200

1
从评论中我了解到你想运行多个算法实例,一个线程对应一个实例。而且由于每个实例的种子几乎同时生成,你希望确保这些种子是不同的。如果这确实是你要解决的问题,那么你的genSeed函数并不能保证这一点。
在我看来,你需要的是可并行化的随机数生成器(RNG)。这意味着,你只需要一个RNG,你可以用你的genSeed实例化它,然后在非重叠的X个序列中分割通常在顺序环境中生成的随机数序列;其中X是线程数。有一个非常好的库可以提供这些类型的RNGs,在C++中遵循C++标准,并称为TRNG(http://numbercrunch.de/trng)。
这里有更多信息。您可以通过两种方式实现每个线程的非重叠序列。假设来自单个随机数生成器的随机数序列为r = {r(1),r(2),r(3),...},并且您只有两个线程。如果您事先知道每个线程需要多少个随机数,例如M,则可以将r序列的前M个数字分配给第一个线程,即{r(1),r(2),...,r(M)},并将第二个M分配给第二个线程,即{r(M + 1),r(M + 2),... r(2M)}。这种技术称为块分割,因为您将序列分成两个连续块。
第二种方法是将第一个线程的序列创建为{r(1),r(3),r(5),...},将第二个线程的序列创建为{r(2),r(4),r(6),...},其优点在于您不需要事先知道每个线程需要多少个随机数。这被称为跳跃式。
请注意,两种方法都保证每个线程的序列确实不重叠。我上面发布的链接有很多示例,并且库本身非常易于使用。希望我的帖子能够帮到您。

谢谢,但是我们也需要在多台计算机上运行此程序,这意味着随机数必须在另一个程序中生成并通过网络发送(而且这还假设计算机当前保持联网状态)。 - Mathhead200

1

您好...

有一些主要的代码启动了线程,并且有多个函数副本在这些线程中运行,每个副本都有自己的Marsenne Twister。我理解的对吗?如果是这样,为什么不在主代码中使用另一个随机生成器呢?它可以用时间戳进行种子化,并将连续的伪随机数发送给函数实例作为它们的种子。


我可以做到这一点,(虽然我更希望序列是独立的);但是,这并没有解决在多台机器上运行模拟的问题。而且,这似乎比只是调整种子生成要麻烦得多。(该程序相当大。) - Mathhead200
我想在之前回复时我误解了你的答案。它仍然没有解决多台机器(或程序同时运行的实例)的问题,但序列将是独立的,并且实现起来不难。但是,有一点变化,PRNG种子可能会发送重复的种子。 - Mathhead200

0

POSIX函数gettimeofday(2)可以提供微秒级别的时间。

POSIX线程函数gettid(2)返回当前线程的ID号码。

您应该能够将自纪元以来的秒数(已经在使用中),微秒数和线程ID组合起来,以获得在一台机器上始终唯一的种子。

如果您还需要它在多台机器上是唯一的,您可以考虑获取主机名、IP地址或MAC地址。

我猜32位可能已经足够了,因为有超过40亿个唯一的种子可用。除非您运行数十亿个进程,这似乎不太可能,否则不需要使用64位种子。


POSIX函数在Windows上可用吗?因为我不想让程序变得依赖于操作系统。另外,我们没有运行四十亿个进程,但是你不需要运行四十亿个进程来获得一个重复的32位种子(参见生日悖论)。 - Mathhead200
因此,使用生日问题概率的近似值来看,在使用32位种子的情况下,使用10k个进程时,您有约99%的机会没有发生碰撞。但是,只有在发生碰撞的进程使用相同参数时,碰撞才会产生影响。如果您具有不同的运行参数,则碰撞将不会产生影响。我只是试图帮你避免不必要的麻烦(除非它确实是必要的)。 - jsw
1
POSIX可以通过Cygwin在Windows上使用,但您是否已经有Windows用户了? 有一个仅适用于某个操作系统但能正常工作的程序总比纠结于操作系统依赖性而没有任何程序要好。 如果您发现这会妨碍采用率,那么您可以回到这个问题,并招募其他开发人员来帮助解决。 这类似于避免过早优化的想法。最终,您可能需要具有开关以选择适当方法的操作系统相关代码。此外,如果操作系统无关性很重要,请将其添加到原始问题(和悬赏)中! - jsw
该程序目前仅在Windows上运行,我没有能力“招募”其他人来帮助(我希望我有;我已经尝试过了)。该程序可能最终会在云或网格系统,甚至是超级计算机上运行,我不想添加任何可以避免的依赖项。感谢您告诉我Cygwin! - Mathhead200
你说得对,只有在存在相同的参数时(这种情况并不常见)才需要使用不同的种子。我只是想澄清一下你的帖子。另外,如果我有64位,我也可以使用它们。 - Mathhead200
都是正确的。可能仍然值得将操作系统信息添加到原始问题中。 - jsw

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接