径向基函数拟合

最近读了几篇论文都用到了径向基函数拟合(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), 并进一步写成矩阵乘法的形式：
[ φ 11 φ 12 . . . φ 1 n φ 21 φ 22 . . . φ 2 n . . . . φ n 1 φ n 2 . . . φ n n ] [ w 1 w 2 . w n ] = [ y 1 y 2 . y n ] (3)

= \tag3 ​φ11​φ21​.φn1​​φ12​φ22​.φn2​​..........​φ1n​φ2n​.φnn​​ ​ ​w1​w2​.wn​​ ​= ​y1​y2​.yn​​ ​(3)
其中，${\phi }_{ij}=\phi (||x_j-x_i||)={\phi }_{ji}$. 为此，权重参数可直接求解：
[ w 1 w 2 . w n ] = [ φ 11 φ 12 . . . φ 1 n φ 21 φ 22 . . . φ 2 n . . . . φ n 1 φ n 2 . . . φ n n ] − 1 [ y 1 y 2 . y n ] (4)

[w1w2.wn]" role="presentation" style="position: relative;">⎡⎣⎢⎢⎢w1w2.wn⎤⎦⎥⎥⎥$$\left[\begin{array}{c}{w}_1\\ w_2\\ .\\ w_n\end{array}\right]$$

=

[φ11φ12...φ1nφ21φ22...φ2n....φn1φn2...φnn]" role="presentation" style="position: relative;">⎡⎣⎢⎢⎢φ11φ21.φn1φ12φ22.φn2..........φ1nφ2n.φnn⎤⎦⎥⎥⎥$$\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}

[y1y2.yn]" role="presentation" style="position: relative;">⎡⎣⎢⎢⎢y1y2.yn⎤⎦⎥⎥⎥$$\left[\begin{array}{c}{y}_1\\ y_2\\ .\\ y_n\end{array}\right]$$

\tag4 ​w1​w2​.wn​​ ​= ​φ11​φ21​.φn1​​φ12​φ22​.φn2​​..........​φ1n​φ2n​.φnn​​ ​−1 ​y1​y2​.yn​​ ​(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()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63

下图绘制了这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