2021SC@SDUSC-PALISADE（十二）具体分析PALISADE中的拒绝采样法

2021SC@SDUSC
PALISADE中离散高斯取样实现的分析专题：

1. 高斯分布采样原理——拒绝采样法
2. PALISADE中的采样方法
3. 具体分析PALISADE中的拒绝采样法

---

具体分析PALISADE中的拒绝采样法

拒绝采样法定义在discretegaussiangenerator.cpp的GenerateInteger方法中。该方法的申明如下：

template <typename VecType>
int32_t DiscreteGaussianGeneratorImpl<VecType>::GenerateInteger(
    double mean, double stddev, size_t n) const 

该采样法返回一个int32_t的样本，传入均值与采样分布的参数。接着是一连串的DEBUG信息：

DEBUG("DiscreteGaussianGeneratorImpl =========");
DEBUG("mean " << mean);
DEBUG("stddev " << stddev);
DEBUG("n " << n);
DEBUG("t " << t);

判断输入的平均值与分布参数是否是无穷大的数，如果是则抛出错误：

if (std::isinf(mean)) {
    PALISADE_THROW(not_available_error,
                   "DiscreteGaussianGeneratorImpl called with mean == +-inf");
  }
if (std::isinf(stddev)) {
    PALISADE_THROW(not_available_error,
                   "DiscreteGaussianGeneratorImpl called with stddev == +-inf");
  }

定义作为返回值的变量result：

typename VecType::Integer result;

在均匀分布中采样，采样一个整数与一个浮点数：

double t = log2(n) * stddev;
std::uniform_int_distribution<int32_t> uniform_int(floor(mean - t),
                                                     ceil(mean + t));
DEBUG("uniform( " << floor(mean - t) << ", " << ceil(mean + t) << ")");
std::uniform_real_distribution<double> uniform_real(0.0, 1.0);

定义一个表示采样成功变量flagSuccess与统计采样次数count和采样次数的上限limit

bool flagSuccess = false;
usint count = 0;
const usint limit = 10000;

采样的核心部分：只要没有成功采样就一直循环，直到超过采样次数的上限limit，每一次循环：先在均匀分布中随机取一个整数x，再从均匀分布中随机取一个浮点数dice，比较dice与非正态高斯分布的最优值的大小，如果小于则保留采样x，作为结果输出

while (!flagSuccess) {
    //  pick random int
    x = uniform_int(PseudoRandomNumberGenerator::GetPRNG());

    //  roll the uniform dice
    double dice = uniform_real(PseudoRandomNumberGenerator::GetPRNG());
    //  check if dice land below pdf
    if (dice <= UnnormalizedGaussianPDFOptimized(mean, sigmaFactor, x)) {
      flagSuccess = true;
    }
    // DEBUG("x "<<x<<" dice "<<dice);
    count++;
    if (count > limit) {
      DEBUG("x " << x << " dice " << dice);
      PALISADE_THROW(
          not_available_error,
          "GenerateInteger could not find success after repeated attempts");
    }
  }

return x;