Benchmarking linear N-dimensional interpolators

Rationale

A fast and simple interpolation routine is needed in the Stochastic Dynamic Programming algorithm. The "fast and simple" requirement leads to choose the linear interpolation method.

The method needs to be applicable to N-dimensional input, i.e. to interpolate a scalar function :

$$ x\mapsto f(x)\in\mathbb{R} ,\quad x\in D \subset \mathbb{R}^n$$

where the dimension $n$ will be in practice less than 3 or 4.

We here try to evaluate what routines can be found "off-the-shelf".

References

Update February 2015: add scipy.interpolate.RegularGridInterpolator in the benchmark (added in scipy 0.14, released in May 2014)

In [8]:
import numpy as np
from numpy import pi
import matplotlib.pyplot as plt
import matplotlib as mpl
%matplotlib inline
In [90]:
from scipy.interpolate import (
    LinearNDInterpolator, RectBivariateSpline,
    RegularGridInterpolator)
from scipy.ndimage import map_coordinates
from dolointerpolation import MultilinearInterpolator
In [64]:
# Tweak how images are plotted with imshow
mpl.rcParams['image.interpolation'] = 'none' # no interpolation
mpl.rcParams['image.origin'] = 'lower' # origin at lower left corner
mpl.rcParams['image.cmap'] = 'RdBu_r'

1) Define a simple interpolation problem

In [9]:
def f_2d(x,y):
    '''a function with 2D input to interpolate on [0,1]'''
    twopi = 2*pi
    return np.exp(-x)*np.cos(x*2*pi)*np.sin(y*2*pi)

# quick check :
f_2d(0,0.25)
Out[9]:
1.0
In [10]:
def f_3d(x,y,z):
    '''a function with 3D input to interpolate on [0,1]'''
    twopi = 2*pi
    return np.sin(x*2*pi)*np.sin(y*2*pi)*np.sin(z*2*pi)

# quick check :
f_3d(0.25,0.25,0.25)
Out[10]:
1.0

Build the 2D and 3D data grids

In [24]:
Ndata = 50
xgrid = np.linspace(0,1, Ndata)
ygrid = np.linspace(0,1, Ndata+1) # use a slighly different size to check differences
zgrid = np.linspace(0,1, Ndata+2)

f_2d_grid = f_2d(xgrid.reshape(-1,1), ygrid)

plt.imshow(f_2d_grid.T)
plt.title(u'image of a 2D function ({}² pts)'.format(Ndata));

f_2d_grid.shape
Out[24]:
(50, 51)

The 3D case

In [25]:
f_3d_grid = f_3d(xgrid.reshape(-1,1,1), ygrid.reshape(1,-1,1), zgrid)
f_3d_grid.shape
Out[25]:
(50, 51, 52)

2) Use interpolators

Try several interpolation routines.

Notice how each routine has its special way (API) to

  1. build the interpolator (instanciation)
  2. call the interpolator (evaluation)
In [26]:
# Define the grid to interpolate on :
Ninterp = 1000
xinterp = np.linspace(0,1, Ninterp)
yinterp = np.linspace(0,1, Ninterp+1) # use a slighly different size to check differences
zinterp = np.linspace(0,1, 5) # small dimension to avoid size explosion

a) LinearNDInterpolator

LinearNDInterpolator uses an unstructured data which is provided as a list of points. There is also a cousin : NearestNDInterpolator

LinearNDInterpolator(points, values) (documentation)

Performance :

  • instanciation :

    • 18 ms for 50x50 pts,
    • but 19.1 s for 50x50x50 pts
  • evaluation : 45 ms for 1 Mpts.

In [27]:
# Build data for the interpolator
points_x, points_y = np.broadcast_arrays(xgrid.reshape(-1,1), ygrid)
points = np.vstack((points_x.flatten(),
                    points_y.flatten())).T
