Fast sampler of exponential distribution （高速指数分布サンプラー・擬似焼きなまし法の遷移判定）
(heuristic/exponential_dist_sampler.hpp)

• View this file on GitHub
• Last update: 2024-03-01 09:07:31+09:00
• Include: #include "heuristic/exponential_dist_sampler.hpp"

母数 $\lambda = 1$ の指数分布を模擬するサンプラー．

典型的な用途として，擬似焼きなまし法 (simulated annealing, SA) の遷移判定を（指数関数や対数関数の計算を毎回行うことなく）高速に行える．

使用方法

constexpr int D = 18;
const ExponentialDistSampler<D> eds;  // 2^D 個のサンプルを前計算

// 下位 D bit がランダムに分布した mask を与えると x ~ Ex(1) をサンプル
uint32_t mask;
double x = eds.sample(mask);

double dx;
double T;

// コストの変化が abs(dx), 温度 T のとき焼きなまし法の遷移を受理するか確率的に判定
bool upd = eds.check_sa(abs(dx), T, mask);

Code

#pragma once

#include <cmath>
#include <cstdint>
#include <array>

template <int D> struct ExponentialDistSampler {
    std::array<double, (1 << D)> minuslogps;

    constexpr ExponentialDistSampler() {
        for (int i = 0; i < (1 << D); ++i) minuslogps.at(i) = -log((0.5 + i) / (1 << D));
    }

    double sample(uint32_t random_mask) const {
        return minuslogps.at(random_mask & ((1 << D) - 1));
    }

    // p ~ U(0, 1) => -log(p) ~ Ex(1)
    // P[exp(-|dx| / T) >= p] = P[|dx| <= -log(p) * T]
    bool check_sa(double abs_dx, double T, uint32_t random_mask) const {
        return abs_dx <= minuslogps.at(random_mask & ((1 << D) - 1)) * T;
    }
};
// const ExponentialDistSampler<16> log_ps;

#line 2 "heuristic/exponential_dist_sampler.hpp"

#include <cmath>
#include <cstdint>
#include <array>

template <int D> struct ExponentialDistSampler {
    std::array<double, (1 << D)> minuslogps;

    constexpr ExponentialDistSampler() {
        for (int i = 0; i < (1 << D); ++i) minuslogps.at(i) = -log((0.5 + i) / (1 << D));
    }

    double sample(uint32_t random_mask) const {
        return minuslogps.at(random_mask & ((1 << D) - 1));
    }

    // p ~ U(0, 1) => -log(p) ~ Ex(1)
    // P[exp(-|dx| / T) >= p] = P[|dx| <= -log(p) * T]
    bool check_sa(double abs_dx, double T, uint32_t random_mask) const {
        return abs_dx <= minuslogps.at(random_mask & ((1 << D) - 1)) * T;
    }
};
// const ExponentialDistSampler<16> log_ps;

Back to top page