为何Gamma分布采样速度远慢于正态分布采样？

我正在为论文开展模型测试模拟，其中一个可变参数是分布类型，但代码运行速度未达预期。经过性能分析（profiling），我发现Gamma分布采样与正态分布采样的耗时差异极大：从VTune生成的火焰图可见，Gamma采样耗时212秒，而正态采样仅14秒。

我是否应对此差异感到意外？若这并非实现层面的问题，我不想耗费大量时间研究所用RNG（随机数生成器）库的源码。我并非RNG专家，但依稀记得正态采样有Box-Muller等高效专用算法，这是否能解释如此大的耗时差异？

补充说明：我因并行能力及与R生态的集成性，使用C++ RNG库TRNG。主循环由Fortran编写，因此需要编写接口代码：

#include <trng/yarn2.hpp>
#include <trng/normal_dist.hpp>
#include <trng/gamma_dist.hpp>

using namespace trng;

extern "C" {

void* init_rng_f(int seed) {

  void* res = (void*) new yarn2((unsigned long) seed);

  return(res);

}

void* init_lrng_f(void* rng, int thread_total, int thread_num) {

  yarn2* ptr = static_cast<yarn2*>(rng);

  yarn2* lrng = new yarn2(*ptr);

  (*lrng).split((unsigned int) thread_total, (unsigned int) thread_num);

  void* res = (void*) &lrng;

  return(res);
}

double rnorm_par(void* lrng, double mean, double sd) {
  normal_dist<> dist(mean, sd);
  yarn2* ptr = static_cast<yarn2*>(lrng);

  double res = dist(*ptr);
  return(res);
}

double rgamma_par(void* lrng, double shape, double scale) {
  gamma_dist<> dist(shape, scale);
  yarn2* ptr = static_cast<yarn2*>(lrng);
  double res = dist(*ptr);
  return(res);
}

}

Fortran部分：

use iso_c_binding

implicit none

interface

      function rnorm_par(lrng, mean, sd) result(res) bind(c)
         import
         type(c_ptr), value :: lrng
         real(c_double), value :: mean
         real(c_double), value :: sd
         real(c_double) :: res
      end function rnorm_par

      function rgamma_par(lrng, shape, scale) result(res) bind(c)
         import
         type(c_ptr), value :: lrng
         real(c_double), value :: shape
         real(c_double), value :: scale
         real(c_double) :: res
      end function rgamma_par

      function init_rng_f(seed) result(res) bind(c)
         import
         integer(c_int), value :: seed
         type(c_ptr) :: res
      end function init_rng_f

      function init_lrng_f(rng, thread_total, thread_num) result(res) bind(c)
         import
         type(c_ptr), value :: rng
         integer(c_int), value :: thread_total
         integer(c_int), value :: thread_num
         type(c_ptr) :: res
      end function init_lrng_f

end interface

主循环逻辑如下：

rng = init_rng_f(42)

allocate(lrngs(n_threads))

do i = 1, n_threads
   lrngs(i) = init_lrng_f(rng, n_threads, i - 1)
end do

do i_sim = 1, n_sim
  ! Some setup code calculating parameters.
  if (dist == GAMMA) then
    rgamma_par(lrng, shape, scale)
  else
    rnorm_par(lrng, mean, sd)
  end if
  ! Some code to store results.
end do

1. 耗时差异是否正常？ #

是的，Gamma采样比正态采样慢一个数量级完全合理，核心原因是两种分布的采样算法复杂度存在本质差异：

• 正态分布有Box-Muller、Ziggurat这类高效算法，通常仅需2-3个均匀随机数、少量浮点运算（平方根、三角函数或查表）即可生成样本，计算量极低。
• Gamma分布的采样算法复杂得多：
  • 当形状参数shape为整数时，需通过shape个指数分布样本求和实现，而指数分布本身依赖均匀数的对数运算；
  • 非整数形状时，常用Marsaglia-Tsang算法、Acceptance-Rejection类方法，这类算法需要多次生成均匀随机数、额外浮点运算（幂运算、对数、比较判断），甚至可能出现采样失败重试的情况，平均计算量远高于正态采样。

TRNG库的实现符合常规逻辑，这种耗时差异并非库的实现缺陷，属于分布采样的固有复杂度差异。

2. 代码可优化点 #

尽管算法本身有差异，但当前代码存在额外性能损耗，会进一步放大耗时差距：

• 每次采样重复创建分布对象：rnorm_par和rgamma_par中，每次调用都会构造normal_dist<>或gamma_dist<>实例，这些对象的构造涉及参数校验、预计算（如正态分布的标准差倒数、Gamma分布的形状参数相关系数），反复构造会带来不必要的开销。
• init_lrng_f的内存错误：函数返回&lrng（指针的地址），而lrng是局部指针变量，函数返回后该变量会被销毁，导致Fortran端拿到悬空指针，可能引发未定义行为，甚至间接影响性能。

代码优化建议：

优化C++接口：

1. 避免重复构造分布对象：可以将分布对象与RNG绑定，或让Fortran端传递预初始化的分布指针；也可修改接口，将分布参数初始化与采样分离。
2. 修复init_lrng_f的指针问题：直接返回lrng指针即可，无需返回指针的地址。

修改后的C++代码示例：

#include <trng/yarn2.hpp>
#include <trng/normal_dist.hpp>
#include <trng/gamma_dist.hpp>

using namespace trng;

extern "C" {

void* init_rng_f(int seed) {
  return new yarn2((unsigned long)seed);
}

// 修复指针返回问题
void* init_lrng_f(void* rng, int thread_total, int thread_num) {
  yarn2* ptr = static_cast<yarn2*>(rng);
  yarn2* lrng = new yarn2(*ptr);
  lrng->split((unsigned int)thread_total, (unsigned int)thread_num);
  return lrng; // 直接返回指针，而非指针的地址
}

// 预初始化分布对象的版本（示例）
void* init_normal_dist(double mean, double sd) {
  return new normal_dist<>(mean, sd);
}

void* init_gamma_dist(double shape, double scale) {
  return new gamma_dist<>(shape, scale);
}

double sample_normal(void* rng, void* dist) {
  return (*static_cast<normal_dist<>*>(dist))(*static_cast<yarn2*>(rng));
}

double sample_gamma(void* rng, void* dist) {
  return (*static_cast<gamma_dist<>*>(dist))(*static_cast<yarn2*>(rng));
}

// 保留原接口但优化构造问题（如果不想修改Fortran端）
double rnorm_par(void* lrng, double mean, double sd) {
  // 使用thread_local静态变量避免重复构造（TRNG分布对象线程安全，多线程独立使用无问题）
  static thread_local normal_dist<> dist;
  dist.set(mean, sd);
  yarn2* ptr = static_cast<yarn2*>(lrng);
  return dist(*ptr);
}

double rgamma_par(void* lrng, double shape, double scale) {
  static thread_local gamma_dist<> dist;
  dist.set(shape, scale);
  yarn2* ptr = static_cast<yarn2*>(lrng);
  return dist(*ptr);
}

}

Fortran端对应调整：

如果使用预初始化分布的版本，需新增接口函数，在模拟开始前初始化分布对象，循环中仅调用采样函数，减少每次传递参数的开销。

3. 其他性能建议 #

• 批量采样：如果模拟允许，尽量一次性生成多个样本而非单个样本，TRNG支持批量生成，可减少函数调用和循环开销。
• 检查形状参数：如果Gamma分布的形状参数接近整数，可手动实现整数Gamma采样（指数分布求和），可能比通用算法更快。

内容的提问来源于stack exchange，提问作者Othman El Hammouchi