如何提高Python中估算π的性能

Question

我用Python编写了以下代码，以估算Pi的值。它被称为蒙特卡罗方法。显然，通过增加样本数量，代码变得更慢，我认为代码中最慢的部分在于采样部分。 如何使它更快？

from __future__ import division
import numpy as np

a = 1
n = 1000000

s1 = np.random.uniform(0,a,n)
s2 = np.random.uniform(0,a,n)

ii=0
jj=0

for item in range(n):
    if ((s1[item])**2 + (s2[item])**2) < 1:
        ii = ii + 1

print float(ii*4/(n))

你是否建议使用其他（更快）的代码？

- NKN

---

你需要计算这个的原因是什么？这已经做过无数次了（例如从文件中读取，或硬编码）。 - Caramiriel

算法本身的缓慢是不可避免的 - 每次你想要更多位数的精度，运行时间就会增加10倍。请查看“圆周率-追求更多位数的现代探索”了解更快的方法。 - Kevin

有没有理由提前构建两个数组？为什么不在循环内部获取两个 random.random() 数字呢？但是，正如 @Kevin 指出的那样，这个算法基本上是 O(n) 的，因此对于大的 n，对精确实现的任何更改都只会对总运行时间产生最小的影响。 - jonrsharpe

@Caramiriel 不要问“为什么”。问“为什么不呢？”学习算法的工作原理、学校、个人经验和乐趣都是原因。 - mbomb007

@Caramiriel 我知道这已经做了很多年，但我只是出于教育目的而使用它。谢谢！ - NKN

Answer 1

瓶颈实际上在于您的for循环。Python的for循环相对较慢，因此如果您需要迭代一百万个项目，则完全避免它们可以获得很多速度。在这种情况下，这很容易做到。取而代之的是：

for item in range(n):
    if ((s1[item])**2 + (s2[item])**2) < 1:
        ii = ii + 1

请执行以下操作：

ii = ((s1 ** 2 + s2 ** 2) < 1).sum()

这能够实现是因为numpy具有内置的支持优化数组操作的功能。循环发生在c中而不是Python中，所以它会更快。我进行了一个简单的测试，让你可以看到差异：

>>> def estimate_pi_loop(x, y):
...     total = 0
...     for i in xrange(len(x)):
...         if x[i] ** 2 + y[i] ** 2 < 1:
...             total += 1
...     return total * 4.0 / len(x)
... 
>>> def estimate_pi_numpy(x, y):
...     return ((x ** 2 + y ** 2) < 1).sum()
... 
>>> %timeit estimate_pi_loop(x, y)
1 loops, best of 3: 3.33 s per loop
>>> %timeit estimate_pi_numpy(x, y)
100 loops, best of 3: 10.4 ms per loop

这里有一些可行的操作示例，让您了解它是如何工作的。

将数组平方：

>>> a = numpy.arange(5)
>>> a ** 2
array([ 0,  1,  4,  9, 16])

添加数组：

>>> a + a
array([0, 2, 4, 6, 8])

比较数组：

>>> a > 2
array([False, False, False,  True,  True], dtype=bool)

布尔值求和：

>>> (a > 2).sum()
2

也许你已经意识到，有更快的方法来估算圆周率，但我必须承认，我一直很欣赏这种方法的简单和有效性。

- senderle

Answer 2

你已经分配了NumPy数组，因此你应该充分利用它们。

for item in range(n):
    if ((s1[item])**2 + (s2[item])**2) < 1:
        ii = ii + 1

变成

s1sqr = s1*s1
s2sqr = s2*s2
s_sum = s1sqr + s2sqr
s_sum_bool = s_sum < 1
ii = s_sum_bool.sum()
print float(ii*4/(n))

你需要对数组进行平方运算，然后将它们相加，检查总和是否小于1，接着将布尔值的和（false = 0, true = 1）相加以得到符合条件的总数。

- James

Answer 3

我赞成senderle的答案，但如果您不想太大改变您的代码：

numba是专门为此目的设计的库。

只需将算法定义为函数，并添加@jit装饰器：

 from __future__ import division
 import numpy as np
 from numba import jit

 a = 1
 n = 1000000

 s1 = np.random.uniform(0,a,n)
 s2 = np.random.uniform(0,a,n)

 @jit
 def estimate_pi(s1, s2):
     ii = 0
     for item in range(n):
         if ((s1[item])**2 + (s2[item])**2) < 1:
             ii = ii + 1
     return float(ii*4/(n))

 print estimate_pi(s1, s2)

在我的笔记本电脑上，n = 100000000 的速度提高了约20倍，n = 1000000 的速度提高了约3倍。

- gbrener