values = f_2d_grid.flatten()
In [29]:
# Build
%timeit f_2d_interp = LinearNDInterpolator(points, values)
100 loops, best of 3: 17.7 ms per loop
In [30]:
f_2d_interp = LinearNDInterpolator(points, values)
In [32]:
# Evaluate
%timeit f_2d_interp(xinterp.reshape(-1,1), yinterp)
10 loops, best of 3: 44.7 ms per loop
In [33]:
# Display
plt.imshow(f_2d_interp(xinterp.reshape(-1,1), yinterp).T)
plt.title(u'interpolation of a 2D function ({}² pts)'.format(Ninterp));

The 3D case : terrible performance !

In [34]:
# Build data for the interpolator
points_x, points_y, points_z = np.broadcast_arrays(xgrid.reshape(-1,1,1), ygrid.reshape(1,-1,1), zgrid)
points = np.vstack((points_x.flatten(),
                    points_y.flatten(),
                    points_z.flatten()
                  )).T
values = f_3d_grid.flatten()
points.shape
Out[34]:
(132600, 3)
In [35]:
%time f_3d_interp = LinearNDInterpolator(points, values) 
CPU times: user 19.1 s, sys: 88 ms, total: 19.2 s
Wall time: 19.2 s
In [36]:
f_3d_interp = LinearNDInterpolator(points, values)
In [37]:
# Evaluate:
%time f_3d_interp(xinterp.reshape(-1,1), np.linspace(0,1,5), 0.25);
CPU times: user 2min 44s, sys: 0 ns, total: 2min 44s
Wall time: 2min 44s
Out[37]:
array([[  0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00],
       [  0.00000000e+00,   6.23317279e-03,   7.64852773e-19,
         -6.23317279e-03,  -1.53553034e-18],
       [  0.00000000e+00,   1.24663456e-02,   1.52970555e-18,
         -1.24663456e-02,  -3.07106068e-18],
       ..., 
       [  0.00000000e+00,  -1.25138148e-02,  -1.53553034e-18,
          1.25138148e-02,   3.07106068e-18],
       [  0.00000000e+00,  -6.25690738e-03,  -7.67765169e-19,
          6.25690738e-03,   1.53553034e-18],
       [  0.00000000e+00,  -2.44098405e-16,  -2.99525375e-32,
          2.44098405e-16,   5.99050750e-32]])

Interpolation on a slice along the x-axis, just to check that it works.

In [38]:
f_3d_interp_grid = f_3d_interp(xinterp.reshape(-1,1,1), 0.25, 0.25)
In [40]:
plt.plot(f_3d_interp_grid.flatten());

b) RectBivariateSpline (scipy.interpolate)

for 2D interpolation only !

RectBivariateSpline(x, y, z, bbox=[None, None, None, None], kx=3, ky=3, s=0) (documentation)

  • x,y : 1-D arrays of coordinates in strictly ascending order.
  • z : 2-D array of data with shape (x.size,y.size).

Performance

  • instanciation : 0.2 ms for 50x50 pts
  • evaluation : 21 ms for 1 Mpts
In [58]:
%%timeit # Build
f_2d_interp = RectBivariateSpline(xgrid, ygrid, f_2d_grid)
10000 loops, best of 3: 192 µs per loop
In [59]:
f_2d_interp = RectBivariateSpline(xgrid, ygrid, f_2d_grid)
In [60]:
%%timeit # Evaluate
f_2d_interp(xinterp, yinterp)
10 loops, best of 3: 21.2 ms per loop
In [65]:
# Display
plt.imshow(f_2d_interp(xinterp, yinterp).T)
plt.title(u'interpolation of a 2D function ({}² pts)'.format(Ninterp));

c) map_coordinates (scipy.ndimage)

notice that this method has no separation between instanciation and evaluation !

map_coordinates(input, coordinates, output=None, order=3, mode='constant', cval=0.0, prefilter=True) (documentation)

Performance :

  • 56 ms (instanciation + evaluation 1Mpts)
  • 511 ms (5 Mpts in 3D)
In [48]:
# Prepare the coordinates to evaluate the array on :
points_x, points_y = np.broadcast_arrays(xinterp.reshape(-1,1), yinterp)
coord = np.vstack((points_x.flatten()*(len(xgrid)-1) , # a weird formula !
                   points_y.flatten()*(len(ygrid)-1)))
