python - scipy 插值速度 : Is there anything faster than RectBivariateSpline?

Question

我正在从事一个严重依赖 2D 插值的项目。经过一些初步搜索后，我发现很多东西都说只要我对结构化网格没问题， map 坐标或可能的 RegularGridInterpolator 就应该提供最佳速度。我使用 scipy 中的所有多元插值函数创建了以下代码作为试验。似乎 RectBivariateSpline 是明显的赢家。

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)

产生以下结果:

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]

我在我的实现中遗漏了什么吗？当 RectBivariateSpline 是高阶拟合并且与 map 坐标相比还允许不均匀的网格间距时，这似乎很难相信。是否有其他一些更快的 2D 插值？请原谅这段代码的任何明显问题，Python 和我才刚刚认识。

谢谢，

Answer 1

这种行为的解释似乎是不同的函数在每次调用时都有不同的 Python 开销。我们可以对调用进行矢量化以减少这种开销，如以下代码所示:

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

这是示例输出:

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 是最快的，正如您所料，尽管 RectBivariateSpline 紧随其后(请注意，我将其阶数减少到 1)。