贝塞尔函数（来自 scipy.special）可以与 Numba 一起使用吗答案

【问题标题】：Can bessel functions (from scipy.special) be used with Numba贝塞尔函数（来自 scipy.special）可以与 Numba 一起使用吗
【发布时间】：2018-03-29 09:51:46
【问题描述】：

我正在尝试使用 Numba 优化对包含 Bessel 函数的函数的积分 (scipy.integrate.quad) 的评估。

虽然 Numba 似乎适用于“常见”numpy 函数，但当我尝试包含 Bessel 函数时它会引发错误：

Untyped global name 'jn': cannot determine Numba type of <class 'numpy.ufunc'>

通过谷歌搜索，我从 Numba 存储库中找到了一个 Jupyter 笔记本，其中讨论了创建 j0 函数 (https://github.com/numba/numba/blob/08d5c889491213288be0d5c7d726c4c34221c35b/examples/notebooks/j0%20in%20Numba.ipynb)。

在 numba 中创建该功能的笔记本电脑会很快，但它们最后显示的计时结果表明使用 numba 的性能降低了约 100 倍。我在这里遗漏了什么明显的东西吗？

更一般地说，是否有可能从 Numba 编译中受益于 scipy Bessel 函数？

【问题讨论】：

不是我的专业领域，但这些函数可能是相当优化的 fortran 代码，这可能是 numba 的麻烦来源（不确定）并且肯定会限制函数本身的优化。
看起来 scipy 函数在 C (github.com/scipy/scipy/tree/…) 中。我想如果 Numba 可以只使用这些文件中的代码会很棒，但是从网上大量阅读来看，我还没有找到解决方案！

标签： python optimization jit numba numpy-ufunc

【解决方案1】：

通常 NumPy 和 SciPy 提供非常快的实现。另一方面，Numba 会基于 Python 函数自动生成代码。

因此，除了 Python 函数之外，您不能将 numba 应用于任何东西，如果您想要 nopython 模式（如果您对速度感兴趣，则需要它）甚至不是每个 Python 函数。 Numba 仅支持非常有限的一组函数和类型。而且这些函数都是在 Numba 中重新实现的，它根本不使用 Python 或 NumPy 函数，即使看起来像！

因此，您拥有自动生成的 LLVM 代码与高度优化的定制 C/Fortran 代码。所以你不应该期望获得任何东西（尽管与 NumPy/SciPy 函数相比，numba 通常对非常小的数组表现出色 - 但即使是中等大小的数组，numba 也会更慢）。

但如果您想在 numba jitted 函数中使用某些（当前不受支持的）函数，您必须自己重新实现它。除非在长的紧密循环中调用它，否则不值得麻烦，只需使用普通函数即可。开发人员时间通常远比运行时间重要。

这并不意味着 numba 不好。它非常适合需要大量计算/循环且无法使用现有 NumPy 或 SciPy 函数实现的任务。

【讨论】：

感谢您的洞察力。通过使用 Numba，我觉得通过在编译代码中将我的 intragrand 传递给 scipy 的 quad（使用其他非 Bessel 函数进行测试确认 > 10 倍的加速）而不是 Python 代码，我可以获得显着的性能提升。但是，这只有在 Numba 可以编译目标函数时才有效——如果存在 Bessel 函数，这现在看起来是不可能的。这里有没有可以提高性能的替代方案？
这取决于瓶颈在哪里。如果瓶颈是贝塞尔函数，你可能什么也做不了（除了寻找更快的实现——可能存在也可能不存在）。但如果瓶颈在其他地方，您可以使用 Cython。 Cython 比 numba 稍微复杂一些，需要编译它（但是，如果您使用 IPython 或 Jupyter 笔记本，这真的很容易）。也可以使用 numbas CFFI 挂钩到 bessel 函数的 SciPy 实现，但这有点超出我的能力（我不确定它是否会成功）。

【解决方案2】：

这是您要查找的代码：

# Bessel function of order 1 - note that smaller arguments are added first
@nb.jit(nopython = True, nogil = True, cache = False)
def Bessel1(z):
    if z.real <= 8.:
        t = z / 8.
        fz = z * (-0.2666949632 * t**14 + 1.7629415168000002 * t**12 + -5.6392305344 * t**10 + 11.1861160576 * t**8 + -14.1749644604 * t**6 + 10.6608917307 * t**4 + -3.9997296130249995 * t**2 + 0.49999791505)
    else:
        t = 8. / z
        eta = z - 0.75 * cmath.pi
        fz = cmath.sqrt(2 / (cmath.pi * z)) * ((1.9776e-06 * t**6 + -3.48648e-05 * t**4 + 0.0018309904000000001 * t**2 + 1.00000000195) * cmath.cos(eta) - (-6.688e-07 * t**7 + 8.313600000000001e-06 * t**5 + -0.000200241 * t**3 + 0.04687499895 * t) * cmath.sin(eta))
    return fz

不过，我不记得这个特定实现的精度。请注意，我已经为复数执行了此操作，您可能需要也可能不需要。

大约一年前，我从一本旧的计算教科书中提取了这个。这有点像泰勒级数展开。具体方法叫什么我忘记了。

请注意，代码中必须先添加“t”的高次幂。这是因为 t 总是小于 1，如果将较小的浮点数添加到较大的浮点数中，则舍入误差会大得多。

【讨论】：