径向基函数拟合(RBF Fitting)

最近读了几篇论文都用到了径向基函数拟合(RBF Fitting)，感觉功能很强大，因此学习一下。

一、径向基函数

径向基函数跟高斯分布的概率密度函数类似，因此也叫高斯核函数，一般定义为：
φ(x)=e−(x−c)22σ2(1)\varphi(x)=e^{-\frac{(x-c)^2}{2\sigma^2}} \tag1φ(x)=e−2σ2(x−c)2(1)
对于给定的σ\sigmaσ, 径向基函数的取值仅仅跟 xxx 离中心点ccc的距离相关，当c=0c=0c=0时，图像如下所示:

二、径向基函数拟合

给定二维平面上的点集：P={(x1,y1),(x2,y2),...(xn,yn)}P =\{(x_1,y_1),(x_2,y_2),...(x_n,y_n)\}P={(x1,y1),(x2,y2),...(xn,yn)}, 我们期望拟合一个函数：
f(x)=∑i=1nwiφ(∣∣x−xi∣∣)(2)f(x) = \sum^n_{i=1}w_i\varphi(||x-x_i||)\tag2f(x)=i=1∑nwiφ(∣∣x−xi∣∣)(2)
也就是说，针对PPP中每一个点(xi,yi)(x_i,y_i)(xi,yi)，构造一个以xix_ixi为中心的径向基函数。给定任意横坐标 xxx, 其 yyy 值可以通过上式预测.

那么，如何计算这里的插值权重wiw_iwi呢？

针对拟合任务，首先要保证拟合的值跟真实值尽可能接近。因此，针对任意点 Pj(xj,yj)P_j(x_j,y_j)Pj(xj,yj), 我们期望完美拟合：

f(xj)=∑i=1nwiφ(∣∣xj−xi∣∣)=yjf(x_j) = \sum^n_{i=1}w_i\varphi(||x_j-x_i||)=y_jf(xj)=i=1∑nwiφ(∣∣xj−xi∣∣)=yj
将PPP中所有点代入公式(2), 并进一步写成矩阵乘法的形式：
[φ11φ12...φ1nφ21φ22...φ2n....φn1φn2...φnn][w1w2.wn]=[y1y2.yn](3)\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} y_1 \\ y_2 \\ . \\ y_n \end{bmatrix}\tag3 ⎣⎡φ11φ21.φn1φ12φ22.φn2..........φ1nφ2n.φnn⎦⎤⎣⎡w1w2.wn⎦⎤=⎣⎡y1y2.yn⎦⎤(3)
其中，φij=φ(∣∣xj−xi∣∣)=φji\varphi_{ij} = \varphi(||x_j-x_i||)=\varphi_{ji}φij=φ(∣∣xj−xi∣∣)=φji. 为此，权重参数可直接求解：
[w1w2.wn]=[φ11φ12...φ1nφ21φ22...φ2n....φn1φn2...φnn]−1[y1y2.yn](4)\begin{bmatrix} w_1 \\ w_2 \\ . \\ w_n \end{bmatrix}= \begin{bmatrix} \varphi_{11} & \varphi_{12}&... & \varphi_{1n}\\ \varphi_{21} & \varphi_{22}&... & \varphi_{2n}\\ . & .& . & .\\ \varphi_{n1} & \varphi_{n2}&... & \varphi_{nn} \end{bmatrix}^{-1} \begin{bmatrix} y_1 \\ y_2 \\ . \\ y_n \end{bmatrix}\tag4 ⎣⎡w1w2.wn⎦⎤=⎣⎡φ11φ21.φn1φ12φ22.φn2..........φ1nφ2n.φnn⎦⎤−1⎣⎡y1y2.yn⎦⎤(4)

三、实际案例

给定函数：
y=−x2+sin(3x)+20cos(2x)y = -x^2+sin(3x)+20cos(2x)y=−x2+sin(3x)+20cos(2x)
函数图像大致如下：

下面的圆圈是我们的均匀采样的20个顶点，我们期望通过这些采样点拟合出该函数。

拟合的结果如下所示，可以看出基本恢复了原始函数。

代码 (参考并修改：https://blog.csdn.net/xfijun/article/details/105670892)：

import numpy as np
import matplotlib.pyplot as plt# y = -x^2 + sin(3x)+20cos(2x)
def f(x):return -x*x + np.sin(3*x) + 20*np.cos(2*x)# y = e^(-(x-xc)^2/2)
def gaussian(x,xc):return np.exp(-(x-xc)**2/(2*(1**2)))# get w from y' = \sum w_i\phi |x'-xc|
def calculate_weights(x_sample,y_sample):num = len(y_sample)int_matrix = np.asmatrix(np.zeros((num,num)))for i in range(num):int_matrix[i,:] = gaussian(x_sample, x_sample[i])w = int_matrix.I * np.asmatrix(y_sample).T      return w# y' = \sum w_i\phi |x'-xc|
def calculate_rbf_value(w,x_sample,x):num = len(x)y_= np.zeros(num)for i in range(num):for k in range(len(w)):y_[i] = y_[i] + w[k]*gaussian(x[i],x_sample[k])return y_def draw_kernels(x_samples,w):for i in range(len(x_samples)):x = x_samples[i]x_ = np.linspace(x-2,x+2,100)y_ = []for a in x_:y_.append(w[i]*gaussian(a,x))y_ = np.array(y_).reshape(100,-1)plt.plot(x_,y_,'m:')if __name__ == '__main__':sample_cnt = 20x_start, x_end = -8, 8x = np.linspace(x_start,x_end,500)y = f(x)x_sample = np.linspace(x_start,x_end,sample_cnt)y_sample = f(x_sample)w = calculate_weights(x_sample,y_sample)y_fit = calculate_rbf_value(w,x_sample,x)plt.figure(1)plt.plot(x,y_fit,'k')plt.plot(x,y,'r:')plt.ylabel('y')plt.xlabel('x')draw_kernels(x_sample,w)for i in range(len(x_sample)):plt.plot(x_sample[i],y_sample[i],'go',markerfacecolor='none')plt.legend(labels=['fitted','original','sample'],loc='lower left')plt.title('kernel interpolation:$y=sin(\pi x/2)+cos(\pi x/3)$')       plt.show()

下图绘制了这20个带权高斯核函数(径向基函数) wiφ(∣∣xj−xi∣∣)w_i\varphi(||x_j-x_i||)wiφ(∣∣xj−xi∣∣)的图像，需要注意这里的权值可正可负且∑wi≠1\sum w_i \neq 1∑wi=1. 这跟高斯混合模型，广义重心坐标插值(MVC)等插值方式的显著区别。

y=−x2y = -x^2y=−x2 的拟合结果：

y=xy = xy=x 的拟合结果：

为什么拟合能力这么强呢？其实从某种角度可以认为是高斯混合模型吧.

参考博客：https://blog.csdn.net/xfijun/article/details/105670892