基于径向基函数(RBF)的函数插值

1. 函数插值

函数插值问题：　用形式简单的插值函数 $\hat{f}(x)$ 近似原函数

$(1)$ 设函数 $y=f(x)$ 在某个区间上有定义，并且已知该区间上的一些数据点 { x i , y i } \{x_i,y_i\} {
xi​,yi​} 严格满足 $y_i=f(x_i),i=1,\cdots ,N$，这些数据点称为“控制节点”或“插值节点”

$(2)$ 如果存在一个形式上比较简单（比如 $n$ 次多项式）的函数 $\hat{f}(x)$，使得 $\hat{f}(x_i)=y_i,i=1,\cdots ,N$ 都成立，就称 $\hat{f}(x)$ 为 $f(x)$ 的插值函数。

\uad 典型的函数插值方法：拉格朗日插值和牛顿插值、$Hermite$插值、样条插值等。

2. RBF函数插值

\uad 与拉格朗日插值之类的常规函数插值不同，基于核函数的函数插值“通过引入核函数”来刻画数据的局部化特征。

\uad 径向基函数 $\text{(Radial Basis Function,RBF)}$ 就是一类特殊的基函数，最常用的就是“高斯基函数”，定义为：

$\phi (x)=e^{-\frac{x^2}{2{\sigma }^2}}$　　(以一维情况为例)

\uad
RBF函数插值：　$\hat{f}(x)=\sum _{i=1}^N{w}_i\phi (\parallel x-x_i\parallel )$

\uad 假设有 $N$ 个插值节点，也就是已知 { x j , y j } ∣ j = 1 N \{x_j,y_j\}\big |_{j=1}^N {
xj​,yj​}∣
∣​j=1N​，其中 $\hat{f}(x_j)=y_j=f(x_j)$，如下图所示。

\uad

图中，红色实线为真实函数曲线，绿色空心圆圈代表插值节点 $(x_j,y_j)$，蓝色实心点为RBF插值所求得的权值 $w_j$

\uad 将 { x j , y j } ∣ j = 1 N \{x_j,y_j\}\big |_{j=1}^N {
xj​,yj​}∣
∣​j=1N​ 带入方程 $\hat{f}(x)=\sum _{i=1}^N{w}_i\phi (\parallel x-x_i\parallel )$，可得到：

$\underset{\Phi }{\left[\begin{array}{cccc}{\phi }_{11}& {\phi }_{12}& \cdots & {\phi }_{1N}\\ {\phi }_{21}& {\phi }_{22}& \cdots & {\phi }_{2N}\\ ⋮& ⋮& & ⋮\\ {\phi }_{11}& {\phi }_{12}& \cdots & {\phi }_{1N}\end{array}\right]}\underset{\mathbf{W}}{\left[\begin{array}{c}{w}_1\\ w_2\\ ⋮\\ w_N\end{array}\right]}=\underset{\mathbf{y}}{\left[\begin{array}{c}{y}_1\\ y_2\\ ⋮\\ y_N\end{array}\right]}$，其中 ${\phi }_{ji}=\phi (\parallel x_j-x_i\parallel )$

\uad
\uad 其中，$\Phi =[{\phi }_{ji}]$ 为插值矩阵。因为 ${\phi }_{ji}=\phi (\parallel x_j-x_i\parallel )={\phi }_{ij}$，因此插值矩阵是对称的。对于高斯核函数而言，插值矩阵的对角线元素的值为 $1$。

\uad 将线性方程组记为 $\Phi \mathbf{W}=\mathbf{y}$，该方程组的第 $j$ 行为：

$\hat{f}(x_j)=y_j=w_1\phi (\parallel x_j-x_1\parallel )+w_2\phi (\parallel x_j-x_2\parallel )+\cdots +w_N\phi (\parallel x_j-x_N\parallel )$

\uad 因此，可求出 $\text{RBF}$ 插值的系数为：$\mathbf{W}={\Phi }^{-1}\mathbf{y}$，其示意图如下图所示。

Micchelli定理可以保证采用高斯函数时，插值矩阵 $\Phi$ 是可逆的（只要插值节点互不相同）。

\uad

代码实现

import numpy as np import matplotlib.pyplot as plt def gen_data(x1,x2): y_sample = np.sin(np.pi*x1/2)+np.cos(np.pi*x1/3) y_all = np.sin(np.pi*x2/2)+np.cos(np.pi*x2/3) return y_sample, y_all def kernel_interpolation(y_sample,x1,sig): gaussian_kernel = lambda x,c,h: np.exp(-(x-x[c])2/(2*(h2))) num = len(y_sample) w = np.zeros(num) int_matrix = np.asmatrix(np.zeros((num,num))) for i in range(num): int_matrix[i,:] = gaussian_kernel(x1,i,sig) w = int_matrix.I * np.asmatrix(y_sample).T return w def kernel_interpolation_rec(w,x1,x2,sig): gkernel = lambda x,xc,h: np.exp(-(x-xc)2/(2*(h2))) num = len(x2) y_rec = np.zeros(num) for i in range(num): for k in range(len(w)): y_rec[i] = y_rec[i] + w[k]*gkernel(x2[i],x1[k],sig) return y_rec if __name__ == '__main__': snum = 20 # control point数量 ratio = 20 # 总数据点数量：snum*ratio sig = 1 # 核函数宽度 xs = -8 xe = 8 x1 = np.linspace(xs,xe,snum) x2 = np.linspace(xs,xe,(snum-1)*ratio+1) y_sample, y_all = gen_data(x1,x2) plt.figure(1) w = kernel_interpolation(y_sample,x1,sig) y_rec = kernel_interpolation_rec(w,x1,x2,sig) plt.plot(x2,y_rec,'k') plt.plot(x2,y_all,'r:') plt.ylabel('y') plt.xlabel('x') for i in range(len(x1)): plt.plot(x1[i],y_sample[i],'go',markerfacecolor='none') plt.legend(labels=['reconstruction','original','control point'],loc='lower left') plt.title('kernel interpolation:$y=sin(\pi x/2)+cos(\pi x/3)$') plt.show() 

在相同区间、分别采用 $8,12,16,20$ 个控制节点 $\text{(control point)}$ 进行函数插值的结果
显然，插值节点过少，无法体现整个函数的特征；插值节点越多，函数插值的结果越精确

扩大插值区间范围，控制节点 $\text{(control point)}$ 也需要增加数量，才能保持函数插值的准确性

\quad
另外，$\text{Scipy}$ 的插值模块也提供了 $\text{RBF}$ 插值，其实现代码如下：

import numpy as np import matplotlib.pyplot as plt from scipy.interpolate import Rbf f = lambda x: np.sin(np.pi*x/2)+np.cos(np.pi*x/3) snum = 20 # control point数量 ratio = 20 # 总数据点数量：snum*ratio xs = -8 xe = 8 x1 = np.linspace(xs,xe,snum) # control points x2 = np.linspace(xs,xe,(snum-1)*ratio+1) # 作图总数据点 y1 = f(x1) # control points rf = Rbf(x1, y1) # reconstructed Rbf function y2 = rf(x2) # Rbf reconstruction plt.plot(x2, y2, 'k-', x2, f(x2),'r-', x1, y1, 'go', markerfacecolor='none') plt.legend(["Radial basis functions", "Orignal", "control point"],loc='best') plt.show() 

\quad
此外，$\text{RBF}$函数插值还可以通过径向基函数网络来实现。

\quad
参考：函数插值的python实现——拉格朗日、牛顿插值