径向基函数拟合(RBF Fitting)

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

一、径向基函数

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

二、径向基函数拟合

给定二维平面上的点集：$P=\{(x_1,y_1),(x_2,y_2),...(x_n,y_n)\}$, 我们期望拟合一个函数：
$$f(x)=\sum _{i=1}^n{w}_i\phi (||x-x_i||)\text{\hspace{1em}(2)}$$
也就是说，针对$P$中每一个点$(x_i,y_i)$，构造一个以$x_i$为中心的径向基函数。给定任意横坐标 $x$, 其 $y$ 值可以通过上式预测.

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

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

$$f(x_j)=\sum _{i=1}^n{w}_i\phi (||x_j-x_i||)=y_j$$
其中，$\phi (||x_j-x_i||)=e^{-\frac{(x_j-x_i{)}^2}{2{\sigma }^2}}$。将$P$中所有点代入公式(2), 并进一步写成矩阵乘法的形式：
$$\left[\begin{array}{cccc}{\phi }_{11}& {\phi }_{12}& ...& {\phi }_{1n}\\ {\phi }_{21}& {\phi }_{22}& ...& {\phi }_{2n}\\ .& .& .& .\\ {\phi }_{n1}& {\phi }_{n2}& ...& {\phi }_{nn}\end{array}\right]\left[\begin{array}{c}{w}_1\\ w_2\\ .\\ w_n\end{array}\right]=\left[\begin{array}{c}{y}_1\\ y_2\\ .\\ y_n\end{array}\right]\text{\hspace{1em}(3)}$$
其中，${\phi }_{ij}=\phi (||x_j-x_i||)={\phi }_{ji}$. 为此，权重参数可直接求解：
$$\left[\begin{array}{c}{w}_1\\ w_2\\ .\\ w_n\end{array}\right]={\left[\begin{array}{cccc}{\phi }_{11}& {\phi }_{12}& ...& {\phi }_{1n}\\ {\phi }_{21}& {\phi }_{22}& ...& {\phi }_{2n}\\ .& .& .& .\\ {\phi }_{n1}& {\phi }_{n2}& ...& {\phi }_{nn}\end{array}\right]}^{-1}\left[\begin{array}{c}{y}_1\\ y_2\\ .\\ y_n\end{array}\right]\text{\hspace{1em}(4)}$$

三、实际案例

给定函数：
$$y=-x^2+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 = 20
    x_start, x_end = -8, 8

    x = 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个带权高斯核函数(径向基函数) $w_i\phi (||x_j-x_i||)$的图像，需要注意这里的权值可正可负且$\sum w_i\ne 1$. 这跟高斯混合模型，广义重心坐标插值(MVC)等插值方式的显著区别。

$y=-x^2$ 的拟合结果：

$y=x$ 的拟合结果：

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

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