大型数组的插值和外插答案

【问题标题】：Interpolation and extrapolation for large arrays大型数组的插值和外插
【发布时间】：2014-07-17 12:54:24
【问题描述】：

我有一个大型数组y 定义在一个非均匀有序网格x 上。数组的长度通常为 N~2^14 到 N~2^18。我想获得数组的样条插值（或二次）。我面临的问题是，即使 N 的值较低，插值也需要很长时间。

import numpy as np
from scipy.interpolate import interp1d
N = 2 ** 12 # = 4096
x = np.linspace(0, 2*np.pi, N)
y = np.sin(x)
%time f = interp1d(x, y, 'cubic', )

CPU times: user 8min 5s, sys: 1.39 s, total: 8min 7s
Wall time: 8min 7s

我看到的一个选项是，我只需要非常有限的一组数据点的插值值。有没有办法只在被要求时计算插值？

您能否提出一个替代方案，该替代方案还具有对低于x.min() 和高于x.max() 的值进行外推的功能？

谢谢！

【问题讨论】：

InterpolatedUnivariateSpline 怎么样？

标签： python numpy scipy interpolation

【解决方案1】：

作为对@HYRY 评论的补充，并支持他使用InterpolatedUnivariateSpline 的建议，这里是我使用一组不同的数组长度进行的一些基准测试。

interp1d 似乎没有很好地缩放，如下所示（y 轴是 每个点的对数时间 [大多数负值对应于每个插值点的最快计算]，x 轴 power 2 在N)。

即使interp1d 表现最好（接近N=2**4 或2**5）InterpolatedUnivariateSpline 也快了大约 2.5 个数量级。绘制代码如下图所示。

import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import interp1d,InterpolatedUnivariateSpline

t=[]

for i in range(2,24):
  N = 2 ** i
  x = np.linspace(0, 2*np.pi, N)
  y = np.sin(x)
  t_=time.time()
  for j in range(20):#to make results more robust
    #f=interp1d(x,y,kind=3)
    f = InterpolatedUnivariateSpline(x, y, k=3)
  t_=time.time()-t_
  t.append(np.log(t_/N))
plt.plot(np.arange(22)+2,t)

请注意，InterpolatedUnivariateSpline 会为大型输入数组消耗更多内存，因此如果考虑到这一点，interp1d 可能是更好的选择。

【讨论】：

【解决方案2】：

对于非均匀分布的横坐标，您可能需要考虑一种广义插值技术（例如 B 样条）。

将您的数据近似为多个系数乘以基函数的总和（例如，具有非均匀选择节点的 B 样条 - 或位置良好的高斯函数的径向基函数网络）。这些函数必须跨越感兴趣的空间。

现在您可以使用最小二乘法来近似加权系数 - 然后以您需要的任何分辨率重新采样。如果您采用这种方法，您还可以根据平滑度对系统进行正则化，以在 x.min() 和 x.max() 之外提供更合理的值。

这是搭配方法：假设您的样本值在向量 x,y 中。将基向量设置为 phi_k(x) 的采样版本

然后设置基础 B = c_[phi_1,phi_2,...,phi_M] 并使用最小二乘法：c,res,rnk,sv = lstsq(B,y)。

如果基多项式的数量很少 - 那么这可以很快。

现在你的向量 c 包含了系数。您可以通过构建在那里采样的新基向量来计算兴趣点的新值：Bnew = c_[phi_1_new,phi_2_new,...,phi_M_new]

和投影 y_new = dot(Bnew,c)

此方法可让您轻松控制以增加您选择的任何类型的正则化
并在任意点对系统重新采样
使用任何对您的问题有意义的基函数
如果 M 足够小，则系统可以快速求解

【讨论】：