scipy.interpolate.LinearNDInterpolator 在大型数据集上无限期挂起

我正在 Python 中插入一些数据，以将其重新网格化到常规网格上，以便我可以部分集成它。数据表示高维参数空间的函数（目前为 3，将扩展到至少 5），并返回可观测值的多值函数（目前为 2，将扩展到 3，然后可能扩展到数十个）。

我正在通过执行插值scipy.interpolate.LinearNDInterpolator http://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.LinearNDInterpolator.html由于缺乏任何其他明显的选择（并且因为我明白griddata无论如何，只是调用它）。在较小的数据集（15,000 行列数据）上它工作得很好。在较大的集合（60,000+）上，该命令似乎无限期地运行。top表示 iPython 正在使用 100% CPU 并且终端完全无响应，包括C-c。到目前为止，我已经等了几个小时了，但毫无结果，最终我想传递几百万个条目。

我怀疑这个问题与这张票 http://mail.scipy.org/pipermail/scipy-user/2011-March/028844.html但这据说是在我昨天升级到的 SciPy 0.10.0 中修补的。

我的问题基本上是如何对大型数据集执行多维插值？根据我的尝试，有几个可能的地方可以提供解决方案，但我没有找到它们。 （我的搜索没有得到帮助，因为 scipy 的几个子域似乎很沮丧 http://www.downforeveryoneorjustme.com/www.scipy.org...)

• 出了什么问题LinearNDInterpolator？或者，至少，我怎样才能找出问题所在并尝试避免挂起？
• 有没有办法重新制定插值以便LinearNDInterpolator将工作？也许通过谨慎地将数据分块来重新网格化？
• 是否还有其他更适合该问题的高维插值器？ （我注意到大多数 SciPy 的替代方案都仅限于
• 是否有其他方法可以将多维数据放到常规的用户定义网格上？这就是我试图通过插值来做的所有事情......

问题很可能是您的数据集太大，因此计算其 Delaunay 三角剖分无法在合理的时间内完成。检查时间缩放scipy.spatial.Delaunay使用从完整数据集中随机挑选的较小数据子集来估计完整数据集计算是否在宇宙结束之前完成。

如果你的原始数据在矩形网格上，例如

v[i,j,k,l] = f(x[i], y[j], z[k], u[l])

那么使用基于三角测量的插值是very效率低下。最好使用张量积插值，即通过一维插值方法依次对每个维度进行插值：

import numpy as np
from scipy.interpolate import interp1d

def interp3(x, y, z, v, xi, yi, zi, method='cubic'):
    """Interpolation on 3-D. x, y, xi, yi should be 1-D
    and z.shape == (len(x), len(y), len(z))"""
    q = (x, y, z)
    qi = (xi, yi, zi)
    for j in range(3):
        v = interp1d(q[j], v, axis=j, kind=method)(qi[j])
    return v

def somefunc(x, y, z):
    return x**2 + y**2 - z**2 + x*y*z

# some input data
x = np.linspace(0, 1, 5)
y = np.linspace(0, 2, 6)
z = np.linspace(0, 3, 7)
v = somefunc(x[:,None,None], y[None,:,None], z[None,None,:])

# interpolate
xi = np.linspace(0, 1, 45)
yi = np.linspace(0, 2, 46)
zi = np.linspace(0, 3, 47)
vi = interp3(x, y, z, v, xi, yi, zi)

import matplotlib.pyplot as plt
plt.subplot(121)
plt.pcolor(xi, yi, vi[:,:,12])
plt.title('interpolated')
plt.subplot(122)
plt.pcolor(xi, yi, somefunc(xi[:,None], yi[None,:], zi[12]))
plt.title('exact')
plt.show()

如果您的数据集分散且对于基于三角测量的方法来说太大，那么您需要切换到其他方法。有些选项是同时处理少量最近邻居的插值方法（可以使用 k-d-树快速检索此信息）。反距离称重是其中之一，但它可能是最糟糕的一种——可能有更好的选择（如果没有进一步研究，我不知道）。