Andrew Kensler 的 permute() 函数：无状态常数时间随机遍历

这个副标题有点拗口，但我早就想聊聊 Andrew Kensler 的 permute() 函数了，而且我想一开始就说清楚它为什么有用。permute() 在图形学圈子里有些应用，尤其是在离线渲染中，相关技术在密码学领域也更为人熟知*，但它的适用范围远不止图形学和密码学。

先来说说它解决的问题。你想从一个长度为 n 的数组中，以随机顺序遍历 k 个元素，其中 k ≤ n。这 k 个元素不能有重复（每个必须对应唯一的索引）。或者换个说法，你想无放回地随机抽取 k 个元素。

说到“随机”，这显然是伪随机，但伪随机也有程度之分。这毫无疑问不是密码学安全的，那也不是它的设计用途。我对伪随机的了解不足以判断它在统计上与真随机的不可区分程度，也无法断言对哪些应用来说它足够随机。但用于计算渲染积分已经够用了，做很多日常统计分析我也放心用。

2021 年 5 月 27 日更新：我看到一个看起来能显著改善其随机性、并扩展到支持 64 位整数的改进方案。建议读完本文后去看看！

总之，随机遍历数组的经典做法是先用 Fisher-Yates 洗牌打乱，然后遍历打乱后的数组。

*我找到了一份 2009 年的 Stack Overflow 回答，给出了完全相同的思路。所以要说清楚，我在这里讲的方法并非 Andrew Kensler 首创，但我认为他在图形学领域推广了它的使用，而且本文末尾展示的那个特定哈希函数确实是他创建的。

1 2 3 4 5 6 7 8 template<typename T> void fisher_yates_shuffle(std::vector<T> *vec) { int n = vec->size(); for (int i = 0; i < n-1; i++) { int swap_idx = uniform_rand_int(i, n-1); std::swap((*vec)[i], (*vec)[swap_idx]);

uniform_rand_int(min, max) 对 min 和 max 都包含，即 (\textnormal{swap_idx} \in {i, i+1, \dots, n-1})

显然这会原地修改数组，出于种种原因，这可能并不理想。实践中，你可能会传入一个索引数组或指针数组来调用它。

有很多场景下 Fisher-Yates 洗牌不太好用，因为你需要存储打乱后的数组。

在后两种情况下，显然还有更多需要考虑的因素，比如从大数组中读取不连续值的开销。解决第 2 和第 3 种情况还有其他好方法，尤其是在分布式计算场景下。但我所知的很多技术都有一个共同问题：你至少需要扫描整个数组一次。如果你能以常数时间访问数组，Kensler 的方法耗时始终只与随机样本的大小成线性关系，而与原始数组无关。

换句话说，如果你想从一个长度为 n 的数组中随机遍历 k 个不重复元素，Fisher-Yates 洗牌需要 O(n) 时间和 O(n) 内存。我接下来会展示，这可以降到 O(k) 时间和 O(k) 内存。但 permute() 只需 O(k) 时间和 O(1) 内存，这已经是最优了。不仅如此，用 permute() 你可以在常数时间内任意查询随机索引，这意味着你可以在多个线程或计算机上并行使用它。

我觉得值得一提的是，你可以扩展 Fisher-Yates 洗牌来处理第 2 种情况。用一个稀疏映射来记录已交换的索引，而不是把它们写进数组。交换下一个索引时，检查映射中该数组位置（原始位置或交换后位置）是否已经被交换过。如果是，就用映射中的值；否则用原始索引，然后更新映射。下面是一段代码示例。它会从 ({0,1,\dots,n-1}) 中选出 k 个随机不重复整数。

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 std::vector<int> get_k_indices(int n, int k) { std::vector<int> indices(k); std::unordered_map<int, int> index_map; for (int i = 0; i < k; i++) { int swap_idx = uniform_rand_int(i, n-1); // 获取新索引和旧值。如果之前交换过， // 它们会在映射中。 auto swap_val_it = index_map.find(swap_idx); auto prev_val_it = index_map.find(i); int swap_val = (swap_val_it == index_map.end()) ? swap_idx : swap_val_it->second; int prev_val = (prev_val_it == index_map.end()) ? i : prev_val_it->second; indices[i] = swap_val; index_map[swap_idx] = prev_val; // 不需要在映射中设置 prev_idx， // 它之后再也不会被检查了。 return indices;

但即便如此，这种“增量式”Fisher-Yates 洗牌仍需要 O(k) 内存，因此不适用于第 1 或第 3 种情况。

更新：我认为这等价于 Floyd 的无放回随机抽样算法，但不如后者优雅。所以如果你想用那个算法（比如你处于第 2 种情况且确实在意良好的随机性），请实现 Floyd 算法，而不是我在这里描述的方式。

Kensler 在他的论文 Correlated Multi-Jittered Sampling 中提出了 permute()。我第一次读的时候不太理解它的原理，所以也想写这篇博客来拆解它，希望能帮助其他人理解。不过如果你只想要代码，直接跳到文末！

这个函数充当数组上的“无状态”伪随机迭代器。它接收一个置换的种子值、数组长度，以及你想映射到另一个索引的索引。它长这样：

uint32_t permute(uint32_t idx, uint32_t len, uint32_t seed) { // 以某种方式返回 [0, len-1] 区间内随机打乱后的整数。 ...

（原文较长，此处译出前 6000 字，完整内容请点「打开原文」。）