如何在 C++ 中有效地生成排序的均匀分布的随机数？

Question

我想在 C++ 中生成大量有序且均匀分布的随机数，n, （即n >= 1,000,000,000）。

一首我认为，简单的方法是

1. n使用std::uniform_real_distribution<double>,顺序生成均匀分布的数字，
2. 然后使用std::sort.

但是，这需要几分钟时间。

一个第二和更先进的方法是做并行的两个步骤为：

template <typename T>
void computeUniformDistribution(std::vector<T>& elements)
{
    #pragma omp parallel
    {
        std::seed_seq seed{distribution_seed, static_cast<size_t>(omp_get_thread_num())};
        std::mt19937 prng = std::mt19937(seed);
        std::uniform_real_distribution<double> uniform_dist(0, std::numeric_limits<T>::max());

        #pragma omp for
        for (size_t i = 0; i < elements.size(); ++i)
        {
            elements[i] = static_cast<T>(uniform_dist(prng));
        }
    }

    std::sort(std::execution::par_unseq, elements.begin(), elements.end());
}

但是，即使这样也需要大约30秒。鉴于均匀分布数字的生成只需要大约1.5秒，瓶颈仍然是排序阶段。

因此，我想问以下问题：如何以排序的方式有效地生成均匀分布的数据？

Answer 1

有一些方法可以生成已经排序的样本，但我认为生成部分排序的样本可能会更好。

将输出范围划分为 k 个等宽的桶。每个桶中的样本数将具有相等概率的多项分布。对多项式分布进行采样的慢速方法是在 [0, k) 中生成 n 个整数。一种更有效的方法是绘制 k 个泊松样本，其速率 n/k 以它们的总和不超过 n 为条件，然后使用慢速方式添加另一个 n - sum 样本。对泊松分布进行采样很难做到完美，但是当 n/k 非常大时（就像这里一样），泊松分布可以通过对均值和方差为 n/k 的正态分布四舍五入来很好地近似。如果这是不可接受的，那么慢速方法确实可以很好地并行化。

给定桶数，计算前缀总和以找到桶边界。对于并行的每个桶，在桶化范围内生成给定数量的样本并对其进行排序。如果我们选择好 n/k，每个桶几乎肯定会适合 L1 缓存。对于 n = 1e9，我想我会尝试 k = 1e5 或 k = 1e6。

这是一个顺序实现。有点粗糙，因为我们确实需要避免对封闭的桶边界进行 2 倍过采样，但我会把它留给你。我不熟悉 OMP，但我认为您可以通过在SortedUniformSamples.

#include <algorithm>
#include <cmath>
#include <iostream>
#include <numeric>
#include <random>
#include <span>
#include <vector>

template <typename Dist, typename Gen>
void SortedSamples(std::span<double> samples, Dist dist, Gen& gen) {
  for (double& sample : samples) {
    sample = dist(gen);
  }
  std::sort(samples.begin(), samples.end());
}

template <typename Gen>
void ApproxMultinomialSample(std::span<std::size_t> samples, std::size_t n,
                             Gen& gen) {
  double lambda = static_cast<double>(n) / samples.size();
  std::normal_distribution<double> approx_poisson{lambda, std::sqrt(lambda)};
  std::size_t sum;
  do {
    for (std::size_t& sample : samples) {
      sample = std::lrint(approx_poisson(gen));
    }
    sum = std::accumulate(samples.begin(), samples.end(), std::size_t{0});
  } while (sum > n);
  std::uniform_int_distribution<std::size_t> uniform{0, samples.size() - 1};
  for (; sum < n; sum++) {
    samples[uniform(gen)]++;
  }
}

template <typename Gen>
void SortedUniformSamples(std::span<double> samples, Gen& gen) {
  static constexpr std::size_t kTargetBucketSize = 1024;
  if (samples.size() < kTargetBucketSize) {
    SortedSamples(samples, std::uniform_real_distribution<double>{0, 1}, gen);
    return;
  }
  std::size_t num_buckets = samples.size() / kTargetBucketSize;
  std::vector<std::size_t> bucket_counts(num_buckets);
  ApproxMultinomialSample(bucket_counts, samples.size(), gen);
  std::vector<std::size_t> prefix_sums(num_buckets + 1);
  std::partial_sum(bucket_counts.begin(), bucket_counts.end(),
                   ++prefix_sums.begin());
  for (std::size_t i = 0; i < num_buckets; i++) {
    SortedSamples(std::span<double>{&samples[prefix_sums[i]],
                                    &samples[prefix_sums[i + 1]]},
                  std::uniform_real_distribution<double>{
                      static_cast<double>(i) / num_buckets,
                      static_cast<double>(i + 1) / num_buckets},
                  gen);
  }
}

int main() {
  std::vector<double> samples(100000000);
  std::default_random_engine gen;
  SortedUniformSamples(samples, gen);
  if (std::is_sorted(samples.begin(), samples.end())) {
    std::cout << "sorted\n";
  }
}

如果您的标准库具有高质量的 实现poisson_distribution，您也可以这样做：

template <typename Gen>
void MultinomialSample(std::span<std::size_t> samples, std::size_t n,
                       Gen& gen) {
  double lambda = static_cast<double>(n) / samples.size();
  std::poisson_distribution<std::size_t> poisson{lambda};
  std::size_t sum;
  do {
    for (std::size_t& sample : samples) {
      sample = poisson(gen);
    }
    sum = std::accumulate(samples.begin(), samples.end(), std::size_t{0});
  } while (sum > n);
  std::uniform_int_distribution<std::size_t> uniform{0, samples.size() - 1};
  for (; sum < n; sum++) {
    samples[uniform(gen)]++;
  }
}

Answer 2

我倾向于依赖这样一个事实，即一组有序分布的变量的连续元素之间的差异呈指数分布。这可以被利用来及时运行O(N)而不是O(N*log N).

快速实现将执行以下操作：

template<typename T> void
computeSorteUniform2(std::vector<T>& elements)
{
    std::random_device rd;
    std::mt19937 prng(rd());

    std::exponential_distribution<T> dist(static_cast<T>(1));

    auto sum = dist(prng);

    for (auto& elem : elements) {
        elem = sum += dist(prng);
    }

    sum += dist(prng);

    for (auto& elem : elements) {
        elem /= sum;
    }
}

假设您想要 Uniform(0, 1) 中的值，此示例已简化，但它应该很容易概括。使用 OMP 完成这项工作并非易事，但应该不会太难。

如果您关心最后大约 50% 的性能，那么有一些数字技巧可能会加速生成随机偏差（例如，有比 MT 更快更好的 PRNG）以及将它们转换为doubles（但最近的编译器可能知道这些技巧）。一些参考资料：Daniel Lemire 的博客和Melissa O'Neill 的 PCG 网站。

我刚刚对此进行了基准测试，发现 clang std::uniform_real_distribution和std::exponential_distribution都很慢。 numpy基于 Ziggurat 的实现快 8 倍，因此我可以使用上述算法double在笔记本电脑上的单个线程（即std实现需要约 80 秒）内在约 10 秒内生成 1e9 。我没有在 1e9 元素上尝试过 OP 的实现，但是我的 1e8 元素快了大约 15 倍。