径向基函数（RBF）插值

RBF函数插值

径向基函数（Radial Basis Function, RBF）插值的基本形式为
F(r)=∑i=1Nwiφ(∥r−ri∥)F(\boldsymbol{r})=\sum_{i=1}^N w_i \varphi\left(\left\|\boldsymbol{r}-\boldsymbol{r}_i\right\|\right) F(r)=i=1∑Nwiφ(∥r−ri∥)
式中， F(r)F(\boldsymbol{r})F(r) 是插值函数， NNN 为插值问题所使用的径向基函数总数目（控制点总数目）， φ(∥r−ri∥)\varphi\left(\left\|\boldsymbol{r}-\boldsymbol{r}_i\right\|\right)φ(∥r−ri∥) 是采用的径向基函数的通用形式， ∥r−ri∥\left\|\boldsymbol{r}-\boldsymbol{r}_i\right\|∥r−ri∥ 是两个位置矢量的欧氏距离， ri\boldsymbol{r}_iri 是第 iii 号径向基函数的控制点位置， wiw_iwi 是第 iii 号径向基函数对应的权重系数[2]^{[2]}[2]。

径向基函数类型很多，总结有如下六种：

Gaussian（高斯曲面函数）：φ(x)=exp⁡(−x22σ2)\quad \varphi(x)=\exp \left(-\frac{x^2}{2 \sigma^2}\right)φ(x)=exp(−2σ2x2)
Multiquadrics（多项式函数）：φ(x)=1+x2σ2\quad \varphi(x)=\sqrt{1+\frac{x^2}{\sigma^2}}φ(x)=1+σ2x2
Linear（线性函数）：φ(x)=x\varphi(x)=xφ(x)=x
Cubic（立方体曲面函数）：φ(x)=x3\varphi(x)=x^3φ(x)=x3
Thinplate（薄板曲面函数）：φ(x)=x2ln⁡(x+1)\quad \varphi(x)=x^2 \ln (x+1)φ(x)=x2ln(x+1)
Wendland’s C2C^2C2 函数（网格变形常用）：φ(x)=(1−x)4(4x+1)\varphi(x)=(1-x)^4(4 x+1)φ(x)=(1−x)4(4x+1)

现有一系列控制（插值）节点 {rj,F(r)j}∣j=1N\left.\left\{r_j, F(\boldsymbol{r})_j\right\}\right|_{j=1} ^N{rj,F(r)j}∣j=1N（插值函数 F(r)F(\boldsymbol{r})F(r) 必须经过控制节点），将其带入方程 F(r)=∑i=1Nwiφ(∥r−ri∥)F(\boldsymbol{r})=\sum_{i=1}^N w_i \varphi\left(\left\|\boldsymbol{r}-\boldsymbol{r}_i\right\|\right)F(r)=∑i=1Nwiφ(∥r−ri∥) ，可得到:
[φ11φ12⋯φ1Nφ21φ22⋯φ2N⋮⋮⋮φ11φ12⋯φ1N]⏟Φ[w1w2⋮wN]⏟W=[y1y2⋮yN]⏟y,其中φji=φ(∥rj−ri∥)\underbrace{\left[\begin{array}{cccc}\varphi_{11} & \varphi_{12} & \cdots & \varphi_{1 N} \\ \varphi_{21} & \varphi_{22} & \cdots & \varphi_{2 N} \\ \vdots & \vdots & & \vdots \\ \varphi_{11} & \varphi_{12} & \cdots & \varphi_{1 N}\end{array}\right]}_{\Phi} \underbrace{\left[\begin{array}{c}w_1 \\ w_2 \\ \vdots \\ w_N\end{array}\right]}_{\mathbf{W}}=\underbrace{\left[\begin{array}{c}y_1 \\ y_2 \\ \vdots \\ y_N\end{array}\right]}_{\mathbf{y}} ,其中 \varphi_{j i}=\varphi\left(\left\|r_j-r_i\right\|\right) Φ⎣⎡φ11φ21⋮φ11φ12φ22⋮φ12⋯⋯⋯φ1Nφ2N⋮φ1N⎦⎤W⎣⎡w1w2⋮wN⎦⎤=y⎣⎡y1y2⋮yN⎦⎤,其中φji=φ(∥rj−ri∥)
上式中，Φ=[φji]\Phi=\left[\varphi_{j i}\right]Φ=[φji] 为插值矩阵。将线性方程组记为 ΦW=y\Phi \mathbf{W}=\mathbf{y}ΦW=y ，可求出 RBF\mathrm{RBF}RBF 插值的权重系数为: W=Φ−1y\mathbf{W}=\Phi^{-1} \mathbf{y}W=Φ−1y 。在得到每个控制点的权重系数后，就能求出定义域内任意插值点所对应的 F(r)F(\boldsymbol{r})F(r) ，实现插值的功能。

