题 为蒙特卡洛模拟种子mt19937_64的最佳方法


我正在研究一个运行蒙特卡罗模拟的程序;具体来说,我正在使用Metropolis算法。该程序需要产生数十亿的“随机”数字。我知道梅森捻线机在蒙特卡罗模拟中非常受欢迎,但我想确保以最佳方式播种发电机。

目前我正在使用以下方法计算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位。有没有人有什么建议?


29
2018-06-20 19:02


起源


为什么这有关系?为什么你需要确保模拟的不同运行得到不同的种子?你为什么要这么做呢?它不会给你“更好”的随机数。 - jalf
因为我们可以使用相同的参数集运行模拟,在这种情况下,我们不一定期望完全相同的结果(如果我们使用相同的种子会发生这种情况。) - Mathhead200
当然,但是用时间戳这么简单的方式播种可以确保这一点。为什么你绝对需要NASA级别的复杂性 保证 那...我甚至都不知道你要保证的是什么。这听起来有点荒谬过度。 - jalf
@jalf <ctime>中time()的时间戳只有一秒的精度。但即使我使用毫秒精度(或其他),很可能许多模拟都是从相同的种子开始的。我同时运行了几个这样的模拟,(通常在不同的线程中以编程方式启动。) - Mathhead200
然而你在问题中没有说出这一点。你要求“最好的种子”,这是一个荒谬的问题。你显然是什么 通缉 回答是“我如何选择种子,以便不同的线程(或进程?),即使它们同时启动,也只有很少的机会选择相同的种子”。这是一个合理的问题。但它与选择“最好的种子”无关。您应该更新您的问题以询问您想要回答的实际问题。 - jalf


答案:


让我们回顾一下(评论也是如此),我们想要生成不同的种子,以便在以下每个事件中获得独立的随机数序列:

  1. 该程序稍后在同一台机器上重新启动,
  2. 两个线程同时在同一台机器上启动,
  3. 该程序同时在两台不同的机器上启动。

