使用 numpy 向量化或映射来加速循环 - Python NumPy 3D 矩阵“摆脱循环” Python 问题，蒙特卡洛答案

【问题标题】：Use numpy vectorize or map to speed up a loop - Python NumPy 3D matrix "get rid of a loop" Python question, Monte Carlo使用 numpy 向量化或映射来加速循环 - Python NumPy 3D 矩阵“摆脱循环” Python 问题，蒙特卡洛
【发布时间】：2021-12-13 06:16:03
【问题描述】：

现在我有 1 个循环来填充 3D NumPy 矩阵。尽管我知道它实际上只是我习惯于在（2D）中思考的正常 XxY 的 XxYxZ 表示，但我并不是最擅长理解 3D 数组结构。因此，如果您想知道这是什么，它是用于财务问题的蒙特卡洛模拟中的布朗桥 (BB) 结构。原始代码的信用（源自修复作者 Kenta Oono 的原始帖子的评论，位于此处）：https://gist.github.com/delta2323/6bb572d9473f3b523e6e。你真的不需要知道它背后的数学。它基本上只是切断了一个步骤路径（在本例中为 21），从 0 开始，应用正态分布的冲击（因此 np.random.randn）直到它到达终点，这也是 0。每条路径都应用于模拟价格随时间随机“冲击”它，生成资产在到期时可能遵循的潜在路径。虽然这些是完全不相关的，所以我想我也会传入一个 V 矩阵来关联正确的路径，但是，让我们保持简单：

import numpy as np
from matplotlib import pyplot
import timeit

steps = 21
underlyings = 3
sims = 131072

seed = 0 # fix the seed for replicating results
np.random.seed(seed)

def sample_path_batches(underlyings, steps, sims):
    dt = 1.0 / (steps-1)
    dt_sqrt = np.sqrt(dt)
    B = np.empty((underlyings, steps, sims), dtype=float)
    B[:,0, :] = 0 # set first step to 0
    for n in range(steps - 2):
        t = n * dt
        xi = np.random.randn(underlyings, sims) * dt_sqrt
        B[:, n + 1, :] = B[:, n, :] * (1 - dt / (1 - t)) + xi
        B[:, -1, :] = 0 # set last step to 0
    return B

start_time = timeit.default_timer()
B = sample_path_batches(underlyings, steps, sims)
print('\n' + 'Run time for ', sims, ' simulation steps * underlyings: ', 
np.round((timeit.default_timer() - start_time),3), ' seconds')

pyplot.plot(B[:,:,np.random.randint(0,sims)].T); # plot a random simulation set of paths
pyplot.show()

131072 个模拟步骤的运行时间 * 底层代码：2.014 秒

无论如何，这对我的应用程序来说太慢了，尽管我最初的带有第二个内循环的版本大约需要 15 秒。所以我已经看到人们通过 np.vectorize 对 NumPy 进行矢量化或使用映射来“展平”循环，但我无法想象自己如何实际做到这一点。我正在寻找将产生相同数字的最佳“本机 Python”实现。 B 是 3D NumPy 数组。如果需要，您可以复制并粘贴它并在线运行它：https://mybinder.org/v2/gh/jupyterlab/jupyterlab-demo/HEAD?urlpath=lab/tree/demo

欢迎提出任何建议！！！即使它只是“像这样重组循环，然后应用 np.vectorize”或其他什么，我也很擅长提出建议并使其从一个简单的“新观点”中解决如何可视化问题。我通常只会在 Cython（nogil / OpenMP / prange）中做这种类型的事情，但我想知道在一般情况下“扁平化”一个循环，在 NumPy 或 Pandas 中内置普通数学库或其他任何工作。

【问题讨论】：

