用于 PyOpenCl 的具有高斯分布的随机数生成器

Question

在我大量使用 PyOpenCl 的 monte carlo 模拟中，使用 numpy 生成高斯随机数被证明是瓶颈。

np.random.randn(int(1e9))

因此我也在寻找一种使用 PyopenCl 生成高斯分布随机数的方法。

我发现一个 6 年前的帖子问过类似的问题。但我不确定如何将 VexCL 库与 PyOpenCl 一起使用： Gaussian distributed random numbers in OpenCL

有什么想法可以实现与 PyOpenCl 中的 np.random.randn 类似的良好 RNG 吗？

Answer 1

似乎 pyopencl 已经包含了一个随机数生成器：

https://github.com/inducer/pyopencl/blob/master/pyopencl/clrandom.py

运行一个简单的测试表明，均值和标准差与 numpy 的实现相当。直方图也对应于均方误差可忽略不计的正态分布。

有人知道进一步的测试来检查随机数生成器的质量吗？

编辑：根据https://documen.tician.de/pyopencl/array.html

import pytest
from pyopencl.clrandom import RanluxGenerator, PhiloxGenerator
import pyopencl as cl
import pyopencl.array as cl_array
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm


def test_compare_to_numpy_rand():
    context = cl.create_some_context()
    queue = cl.CommandQueue(context)
    rng =  PhiloxGenerator(context=context)
    mu=0
    sigma=3
    size = (int(1e8))
    dtype = np.float32
    numbers_cl = cl_array.zeros(queue, size, dtype=dtype)
    rng.fill_normal(ary=numbers_cl, mu=mu, sigma=sigma)
    numbers_cl_host = numbers_cl.get()
    mu_meas = np.mean(numbers_cl_host)
    sigma_meas = np.std(numbers_cl_host)

    hist_cl, bins_edges = np.histogram(numbers_cl_host, bins=1000, normed=True)
    delta = bins_edges[1] - bins_edges[0]
    bins_center = bins_edges[:-1]+delta/2
    plt.bar(bins_center, hist_cl,width=delta,label='cl')

    numbers_np = mu + sigma * np.random.randn(size)
    hist_np, bins_edges = np.histogram(numbers_np, bins=bins_edges, normed=True)
    plt.bar(bins_center, hist_np, width=delta,alpha=0.8,label='numpy')

    p = norm.pdf(bins_center, mu, sigma)
    plt.plot(bins_center, p, color='red',label='Exact')
    plt.title('Mean squared error: cl={:.5E} np={:.5E}'.format(np.mean(np.abs(p-hist_cl)**2),np.mean(np.abs(p-hist_np)**2)))
    plt.legend()
    plt.show()

    assert pytest.approx(mu_meas, 1e-2) == mu
    assert pytest.approx(sigma_meas, 1e-2) == sigma
    assert pytest.approx(mu_meas, 1e-2) == numbers_np.mean()
    assert pytest.approx(sigma_meas, 1e-2) == numbers_np.std()

用于 PyOpenCl 的具有高斯分布的随机数生成器

Random number generator with Gaussian distribution for PyOpenCl

python

numpy

opencl

pyopencl