1使用时间解决,因为时期,2用全局原子计数器求解,3用平台相关id求解(见 如何以跨平台的方式获取(几乎)唯一的系统标识符?

现在重点是将它们组合起来得到一个最好的方法 uint_fast64_t (种子类型 std::mt19937_64)?我在这里假设我们不知道每个参数的范围或它们是否太大,所以我们不能只是通过位移来获得一个独特的种子以一种微不足道的方式。

一个 std::seed_seq 这将是一个简单的方法,但它的返回类型 uint_least32_t 不是我们最好的选择。

一个好的64位hasher是一个更好的选择。 STL提供 std::hash 在下面 functional 标题,一种可能性是将上面的三个数字连接成一个字符串,然后将它传递给hasher。返回类型是 size_t 在64台机器上很可能符合我们的要求。

冲突是不可能的,但当然可能,如果你想确保不会多次建立包含序列的统计数据,你只能存储种子并丢弃重复的运行。

一个 std::random_device 也可以用来生成种子(碰撞可能仍然发生,很难说或多或少经常发生),但是由于实现是依赖于库的,并且可能归结为伪随机生成器,因此必须检查熵的熵。设备并避免使用零熵设备为此目的,因为你可能会破坏上面的点(特别是第3点)。遗憾的是,只有将程序带到特定计算机并使用已安装的库进行测试时,才能发现熵。


4
2017-07-09 10:37



谢谢,完美的总结。我将在获取唯一系统标识符时查看该SO链接。我主要担心的是,我不知道将来该计划将运行哪些机器。 (同样的问题 random_device。) - 我想,一般来说,Mersenne Twister的种子数超过64位......? - Mathhead200
@ Mathhead200 STL实现采用64位种子,但原始C实现应采用任意长种子,请参阅 math.sci.hiroshima-u.ac.jp/~m-mat/MT/efaq.html (寻找 init_by_array这可以节省散列函数的使用,并可能降低碰撞的概率。然而,如果已经煮熟的STL实现,我不确定它是否值得。 - DarioP
我认为 seed_seq 已经提供任意长种子。看一眼 seed_seq::generate。 generate 生成多个32位值,但是 mt19937_64 不需要使用基于64位数据类型的状态;即使你仍然可以使用32位值来填充它(通过发行版/适配器)。 - dyp
该 州 一个MT19937有19937位。一个构造函数 mt19937_64 获取单个64位种子并从此种子计算初始状态。另一个构造函数接受SeedSequence及其 generate 填补其状态的方法。 一些实现 从种子序列中提取19968位;足以填满整个州。 - dyp
算法 seed_seq::generate 使用在标准中完全指定。我不知道该算法的确切数学属性。已经专门发明了从任意长度的(偏置)阵列中种植MT。那个算法 std::mersenne_twister_engine 用于从生成的值初始化其状态也是完全指定的(我刚学到的东西)。据我所知,没有高碰撞率:这些工具专门为此目的而设计。 - dyp


使用 std::random_device 生成种子。它将提供非确定性随机数,只要您的实现支持它。否则允许使用其他随机数引擎。

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

operator() 的 std::random_device 返回一个 unsigned int,所以如果你的平台有32位 ints,你想要一个64位种子,你需要调用它两次。

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。但是,就初始序列中元素的长度或来源而言,我不确定最佳做法是什么。


23
2018-06-20 19:08



AFAIK,它只使用64位种子。我不确定这是否足以满足各种目的。我仍然想知道如何优雅地使用 random_device 提供与PRNG可以使用的种子一样多的种子。 - dyp
@dyp有 std::seed_seq 您可以通过多次调用来提供 random_device,然后传递给 mt19937_64::seed。根据 cppreference  seed_seq 将生成分布在[0≤i<2 ^ 32)的结果,但我不知道这样做是否比位移更好,或者你需要多少次调用 random_device 用于构造输入范围 seed_seq (意味着输入应该考虑多少元素 好)。 - Praetorian
有关种子序列的有趣部分是MT可以根据需要请求尽可能多的种子,这可能不仅仅是64位。我不确定是否使用 seed_seq 与...结合 random_device 是有用的,这取决于MT的行为/要求(seed_seq::generate 消除偏见)。 - dyp
像这样的东西 (可能实现了更多的成员函数);但正如我所说,我不知道是否输出 uniform_int_distribution + random_device对MT有好处。 - dyp
咦?我不太明白你的意思。示例代码我发布了624个32位数量的种子(参见 coliru.stacked-crooked.com/a/5e754fc72c89ed59 ),这是19968位,因此可能填补整个19937位的状态。 - dyp


根据你的评论我可以看出,你感兴趣的似乎是确保如果一个过程在几个同一时间开始你的几个模拟,他们会得到不同的种子。

我用你当前的方法看到的唯一重要问题是竞争条件:如果要同时启动多个模拟,则必须从单独的线程完成。如果是从单独的线程完成的,则需要更新 seed_count 以线程安全的方式,或多个模拟可能最终相同 seed_count。你可以简单地做到 std::atomic<int> 解决这个问题。

除此之外,它似乎比它必须更复杂。使用两个独立的计时器可以获得什么?你可以做一些简单的事情:

  1. 在程序启动时,抓取当前系统时间(使用高分辨率计时器)一次,并存储它。
  2. 为每个模拟分配一个唯一的ID(这可能只是一个初始化为0的整数,(应该在没有任何竞争条件的情况下生成,如上所述),每次模拟开始时都会增加,有效地类似于 seed_count
  3. 播种模拟时,只需使用最初生成的时间戳+唯一ID。如果这样做,过程中的每个模拟都会确保一个独特的种子。

4
2017-07-08 14:33



这就是我试图用上面的代码做的事情。但是,代码是一个来自类的剪辑, MSD,这些实际上是成员变量。 (也许是在制作 seed_count 保证静态。)在多线程的情况下,每个线程都有自己的实例 MSD 因此,对于这种情况,seed_count没有帮助。 (再次制作 seed_count 静态应该解决这个问题。)另一个问题是这些任务可以在多台计算机之间分配,在这种情况下,谈论的解决方案不适用。 - Mathhead200


怎么样...

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


1
2017-07-09 08:53



我能做到这一点,(虽然我宁愿这些序列是独立的);但是,这并没有解决在多台机器上运行模拟的问题。此外,似乎还有更多的工作,只是搞乱种子生成。 (该计划相当大。) - Mathhead200
我想我在回答之前就误解了你的答案。它仍然没有解决多个机器(或程序并发运行的实例)的问题,而是序列 将 独立,并不难实施。然而,PRNG种子可能会发送重复种子。 - Mathhead200


POSIX功能 gettimeofday(2) 给出微秒精度的时间。

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

您应该能够组合自纪元(您已经使用)以来的秒数,以微秒为单位的时间以及获取在一台计算机上始终唯一的种子的线程ID。

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

我猜这32位可能就足够了,因为有超过40亿个独特的种子可用。除非你运行数十亿的进程,这似乎不太可能,否则你应该没有64位种子。


0
2017-07-10 20:24



POSIX功能是否在Windows上可用,因为我不希望程序依赖于操作系统?此外,我们没有运行四十亿个流程,没有;但是,您不需要运行40亿个进程来获取重复的32位种子(请参阅 生日问题。) - Mathhead200
因此,使用生日问题概率的近似值,当使用具有32位种子的10k进程时,看起来你有大约99%的机会没有碰撞。但是,只有碰撞过程使用相同的参数时才发生碰撞。如果您有不同的运行参数,那么碰撞将不会产生任何影响。我只是想为你省去不必要的麻烦(除非事实上是必要的)。 - jsw
Windows上通过Cygwin支持POSIX,但是你在Windows上有用户吗?最好是有一个程序可以工作,但是依赖于操作系统,而不是过分依赖操作系统依赖,根本没有程序。如果您发现自己阻碍了采用,那么您可以回过头来解决这个问题,甚至可以招募更多的开发人员来帮助您。这类似于避免过早优化的想法。最终,您可能需要使用与操作系统相关的代码和开关来选择适当的方法。此外,如果操作系统独立性很重要,请将其添加到原始问题(和赏金)! - jsw
该程序目前只在Windows上运行,我没有能力“招募”其他任何人来帮助(我希望我做过;我已经尝试过了。)该程序可能最终在云或网格系统上运行,或者甚至是某台超级计算机,我不想添加任何我可以避免的依赖项。谢谢你告诉我关于Cygwin的事! - Mathhead200
你是正确的,因为有不同的种子;如果有相同的参数(通常不常见),则非常重要。我只是在澄清你的帖子。另外,如果我有它,我不妨使用64位。 - Mathhead200


从我理解的评论中你想要运行算法的几个实例,每个线程一个实例。并且考虑到每个实例的种子将同时生成,您需要确保这些种子不同。如果这确实是您要解决的问题,那么您的genSeed功能不一定能保证。

在我看来,你需要的是一个并行的随机数发生器(RNG)。这意味着,你只需要  只与您实例化的RNG  种子(您可以使用genSeed生成),然后通常在顺序环境中生成的随机数序列被分成X个非重叠序列;其中X是线程数。有一个很好的库,用C ++提供这些类型的RNG,遵循RNG的C ++标准,称为TRNG(http://numbercrunch.de/trng)。

这里有更多信息。每种线程有两种方法可以实现非重叠序列。让我们假设a的随机数序列  RNG是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)}。这种技术称为blocksplitting,因为您将序列拆分为两个连续的块。

第二种方法是将第一个线程的序列创建为{r(1),r(3),r(5),...},将第二个线程创建为{r(2),r(4), r(6),...},其优点是您不需要事先知道每个线程需要多少个随机数。这叫做越级。

注意这两种方法 保证 每个线程的序列确实不重叠。我上面发布的链接有很多例子,库本身非常容易使用。我希望我的帖子有所帮助。


0
2017-07-11 12:28



谢谢,但我们也在多台机器上运行这个程序,这意味着随机数必须在一个单独的程序中生成并通过网络发送(假设计算机保持联网,就像现在一样。) - Mathhead200