RegularGridInterpolator 与 interp2d 相比慢得令人难以忍受

Question

考虑以下代码示例：

# %%
import numpy
from scipy.interpolate import interp2d, RegularGridInterpolator

x = numpy.arange(9000)
y = numpy.arange(9000)
z = numpy.random.randint(-1000, high=1000, size=(9000, 9000))
f = interp2d(x, y, z, kind='linear', copy=False)

f2 = RegularGridInterpolator((x, y), z, "linear")

mx, my = np.meshgrid(x, y)
M = np.stack([mx, my], axis=-1)

# %%
%timeit f(x, y)

# %%
%timeit f2(M)

它使用

scipy.interpolate.interp2d

和

scipy.interpolate.RegularGridInterpolator

设置一些示例插值器。上面两个单元格的输出是

每次循环 1.09 秒 ± 4.38 毫秒（7 次运行的平均值 ± 标准差，每次 1 次循环）

和

每次循环 10 秒 ± 17.6 毫秒（7 次运行的平均值 ± 标准差，每次 1 次循环）

分别。

RegularGridInterpolator

比

interp2d

慢大约10倍。问题是

interp2d

在 scipy 1.10.0 中已被 标记为已弃用。新代码应该使用

RegularGridInterpolator

。这对我来说似乎有点奇怪，因为它是一个非常糟糕的替代品。我上面的代码示例可能有问题吗？如何加快插值过程？

Answer 1

你的代码没有问题，这可能是 scipy 中的一个错误。 我已经在 github

上举报了

Answer 2

如果其他人偶然发现这一点：

事实证明，手动计算双线性插值比 scipy 的 RegularGridInterpolator 快 10 倍

。当我说“手动”时，我的意思是：

1. 从相同的

   points

   （假设长度为 N 的

   x_vals

   和长度为 M 的

   y_vals

   ）和

   values

   （假设形状为 NxM 的

   z_grid

   ）数组开始，传递给

   RegularGridInterpolator

   。

2. 给定您想要的坐标（假设为

   x,y

   ），使用

   np.searchsorted()

   查找四个最近点的索引（如果

   x_vals

   /

   y_vals

   为 增加 - 如果其中一个为 减少，则此方法本身有效） ，请参阅下面的代码）。

3. 使用标准公式计算双线性插值，我使用Raymond Hettinger 的代码。

这是我自己的实现的稍微修改的摘录：

# since `x_vals` and `y_vals` are decreasing rather than increasing, we do some trickery on top of `np.searchsorted()` to get the desired indices.
i_x = x_vals.shape[0] - np.searchsorted(np.flip(x_vals), x)
i_y = y_vals.shape[0] - np.searchsorted(np.flip(y_vals), y)

points = (
    (
        x_vals[i_x - 1],
        y_vals[i_y - 1],
        z_grid[i_x - 1, i_y - 1]
    ),
    (
        x_vals[i_x],
        y_vals[i_y - 1],
        z_grid[i_x, i_y - 1]
    ),
    (
        x_vals[i_x - 1],
        y_vals[i_y],
        z_grid[i_x - 1, i_y]
    ),
    (
        x_vals[i_x],
        y_vals[i_y],
        z_grid[i_x, i_y]
    )
)

'''here's an option with list comprehension
points = [
    (
        x_vals[i_x-1+i],
        y_vals[i_y-1+j],
        z_grid[i_x-1+i, i_y-1+j]
    )
    for i, j in [(i, j) for i in range(2) for j in range(2)]
]
'''

return bilinear_interpolation(x, y, points) # see https://stackoverflow.com/a/8662355/22122546

就像我之前说过的，根据我自己的测试，这比

RegularGridInterpolator

快大约 10 倍。我不知道为什么 scipy 不在文档中提供有关此的警告，我在这个问题上浪费了几个小时，哈哈。但至少有一个干净的解决方案:>

保重！

Answer 3

从 SciPy 1.13 开始，RGI 张量积样条模式（线性、三次和五次）应该更快。虽然 1.13 尚未发布，但您可以使用 scipy nightly builds。