使用 numba 进行 PI 的蒙特卡洛模拟对于最少的点数来说是最慢的？

Question

作为家庭作业练习的一部分，我正在用 Python 实现 pi 的蒙特卡罗模拟。我正在使用 Numba 来加速和并行化计算。从我之前执行的测试中，我发现并行化

numba

比非并行化

numba

版本和纯

numpy

版本运行得更快，所以这就是我选择的。问题要求测量

n = 10**k

的

k = [1, 2, 3, 4, 5, 6, 7]

算法的性能，当我实现这个问题时，我得到了一个非常奇怪的结果：具有最低

n = 10

值的第一次运行运行了所有算法中的 slowest 。

n = 10

的运行始终比

n = 10,000,000

慢 3 倍：

起初，我认为第一次运行很慢，因为

numba

需要编译我的函数，但我正在使用缓存，所以这个借口只能在我第一次运行脚本时起作用。从上图中可以看到，我第一次运行脚本时，

n = 10

的运行速度比其他时候慢得多，因为它需要编译函数。然而，下一次使用缓存版本仍然比使用更大的

n

运行慢得多。

另一个奇怪的结果是，有时，

n = 1,000

的运行速度比

n = 100

快。我认为这个算法应该是时间线性的。我不明白为什么

n

值最低的第一次运行结果是最慢的，以及为什么某些运行通常比

n

值较小的先前运行更快。有人可以解释我得到的结果吗？

这是我使用的代码：

from datetime import timedelta
from time import perf_counter

from numba import jit, prange
import numpy as np
import numpy.typing as npt

jit_opts = dict(
    nopython=True, nogil=True, cache=True, error_model="numpy", fastmath=True
)

rng = np.random.default_rng()


@jit(**jit_opts, parallel=True)
def count_points_in_circle(points: npt.NDArray[float]) -> tuple[npt.NDArray[bool], int]:
    in_circle = np.empty(points.shape[0], dtype=np.bool_)
    in_circle_count = 0
    for i in prange(points.shape[0]):
        in_ = in_circle[i] = points[i, 0] ** 2 + points[i, 1] ** 2 < 1
        in_circle_count += in_
    return in_circle, in_circle_count


def monte_carlo_pi(n: int) -> tuple[npt.NDArray[float], npt.NDArray[bool], float]:
    points = rng.random((n, 2))
    in_circle, count = count_points_in_circle(points)
    return points, in_circle, 4 * count / n


def main() -> None:
    n_values = 10 ** np.arange(1, 8)
    for n in n_values:
        start = perf_counter()
        points, in_circle, pi_approx = monte_carlo_pi(n)
        end = perf_counter()
        duration = end - start
        delta = timedelta(seconds=duration)
        elapsed_msg = (
            f"[{delta} (Raw time: {duration} s)]"
            if delta
            else f"[Raw time: {duration} s]"
        )
        print(
            f"n = {n:,}:".ljust(20),
            f"\N{GREEK SMALL LETTER PI} \N{ALMOST EQUAL TO} {pi_approx}".ljust(20),
            elapsed_msg,
        )


if __name__ == "__main__":
    main()

Answer 1

n = 10 的最低值的第一次运行运行速度最慢。

这是由于惰性编译造成的。事实上，我可以通过“急切地编译函数”来重现它并解决问题。您可以通过手动指定函数的签名来做到这一点： @jit('(float64[:,:],)', **jit_opts, parallel=True)

我认为第一次运行很慢，因为 numba 需要编译我的函数，但我正在使用缓存，所以这个借口只能在我第一次运行脚本时起作用。

在实践中，缓存是脆弱的。例如，如果环境发生变化或函数中使用的全局变量或者您不从文件运行脚本，则函数将被重新编译。此外，重新加载存储设备的功能并不是免费的，但通常快于0.3秒。

另一个奇怪的结果是，有时，n = 1,000 的运行速度比 n = 100 更快。我认为这个算法应该在时间上是线性的。

它们的运行时间都只有几十微秒。这是一个非常小的时间。在这种规模下，时间安排往往不太稳定。一个问题是代码使用

多线程

和生成线程，共享工作并等待结果很慢：在大多数 PC 机器上通常至少需要几十微秒。事实上，在大型服务器上可能需要数百微秒（生成的线程数量越多，开销就越高）。可变性的来源之一是调度程序、内核和 C 状态（内核可以休眠以消耗更少的能量，但它们需要更多时间来唤醒）以及 P 状态（即频率缩放）的可用性。 在我的机器上，对于 n=100 和 n=1000，不使用线程比使用线程更快。大部分开销似乎来自数组分配。与两者的差异相比，两次运行的

变异性相当大

。在基准测试中计算“标准偏差”（甚至时间分布）对于衡量其准确性非常重要。它使您能够通过计算 p 值来判断一种情况是否“在统计上比另一种情况更快”。仅靠时机是不够的。较长的执行时间可能是由于噪音，例如 CPU 密集型后台程序唤醒一小段时间（运气不好）、冷缓存效应等。