coord.shape
Out[48]:
(2, 1001000)
In [66]:
%%timeit # Build and Evaluate
f_2d_interp = map_coordinates(f_2d_grid, coord, order=1)
# Reshape
f_2d_interp = f_2d_interp.reshape(len(xinterp), len(yinterp))
10 loops, best of 3: 56.4 ms per loop
In [67]:
# Display 
f_2d_interp = map_coordinates(f_2d_grid, coord, order=1).reshape(len(xinterp), len(yinterp))

plt.imshow(f_2d_interp.T)
plt.title(u'interpolation of a 2D function ({}² pts)'.format(Ninterp));

The 3D case

In [68]:
# Prepare the coordinates to evaluate the array on :
points_x, points_y, points_z = np.broadcast_arrays(xinterp.reshape(-1,1,1), yinterp.reshape(1,-1,1), zinterp)
coord = np.vstack((points_x.flatten()*(len(xgrid)-1), # a weird formula !
                   points_y.flatten()*(len(ygrid)-1),
                   points_z.flatten()*(len(zgrid)-1)
                   ))
coord.shape
Out[68]:
(3, 5005000)
In [69]:
%%timeit # Build and Evaluate
f_3d_interp = map_coordinates(f_3d_grid, coord, order=1)
# Reshape
f_3d_interp = f_3d_interp.reshape(len(xinterp), len(yinterp), len(zinterp))
1 loops, best of 3: 511 ms per loop
In [70]:
f_3d_interp = map_coordinates(f_3d_grid, coord, order=1)
f_3d_interp = f_3d_interp.reshape(len(xinterp), len(yinterp), len(zinterp))
f_3d_interp.shape
Out[70]:
(1000, 1001, 5)
In [72]:
n_z = 1
plt.imshow(f_3d_interp[:,:,n_z])
plt.title('f_3d(x,y, z={:.2f})'.format(zinterp[n_z]))
plt.colorbar();

d) MultilinearInterpolator from Pablo Winant

Performance :

  • 0.010 ms for instanciation. 13 ms for evaluation 1Mpts
  • 105 ms (5 Mpts in 3D)
In [73]:
smin = [xgrid[0], ygrid[0]]
smax = [xgrid[-1], ygrid[-1]]
orders = [len(xgrid), len(ygrid)]

print(smin, smax, orders)
([0.0, 0.0], [1.0, 1.0], [50, 51])
In [75]:
%%timeit
minterp = MultilinearInterpolator(smin,smax,orders)
minterp.set_values(np.atleast_2d(f_2d_grid.flatten()))
100000 loops, best of 3: 9.47 µs per loop
In [76]:
minterp.grid.shape
Out[76]:
(2, 2550)
In [77]:
f_2d_grid.shape
Out[77]:
(50, 51)
In [78]:
# Prepare the coordinates to evaluate the array on :
points_x, points_y = np.broadcast_arrays(xinterp.reshape(-1,1), yinterp)
coord = np.vstack((points_x.flatten(),
                   points_y.flatten()))
coord.shape
Out[78]:
(2, 1001000)
In [79]:
%timeit f_2d_interp = minterp(coord).reshape(len(xinterp), len(yinterp))
100 loops, best of 3: 12.6 ms per loop
In [80]:
# Display 
f_2d_interp = minterp(coord).reshape(len(xinterp), len(yinterp))

plt.imshow(f_2d_interp.T)
plt.title(u'interpolation of a 2D function ({}² pts)'.format(Ninterp));
plt.colorbar();
Out[80]:
<matplotlib.colorbar.Colorbar instance at 0x7fe0c49d0b00>

The 3D case

In [81]:
smin = [xgrid[0], ygrid[0], zgrid[0]]
smax = [xgrid[-1], ygrid[-1], zgrid[-1]]
orders = [len(xgrid), len(ygrid), len(zgrid)]

print(smin, smax, orders)
([0.0, 0.0, 0.0], [1.0, 1.0, 1.0], [50, 51, 52])
In [82]:
minterp = MultilinearInterpolator(smin,smax,orders)
minterp.set_values(np.atleast_2d(f_3d_grid.flatten()))
In [83]:
# Prepare the coordinates to evaluate the array on :
points_x, points_y, points_z = np.broadcast_arrays(xinterp.reshape(-1,1,1), yinterp.reshape(1,-1,1), zinterp)
coord = np.vstack((points_x.flatten(), # a weird formula !
                   points_y.flatten(),
                   points_z.flatten()
                   ))
