使用 python/numpy 进行时间和内存高效的随机采样答案

【问题标题】：Time and memory-efficient random sampling with python/numpy使用 python/numpy 进行时间和内存高效的随机采样
【发布时间】：2016-08-28 05:04:30
【问题描述】：

我需要计算一个单位正方形内 2 个随机均匀分布的点之间的距离的期望值。这个问题可以解析求解，答案约为 0.521405。我想要完成的是使用随机数生成器和 Python （有和没有 numpy）来得出这个数字。使用下面的代码

dist=list()
random.seed()
for each in range(100000000):
    x1=random.random()
    x2=random.random()
    y1=random.random()
    y2=random.random()
    dist.append(math.sqrt((x1-x2)**2+(y1-y2)**2))
print(round(sum(dist)/float(len(dist)),6))

它在我的机器上在 125 秒内迭代 1 亿次，但只有 4 个十进制数字是正确的。现在，我使用 numpy 创建了以下代码

dist=list()
start_time = time.time()
np.random.seed()
for each in range(100000000):
    x = np.array((np.random.random(),np.random.random()))
    y = np.array((np.random.random(),np.random.random()))
    dist.append(np.linalg.norm(x-y))
print(round(np.mean(dist),6))

同样的 1 亿次迭代需要 1111 秒！

由于结果最多只有 4 位十进制数字是正确的，我尝试将迭代次数增加到 10 亿使用以前的版本，没有 numpy。我想，因为每个浮点数都在大多数 64 位（我使用的是 64 位 Python）列表大约需要 8 GB。但是，当列表有 7.9 亿项

时，程序用尽了 26GB 的内存并抛出异常

所以我正在就以下问题征求您的意见：

有没有办法利用各种 numpy 优化并真正让它工作得更快？
有没有办法让程序更节省内存？我意识到列表是一种比我需要的更复杂的数据结构
我是否正确假设为了正确获得 6 位十进制数字，我需要接近 10^12 的迭代次数？（因为 N 次测量的标准误差随着 1/sqrt(N) 递减）

提前致谢！

【问题讨论】：

你应该创建一个 x 和 y，每个都有 10m 的值 - 查看 random 的文档。然后尝试找到一种方法来一次性计算距离（即没有循环）。

标签： python list numpy random out-of-memory

【解决方案1】：

为了回答前两个子问题，这里有一种方法可以一次性生成所有随机数，然后使用np.einsum 替换np.linalg.norm 的大部分工作（将微分平方和沿rows) ，保持其余代码不变。实现看起来像这样 -

N = 100000 # Edit this to change number of random points
x,y = np.random.random((2,N,2))
diffs = x-y
out = round(np.mean(np.sqrt(np.einsum('ij,ij->i',diffs,diffs))),6)

【讨论】：

你真的需要在这里使用einsum吗？ hypot 不是更好更简单吗？
@dan-man 从性能的角度来看，我发现einsum 非常有效，尤其是在替换平方时，我们也在减少和。所以，我的直觉随之而来:)
@Divakar - 但这里的总和刚好超过长度为 2 的轴，hypot 将在内部执行，并在同一通道中执行 sqrt。另请注意，hypot 重新排列方程以获得更好的精度。
@dan-man 坦白说，我用的不多。或许可以将其发布为答案？

【解决方案2】：

至于内存问题，您关于“浮点数需要 64 位”的假设并不完全正确。每个浮点对象（实际上是 python 中的盒装对象）将占用 24 个字节。您可以使用 sys.getsizeof(0.0) 检查它。因此，您的程序需要的空间应该是您估计的 3 倍，这与您实际体验的差不多。

【讨论】：

【解决方案3】：

如果您使用的是 Python 3，我建议使用此变体不带 numpy：

random.seed()
dist_count = 100000000
dist_sum = 0

for _ in range(dist_count):
    dx = random.random() - random.random()  # x1 - x2
    dy = random.random() - random.random()  # y1 - y2

    dist += math.sqrt( dx*dx + dy*dy )

dist_average = dist_sum / dist_count
print(round(dist_average, 6))

首先，如果我们只是要计算和总结它们，为什么要将所有距离存储在一个列表中？将每个随机距离直接添加到整数变量会更快。此外，我们已经知道我们创建了多少随机距离，因为这是我们在 for 循环的范围中指定的，所以我们也不需要像 len(dist) 这样的东西或单独的计数器。

此外，我们不必为每个坐标指定名称，我们可以即时计算dx 和dy 的差异。这也有助于我们进行下一步。

将相同的值与自身相乘比将其提高到 2 的幂 (further reading...) 快一个数量级。因此，我们当然将a**2 替换为a*a。现在我们直接存储上面的差异变得有用了。

最后，我们将距离总和除以计数并显示一次结果。

【讨论】：

这在内存方面确实更有效，并且它没有在 10 亿次迭代中始终保持在 30 MB 内存以下！但是，与 Divakar 建议的 numpy 示例相反，执行 1 亿次迭代需要 75 秒 的代码。我想这是我们能得到的最快速度
是的，我没有使用numpy的经验，也不知道它有哪些优化功能。它们可能是直接用快速 C 代码编写的，没有 Python 开销，所以我无法与之竞争。出于好奇，您能否告诉我接受的 numpy 答案在内存消耗和可扩展性方面的得分如何？
好吧，由于 numpy 的答案会生成 2 个大小为 N 的数组，并且每个浮点数占用 24 字节（如此处发布的），因此它根本无法扩展，并且在 N=5 亿时内存不足。下面的所有内容都是在不到 30 秒的时间内计算出来的 :) 我想我在一个问题中提出了 3 个问题，这是一个错误：现在我有 2 个答案：时间和记忆方面，我有点迷失了接受哪一个。我会和你一起去，因为 Divakar 有更多的选票，所以访问这个页面的人会看到你的被接受
@NickSlavsky 非常感谢，我很高兴能为您提供帮助。我认为我们可以总结一下，对于中等长度的测试系列，人们应该更喜欢 numpy 的性能，但大型系列需要我的方法的灵活性和低资源使用率。

【解决方案4】：

嗯，hypot 在我的笔记本电脑上要快一点，13.5 秒，einsum 为 15.2 秒内存方面可能无法使用 10 亿，当前版本需要 ~6G

代码：

import numpy as np

np.random.seed(12345)

N = 100000000 # Edit this to change number of random points
x,y = np.random.random((2,N,2))
diffs = x - y

out = round(np.mean( np.hypot(diffs[:,0], diffs[:,1]) ),6)

#out = round(np.mean(np.sqrt(np.einsum('ij,ij->i',diffs,diffs))),6)

print(out)

【讨论】：