随机引擎差异答案

【问题标题】：Random Engine Differences随机引擎差异
【发布时间】：2013-05-08 08:19:15
【问题描述】：

C++11 标准规定了许多不同的随机数生成引擎：linear_congruential_engine、mersenne_twister_engine、subtract_with_carry_engine 等等。显然，这与std::rand 的旧用法相比有了很大的变化。

显然，这些引擎（至少有一些）的主要好处之一是大大增加了周期长度（它内置在 std::mt19937 的名称中）。

但是，引擎之间的差异不太明显。不同引擎的优缺点是什么？什么时候应该使用一个而不是另一个？是否存在通常应该首选的合理默认值？

【问题讨论】：

+1 好问题。我会说std::mt19937 似乎是最常用的通用 RNG，但我绝对没有任何事实可以以此为基础......
@anthony-arnold 我同意。几乎所有的例子都使用了我见过的std::mt19937，而且我现在似乎在每次进行随机数生成时都会使用它，但除了习惯之外我没有这样做的具体原因。
原mersenne twister论文大概有详细解释为什么好：math.sci.hiroshima-u.ac.jp/~m-mat/MT/ARTICLES/mt.pdf
@DavidBrown 感谢您的链接。我同意带进位减法的信息似乎是最难获得的。

标签： c++ random c++11

【解决方案1】：

从下面的解释来看，线性引擎似乎更快但随机性更小，而 Marsenne Twister 具有更高的复杂性和随机性。进位减法随机数引擎是对线性引擎的改进，它肯定更加随机。在最后的参考文献中，指出 Mersenne Twister 比带进位减法的随机数引擎具有更高的复杂度

线性同余随机数引擎

一个产生无符号整数的伪随机数生成器引擎。

这是标准库中最简单的生成器引擎。它的状态是单个整数值，具有以下转换算法：

x = (ax+c) mod m

其中 x 是当前状态值，a 和 c 是它们各自的模板参数，如果大于 0，则 m 是其各自的模板参数，否则为 numerics_limits::max() 加 1。

它的生成算法是状态值的直接拷贝。

这使得它在处理和内存消耗方面非常有效，但生成的数字具有不同程度的序列相关性，具体取决于所使用的具体参数。

linear_congruential_engine 生成的随机数的周期为 m。 http://www.cplusplus.com/reference/random/linear_congruential_engine/

梅森捻线机随机数引擎

一个伪随机数生成器引擎，在闭区间 [0,2^w-1] 中生成无符号整数。

此引擎使用的算法经过优化，可计算范围内几乎均匀分布的大量数字（例如在蒙特卡罗实验中）。

引擎有一个由n个整数元素组成的内部状态序列，其中填充了构造或调用成员函数seed生成的伪随机序列。

内部状态序列成为 n 个元素的来源：当状态被推进时（例如，为了产生一个新的随机数），引擎通过在 a 上使用 xor mask a 扭曲当前值来改变状态序列由参数 r 确定的来自该值和 m 个元素之外的值的混合位（有关详细信息，请参见 operator()）。

产生的随机数是这些扭曲值的缓和版本。回火是一系列移位和异或操作，由应用于所选状态值的参数 u、d、s、b、t、c 和 l 定义（参见 operator()）。

mersenne_twister_engine 生成的随机数的周期相当于 mersenne 数 2^((n-1)*w)-1.http://www.cplusplus.com/reference/random/mersenne_twister_engine/

进位减法随机数引擎

一个产生无符号整数的伪随机数生成器引擎。

该引擎使用的算法是一个滞后斐波那契生成器，具有 r 个整数元素的状态序列，加上一个进位值。 http://www.cplusplus.com/reference/random/subtract_with_carry_engine/

如果使用加法或减法，滞后斐波那契生成器的最大周期为 (2k - 1)*^(2M-1)。 LFG 的初始化是一个非常复杂的问题。 LFG 的输出对初始条件非常敏感，统计缺陷可能会在最初出现，但也会周期性地出现在输出序列中，除非特别小心。 LFG 的另一个潜在问题是它们背后的数学理论不完整，因此有必要依赖统计测试而不是理论性能。 http://en.wikipedia.org/wiki/Lagged_Fibonacci_generator

