大数组的插值与外推

Question

我有一个很大的数组 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() 以上的值进行外推。

谢谢！

Answer 1

为了补充@HYRY的评论，并支持他建议使用 InterpolatedUnivariateSpline 的观点，我做了一些不同数组长度的基准测试。

从下面的结果来看，interp1d 的表现似乎不太理想（纵轴是每个点的计算时间的对数，越负越快，横轴是 N 的2的幂）。

即使在 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 可能是更好的选择。

Answer 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足够小，那么这个系统可以快速解决。