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

zfoox

已于 2022-10-12 15:49:44 修改

阅读量2.3w

点赞数 67

分类专栏： CV 文章标签：算法机器学习人工智能

于 2020-07-07 01:16:19 首次发布

本文链接：https://blog.csdn.net/xfijun/article/details/105670892

版权

CV 专栏收录该内容

31 篇文章 12 订阅

订阅专栏

基于径向基函数的函数插值

1. 函数插值
2. RBF函数插值
代码实现

1. 函数插值

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

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

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

$\qquad$ 典型的函数插值方法：拉格朗日插值和牛顿插值、 $Her mi t e$ 插值、样条插值等。

与“函数逼近”的主要区别：
插值函数 $\hat f(x)$ 必须经过“插值节点”，也就是要满足 $\hat f(x_i)=y_i,i=1,\cdots,N$

2. RBF函数插值

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

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

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

$\qquad$ 在这里插入图片描述
RBF函数插值：　 $\hat{f}(x)=\displaystyle\sum_{i=1}^Nw_i\varphi(\parallel x-x_i\parallel)$

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

$\qquad$ 在这里插入图片描述

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

$\qquad$ 将 $\{x_j,y_j\}\big |_{j=1}^N$ 带入方程 $\hat{f}(x)=\displaystyle\sum_{i=1}^Nw_i\varphi(\parallel x-x_i\parallel)$ ，可得到：

$\qquad\qquad\underbrace{\left[ \begin{matrix} \varphi_{11} & \varphi_{12} & \cdots & \varphi_{1N} \\ \varphi_{21} & \varphi_{22} & \cdots & \varphi_{2N} \\ \vdots & \vdots & &\vdots \\ \varphi_{11} & \varphi_{12} & \cdots & \varphi_{1N} \end{matrix} \right] }_{\Phi}\underbrace{ \left[ \begin{matrix} w_1 \\ w_2 \\ \vdots \\ w_N \end{matrix} \right]}_{\bold W}=\underbrace{\left[ \begin{matrix} y_1 \\ y_2 \\ \vdots \\ y_N \end{matrix} \right]}_{\bold y}$ ，其中 $\varphi_{ji}=\varphi(\parallel x_j-x_i\parallel)$

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

$\qquad$ 将线性方程组记为 $\Phi\bold W=\bold y$ ，该方程组的第 $j$ 行为：

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

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

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

$\qquad$ 在这里插入图片描述

代码实现

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*(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,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*(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_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()