不，np.vectorize 和 map 都不会显着提升与普通 for 循环相比的性能。
map 或 np.vectorize 都不会加快您的代码速度。
感谢 cmets - 我想 Cython 是最好的选择。
您可以尝试@numba.vectorize 以尽量减少重写，但您可能会遇到 Numba 限制。 Cython 需要大量类型注释才能工作，但总体而言这可能是一件好事。
是的，numba @jit 从来没有为我编写代码的方式带来任何加速；称它为我不了解编译更快例程所需的格式，或者我的例程无法使用该软件加速。过去，真正加快我的代码速度的唯一因素是 Cython / nogil / inline / OpenMP / prange 语句，尽管转换它的努力和时间很长；全部使用 NumPy 内存视图。是的，有一个函数可以在 Cython 中执行 Mersenne Twister，但它不使用任何这些优化：github.com/ananswam/cython_random

标签： python numpy vectorization montecarlo

【解决方案1】：

加速此代码的一个简单解决方案是使用 Numba 并行化它。您只需将装饰器@nb.njit('float64[:,:,::1](int64, int64, int64)', parallel=True) 用于函数sample_path_batches（其中nb 是Numba 模块）。注意dtype=float必须在函数中替换为dtype=np.float64，这样Numba才能正确编译代码。请注意，parallel=True 应该自动并行化 np.random.randn 调用以及循环中的基本后续操作。在 10 核机器上，速度要快 7 倍（使用 Numpy 需要 0.253 秒，使用 Numba 的并行实现需要 0.036 秒）。如果您没有看到任何改进，您也可以尝试使用 prange 手动并行化它。

此外，您可以使用 np.float32 类型来显着提高性能（理论上最多快 2 倍）。但是，Numpy 目前不支持 np.random.randn 的此类类型。相反，np.random.default_rng().random(size=underlyings*sims, dtype=np.float32).reshape(underlyings, sims)should be used。不幸的是，Numba 可能还不支持它，因为 Numpy 最近添加了这个......

如果您有 Nvidia GPU，另一种解决方案是使用 CUDA 在 GPU 上执行功能。这应该快得多。请注意，Numba 具有特定的优化函数，可以使用 CUDA 在 GPU 上生成随机的np.float32 值（请参阅here）。

【讨论】：

谢谢老兄，这正是我想要的。我从来没有得到适合 numba 的格式，因此从未见过实际的加速，但在我的 PC 上，它运行在 0.25-0.28 秒！看起来 Numba 并行格式遵循 NumPy Cython memoryview 格式。我尝试切换到 np.float32（对于内部的所有变量，可以将其更改为 float64 以及 B），但速度与 Numba 相同。所以我有足够的内存来处理 float64，所以这个解决方案效果很好！
你说得对，顺便说一句，Numba 还不支持 NumPy 的新生成器对象——正确的是 PCG64 用于标准正态分布。当我添加 Numba 以生成 32 位浮点数时，Numba 就会崩溃，错误消息说“我不知道如何处理这个函数”这样的错误信息被松散地翻译......
猜测使用 float32 类型的速度不会更好，因为随机值是使用 float64 生成的，而计算随机值的时间可能是这里的瓶颈。 Numba 还很年轻，很多高级功能还没有（完全）实现。这仍然是一个很好的加速内核的工具，而 Numpy 并不简单。对于使用 Numpy 函数的顺序代码，Numba 通常不会（太多）快，因为 Numpy 已经得到了很好的优化。但是，它在基本循环和可以并行化的代码上大放异彩。
是的，Jérôme Richard 似乎 Numba 的想法是将其应用于 NumPy 特定例程，而我在源代码中看到的大部分内容无论如何都是在 float64 上运行。尚不支持 SciPy 的新生成器，但“似乎”它们可能是由附加项目 github.com/numba/numba-scipy 提供的，但不幸的是，这不支持包含 QMC (Sobol) rng 生成器例程的 Scipy 1.7.1。这实际上是我更喜欢使用而不是 Mersenne twister 的方法，但对于可用的方法，Numba 可以解决问题。没有 Cython 的快速循环。总是很方便:) 非常感谢。