最后：选择使用哪个引擎涉及许多权衡：线性同余引擎速度适中，并且对状态的存储需求非常小。即使在没有高级算术指令集的处理器上，滞后的斐波那契生成器也非常快，但代价是更大的状态存储和有时不太理想的频谱特性。 Mersenne twister 速度较慢，并且具有更高的状态存储要求，但具有正确的参数具有最长的非重复序列和最理想的光谱特性（对于给定的理想定义）。在http://en.cppreference.com/w/cpp/numeric/random

【讨论】：

这类回答了我的问题。 LCG 的周期长度比 mersenne twister 的周期长度要短得多，但是，它并没有真正揭示出subtract_with_carry_engine。
Boost 有一些关于循环长度和速度的有用信息，link

【解决方案2】：

我认为关键在于随机生成器具有不同的属性，这可以使它们更适合或不适合给定问题。

周期长度是属性之一。
随机数的质量也很重要。
生成器的性能也可能是个问题。

根据您的需要，您可以使用一台或另一台发电机。例如，如果您需要快速随机数但并不真正关心质量，LCG 可能是一个不错的选择。如果您想要更高质量的随机数，Mersenne Twister 可能是更好的选择。

为了帮助您做出选择，有一些标准测试和结果（我非常喜欢this paper 的第 29 页表格）。

编辑：来自论文，

LCG（论文中的LCG(***)）系列是最快的生成器，但质量最差。
Mersenne Twister (MT19937) 速度稍慢，但产生的随机数更好。
带有进位的减法（我认为是SWB(***)）要慢得多，但如果调整得当，可以产生更好的随机属性。

【讨论】：

那篇论文正是我想要的：性能数据。

【解决方案3】：

由于其他答案忘记了ranlux，这是最近将其移植到 OpenCL 的 AMD 开发人员的一个小笔记：

https://community.amd.com/thread/139236

RANLUX 也是极少数（我实际上知道的唯一一个）PRNG 之一，它有一个基本理论来解释它为什么会生成“随机”数字，以及为什么它们是好的。事实上，如果这个理论是正确的（我不知道有谁对此提出异议），最高奢华级别的 RANLUX 会产生完全去相关的数字，直到最后一点，只要我们保持良好状态，就没有长期相关性低于周期 (10^171)。大多数其他生成器对其质量只能说很少（如 Mersenne Twister、KISS 等）。它们必须依靠通过统计测试。

欧洲核子研究中心的物理学家是这个 PRNG 的粉丝。 'nuff 说。

【讨论】：

您好，链接已损坏。较新的 URL 是 community.amd.com/thread/139236。页面标题为“OpenCL中实现的RANLUX伪随机数生成器”。

【解决方案4】：

这些其他答案中的某些信息与我的发现相冲突。我已经使用 Visual Studio 2013 在 Windows 8.1 上运行了测试，并且始终发现mersenne_twister_engine 的质量比linear_congruential_engine 或subtract_with_carry_engine 的质量更高且速度明显更快。这让我相信，当考虑到其他答案中的信息时，引擎的具体实现会对性能产生重大影响。

我敢肯定，这对任何人来说都令人惊讶，但在其他答案中没有提到 mersenne_twister_engine 据说速度较慢。我没有针对其他平台和编译器的测试结果，但在我的配置下，mersenne_twister_engine 在考虑周期、质量和速度性能时显然是更好的选择。我没有分析内存使用情况，所以我不能谈论空间需求属性。

这是我用来测试的代码（为了便于移植，您应该只需要用适当的计时机制替换 windows.h QueryPerformanceXxx() API 调用）：

// compile with: cl.exe /EHsc
#include <random> 
#include <iostream>
#include <windows.h>

using namespace std;

void test_lc(const int a, const int b, const int s) {
    /*
    typedef linear_congruential_engine<unsigned int, 48271, 0, 2147483647> minstd_rand;
    */
    minstd_rand gen(1729);

    uniform_int_distribution<> distr(a, b);

    for (int i = 0; i < s; ++i) {
        distr(gen);
    }
}

