在二维数组上使用numpy.interp的最快方法

Question

我有以下问题。我正试图找到在一个二维 x 坐标数组上使用 numpy 插值方法的最快方式。

import numpy as np

xp = [0.0, 0.25, 0.5, 0.75, 1.0]

np.random.seed(100)
x = np.random.rand(10)
fp = np.random.rand(10, 5)

基本上，xp 是数据点的 x 坐标，x 是包含我想要插值的值的 x 坐标的数组，fp 是包含数据点 y 坐标的 2-D 数组。

xp
[0.0, 0.25, 0.5, 0.75, 1.0]

x
array([ 0.54340494,  0.27836939,  0.42451759,  0.84477613,  0.00471886,
        0.12156912,  0.67074908,  0.82585276,  0.13670659,  0.57509333])

fp
array([[ 0.89132195,  0.20920212,  0.18532822,  0.10837689,  0.21969749],
       [ 0.97862378,  0.81168315,  0.17194101,  0.81622475,  0.27407375],
       [ 0.43170418,  0.94002982,  0.81764938,  0.33611195,  0.17541045],
       [ 0.37283205,  0.00568851,  0.25242635,  0.79566251,  0.01525497],
       [ 0.59884338,  0.60380454,  0.10514769,  0.38194344,  0.03647606],
       [ 0.89041156,  0.98092086,  0.05994199,  0.89054594,  0.5769015 ],
       [ 0.74247969,  0.63018394,  0.58184219,  0.02043913,  0.21002658],
       [ 0.54468488,  0.76911517,  0.25069523,  0.28589569,  0.85239509],
       [ 0.97500649,  0.88485329,  0.35950784,  0.59885895,  0.35479561],
       [ 0.34019022,  0.17808099,  0.23769421,  0.04486228,  0.50543143]])

所期望的结果应该是这样的：

array([ 0.17196795,  0.73908678,  0.85459966,  0.49980648,  0.59893702,
        0.9344241 ,  0.19840596,  0.45777785,  0.92570835,  0.17977264])

再次寻找最快的方法来处理这个问题，因为这只是一个长度约为1百万与10之间简化版本的问题。

谢谢

- Eric B

---

现在你如何获得所需的输出？看起来你想让x中的每个元素插值到fp中相应的行，但我可能会误解，因为我不得不通过反向工程来确定你的值。 - MB-F

最初，我一直在使用 pandas 进行工作，但现在我正在尝试脱离 pandas，主要是因为速度问题。所以我想到了一个通过欺骗来得出答案的办法。但是，是的，你的理解是正确的。今天早上我会测试你们的答案，谢谢。附带说明，这是使用 pandas 的我的问题链接（稍微修改了问题）：根据列值插值数据帧中的值 - Eric B

Answer 1

所以基本上您想要的输出等效于

np.array([np.interp(x[i], xp, fp[i]) for i in range(x.size)])

但是对于大的 x.size，那个 for 循环会使它变得非常缓慢。

这样做应该可以：

def multiInterp(x, xp, fp):
    i, j = np.nonzero(np.diff(np.array(xp)[None,:] < x[:,None]))
    d = (x - xp[j]) / np.diff(xp)[j]
    return fp[i, j] + np.diff(fp)[i, j] * d

编辑：这个方案甚至更好，能够处理更大的数组：

def multiInterp2(x, xp, fp):
    i = np.arange(x.size)
    j = np.searchsorted(xp, x) - 1
    d = (x - xp[j]) / (xp[j + 1] - xp[j])
    return (1 - d) * fp[i, j] + fp[i, j + 1] * d

测试：

multiInterp2(x, xp, fp)
Out: 
array([ 0.17196795,  0.73908678,  0.85459966,  0.49980648,  0.59893702,
        0.9344241 ,  0.19840596,  0.45777785,  0.92570835,  0.17977264])

使用原始数据进行计时测试：

    %timeit multiInterp2(x, xp, fp)
The slowest run took 6.87 times longer than the fastest. This could mean that an intermediate result is being cached.
10000 loops, best of 3: 25.5 µs per loop

%timeit np.concatenate([compiled_interp(x[[i]], xp, fp[i]) for i in range(fp.shape[0])])
The slowest run took 4.03 times longer than the fastest. This could mean that an intermediate result is being cached.
10000 loops, best of 3: 39.3 µs per loop

即使对于较小的 x 大小，似乎也更快。

让我们尝试一些更大的东西：

n = 10000
m = 10000

xp = np.linspace(0, 1, n)
x = np.random.rand(m)
fp = np.random.rand(m, n)

%timeit b()  # kazemakase's above
10 loops, best of 3: 38.4 ms per loop

%timeit multiInterp2(x, xp, fp)
100 loops, best of 3: 2.4 ms per loop

优势比np.interp的编译版本更易于扩展

- Daniel F

Answer 2

np.interp基本上是编译好的numpy.core.multiarray.interp的包装器。我们可以直接使用它来提高一些性能：

from numpy.core.multiarray import interp as compiled_interp

def a(x=x, xp=xp, fp=fp):
    return np.array([np.interp(x[i], xp, fp[i]) for i in range(fp.shape[0])])

def b(x=x, xp=xp, fp=fp):
    return np.concatenate([compiled_interp(x[[i]], xp, fp[i]) for i in range(fp.shape[0])])

def multiInterp(x=x, xp=xp, fp=fp):
    i, j = np.nonzero(np.diff(xp[None,:] < x[:,None]))
    d = (x - xp[j]) / np.diff(xp)[j]
    return fp[i, j] + np.diff(fp)[i, j] * d

时间测试显示，对于这些示例数组，这与Daniel Forsman的优雅解决方案相当:

%timeit a()
10000 loops, best of 3: 44.7 µs per loop

%timeit b()
10000 loops, best of 3: 32 µs per loop

%timeit multiInterp()
10000 loops, best of 3: 33.3 µs per loop

更新

对于较大的数组，multiInterp是最优选择:

n = 100
m = 1000

xp = np.linspace(0, 1, n)
x = np.random.rand(m)
fp = np.random.rand(m, n)

%timeit a()
100 loops, best of 3: 4.14 ms per loop

%timeit b()
100 loops, best of 3: 2.97 ms per loop

%timeit multiInterp()
1000 loops, best of 3: 1.42 ms per loop

但对于更大的文件，它就会落后:

n = 1000
m = 10000

%timeit a()
10 loops, best of 3: 43.3 ms per loop

%timeit b()
10 loops, best of 3: 32.9 ms per loop

%timeit multiInterp()
10 loops, best of 3: 132 ms per loop

最后，对于非常大的数组（我使用的是32位），临时数组会成为一个问题：

n = 10000
m = 10000

%timeit a()
10 loops, best of 3: 46.2 ms per loop

%timeit b()
10 loops, best of 3: 32.1 ms per loop

%timeit multiInterp()
# MemoryError

- MB-F

---

嗯，是的。我认为很多事情都取决于xp有多大。如果太大了，我的布尔矩阵会出现“内存错误”。也许我可以找到更好的方法…… - Daniel F

看看你能否打破我的新版本，我认为它应该可以处理更大的数组。 - Daniel F

@DanielForsman 很好地使用了 searchsorted。我不想破坏你的版本，但如果我现在有时间，也许可以帮助进一步优化它 :) - MB-F