欢迎关注更多精彩
关注我，学习常用算法与数据结构，一题多解，降维打击。
问题描述
在平面内给定n个点 {Pi=(xi, yi)}, i = [1, n]。
分别使用多项式插值, 拉格朗日插值, 高斯基函数插值, 最小二乘多项式(最多m次方)逼近, 最小二乘多项式(最多m次方)岭回归逼近来拟合函数f(x)
多项式插值
设多项式 f ( x ) = ∑ i = 0 n − 1 a i x i 设多项式f(x) = \displaystyle\sum_{i=0}^{n-1}a_ix^i 设多项式f(x)=i=0∑n−1aixi
对于n个点就有n个方程
{ y 1 = a 0 + a 1 x 1 + a 2 x 1 2 + . . . + a n − 1 x 1 n − 1 y 2 = a 0 + a 1 x 2 + a 2 x 2 2 + . . . + a n − 1 x 2 n − 1 . . . y n = a 0 + a 1 x n + a 2 x n 2 + . . . + a n − 1 x n n − 1 \left \{\begin{array}{lc} y_1=a_0+a_1x_1+a_2x_1^2+...+a_{n-1}x_1^{n-1} \\ y_2=a_0+a_1x_2+a_2x_2^2+...+a_{n-1}x_2^{n-1} \\ ...\\y_n=a_0+a_1x_n+a_2x_n^2+...+a_{n-1}x_n^{n-1} \end{array} \right. ⎩ ⎨ ⎧y1=a0+a1x1+a2x12+...+an−1x1n−1y2=a0+a1x2+a2x22+...+an−1x2n−1...yn=a0+a1xn+a2xn2+...+an−1xnn−1
令
Y = [ y 1 y 2 . . . y n ] Y=\begin{bmatrix} y_1\\y_2\\...\\y_n \end{bmatrix} Y= y1y2...yn
B = [ 1 x 1 x 1 2 . . . x 1 n − 1 1 x 2 x 2 2 . . . x 2 n − 1 . . . . . . . . . . . . . . . 1 x n x n 2 . . . x n n − 1 ] B=\begin{bmatrix} 1&x_1&x_1^2&...&x_1^{n-1}\\1&x_2&x_2^2&...&x_2^{n-1}\\...&...&...&...&...\\1&x_n&x_n^2&...&x_n^{n-1} \end{bmatrix} B= 11...1x1x2...xnx12x22...xn2............x1n−1x2n−1...xnn−1
a = [ a 0 a 1 . . . a n − 1 ] a=\begin{bmatrix} a_0\\a_1\\...\\a_{n-1} \end{bmatrix} a= a0a1...an−1
求解方程组Ba=Y, 即可得到a向量,
拉格朗日插值
f ( x ) = ∑ k = 0 n − 1 y k l k ( x ) f(x) = \displaystyle\sum_{k=0}^{n-1}y_kl_k(x) f(x)=k=0∑n−1yklk(x)
上式中 l k ( x ) = ∏ j = 0 , j ≠ k n − 1 x − x j x i − x j 上式中l_k(x) = \displaystyle\prod_{j=0,j\ne k}^{n-1}\frac{x-x_j}{x_i-x_j} 上式中lk(x)=j=0,j=k∏n−1xi−xjx−xj
拉格朗日插值结果与多项式插值法得出的曲线是相同, 计算量比直接解方程少.
牛顿插值也有相同的功能.
高斯基函数插值
设 f ( x ) = a 0 + ∑ i = 1 n a i g i ( x ) 设f(x) = a_0 + \displaystyle\sum_{i=1}^{n}a_ig_i(x) 设f(x)=a0+i=1∑naigi(x)
g i ( x ) = e − ( x − x i ) 2 2 σ 2 , σ 为用户给定参数 g_i(x) = e^{-\frac {(x-x_i)^2}{2\sigma^2}}, \sigma为用户给定参数 gi(x)=e−2σ2(x−xi)2,σ为用户给定参数
a = [ a 0 a 1 a 2 . . . a n ] 为未知数 a=\begin{bmatrix} a_0\\a_1\\a_2\\...\\a_n \end{bmatrix} 为未知数 a= a0a1a2...an 为未知数
点代入f(x)
得到方程组
{ y 1 = a 0 + a 1 g 1 ( x 1 ) + a 2 g 2 ( x 1 ) + . . . + a n g n ( x 1 ) y 2 = a 0 + a 1 g 1 ( x 2 ) + a 2 g 2 ( x 2 ) + . . . + a n g n ( x 2 ) . . . y n = a 0 + a 1 g 1 ( x n ) + a 2 g 2 ( x n ) + . . . + a n g n ( x n ) \left \{\begin{array}{lc} y_1=a_0+a_1g_1(x_1)+a_2g_2(x_1)+...+a_ng_n(x_1) \\ y_2=a_0+a_1g_1(x_2)+a_2g_2(x_2)+...+a_ng_n(x_2) \\ ...\\ y_n=a_0+a_1g_1(x_n)+a_2g_2(x_n)+...+a_ng_n(x_n) \end{array} \right. ⎩ ⎨ ⎧y1=a0+a1g1(x1)+a2g2(x1)+...+angn(x1)y2=a0+a1g1(x2)+a2g2(x2)+...+angn(x2)...yn=a0+a1g1(xn)+a2g2(xn)+...+angn(xn)
上面只有n条方程, 有n+1个未知数, 我们加入一条构造方程
f ( x 1 + x 2 2 ) = y 1 + y 2 2 f(\frac {x_1+x_2}2) = \frac {y_1+y_2}2 f(2x1+x2)=2y1+y2
求解方程组即可得到a向量.
最小二乘多项式逼近
在多项式插值中, 随着点数增多,曲线阶会升高, 阶数过高时会带来不稳定因素。
我们可以通过固定幂基函数来作拟合逼近。
设 f ( x ) = ∑ i = 0 m a i x i , m < n 设f(x) = \sum_{i=0}^ma_ix_i, m<n 设f(x)=∑i=0maixi,m<n
构建方程组
{ y 1 = a 0 + a 1 x 1 + a 2 x 1 2 + . . . + a m x 1 m y 2 = a 0 + a 1 x 2 + a 2 x 2 2 + . . . + a m x 2 m . . . y n = a 0 + a 1 x n + a 2 x n 2 + . . . + a m x n m \left \{\begin{array}{lc} y_1=a_0+a_1x_1+a_2x_1^2+...+a_{m}x_1^{m} \\ y_2=a_0+a_1x_2+a_2x_2^2+...+a_{m}x_2^{m} \\ ...\\y_n=a_0+a_1x_n+a_2x_n^2+...+a_{m}x_n^{m} \end{array} \right. ⎩ ⎨ ⎧y1=a0+a1x1+a2x12+...+amx1my2=a0+a1x2+a2x22+...+amx2m...yn=a0+a1xn+a2xn2+...+amxnm
上述方程数多于未知数，可以利用最小二乘来得到最优解。
Y = [ y 1 y 2 . . . y n ] Y=\begin{bmatrix} y_1\\y_2\\...\\y_n \end{bmatrix} Y= y1y2...yn
B = [ 1 x 1 x 1 2 . . . x 1 m 1 x 2 x 2 2 . . . x 2 m . . . . . . . . . . . . . . . 1 x n x n 2 . . . x n m ] B=\begin{bmatrix} 1&x_1&x_1^2&...&x_1^{m}\\1&x_2&x_2^2&...&x_2^{m}\\...&...&...&...&...\\1&x_n&x_n^2&...&x_n^{m} \end{bmatrix} B= 11...1x1x2...xnx12x22...xn2............x1mx2m...xnm
a = [ a 0 a 1 . . . a n − 1 ] a=\begin{bmatrix} a_0\\a_1\\...\\a_{n-1} \end{bmatrix} a= a0a1...an−1
最小二乘求解
最小二乘优化 m i n ∣ ∣ B a − Y ∣ ∣ 2 ， a = ( B T B ) − 1 ( B T Y ) 最小二乘优化min||Ba-Y||^2，a = (B^{T}B)^{-1}(B^TY) 最小二乘优化min∣∣Ba−Y∣∣2，a=(BTB)−1(BTY)
岭回归逼近
在上述逼近算法中, 会存在过拟合的情况，可以加入2模正则项作为惩罚函数，希望系数的模越小越好。
最小二乘优化 m i n ∣ ∣ B a − Y ∣ ∣ 2 + λ ∣ ∣ a ∣ ∣ 2 ， λ = 0.5 最小二乘优化min ||Ba-Y||^2+\lambda||a||^2，\lambda=0.5 最小二乘优化min∣∣Ba−Y∣∣2+λ∣∣a∣∣2，λ=0.5
对 ∣ ∣ B a − Y ∣ ∣ 2 + λ ∣ ∣ a ∣ ∣ 2 求导 , 使得导数等于 0 ，得到 对 ||Ba-Y||^2+\lambda||a||^2求导, 使得导数等于0，得到 对∣∣Ba−Y∣∣2+λ∣∣a∣∣2求导,使得导数等于0，得到
2 B T Y − 2 B T B a − 2 λ I a = 0 2B^TY-2B^TBa-2\lambda I a=0 2BTY−2BTBa−2λIa=0
a = ( B T B + λ I ) − 1 B T Y a=(B^TB+\lambda I)^{-1}B^TY a=(BTB+λI)−1BTY
运行结果
代码链接：
https://gitcode.com/chenbb1989/geometric_model/tree/master/gohw1
sigma=50
sigma = 20
高斯基函数拟合如果没有合适的sigma，效果会很不好，可以通过RBF神经网络对sigma和miu进行训练。
本人码农，希望通过自己的分享，让大家更容易学懂计算机知识。创作不易，帮忙点击公众号的链接，帮忙转发，感激不尽。
