没有任何循环的 Scipy 快速一维插值答案

【问题标题】：Scipy Fast 1-D interpolation without any loop没有任何循环的 Scipy 快速一维插值
【发布时间】：2013-01-28 10:09:13
【问题描述】：

我有两个二维数组 x(ni, nj) 和 y(ni,nj)，我需要在一个轴上进行插值。我想为每个 ni 沿最后一个轴进行插值。

我写的

import numpy as np
from scipy.interpolate import interp1d

z = np.asarray([200,300,400,500,600])
out = []
for i in range(ni):
    f = interp1d(x[i,:], y[i,:], kind='linear')
    out.append(f(z))
out = np.asarray(out)

但是，如果数组大小太大，我认为这种方法效率低下且速度慢。像这样插入多维数组的最快方法是什么？有没有办法在没有循环的情况下执行线性和三次插值？谢谢。

【问题讨论】：

您是否对您的代码进行了任何分析？

标签： python numpy matplotlib scipy

【解决方案1】：

您提出的方法确实有一个 python 循环，因此对于较大的 ni 值，它会变慢。也就是说，除非您要拥有大的ni，否则您不必担心太多。

我使用以下代码创建了示例输入数据：

def sample_data(n_i, n_j, z_shape) :
    x = np.random.rand(n_i, n_j) * 1000
    x.sort()
    x[:,0] = 0
    x[:, -1] = 1000
    y = np.random.rand(n_i, n_j)
    z = np.random.rand(*z_shape) * 1000
    return x, y, z

并用这两个版本的线性插值对它们进行了测试：

def interp_1(x, y, z) :
    rows, cols = x.shape
    out = np.empty((rows,) + z.shape, dtype=y.dtype)
    for j in xrange(rows) :
        out[j] =interp1d(x[j], y[j], kind='linear', copy=False)(z)
    return out

def interp_2(x, y, z) :
    rows, cols = x.shape
    row_idx = np.arange(rows).reshape((rows,) + (1,) * z.ndim)
    col_idx = np.argmax(x.reshape(x.shape + (1,) * z.ndim) > z, axis=1) - 1
    ret = y[row_idx, col_idx + 1] - y[row_idx, col_idx]
    ret /= x[row_idx, col_idx + 1] - x[row_idx, col_idx]
    ret *= z - x[row_idx, col_idx]
    ret += y[row_idx, col_idx]
    return ret

interp_1 是您的代码的优化版本，遵循 Dave 的回答。 interp_2 是线性插值的矢量化实现，它避免了任何 python 循环。编写这样的代码需要对 numpy 中的广播和索引有充分的了解，并且有些事情的优化程度将低于 interp1d 所做的事情。一个典型的例子是找到插入值的 bin：interp1d 一旦找到 bin，肯定会提前跳出循环，上面的函数将值与所有 bin 进行比较。

所以结果将非常依赖于 n_i 和 n_j 是什么，甚至你的数组 z 要插入的值有多长。如果n_j 很小而n_i 很大，那么您应该期望interp_2 和interp_1 有优势，如果相反的话。较小的 z 应该是 interp_2 的优势，较长的 interp_1。

我实际上已经用各种n_i 和n_j 对这两种方法进行了计时，对于形状为(5,) 和(50,) 的z，这是图表：

因此，对于形状为(5,) 的z，无论何时n_j < 1000，您都应该使用interp_2，并在其他地方使用interp_1。毫不奇怪，形状为(50,) 的z 的阈值不同，现在在n_j < 100 附近。似乎很容易得出结论，如果n_j * len(z) > 5000，您应该坚持使用您的代码，但如果不是，则将其更改为上面的interp_2，但该声明中有大量推断！如果您想进一步自己试验，这里是我用来生成图表的代码。

n_s = np.logspace(1, 3.3, 25)
int_1 = np.empty((len(n_s),) * 2)
int_2 = np.empty((len(n_s),) * 2)
z_shape = (5,)

for i, n_i in enumerate(n_s) :
    print int(n_i)
    for j, n_j in enumerate(n_s) :
        x, y, z = sample_data(int(n_i), int(n_j), z_shape)
        int_1[i, j] = min(timeit.repeat('interp_1(x, y, z)',
                                        'from __main__ import interp_1, x, y, z',
                                        repeat=10, number=1))
        int_2[i, j] = min(timeit.repeat('interp_2(x, y, z)',
                                        'from __main__ import interp_2, x, y, z',
                                        repeat=10, number=1))

cs = plt.contour(n_s, n_s, np.transpose(int_1-int_2))
plt.clabel(cs, inline=1, fontsize=10)
plt.xlabel('n_i')
plt.ylabel('n_j')
plt.title('timeit(interp_2) - timeit(interp_1), z.shape=' + str(z_shape))
plt.show()

【讨论】：

@Tetsuro 如果您发现答案有用，您可以通过单击关闭答案标题的向上箭头来投票。此外，如果您认为它回答了您的问题，您可以通过单击勾勒出的复选标记来接受它。
这已经是一个很好的答案了。如果对interp_2 的工作原理有一个简短的概述，那就更好了。（也许不是逐行，但总的来说，进行插值的方法是什么）
您的线性插值是否“仅”是一条样条曲线，即通过最近网格点的线性近似值 - 而不是通过整个网格拟合线性函数？

【解决方案2】：

一种优化是像这样分配一次结果数组：

import numpy as np
from scipy.interpolate import interp1d

z = np.asarray([200,300,400,500,600])
out = np.zeros( [ni, len(z)], dtype=np.float32 ) 
for i in range(ni):
    f = interp1d(x[i,:], y[i,:], kind='linear')
    out[i,:]=f(z)

这将为您节省一些在您的实现中发生的内存复制，这发生在对out.append(...)的调用中。

【讨论】：