scipy Interpolation Speed: Is there anything faster than RectBivariateSpline?

Question

I am working on a project that relies heavily on 2D interpolation. After some initial searching I found quite a few things that said that map coordinates or possibly RegularGridInterpolator should provide the best speed as long as I was OK with a structured gird. I created the following code as a trial using all the multivariate interpolation functions in scipy. It seems that RectBivariateSpline is the clear winner.

from __future__ import division
import numpy as np
import time
import math
from scipy.interpolate import RectBivariateSpline, RegularGridInterpolator, interp2d, CloughTocher2DInterpolator, Rbf, griddata, interpn
from scipy.ndimage.interpolation import map_coordinates
import random

grid_size = 100
xrange = np.linspace(0, 1, grid_size)
yrange = np.linspace(0, 1, grid_size)
x2d, y2D = np.meshgrid(xrange, yrange, indexing='ij')

z2D = x2d ** 2.0 + y2D

i2d = interp2d(xrange, yrange, z2D, 'linear')
rbs = RectBivariateSpline(xrange, yrange, z2D)
RGI = RegularGridInterpolator((xrange, yrange), z2D)
CT2d = CloughTocher2DInterpolator(np.array(list(zip(x2d.flatten(), y2D.flatten()))), z2D.flatten())
rbfi = Rbf(x2d, y2D, z2D, function='linear')
grd_inputs1 = np.array(list(zip(x2d.flatten(), y2D.flatten())))
grd_inputs2 = z2D.flatten()

print 'Exact: ',0.5**2.0 + 0.5
print 'interp2d',i2d(0.5, 0.5)[0]
print 'RectBivariateSpline',rbs(0.5, 0.5)[0][0]
print 'RegularGridInterpolator',RGI([0.5, 0.5])[0]

coords = np.asarray(([0.5], [0.5]))
coords = [(c - lo) * (n - 1) / (hi - lo) for (lo, hi), c, n in zip([[0,1],[0,1]], coords, np.shape(z2D.flatten()))]
print 'Map Coordinates:', map_coordinates(z2D.flatten(), coords, order=1)[0]

print 'Rbf', rbfi([0.5, 0.5])[0]
print 'CloughTocher2DInterpolator', CT2d([0.5, 0.5])[0]
print 'griddata', griddata(grd_inputs1, grd_inputs2, [0.5, 0.5], 'linear')[0]
print 'interpn', interpn([xrange, yrange], z2D, [0.5, 0.5])[0]

print ''

samples = 25000

start = time.time()
for n in range(samples):
    i2d(0.5, 0.5)
end = time.time()
print 'interp2d: %0.4f [us]' % (((end-start)/samples)*1e6)

start = time.time()
for n in range(samples):
    rand = random.uniform(0,1)
    rbs(0.5, 0.5)
end = time.time()
print 'RectBivariateSpline: %0.4f [us]' % (((end-start)/samples)*1e6)

start = time.time()
for n in range(samples):
    rand = random.uniform(0,1)
    RGI([0.5, 0.5])
end = time.time()
print 'RegularGridInterpolator: %0.4f [us]' % (((end-start)/samples)*1e6)

start = time.time()
for n in range(samples):
    map_coordinates(z2D.flatten(), coords, order=1)
end = time.time()
print 'Map Coordiantes: %0.4f [us]' % (((end-start)/samples)*1e6)

start = time.time()
for n in range(samples):
    rand = random.uniform(0,1)
    rbfi([0.5, 0.5])[0]
end = time.time()
print 'Rbf: %0.4f [us]' % (((end-start)/samples)*1e6)

start = time.time()
for n in range(samples):
    rand = random.uniform(0,1)
    CT2d([0.5, 0.5])[0]
end = time.time()
print 'CloughTocher2DInterpolator: %0.4f [us]' % (((end-start)/samples)*1e6)

start = time.time()
for n in range(int(samples/100)):
    rand = random.uniform(0,1)
    griddata(grd_inputs1, grd_inputs2, [0.5, 0.5], 'linear')
end = time.time()
print 'griddata: %0.4f [us]' % (((end-start)/int(samples/100))*1e6)

start = time.time()
for n in range(int(samples)):
    rand = random.uniform(0,1)
    #grdd(250e3+rand, bar2Pa(125+rand))[0]
    interpn([xrange, yrange], z2D, [0.5, 0.5])[0]
end = time.time()
print 'interpn: %0.4f [us]' % (((end-start)/int(samples))*1e6)

Yeilds the following:

