RBF（Radial-Basis Function）网路之一：多元插值法_综合

1 插值：在离散数学的基础上补插连续函数，使得这条连续曲线通过全部给定的离散数据点。（来自百度百科）插值除了给定的离散数据点外，连续曲线的其他点都属于“插入”进来的，所以称插值。

在多维空间的插值法定义为：

给定N个不同的点集 $x\epsilon R^{m_{0}}|i=1,2,...N$ 和一组相对应的N个实数集合 $d\epsilon R^{1}|i=1,2,...N$ ，找到一个函数满足如下的条件：

$F(x_{i})=d_{i}$ (1) 插值和拟合的区别：

(1)共同点：两者都是寻求函数曲线（曲线）

(2)不同点：插值要求曲线（或曲面）过全部的离散数据点，而拟合是寻求全局最接近的结果，即要求拟合函数和数据点的最小二乘意义最小。

2插值的常用方法

（1）牛顿法

（2）逐次线性插值

（3）拉格朗日插值法

（4）等距节点插值

（5）Herminte法

（6）分段低次插值

（7）样条插值

本篇主要介绍的rbf插值法

3 rbf插值法

rbf插值法是由Powell 于1985年提出。其基本思想不是直接使用原数据点进行插值，而是使用每个点与中心点的距离，距离一般是欧几里得距离，也就是每个点离中心点的径向距离如下表示：

$||\mathbf{x}-\mathbf{x}_{i}||)$

其中已知的插值节点 $\mathbf{x}_{i}$ 是中心点， $\mathbf{x}_{i}\epsilon R^{m_{0}}|i=1,2,...N$ ， $\mathbf{x}\epsilon R^{m_{0}}$ 是等待插值的节点， $m_{0}$ 是插值节点的维度。而RBF插值的基函数具有如下形式：

$\varphi (||\mathbf{x}-\mathbf{x}_{i}||)| i=1,2,...,N$ ,

此函数称为径向基函数，通常是非线性的，此函数常用的形式见下文。而RBF的插值函数有如下形式：

$F(\mathbf{x})=\sum_{i=1}^{N}w_{i}\varphi (||\mathbf{x}-\mathbf{x}_{i}||)$

其中 $w_{i}$ 是插值矩阵。由于 $F(x_{i})=d_{i}$ 我们可以得到如下变换式：

$\begin{bmatrix} \varphi _{11} & \varphi _{12} &... & \varphi _{1N}\\ \varphi _{21}& \varphi _{22} &... &\varphi _{2N} \\ ... & ... & ... & ...\\ \varphi _{N1}& \varphi _{N2} &... &\varphi _{NN} \end{bmatrix}\begin{bmatrix} w_{1}\\ w_{2}\\ ...\\ w_{N} \end{bmatrix}= \begin{bmatrix} d_{1}\\ d_{2}\\ ...\\ d_{N} \end{bmatrix}$

其中记 $\varphi _{i,j}=\varphi (||\mathbf{x}_{i}-\mathbf{x}_{j}||)|i,j=1,2,...N$

便于表示，我们有如下记号：

记 $\mathbf{w}=\left [ w_{1} ,w_{1} ,... ,w_{1N} \right ]^{T}$

$\mathbf{d}=\left [ d_{1} ,d_{1} ,... ,d_{1N} \right ]^{T}$

记 $\mathbf{\phi }=\left \{ \varphi _{i,j}\right \}_{i=1}^{N}$ (左边是粗体)，

其中 $\mathbf{d,w }$ 均是N行1列的向量，称 $\mathbf{d}$ 为要求的结果向量，称 $\mathbf{w}$ 为线性权重向量，而 $\phi$ 称为插值矩阵，它是N行N列，且是对称矩阵

由此，上述公式可简化为：

$\mathbf{\phi }\mathbf{w}=\mathbf{d}$

如果插值矩阵 $\phi$ 是非奇异的，则其逆矩阵 $\phi ^{-1}$ 存在，此时可求得线性权重向量

$\mathbf{w}=\phi ^{-1}\mathbf{d}$

插值向量是不是非奇异的，可以由Micchelli理论来解释。

Micchelli定理：

$\mathbf{x}_{i}\epsilon R^{m_{0}}|i=1,2,...N$ 是实数域 $R^{m_{0}}$ 里互不相同的节点，则其插值矩阵 $\mathbf{\phi }=\left \{ \varphi _{i,j}\right \}_{i=1}^{N}$ 是非奇异的。

满足Micchelli定理的径向基函数 $\varphi$ 有如下：

(1)多二次函数

