Python SciPy interpolate.RBFInterpolator用法及代码示例

本文简要介绍 python 语言中 scipy.interpolate.RBFInterpolator 的用法。

用法:

class  scipy.interpolate.RBFInterpolator(y, d, neighbors=None, smoothing=0.0, kernel='thin_plate_spline', epsilon=None, degree=None)#

N 维径向基函数 (RBF) 插值。

参数 ：：

y： (npoints, ndims) 类似数组

数据点坐标的二维数组。

d： (npoints, …) 数组

N-D y 处的数据值数组。 d 沿第一轴的长度必须等于 y 的长度。与某些插补器不同，插补轴无法更改。

neighbors： 整数，可选

如果指定，将仅使用这么多最近的数据点计算每个评估点的插值值。默认使用所有数据点。

smoothing： float 或 (npoints, ) 数组，可选

平滑参数。将其设置为 0 时，插值与数据完美拟合。对于较大的值，插值接近具有指定次数的多项式的最小二乘拟合。默认值为 0。

kernel： str，可选

RBF 类型。这应该是其中之一

• ‘linear’ : -r

• ‘thin_plate_spline’ : r**2 * log(r)

• ‘cubic’ : r**3

• ‘quintic’ : -r**5

• ‘multiquadric’ : -sqrt(1 + r**2)

• ‘inverse_multiquadric’ : 1/sqrt(1 + r**2)

• ‘inverse_quadratic’ : 1/(1 + r**2)

• ‘gaussian’ : exp(-r**2)

默认为‘thin_plate_spline’。

epsilon： 浮点数，可选

缩放 RBF 输入的形状参数。如果内核是‘linear’, ‘thin_plate_spline’、‘cubic’或‘quintic’，则默认为1并且可以忽略，因为它与缩放平滑参数具有相同的效果。否则，必须指定这一点。

degree： 整数，可选

添加多项式的次数。对于某些 RBF，如果多项式次数太小，则插值可能不是well-posed。这些 RBF 及其相应的最小度数是

• ‘multiquadric’ : 0

• ‘linear’ : 0

• ‘thin_plate_spline’ : 1

• ‘cubic’ : 1

• ‘quintic’ : 2

默认值为内核的最小度数，如果没有最小度数，则为 0。将此设置为 -1 表示不添加多项式。

注意：

RBF 是 N 维空间中的标量值函数，其在 \(x\) 处的值可以用 \(r=||x - c||\) 表示，其中 \(c\) 是 RBF 的中心。

来自位置 \(y\) 的数据值向量 \(d\) 的 RBF 插值是以 \(y\) 为中心的 RBF 加上具有指定次数的多项式的线性组合。 RBF 插值写为

\[f(x) = K(x, y) a + P(x) b,\]

其中 \(K(x, y)\) 是 RBF 的矩阵，其中心在 \(y\) 在点 \(x\) 处评估，而 \(P(x)\) 是单项式矩阵，它跨越具有指定次数的多项式，在 \(x\) 处评估。系数 \(a\) 和 \(b\) 是线性方程组的解

\[(K(y, y) + \lambda I) a + P(y) b = d\]

和

\[P(y)^T a = 0,\]

其中\(\lambda\) 是一个非负平滑参数，用于控制我们想要拟合数据的程度。当平滑参数为 0 时，数据完全拟合。

如果满足以下要求，则上述系统是唯一可解决的：

• \(P(y)\) must have full column rank. \(P(y)\) always has full column rank when degree is -1 or 0. When degree is 1, \(P(y)\) has full column rank if the data point locations are not all collinear (N=2), coplanar (N=3), etc.

• If kernel is ‘multiquadric’, ‘linear’, ‘thin_plate_spline’, ‘cubic’, or ‘quintic’, then degree must not be lower than the minimum value listed above.

• If smoothing is 0, then each data point location must be distinct.

当使用非尺度不变的 RBF(‘multiquadric’, ‘inverse_multiquadric’、‘inverse_quadratic’ 或 ‘gaussian’)时，必须选择适当的形状参数(例如，通过交叉验证)。较小的形状参数值对应于较宽的 RBF。当形状参数太小时，问题可能会变成ill-conditioned 或单数。

求解 RBF 插值系数所需的内存随着数据点的数量呈二次方增加，当插值超过大约一千个数据点时，这可能变得不切实际。为了克服大型插值问题的内存限制，可以指定近邻参数以仅使用最近的数据点为每个评估点计算 RBF 插值。

参考：

[1]

Fasshauer, G.，2007 年。使用 Matlab 的无网格近似方法。世界科学出版公司

[2]

http://amadeus.math.iit.edu/~fass/603_ch3.pdf

[3]

Wahba, G.，1990。观测数据的样条模型。暹。

[4]

http://pages.stat.wisc.edu/~wahba/stat860public/lect/lect8/lect8.pdf

例子：

演示将分散的数据插值到二维网格中。

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from scipy.interpolate import RBFInterpolator
>>> from scipy.stats.qmc import Halton

>>> rng = np.random.default_rng()
>>> xobs = 2*Halton(2, seed=rng).random(100) - 1
>>> yobs = np.sum(xobs, axis=1)*np.exp(-6*np.sum(xobs**2, axis=1))

>>> xgrid = np.mgrid[-1:1:50j, -1:1:50j]
>>> xflat = xgrid.reshape(2, -1).T
>>> yflat = RBFInterpolator(xobs, yobs)(xflat)
>>> ygrid = yflat.reshape(50, 50)

>>> fig, ax = plt.subplots()
>>> ax.pcolormesh(*xgrid, ygrid, vmin=-0.25, vmax=0.25, shading='gouraud')
>>> p = ax.scatter(*xobs.T, c=yobs, s=50, ec='k', vmin=-0.25, vmax=0.25)
>>> fig.colorbar(p)
>>> plt.show()