代码测试

文章[1]^{[1]}[1]中代码主要实现的是一维下高斯基函数插值，文章[2]^{[2]}[2]中代码实现的主要是多维下的Wendland’s C2C^2C2基函数插值。借鉴上面篇代码，稍微修改精简一下。

RBF.py

from scipy.linalg import solve
import numpy as npdef rbf_coefficient(support_points, support_values, function_name = 'C2', radius = None):"""计算并返回径向基(radical basis function, RBF)插值函数的插值系数:param support_points: 径向基插值的支撑点:param support_values: 支撑点上的物理量，如位移、压力等:param function_name: 使用的径向基函数，默认为 Wendland C2 函数:param radius: 径向基函数的作用半径，默认作用范围包含所有支撑点:return: coefficient_mat, 径向基函数的插值系数矩阵"""num_support_points, dim = np.shape(support_points)phi_mat = np.zeros((num_support_points, num_support_points), dtype = np.float)for i in range(num_support_points):for j in range(num_support_points):eta = np.linalg.norm(support_points[i] - support_points[j])if radius is not None:eta = eta / radiusif eta > 1:continueif function_name == 'C2':phi_mat[i, j] = (1 - eta) ** 4 * (4 * eta + 1)elif function_name == 'cubic':phi_mat[i, j] = eta ** 3elif function_name == 'linear':phi_mat[i, j] = etaelif function_name == 'gaussian':sig = 1  # 高斯核宽度phi_mat[i, j] = np.exp(-1 * eta ** 2 / (2 * (sig ** 2)))else:print('暂不支持此插值函数')returncoefficient_mat = solve(phi_mat, support_values)return coefficient_matdef rbf_interpolation(support_points, coefficient_mat, interpolation_points, function_name = 'C2', radius = None):"""计算并返回RBF插值的结果:param support_points: 支撑点:param coefficient_mat: 插值系数矩阵:param interpolation_points: 插值点:param function_name: 插值函数名，默认为 Wendland C2:param radius: 插值函数作用半径，默认作用范围包含所有支撑点:return: interpolation_values, 插值点的物理量"""num_interpolation_points, dim = np.shape(interpolation_points)num_support_points = np.shape(support_points)[0]interpolation_values = np.zeros((num_interpolation_points, dim), dtype = np.float)for i in range(num_interpolation_points):for j in range(num_support_points):eta = np.linalg.norm(interpolation_points[i] - support_points[j])try:eta = eta / radiusif eta > 1:interpolation_values[i] += coefficient_mat[j] * 0  # phi = 0continueexcept TypeError:passfinally:if function_name == 'C2':phi = (1 - eta) ** 4 * (4 * eta + 1)elif function_name == 'cubic':phi = eta ** 3elif function_name == 'linear':phi = etaelif function_name == 'gaussian':sig = 1  # 高斯核宽度phi = np.exp(-1 * eta ** 2 / (2 * (sig ** 2)))else:print('暂不支持此插值函数')returninterpolation_values[i] += coefficient_mat[j] * phireturn interpolation_valuesdef calculation_rmse(y_hat, y):rmse = np.sqrt(np.mean(np.square(y - y_hat)))return rmse

测试代码如下：

test.py

import numpy as np
import matplotlib.pyplot as plt
import RBFdef 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_allif __name__ == '__main__':function_name = "gaussian"snum = 20  # control point数量ratio = 20  # 总数据点数量：snum*ratioxs = -8xe = 8x1 = np.linspace(xs, xe, snum).reshape(-1, 1)x2 = np.linspace(xs, xe, (snum - 1) * ratio + 1).reshape(-1, 1)y_sample, y_all = gen_data(x1, x2)plt.figure(1)w = RBF.rbf_coefficient(x1, y_sample, function_name)y_rec = RBF.rbf_interpolation(x1, w, x2, function_name)rmse = RBF.calculation_rmse(y_rec, y_all)print(rmse)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(function_name + ' kernel interpolation:$y=sin(\pi x/2)+cos(\pi x/3)$')plt.show()

测试结果如下（一维插值）：

总结

对于rbf插值中，一般控制点数量越多越密集，插值精度也就越高，但是随之而来的插值矩阵 Φ\PhiΦ 增大，计算量增大，模型训练（权值系数计算）效率变低，正所谓“No Free Lunch”。其次，不同的基函数在同样的数据上插值的表现也有所差异，因此需要选择最优的基函数进行插值计算。

参考

[1] https://blog.csdn.net/xfijun/article/details/105670892

[2] https://zhuanlan.zhihu.com/p/339854179