Exact:  0.75
interp2d 0.7500255076012651
RectBivariateSpline 0.7499999999999997
RegularGridInterpolator 0.750025507601265
Map Coordinates: 0.7500255076012652
Rbf 0.7500000434270742
CloughTocher2DInterpolator 0.750003857158422
griddata 0.7500255076012653
interpn 0.750025507601265

interp2d: 14.4000 [us]
RectBivariateSpline: 2.8400 [us]
RegularGridInterpolator: 90.9200 [us]
Map Coordiantes: 8.9600 [us]
Rbf: 105.8000 [us]
CloughTocher2DInterpolator: 18.0000 [us]
griddata: 138824.0004 [us]
interpn: 158.1600 [us]

Am I missing something in my implementation? This seems hard to believe when RectBivariateSpline is a higher order fit and also allows for non uniform grid spacing compared to map coordinates. Is there some other 2D interpolation that would be faster? Please excuse any glaring issues with this code, Python and I are just barely getting acquainted.

Thanks,

Answer 1

The explanation for this behavior seems to be that the different functions have different Python overhead on each call. We can vectorize the calls to reduce this overhead, as in the following code:

import numpy as np
import time
from scipy.interpolate import RectBivariateSpline, RegularGridInterpolator, CloughTocher2DInterpolator, Rbf, interpn
from scipy.ndimage import map_coordinates

grid_size = 100
xrange = np.linspace(0, 1, grid_size)
yrange = np.linspace(0, 1, grid_size)
x2D, y2D = np.meshgrid(xrange, yrange, indexing='ij')

z2D = x2D ** 2.0 + y2D

rbs = RectBivariateSpline(xrange, yrange, z2D,kx=1,ky=1)
RGI = RegularGridInterpolator((xrange, yrange), z2D)
CT2d = CloughTocher2DInterpolator(np.array(list(zip(x2D.flatten(), y2D.flatten()))), z2D.flatten())
rbfi = Rbf(x2D, y2D, z2D, function='linear')

samples = 25000
xpt = np.random.rand(samples)
ypt = np.random.rand(samples)

start = time.time()
rbs_vals = rbs.ev(xpt, ypt)
end = time.time()
print('RectBivariateSpline: %0.4f [us]' % (((end-start)/samples)*1e6))

start = time.time()
rgi_vals = RGI(np.asarray([xpt, ypt]).T)
end = time.time()
print('RegularGridInterpolator: %0.4f [us]' % (((end-start)/samples)*1e6))

start = time.time()
mc_vals = map_coordinates(z2D, np.asarray([xpt, ypt])*(grid_size-1), order=1)
end = time.time()
print('map_coordinates: %0.4f [us]' % (((end-start)/samples)*1e6))

start = time.time()
rbf_vals = rbfi(xpt,ypt)
end = time.time()
print('Rbf: %0.4f [us]' % (((end-start)/samples)*1e6))

start = time.time()
ct_vals = CT2d(np.asarray([xpt, ypt]).T)
end = time.time()
print('CloughTocher2DInterpolator: %0.4f [us]' % (((end-start)/samples)*1e6))

start = time.time()
in_vals = interpn([xrange, yrange], z2D, np.asarray([xpt, ypt]).T)
end = time.time()
print('interpn: %0.4f [us]' % (((end-start)/int(samples))*1e6))

vals = xpt**2 + ypt
print('\nRMS error:')
print('  RectBivariateSpline: %g' % np.sqrt(np.mean((rbs_vals-vals)**2)))
print('  RegularGridInterpolator: %g' % np.sqrt(np.mean((rgi_vals-vals)**2)))
print('  map_coordinates: %g' % np.sqrt(np.mean((mc_vals-vals)**2)))
print('  Rbf: %g' % np.sqrt(np.mean((rbf_vals-vals)**2)))
print('  CloughTocher2DInterpolator: %g' % np.sqrt(np.mean((ct_vals-vals)**2)))
print('  interpn: %g' % np.sqrt(np.mean((in_vals-vals)**2)))

Here is example output:

RectBivariateSpline: 0.1998 [us]
RegularGridInterpolator: 0.3595 [us]
map_coordinates: 0.0799 [us]
Rbf: 114.0414 [us]
CloughTocher2DInterpolator: 7.0727 [us]
interpn: 0.2797 [us]

RMS error:
  RectBivariateSpline: 1.85823e-05
  RegularGridInterpolator: 1.85823e-05
  map_coordinates: 1.85823e-05
  Rbf: 2.00083e-05
  CloughTocher2DInterpolator: 1.63274e-06
  interpn: 1.85823e-05

map_coordinates is the fastest, as you expect, although RectBivariateSpline is not far behind (note that I reduced its order to 1).