performance - 一种快速逼近大尺寸 np.random.dirichlet 的方法

Question

我想尽快评估大尺寸的 np.random.dirichlet。更准确地说，我想要一个至少快 10 倍的函数。根据经验，我观察到这个函数的小维度版本输出一到两个具有 0.1 顺序的条目，而其他所有条目都非常小，以至于它们无关紧要。但这一观察并非基于任何严格的评估。近似值不需要那么准确，但我想要一些不太粗糙的东西，因为我正在将这种噪声用于 MCTS。

def g():
   np.random.dirichlet([0.03]*4840)

>>> timeit.timeit(g,number=1000)
0.35117408499991143

score 1 · Accepted Answer

假设您的 alpha 在组件上是固定的并用于多次迭代，您可以将相应 gamma 分布的 ppf 制成表格。这可能是可用的，scipy.stats.gamma.ppf但我们也可以使用scipy.special.gammaincinv. 此功能似乎相当缓慢，因此这是一笔可观的前期投资。

这是一般想法的粗略实现：

import numpy as np
from scipy import special

class symm_dirichlet:
    def __init__(self, alpha, resolution=2**16):
        self.alpha = alpha
        self.resolution = resolution
        self.range, delta = np.linspace(0, 1, resolution,
                                        endpoint=False, retstep=True)
        self.range += delta / 2
        self.table = special.gammaincinv(self.alpha, self.range)
    def draw(self, n_sampl, n_comp, interp='nearest'):
        if interp != 'nearest':
            raise NotImplementedError
        gamma = self.table[np.random.randint(0, self.resolution,
                                             (n_sampl, n_comp))]
        return gamma / gamma.sum(axis=1, keepdims=True)

import time, timeit

t0 = time.perf_counter()
X = symm_dirichlet(0.03)
t1 = time.perf_counter()
print(f'Upfront cost {t1-t0:.3f} sec')
print('Running cost per 1000 samples of width 4840')
print('tabulated           {:3f} sec'.format(timeit.timeit(
    'X.draw(1, 4840)', number=1000, globals=globals())))
print('np.random.dirichlet {:3f} sec'.format(timeit.timeit(
    'np.random.dirichlet([0.03]*4840)', number=1000, globals=globals())))

样本输出：

Upfront cost 13.067 sec
Running cost per 1000 samples of width 4840
tabulated           0.059365 sec
np.random.dirichlet 0.980067 sec

最好检查它是否大致正确：

performance - 一种快速逼近大尺寸 np.random.dirichlet 的方法

1 回答 1

Related

Reference