使用多元随机生成器提高迭代二维 Numpy 数组的性能

Improving performance of iterative 2D Numpy array with multivariate random generator

UxU 周期域中,我模拟了二维数组的动态,其中的条目表示 x-y 坐标。在每个时间步,"parent" 条目被从它们的正态分布 "offsprings" 中选择的新坐标替换,保持数组大小相同。举例说明:

import numpy as np
import random
np.random.seed(13)

def main(time_step=10):
    def dispersal(self, litter_size_):
        return np.random.multivariate_normal([self[0], self[1]], [[sigma**2*1, 0], [0, 1*sigma**2]], litter_size_) % U

    U = 10
    sigma = 2.
    parent = np.random.random(size=(4,2))*U

    for t in range(time_step):
        offspring = []
        for parent_id in range(len(parent)):
             litter_size = np.random.randint(1,4)   # 1-3 offsprings reproduced per parent  
             offspring.append(dispersal(parent[parent_id], litter_size))
        offspring = np.vstack(offspring)

        indices = np.arange(len(offspring))
        parent = offspring[np.random.choice(indices, 4, replace=False)]  # only 4 survives to parenthood

    return parent

但是,该功能可能效率低下 运行,表示为:

from timeit import timeit
timeit(main, number=10000)

那 returns 40.13353896141052 秒。

对 cProfile 的快速检查似乎将 Numpy 的 multivariate_normal 函数确定为瓶颈。

有没有更高效的方法来实现这个操作?

是的,如果您将 Numpy 中的许多函数用于单个数字,它们的成本相对较高,如本例中的 multivariate_normal 所示。因为后代的数量在 [1, 3] 的狭窄范围内,所以 pre-compute 个随机样本是值得的。我们可以在 mean=(0,0) 周围采样,并在迭代过程中添加 parents.

的实际坐标

内循环也可以矢量化。结果:

def main_2(time_step=10, n_parent=4, max_offspring=3):
    U = 10
    sigma = 2.
    cov = [[sigma**2, 0], [0, sigma**2]]
    size = n_parent * max_offspring * time_step
    samples = np.random.multivariate_normal(np.zeros(2), cov, size)

    parents = np.random.rand(n_parent, 2) * U
    for _ in range(time_step):
        litter_size = np.random.randint(1, max_offspring+1, n_parent)
        n_offspring = litter_size.sum()
        parents = np.repeat(parents, litter_size, axis=0)
        offspring = (parents + samples[:n_offspring]) % U
        samples = samples[n_offspring:]
        parents = np.random.permutation(offspring)[:n_parent]

    return parents

我得到的时间是:

In [153]: timeit(main, number=1000)
Out[153]: 9.255848071099535

In [154]: timeit(main_2, number=1000)
Out[154]: 0.870663221881841