$\varphi(r)=(r^{2}+c^{2})^{1/2},c>0,r\epsilon R$

(2)逆多二次函数

$\varphi(r)=\frac{1}{(r^{2}+c^{2})^{1/2}} ,c>0,r\epsilon R$

(3)高斯函数

$\varphi(r)=exp(-\frac{r^{2}}{2\sigma ^{2}}) ,\sigma >0,r\epsilon R$

对于上述三种函数，如果插值节点是互不相同的，无论输入样本数N是多少，也无论输入样本的维度 $m_{0}$ 是多少，它的插值矩阵 $\phi$ 总是非奇异的。注意：此方法介绍的插值法是完全内插法，要求插值曲线过每一个插值节点。对于逆多二次函数和高斯函数，当 $r\rightarrow \infty$ 时， $\varphi(r)\rightarrow 0$ ,插值矩阵是正定的，这两种函数被称为本地化函数，而多二次函数则 $r\rightarrow \infty$ 时， $\varphi(r)\rightarrow \infty$ ，被称为非本地化函数，插值矩阵具有N-1个负的特征值和1个正的特征值，而且它还不是正定的

插值节点的计算公式如下：

4实现

（1）使用RBF进行一维插值

#rbf插值
import numpy as np
import matplotlib.pyplot as pltdef gen_data(x1,x2):y_sample = np.sin(np.pi*x1/2)+np.cos(np.pi*x1/2)y_all = np.sin(np.pi*x2/2)+np.cos(np.pi*x2/2)return y_sample, y_alldef kernel_interpolation(y_sample,x1,sig):gaussian_kernel = lambda x,c,h: np.exp(-(x-c)**2/(2*(h**2)))#使用高斯核函数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,x1[i],sig)  w = int_matrix.I * np.asmatrix(y_sample).T  #计算权重矩阵    return wdef kernel_interpolation_rec(w,x1,x2,sig):#根据权重矩阵，进行插值gkernel = lambda x,xc,h: np.exp(-(x-xc)**2/(2*(h**2)))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_recif __name__ == '__main__':snum = 20   # 插值节点数量ratio = 20  # 总数据点数量：snum*ratiosig = 1		# 核函数宽度xs = -8xe = 8x1 = 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)plt.rcParams['font.sans-serif']=['SimHei'] #显示中文标签plt.rcParams['axes.unicode_minus']=False   #这两行需要手动设置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=['插值曲线','原函数曲线','插值节点'],loc='lower left')plt.title('kernel interpolation:$y=sin(\pi x/2)+cos(\pi x/2)$')       plt.show()MSE = np.linalg.norm(y_all-y_rec, ord=2)**2/len(y_all)

MSE误差为0.00013145540771572465

(2)使用RBF进行二维插值

#rbf二维插值
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
#定义二维函数
def gen_data(x1,x2):y=20+x1**2-10*np.cos(2*np.pi*x1)+x2**2-10*np.cos(2*np.pi*x2) #插值函数return y
def kernel_interpolation(y_sample,x1,x2,sig):gaussian_kernel = lambda x1,x2,y1,y2,h:np.exp((-(x1-x2)**2-(y1-y2)**2)/(2*(h**2)))#计算二维空间欧式距离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,x1[i],x2,x2[i],sig)  w = int_matrix.I * np.asmatrix(y_sample).T    #权重矩阵  return wdef kernel_interpolation_rec(w,x1,x2,y1,y2,sig):#进行插值gkernel = lambda x1,y1,x2,y2,h: np.exp((-(x1-y1)**2-(x2-y2)**2)/(2*(h**2)))num = len(y1)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(y1[i],x1[k],y2[i],x2[k],sig)return y_rec
if __name__ == '__main__':snum = 20  # 插值节点数量ratio = 20  # 总数据点数量：snum*ratiosig = 1		# 核函数宽度xs = -8xe = 8x1 = np.linspace(xs,xe,snum)x2=x1y1 = np.linspace(xs,xe,(snum-1)*ratio+1)y2=y1y_sample = gen_data(x1,x2)y_all=gen_data(y1,y2)w = kernel_interpolation(y_sample,x1,x2,sig) y_rec = kernel_interpolation_rec(w,x1,x2,y1,y2,sig)MSE = np.linalg.norm(y_all-y_rec, ord=2)**2/len(y_all)print(MSE)

参考资料：

1 Nerual networks and Learning Machines Third Edition ,.神经网络与机器学习 Simon Haykin第三版，第五章。

2 基于径向基函数(RBF)的函数插值.https://blog.csdn.net/xfijun/article/details/105670892

3 数值分析第五版李庆扬、王能超、易大义编