coord.shape
Out[83]:
(3, 5005000)

A very nice interpolation performance !

In [84]:
%timeit minterp(coord)
10 loops, best of 3: 105 ms per loop
In [85]:
f_3d_interp = minterp(coord)
f_3d_interp = f_3d_interp.reshape(len(xinterp), len(yinterp), len(zinterp))
f_3d_interp.shape
Out[85]:
(1000, 1001, 5)
In [87]:
n_z = 1
plt.imshow(f_3d_interp[:,:,n_z])
plt.title('f_3d(x,y, z={:.2f})'.format(zinterp[n_z]))
plt.colorbar();

Interpolation throuput:

In [88]:
print('{:.1f} Mpts/s'.format(5e6 / 0.105 /1e6) )
47.6 Mpts/s

e) RegularGridInterpolator (scipy.interpolate, starting v0.14)

documentation

actual interpolation code: https://github.com/scipy/scipy/blob/master/scipy/interpolate/interpolate.py#L1577 (pure Python, no Cython involved!)

API:

RegularGridInterpolator(points, values, [...])

  • points: tuple of ndarray of float, with shapes (m1, ), ..., (mn, ).
    → The points defining the regular grid in n dimensions.
  • value: array_like, shape (m1, ..., mn, ...).
    → The data on the regular grid in n dimensions.

Performance

  • 0.022 ms for instanciation, 116 ms for evaluation (1Mpts)
  • in 1.8 s in 3D (5 Mpts)
In [93]:
%%timeit # Instanciate the interpolator
f_2d_interp = RegularGridInterpolator((xgrid, ygrid), f_2d_grid)
10000 loops, best of 3: 21.6 µs per loop
In [110]:
f_2d_interp = RegularGridInterpolator((xgrid, ygrid), f_2d_grid)
In [111]:
# Prepare the coordinates to evaluate the array on :
points_x, points_y = np.broadcast_arrays(xinterp.reshape(-1,1), yinterp)
coord = np.vstack((points_x.flatten(),
                   points_y.flatten()))
coord.shape
Out[111]:
(2, 1001000)
In [113]:
f_2d_interp
Out[113]:
<scipy.interpolate.interpolate.RegularGridInterpolator at 0x7fe0c4933ed0>
In [112]:
%%timeit # Interpolate
f_2d_interp_res = f_2d_interp(coord.T).reshape(len(xinterp), len(yinterp))
10 loops, best of 3: 112 ms per loop
In [108]:
# Display 

f_2d_interp_res = f_2d_interp(coord.T).reshape(len(xinterp), len(yinterp))

plt.imshow(f_2d_interp_res.T)
plt.title(u'interpolation of a 2D function ({}² pts)'.format(Ninterp));
plt.colorbar();

The 3D case

In [116]:
f_3d_interp = RegularGridInterpolator((xgrid, ygrid, zgrid), f_3d_grid)
In [117]:
# Prepare the coordinates to evaluate the array on :
points_x, points_y, points_z = np.broadcast_arrays(xinterp.reshape(-1,1,1), yinterp.reshape(1,-1,1), zinterp)
coord = np.vstack((points_x.flatten(), # a weird formula !
                   points_y.flatten(),
                   points_z.flatten()
                   ))
coord.shape
Out[117]:
(3, 5005000)
In [118]:
%%timeit # Interpolate
f_3d_interp_res = f_3d_interp(coord.T).reshape(len(xinterp), len(yinterp), len(zinterp))
1 loops, best of 3: 1.8 s per loop

table

Method Instanciation Evaluation (1Mpts 2D) (5 Mpts, 3D)
dolo 0.010 ms 13 ms 105 ms
RegularGridInterpolator 0.022 ms 116 ms 1.8 s

Conclusion: Dolo's Multilinear interpolation (for an even grid only) is still the best!

In [ ]: