具有任意边界的快速、无偏整数伪随机生成器答案

【问题标题】：Fast, unbiased, integer pseudo random generator with arbitrary bounds具有任意边界的快速、无偏整数伪随机生成器
【发布时间】：2014-10-24 09:30:17
【问题描述】：

对于蒙特卡洛积分过程，我需要从中提取很多个随机样本具有 N 个桶的直方图，其中 N 是任意的（即不是 2 的幂），但在计算过程中完全没有变化。

很多，我的意思是大约 10^10、100 亿，几乎任何面对大量的样本）。

我有一个非常快速的统一伪随机数生成器可供我使用通常产生无符号的 64 位整数（讨论中的所有整数以下未签名）。

提取样本的幼稚方法：histogram[ prng() % histogram.size() ]

天真的方法是非常慢：模运算使用整数除法（IDIV）这是非常昂贵的编译器，不知道histogram.size() 的价值在编译时，不能达到它通常的魔力（即http://www.azillionmonkeys.com/qed/adiv.html）

事实上，我的大部分计算时间都花在了提取那个该死的模数上。

稍微不那么天真的方式：我使用有能力的 libdivide (http://libdivide.com/) 实现非常快的“除以编译时未知的常数”。

这给了我一个很好的胜利（25% 左右），但我有一种唠叨的感觉，我可以做到更好，原因如下：

第一直觉：libdivide 计算除法。我需要的是模数，然后到达那里我必须做一个额外的 mult 和一个 sub ：mod = dividend - divisor*(uint64_t)(dividend/divisor)。我怀疑那里可能会有一个小胜利，使用 libdivide-type 直接产生模数的技术。
第二个直觉：我实际上对模本身不感兴趣。我真正想要的是有效地产生一个均匀分布的整数值，保证严格小于N。

模数是一种相当标准的方法，因为它有两个属性：

A) 如果prng() 是，则mod(prng(), N) 保证均匀分布
B) mod(prgn(), N) 保证属于 [0,N[

但是模是/做更多的只是满足上面的两个约束，事实上它可能做的工作太多了。

所有需要的是一个函数，any 函数遵循约束 A) 和 B) 并且快速。

所以，介绍很长，但这里有两个问题：

是否有与 libdivide 等效的方法直接计算整数模数？
是否有一些整数 X 和 N 的函数 F(X, N) 服从以下两个约束：
- 如果 X 是均匀分布的随机变量，则 F(X,N)也是不均匀分布
- F(X, N) 保证在 [0, N[

(PS : 我知道如果 N 很小，我不需要计算出所有的 64 位 PRNG。事实上，我已经这样做了。但就像我说的，即使是优化与必须计算模数的巨大脂肪损失相比，这是一个小胜利。

编辑：prng() % N 确实不是完全均匀分布的。但是对于足够大的 N，我认为这不是什么大问题（或者是吗？）

编辑 2：prng() % N 确实可能分布非常糟糕。我从来没有意识到它会变得多么糟糕。哎哟。我找到了一篇很好的文章：http://ericlippert.com/2013/12/16/how-much-bias-is-introduced-by-the-remainder-technique

【问题讨论】：

(1) - 在大多数平台上，余数是“免费”计算的，是硬件级别除法的结果。 (2) - 取模和除法都不会给你一个公正的结果。
"A) mod(prng(), N) 保证均匀分布，如果 prng() is" 仅当N 均分M 时才为真，其中prng() 返回数字统一在[0, M[。
你试过std::uniform_int_distribution吗？
“编辑：prng() % N 确实不是完全均匀分布的。但是对于足够大的 N，我认为这不是什么大问题（或者是吗？）”它实际上变得更糟N 变大（即小 N 更好），例如如果 prng() 均匀地生成 0, 1, 2, ..., 15，则 prng() % 4 与 3 有一个小的偏差（3 出现 20% 的时间，0、1、2 出现 26%） , 但prng() % 10 对 0、1、...、5 有很大的偏差（它们的出现频率是 6、7、8、9 的 2 倍）。
你可以试试histogram[ (int)(prng() * (HISTOGRAM_SIZE / (PRNG_MAX + 1.0))) ]。每个直方图预计算一次常数。这编译为一个浮点乘法和一个整数转换。使用 MMX 或 GPU 实现一次做多个。但我同意@OliCharlesworth 的观点，如果你用随机直方图访问来炸毁缓存，这是一个巨大的成本。

标签： c++ performance algorithm random random-sample

【解决方案1】：

在这种情况下，最简单的方法可能效果最好。如果您的 PRNG 足够快，一种非常简单的方法可能会奏效，那就是预先计算比您的 N 的下一个更大的 2 的幂小 1 以用作掩码。即，给定一些看起来像 0001xxxxxxxx 的二进制数（其中 x 表示我们不在乎它是 1 还是 0），我们想要像 000111111111 这样的掩码。

从那里，我们生成如下数字：

生成一个数字
and戴上你的面具
如果结果 > n，转到 1

这种方法的确切效果将取决于 N 与 2 的幂的接近程度。每个连续的 2 幂（显然足够）是其前任的两倍。因此，在最好的情况下，N 正好是 2 的幂次方，并且我们在步骤 3 中的测试总是通过。我们只添加了一个掩码以及与 PRNG 本身所用时间的比较。

在最坏的情况下，N 正好等于 2 的幂。在这种情况下，我们预计会丢弃大约一半生成的数字。

平均而言，N 大约在 2 的幂次方之间。这意味着，平均而言，我们丢弃了大约四分之一的输入。我们几乎可以忽略掩码和自己的比较，因此与“原始”生成器相比，我们的速度损失基本上等于我们丢弃的输出数量，或平均 25%。

【讨论】：

我已经实现了你描述的拒绝方法。它工作得相当好，并且在质量（均匀性）和速度方面都优于模数。但对于大多数 N 来说，它比“通过双精度标准化”方法慢，即使有分支预测提示。

【解决方案2】：

如果您可以快速访问所需的指令，您可以将prng() 乘以N 64 位并返回 128 位结果的高 64 位。这有点像将 [0, 1) 中的统一实数乘以 N 并截断，模数版本的顺序有偏差（即实际上可以忽略不计；这个答案的 32 位版本可能很小但可能明显的偏见）。

另一种探索的可能性是在单个位上运行的无分支模算法上使用字并行性，以批量获取随机数。

【讨论】：

我实现了这个。它运行良好，并且比浮点方法更快（FP 为 20.34 秒，而您的方法为 18.66 秒，因此快了大约 9%）。不错。

【解决方案3】：

Libdivide 或任何其他优化该模数的复杂方法简直是矫枉过正。在您的情况下，唯一明智的方法是

确保您的表格大小是 2 的幂（如果必须添加填充！）

将模运算替换为位掩码运算。像这样：

size_t tableSize = 1 << 16;
size_t tableMask = tableSize - 1;

...

histogram[prng() & tableMask]

位掩码操作是任何物有所值的 CPU 上的单个周期，您无法超越它的速度。

--

注意：
我不知道您的随机数生成器的质量，但使用随机数的最后几位可能不是一个好主意。一些 RNG 在最后一位中产生较差的随机性，而在高位中产生更好的随机性。如果您的 RNG 是这种情况，请使用位移来获取最高有效位：

size_t bitCount = 16;

...

histogram[prng() >> (64 - bitCount)]

这与位掩码一样快，但它使用不同的位。

【讨论】：

如果你不追求 2 的幂怎么办？
@cmaster : 25% 的收益并不是我所说的矫枉过正。就我而言，它非常好。至于四舍五入……嗯。我的直方图很大，所以下一个二的幂是“远”。此外，直方图由来自现实世界的观察值组成，因此您的填充想法是 - 至少可以说 - 不是微不足道的：你会用什么“填充”？
@OliCharlesworth 我明确提到：如果您的用例将您限制为非 2 的表大小，请将其填充为 2 的幂。如果必须的话，丢弃所有会进入填充的随机数，但将物理表设为 2 的幂以避免除法。
填充可以是任何东西，甚至可以是，您将结果地址与正确的最大值进行比较，如果它不在范围内，则将其丢弃。将为您节省空间开销。但是，如果条件对性能也很重要，请尽量避免。在大多数情况下，只需对填充部分进行一些虚假计算就可以逃脱，避免它并没有回报。当然，由于填充，您有开销，但在所有情况下，开销都小于有用工作量，即。 e.在大多数情况下，这是一个明智的权衡。

【解决方案4】：

您可以通过循环将 histogram 扩展为 2 的“大”幂，用一些虚拟值填充尾随空格（保证不会出现在真实数据中）。例如。给定一个直方图

[10, 5, 6]

像这样将其延长到 16（假设 -1 是一个合适的标记）：

[10, 5, 6, 10, 5, 6, 10, 5, 6, 10, 5, 6, 10, 5, 6, -1]

然后可以通过二进制掩码histogram[prng() & mask] 其中mask = (1 << new_length) - 1 进行采样，并检查要重试的标记值，即

int value;
do {
    value = histogram[prng() & mask];
} while (value == SENTINEL);

// use `value` here

通过确保绝大多数元素有效（例如，在上面的示例中，只有 1/16 的查找将“失败”，并且可以通过扩展它进一步降低此比率例如 64)。您甚至可以在检查中使用“分支预测”提示（例如 __builtin_expect in GCC），以便编译器命令代码在 value != SENTINEL 时是最佳的，希望这是常见的情况。

这在很大程度上是内存与速度的权衡。

【讨论】：

现代 x86 实现不支持分支预测提示。哦，没关系，你说的是编译器提示。
不过，我还不清楚（因为我手头没有证据）分支预测未命中率是否超过了使用模数的成本。
@OliCharlesworth，我的意思是__builtin_expect 编译器内在函数，纯粹是为了让编译器可以重新排序代码（如我的回答所述）；我只是花了一点时间来挖掘链接。
@dbaupp ：这与 cmaster 提到的想法相同（填充）。但不幸的是，当 N 不是 2 的幂时，我不知道如何在不扭曲原始直方图的统计特性的情况下“填充”4Meg 直方图。和拒绝......好主意，但可能会缓慢导致数据大小变胖。不过我会试试的，很有趣。编辑 - 我真的很喜欢这个小直方图大小的想法。
@blondiepassesby 我确信在现实世界的直方图中不太可能出现负数。

【解决方案5】：

只是一些想法来补充其他好的答案：

模运算花费了多少时间，您如何知道该百分比是多少？我之所以问，是因为有时人们说某些事情非常慢，而实际上它不到 10% 的时间，他们只是认为它很大，因为他们使用的是一个愚蠢的仅限自拍时间的分析器。（与随机数生成器相比，我很难想象模运算会花费大量时间。）
什么时候知道桶的数量？如果它不经常改变，你可以写一个程序生成器。当桶数发生变化时，自动打印出一个新程序，编译、链接，用于你的海量执行。这样，编译器就会知道桶的数量。
您是否考虑过使用quasi-random number generator，而不是伪随机生成器？它可以在更少的样本中为您提供更高的积分精度。
是否可以减少桶的数量而不会过多地损害积分的准确性？

【讨论】：

1. Monte-Carlo 积分的总时间的 25% 用于取模。通过 perf 和 gprof 验证。当直方图不适合 L2 时，另外 25% 用于访问直方图。我使用的随机数生成器是 Mersenne Twister 的 AVX 优化版本，可以大批量生成随机数。它仅占总计算时间的 7%。模是杀手。
2.桶的数量很少变化，“程序生成”的想法也不错，但在特定情况下 gcc 并没有比 libdivide 做得更好。
3.这是一个有趣的想法，我将进一步对此进行试验（假设它比 AVX-MT 更快，并且我不必对 QRNG 的输出取模来戳我的直方图）。跨度>
4.简短的回答：我不知道。长答案：直方图数据来自现实世界，我没有足够的统计背景来知道在我的计算受到严重影响之前我可以承受多少下降。
@blondiepassesby: 1. 我有点疯狂地告诉人们all the faults with gprof 好像它没有给你行或指令级别的分辨率。（如果 perf 更好，我会感到惊讶。）在我相信百分比之前，我会使用那篇文章中推荐的诊断方法。无论如何，我永远不会认为只有一件事需要解决。如果您可以使用模数获得 30% 的加速，那就太好了，但很可能还有更多。 This explains why. 即另外 68% 的时间是多少？

【解决方案6】：

可以通过拒绝和重绘不小于M*(2^64/M) 的值（在取模之前）来避免不均匀性 dbaupp 的注意事项。
如果M 可以用不超过 32 位表示，则可以通过重复乘法（参见 David Eisenstat 的答案）或 divmod 得到比M 少一个以上的值；或者，您可以使用位操作为M 挑选出足够长的位模式，再次拒绝不小于M 的值。
（我会惊讶于随机数生成的模数在时间/周期/能源消耗方面没有相形见绌。）

【讨论】：

+ 尤其是你的最后一句话。
@greybeard ：模比 RNG 更昂贵。 Mersenne Twister 批处理和矢量化的成本接近于零。它通常可以在 3 个周期内产生 32 位。
@blondiepassesby：由于无法在上面发表评论：一遍又一遍地从直方图中抽取样本有什么用？我看到采样填充直方图的价值。
@greybeard ：直方图是根据来自真实世界（物理）过程的观察值构建的。从某种意义上说，直方图是现实世界过程的“模型”。我正在尝试将此物理过程合并到一个更大系统的模拟中，并且我正在尝试通过蒙特卡洛积分来计算最终状态。为此，我需要从直方图中抽取许多随机样本。

【解决方案7】：

投料桶，可以使用std::binomial_distribution直接投料每个桶，而不是逐个样本地投料桶：

以下可能会有所帮助：

int nrolls = 60; // number of experiments
const std::size_t N = 6;
unsigned int bucket[N] = {};

std::mt19937 generator(time(nullptr));

for (int i = 0; i != N; ++i) {
    double proba = 1. / static_cast<double>(N - i);
    std::binomial_distribution<int> distribution (nrolls, proba);
    bucket[i] = distribution(generator);
    nrolls -= bucket[i];
}

Live example

【讨论】：

我不确定我是否理解您的想法是如何运作的。你能解释一下“喂水桶”是什么意思以及二项分布是如何发挥作用的吗？ bucket[i] 代表什么？我必须承认：我对你的代码有点困惑。
我看错了，我用 feed 而不是 pull，抱歉。据我了解，您循环选择每次选择哪个直方图（桶），我建议以另一种方式计算给定直方图的次数。我不是掷骰子 60 次然后查看分布，而是通过二项分布计算每个数字的掷骰次数。

【解决方案8】：

您可以使用定点数学代替整数除法，即整数乘法和位移。假设您的 prng() 返回 0-65535 范围内的值，并且您希望将其量化到 0-99 范围内，那么您执行 (prng()*100)>>16。只需确保乘法不会溢出您的整数类型，因此您可能必须将 prng() 的结果右移。请注意，此映射比模数更好，因为它保留了均匀分布。

【讨论】：

【解决方案9】：

感谢大家的建议。

首先，我现在完全相信模数确实是邪恶的。
它既非常慢而且在大多数情况下会产生不正确的结果。

在实施和测试了很多建议之后，什么
似乎是最好的速度/质量折衷方案是提出的解决方案
@基因：

预计算 normalizer 为：

auto normalizer = histogram.size() / (1.0+urng.max());
使用以下方法绘制样本：

return histogram[ (uint32_t)floor(urng() * normalizer);

这是迄今为止我尝试过的所有方法中最快的，据我所知，
它产生的分布要好得多，即使它可能不那么完美
作为拒绝方法。

编辑：我实现了 David Eisenstat 的方法，这与 Jarkkol 的建议大致相同：index = (rng() * N) >> 32。它和浮点归一化一样有效，而且速度更快（实际上快了 9%）。所以这是我现在首选的方式。

【讨论】：

来回转换浮动真的比在定点转换更快吗？让我感到惊讶。
@DavidEisenstat ：我还没有尝试过定点方法（准随机积分路线似乎更有希望）。您是指您建议的方法，还是 Jarkkol 的？
区别是（便宜的）移位指令，所以呢？