RegularGridInterpolator excruciatingly slow compared to interp2d

Question

Consider the following code example:

# %%
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)

It sets up some example interpolators using scipy.interpolate.interp2d and scipy.interpolate.RegularGridInterpolator. The output of the two cells above is

1.09 s ± 4.38 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

and

10 s ± 17.6 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

respectively.

The RegularGridInterpolator is about 10 times slower than the interp2d. The problem is that interp2d has been marked as deprecated in scipy 1.10.0. And new code should use RegularGridInterpolator. This seems a bit strange to me since it would be such a bad replacement. Is there maybe a problem in my code example above? How can I speed this interpolation process up?

you could use [RectBivariateSpline](https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.RectBivariateSpline.html)? — Nin17, Feb 12 '23 at 14:23
Good question! I can confirm that performance of RegularGridInterpolator is approx 10x slower than inerp2d. *and* its api is more difficult to use for simple 2d problems with measured points. — Jev, Feb 19 '23 at 21:52
... update - I've tried using RegularGridInterpolator on another example (from [here](https://scipython.com/book/chapter-8-scipy/examples/scipyinterpolateinterp2d/) . It uses cubic interpolation. The difference is dramatic - 175us for interp2d vs 1.2 **seconds** for the new code. I guess we'll just have to stick with scipy 1.10 for now and silence the warnings... a bit sad about this. — Jev, Feb 19 '23 at 22:51
reported as a bug on [github](https://github.com/scipy/scipy/issues/18010) — Jev, Feb 19 '23 at 23:04

score 1 · Accepted Answer · answered Feb 19 '23 at 23:37

1

There is no problem with your code, it's probably a bug in scipy. I've reported it on github

answered Feb 19 '23 at 23:37

Jev

1,163
10
13

1

Commenting here for visibility. It is a known performance regression and we are actively working on it. The new interface was needed in term of maintainability moving forward. – tupui Feb 21 '23 at 18:10

score 1 · Answer 2 · answered Jun 27 '23 at 15:01

In case anyone else stumbles on this:

Turns out that calculating bilinear interpolation manually is about 10x faster than scipy's RegularGridInterpolator. When I say "manually", I mean:

Start with the same points (let's say x_vals with length N and y_vals with length M) and values (let's say z_grid with shape NxM) arrays you'd pass to RegularGridInterpolator.
Given your desired coordinates (let's say x,y), find indices of the four nearest points using np.searchsorted() (this works on its own if x_vals/y_vals is increasing -- if either is decreasing, see code below).
Calculate bilinear interpolation using a standard formula, I use Raymond Hettinger's code.

Here's a slightly modified excerpt from my own implementation:

# 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

Like I said before, this is about 10 times faster than RegularGridInterpolator from my own testing. I have no idea why scipy doesn't provide a warning about this in documentation, I've wasted hours on this lol. But at least there's a clean solution :>

Take care!

RegularGridInterpolator excruciatingly slow compared to interp2d

2 Answers2