径向基函数插值

一、径向基函数的定义

如果 $||x_1||=||x_2||$，那么 $\varphi (x_1)=\varphi (x_2)$ 的函数 $\varphi$ 就是径向函数，即仅由 $r=||x||$ 控制的函数（径向基函数是一个取值仅仅依赖于离原点距离的实值函数，或者还可以是到任意一点 $c$ 的距离）。

给定一个一元函数 $\varphi ：R_{+}\to R$，在定义域 $x\in R^d$ 上，所有形如 $\psi (x)=\varphi (||x-c||)$ 及其线性组合张成的函数空间称为由函数 $\varphi$ 导出的径向基函数空间。

在一定的条件下，只要取 { x j } \{x_j\} {xj​} 两两不同， { ϕ ( x − x j ) } \{\phi(x-x_j)\} {ϕ(x−xj​)} 就是线性无关的，从而形成径向基函数空间中某子空间的一组基。当 { x j } \{x_j\} {xj​} 几乎充满 $R$ 时， { x j } \{x_j\} {xj​} 几乎充满 $R$ 时， { ϕ ( x − x j ) } \{\phi(x-x_j)\} {ϕ(x−xj​)} 及其线性组合可以逼近几乎任何函数。

各类文献中常用的径向基函数有：

1. Kriging 方法的 Gauss 分布函数：$\varphi (r)=e^{-c^2{r}^2}$
2. Kriging 方法的 Markoff 分布函数：$\varphi (r)=e^{-c|r|}$，及其他概率分布函数；
3. Hardy 的 Multi-Quadric 函数：$\varphi (r)=(c^2+r^2{)}^{\beta }$（其中 $\beta$ 是正的实数）；
4. Hardy 的逆 Multi-Quadric 函数：$\varphi (r)=(c^2+r^2{)}^{-\beta }$（其中 $\beta$ 是正的实数）；
5. Duchon 的薄板样条：$d$ 为偶数时，$\varphi (r)=r^{2k-d}\log r$；$d$ 为奇数时，$\varphi (r)=r^{2k-d}$；

二、径向基函数插值

定义：径向基函数插值是对于给定的多元散乱数据 { x j , f j } j = 1 n , x j ∈ R n , f j ∈ R , j = 1 , ⋯ , n \{x_j,f_j\}^n_{j=1},x_j\in R^n,f_j\in R,j=1,\cdots,n {xj​,fj​}j=1n​,xj​∈Rn,fj​∈R,j=1,⋯,n。选取径向函数 $\varphi :R_{+}\to R$ 来构造函数系 { ϕ ( ∣ ∣ x − x j ∣ ∣ ) } j = 1 n \{\phi(||x-x_j||)\}_{j=1}^n {ϕ(∣∣x−xj​∣∣)}j=1n​ 并寻找形如 $S(x)={\sum }_{j=1}^n{\lambda }_j\varphi (||x-x_j||)$ 的插值函数 $S(x)$，使其满足条件 $S(x_j)=f_j,j=1,\cdots ,n$。

为了方便，我们定义
$$\{\begin{array}{l}{f}^T=(f_1,f_2,\cdots ,f_n)\\ {\varphi }^T(x)=(\varphi (||x-x_1||,\varphi (||x-x_2||,\cdots ,\varphi (||x-x_n||))\\ {\lambda }^T=({\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_n)\\ A=(\varphi (||x_j-x_k||){)}_{n\times n}\end{array}$$

上述插值方程对任意两两不同的 $x_j$ 的数据 { x j , f j } \{x_j,f_j\} {xj​,fj​} 有解的充要条件是：对任意两两不同的 $x_j$，对称矩阵 $A$ 都非奇异。

定理：函数 $\varphi :R_{+}\to R$ 是连续的，${\lim }_{r\to \infty }\varphi (r)=0$，那么对于 $n$ 元的径向基函数插值总是存在唯一解的充分条件是：矩阵 $A$ 是正定矩阵。

上面提到的径向基函数中逆 Multi-Quadric 函数和 Gauss 函数在任意维空间上都是正定函数，因此插值是唯一的。

三、用高斯函数进行散乱数据的插值

对于数据量少的情况，径向基函数（尤其是高斯函数）插值的结果较令人满意，而且计算也比较简单。

令径向基函数插值方程为
$$S(x)=\sum _{j=1}^n{\lambda }_j\varphi (||x-x_j||)$$

将已知点 $(x_j,f_j),j=1,\cdots ,n$ 代入方程，可得：
$$\left[\begin{array}{cccc}{\lambda }_1& {\lambda }_2& \cdots & {\lambda }_n\end{array}\right]\left[\begin{array}{cccc}\varphi (||x_1-x_1||)& \varphi (||x_2-x_1||)& \cdots & \varphi (||x_n-x_1||)\\ \varphi (||x_1-x_2||)& \varphi (||x_2-x_2||)& \cdots & \varphi (||x_n-x_2||)\\ ⋮& ⋮& & ⋮\\ \varphi (||x_1-x_n||)& \varphi (||x_2-x_n||)& \cdots & \varphi (||x_n-x_n||)\end{array}\right]=\left[\begin{array}{cccc}{f}_1& f_2& \cdots & f_n\end{array}\right]$$

求解上述方程，可求出 ${\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_n$ 的值，从而求出插值曲线方程。插值曲面方程类似，将 $x$ 替换成向量 $X$ 即可。

具体应用到高斯函数，设高斯函数插值方程为
$$S(x)=\sum _{j=1}^n{\lambda }_j{e}^{-\alpha ||x-x_j|{|}^2}\alpha >0$$

其中，$\alpha$ 为形状调整参数，可根据散乱数据点分布特征选取，当数据点对应的函数值变化较大时，$\alpha$ 可取的稍大些；数据点对应的函数值变化较小时，$\alpha$ 可取得稍小些。

python 代码实现

import numpy as np
import matplotlib.pyplot as pltdef gen_data(x1, x2):# 用于生成插值节点和总数据点 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_alldef kernel_interpolation(y_sample, x1, sig):# 求解插值函数中高斯基函数的系数gaussian_kernel = lambda x, c, h: np.exp(-h*(x - x[c]) ** 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 = np.asmatrix(y_sample) * int_matrix.Iw = w.Treturn wdef kernel_interpolation_rec(w, x1, x2, sig):gkernel = lambda x, xc, h: np.exp(-h*(x - xc) ** 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  # control point数量ratio = 20  # 总数据点数量：snum*ratiosig = 0.5  # 核函数宽度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)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()

运行结果

Matlab 代码实现