void test_mt(const int a, const int b, const int s) {
    /*
    typedef mersenne_twister_engine<unsigned int, 32, 624, 397,
    31, 0x9908b0df,
    11, 0xffffffff,
    7, 0x9d2c5680,
    15, 0xefc60000,
    18, 1812433253> mt19937;
    */
    mt19937 gen(1729);

    uniform_int_distribution<> distr(a, b);

    for (int i = 0; i < s; ++i) {
        distr(gen);
    }
}

void test_swc(const int a, const int b, const int s) {
    /*
    typedef subtract_with_carry_engine<unsigned int, 24, 10, 24> ranlux24_base;
    */
    ranlux24_base gen(1729);

    uniform_int_distribution<> distr(a, b);

    for (int i = 0; i < s; ++i) {
        distr(gen);
    }
}

int main()
{
    int a_dist = 0;
    int b_dist = 1000;

    int samples = 100000000;

    cout << "Testing with " << samples << " samples." << endl;

    LARGE_INTEGER ElapsedTime;
    double        ElapsedSeconds = 0;

    LARGE_INTEGER Frequency;
    QueryPerformanceFrequency(&Frequency);
    double TickInterval = 1.0 / ((double) Frequency.QuadPart);

    LARGE_INTEGER StartingTime;
    LARGE_INTEGER EndingTime;
    QueryPerformanceCounter(&StartingTime);
    test_lc(a_dist, b_dist, samples);
    QueryPerformanceCounter(&EndingTime);
    ElapsedTime.QuadPart = EndingTime.QuadPart - StartingTime.QuadPart;
    ElapsedSeconds = ElapsedTime.QuadPart * TickInterval;
    cout << "linear_congruential_engine time: " << ElapsedSeconds << endl;

    QueryPerformanceCounter(&StartingTime);
    test_mt(a_dist, b_dist, samples);
    QueryPerformanceCounter(&EndingTime);
    ElapsedTime.QuadPart = EndingTime.QuadPart - StartingTime.QuadPart;
    ElapsedSeconds = ElapsedTime.QuadPart * TickInterval;
    cout << "   mersenne_twister_engine time: " << ElapsedSeconds << endl;

    QueryPerformanceCounter(&StartingTime);
    test_swc(a_dist, b_dist, samples);
    QueryPerformanceCounter(&EndingTime);
    ElapsedTime.QuadPart = EndingTime.QuadPart - StartingTime.QuadPart;
    ElapsedSeconds = ElapsedTime.QuadPart * TickInterval;
    cout << "subtract_with_carry_engine time: " << ElapsedSeconds << endl;
}

输出：

使用 100000000 个样本进行测试。 linear_congruential_engine 时间：10.0821 mersenne_twister_engine 时间：6.11615 减载引擎时间：9.26676

【讨论】：

【解决方案5】：

我刚从 Marnos 看到 this answer 并决定自己测试一下。我使用std::chono::high_resolution_clock 计时100000 采样100 次以产生平均值。我测量了std::chrono::nanoseconds 中的所有内容，最终得到了不同的结果：

std::minstd_rand 平均为28991658 纳秒

std::mt19937 平均为 29871710 纳秒

ranlux48_base 平均为 29281677 纳秒

这是在 Windows 7 机器上。编译器是 Mingw-Builds 4.8.1 64bit。这显然是使用了 C++11 标志并且没有优化标志。

当我打开 -O3 优化时，std::minstd_rand 和 ranlux48_base 实际上运行得比 high_precision_clock 的实现可以测量的要快；但是 std::mt19937 仍然需要 730045 纳秒，即 3/4 秒。

因此，正如他所说，它是特定于实现的，但至少在 GCC 中，平均时间似乎坚持接受答案中的描述。 Mersenne Twister 似乎从优化中受益最少，而其他两个实际上只是在考虑编译器优化后以令人难以置信的速度快速抛出随机数。

顺便说一句，我一直在我的噪声生成库中使用 Mersenne Twister 引擎（它不会预先计算梯度），所以我想我会改用其他引擎来真正看到速度的提升。就我而言，“真正的”随机性并不重要。

代码：

#include <iostream>
#include <chrono>
#include <random>

using namespace std;
using namespace std::chrono;

int main()
{
    minstd_rand linearCongruentialEngine;
    mt19937 mersenneTwister;
    ranlux48_base subtractWithCarry;
    uniform_real_distribution<float> distro;

    int numSamples = 100000;
    int repeats = 100;

    long long int avgL = 0;
    long long int avgM = 0;
    long long int avgS = 0;

    cout << "results:" << endl;

    for(int j = 0; j < repeats; ++j)
    {
        cout << "start of sequence: " << j << endl;

        auto start = high_resolution_clock::now();
        for(int i = 0; i < numSamples; ++i)
            distro(linearCongruentialEngine);
        auto stop = high_resolution_clock::now();
        auto L = duration_cast<nanoseconds>(stop-start).count();
        avgL += L;
        cout << "Linear Congruential:\t" << L << endl;

        start = high_resolution_clock::now();
        for(int i = 0; i < numSamples; ++i)
            distro(mersenneTwister);
        stop = high_resolution_clock::now();
        auto M = duration_cast<nanoseconds>(stop-start).count();
        avgM += M;
        cout << "Mersenne Twister:\t" << M << endl;

        start = high_resolution_clock::now();
        for(int i = 0; i < numSamples; ++i)
            distro(subtractWithCarry);
        stop = high_resolution_clock::now();
        auto S = duration_cast<nanoseconds>(stop-start).count();
        avgS += S;
        cout << "Subtract With Carry:\t" << S << endl;
    }

    cout << setprecision(10) << "\naverage:\nLinear Congruential: " << (long double)(avgL/repeats)
    << "\nMersenne Twister: " << (long double)(avgM/repeats)
    << "\nSubtract with Carry: " << (long double)(avgS/repeats) << endl;
}

【讨论】：

对于未来的读者：不要在基准测试中使用high_resolution_clock，使用steady_clock。很可能是相同的，但由于它不是标准化的，high_resolution_clock 可能是system_clock 的别名，具体取决于实现。 system_clock 可能不是单调的，因此不适合进行基准测试。
Gcc (11.1) at -O3 优化级别可以优化对 distro(linearCongruentialEngine); 的调用。该代码需要针对此类影响进行强化。 march=native 标志似乎对于 Mersenne-Twister 在计算效率方面具有竞争力至关重要。

【解决方案6】：

这确实是一种权衡。像Mersenne Twister 这样的 PRNG 更好，因为它具有极大的周期和其他良好的统计特性。

但是大周期的 PRNG 会占用更多的内存（用于维持内部状态），也需要更多的时间来生成随机数（由于复杂的转换和后处理）。

根据您的应用需求选择 PNRG。如有疑问，请使用 Mersenne Twister，这是许多工具中的默认设置。

【讨论】：

【解决方案7】：

一般来说，mersenne twister 是最好的（也是最快的）RNG，但它需要一些空间（大约 2.5 KB）。哪一个适合您的需要取决于您需要实例化生成器对象的次数。（如果您只需要实例化一次或几次，则可以使用 MT。如果您需要实例化数百万次，那么可能会更小一些。）

有些人报告说机器翻译比其他人慢。根据我的实验，这在很大程度上取决于您的编译器优化设置。最重要的是，-march=native 设置可能会产生巨大的差异，具体取决于您的主机架构。

我运行了一个小程序来测试不同发电机的速度，以及它们的大小，得到了这个：

std::mt19937 (2504 bytes): 1.4714 s
std::mt19937_64 (2504 bytes): 1.50923 s
std::ranlux24 (120 bytes): 16.4865 s
std::ranlux48 (120 bytes): 57.7741 s
std::minstd_rand (4 bytes): 1.04819 s
std::minstd_rand0 (4 bytes): 1.33398 s
std::knuth_b (1032 bytes): 1.42746 s

【讨论】：