数值计算-插值与回归

魏尔施特拉斯逼近定理(Stone–Weierstrass theorem、Weierstrass theorem)共有两个:

  1. 闭区间上的连续函数可用多项式级数一致逼近。(泰勒级数特例)
  2. 闭区间上周期为2π2\pi的连续函数可用三角函数级数一致逼近。(傅里叶级数特例)

在一个闭区间内,多项式可以逼近任何一个连续函数。在这种情况下,研究如何使用多项式来进行数据的拟合就成为了一个关键的问题(最小误差)。这就是所谓的多项式的回归算法

在数值分析中,多项式插值是使用多项式对一组给定的数据来进行插值的过程,也就是说,在给定一组数据的情况下,多项式插值的目的就是找到一个多项式,使得它恰好通过这些数据点。拉格朗日插值算法可以对实践中的某些物理量进行观测,然后得到一个多项式,从而表示各个结果之间的内在联系。多项式插值算法也是数值积分和数值常微分方程的算法基础。

给定两个点(x0,y0),(x1,y1)(x_0,y_0),(x_1,y_1),线性插值是一条连接两个点的直线。对于x(x0,x1)x\in(x_0,x_1),y可以通过斜截式给出: yy0xx0=y1y0x1x0    y=y0+y1y0x1x0(xx0),x(x0,x1)\frac{y-y_0}{x-x_0}=\frac{y_1-y_0}{x_1-x_0}\implies y=y_0+\frac{y_1-y_0}{x_1-x_0}(x-x_0),x\in(x_0,x_1)

对于多个点的线性插值来说,x的值之和最近的两点相关,每个采样点都是用折现连接的(如图1)。其公式可以表达为: y=f(x)=yi+yiyi1xixi1(xxi),x(xi1,xi)y=f(x)=y_i+\frac{y_i-y_{i-1}}{x_i-x_{i-1}}(x-x_i),x\in(x_{i-1},x_i)

线性插值

Figure 1: 线性插值

图1 线性插值

双线性插值一般用于z=f(x,y)z=f(x,y)的2+1维模式,在平面图形处理中有广泛应用。其本质上是做2+1次线性插值

典型的示意图如图2:

双线性插值

Figure 2: 双线性插值

图2 双线性插值

在2D网格上,我们给出四个点(x1,y1),(x1,y2),(x2,y1),(x2,y2)(x_1,y_1),(x_1,y_2),(x_2,y_1),(x_2,y_2)和四个点的函数值,要求x(x1,x2),y(y1,y2)x\in(x_1,x_2),y\in(y_1,y_2),双线性插值求f(x,y)f(x,y)则如下:

f(R1)f(Q11)x2xx2x1+f(Q21)xx1x2x1f(R_1)\approx f(Q_{11})\frac{x_2-x}{x_2-x_1}+f(Q_{21})\frac{x-x_1}{x_2-x_1} f(R2)f(Q21)x2xx2x1+f(Q22)xx1x2x1f(R_2)\approx f(Q_{21})\frac{x_2-x}{x_2-x_1}+f(Q_{22})\frac{x-x_1}{x_2-x_1} 有点像加权平均。然后将R1,R2R_1,R_2再做一次线性插值: f(x,y)=f(P)=f(R1)y2yy2y1+f(R2)yy1y2y1f(x,y)=f(P)=f(R_1)\frac{y_2-y}{y_2-y_1}+f(R_2)\frac{y-y_1}{y_2-y_1} 带入f(R1),f(R2)f(R_1),f(R_2)可得: f(x,y)=1(x2x1)(y2y1)[f(Q11)(x2x)(y2y)+f(Q12)(xx1)(y2y)+f(Q21)(x2x)(yy1)+f(Q22)(xx1)(yy1)]f(x,y)=\frac{1}{(x_2-x_1)(y_2-y_1)}[f(Q_{11})(x_2-x)(y_2-y)+f(Q_{12})(x-x_1)(y_2-y) \\ +f(Q_{21})(x_2-x)(y-y_1)+f(Q_{22})(x-x_1)(y-y_1)] 写成矩阵形式: f(x,y)=1(x2x1)(y2y1)[x2xxx1][Q11Q12Q21Q22][y2yyy1]f(x,y)=\frac{1}{(x_2-x_1)(y_2-y_1)}\begin{bmatrix} x_2-x&x-x_1\end{bmatrix}\begin{bmatrix} Q_{11}&Q_{12}\\Q_{21}&Q_{22}\end{bmatrix}\begin{bmatrix} y_2-y\\y-y_1\end{bmatrix} 在CV计算中,x2x1,y2y1x_2-x_1,y_2-y_1经常为1(相邻像素点),因此可直接写成: f(x,y)=[x2xxx1][Q11Q12Q21Q22][y2yyy1]f(x,y)=\begin{bmatrix} x_2-x&x-x_1\end{bmatrix}\begin{bmatrix} Q_{11}&Q_{12}\\Q_{21}&Q_{22}\end{bmatrix}\begin{bmatrix} y_2-y\\y-y_1\end{bmatrix}

给定(n+1)个数据点(xi,yi),i{0,...,n}(x_i,y_i),i\in\{0,...,n\},且任意两个xix_i都不相等。可以用一个p阶多项式(0pn)(0\leq p \leq n)进行插值,使: p(xi)=yi,0inp(x_i)=y_i,0 \leq i \leq n p(x)p(x)具体可以写为: p(x)=anxn+an1xn1++a1x+a0p(x)=a_nx^n+a_{n-1}x^{n-1}+\dotsb+a_1x+a_0 求解系数aia_i这是一个n+1元一次方程组,写成范德蒙矩阵为: [x0nx0n1x01x1nx1n1x11xnnxnn1xn1][anan1a0]=[y0y1yn]\begin{bmatrix}x_0^n&x_0^{n-1}&\dotsb&x_0&1\\ x_1^n&x_1^{n-1}&\dotsb&x_1&1\\ \vdots&\vdots&\ddots&\vdots&\vdots\\ x_n^n&x_n^{n-1}&\dotsb&x_n&1\\ \end{bmatrix} \begin{bmatrix}a_n\\ a_{n-1}\\ \vdots\\ a_0\\ \end{bmatrix}= \begin{bmatrix}y_0\\ y_1\\ \vdots\\ y_n\\ \end{bmatrix} 可以通过标准方式来求解系数aia_i,由于在范德蒙矩阵中,如果xi,xj\forall x_i,x_j两两不同,则范德蒙矩阵的行列式不为0,可以获得唯一解。以后的无论是牛顿法还是Lagrange插值法得到的结果和这种方式解出来是一样的,只是加快了计算。

同样,给点(n+1)个数据点(xi,yi),i{0,...,n}(x_i,y_i),i\in\{0,...,n\},且任意两个xix_i都不相等,可以通过Lagrange插值公式得到多项式,其可以写成: L(x)=i=0nli(x)yiL(x)=\sum_{i=0}^{n}l_i(x)y_i 其中,li(x)l_i(x)为Lagrange插值系数,表达式为: li(x)=0jn,jixxjxixjl_i(x)=\prod_{0 \leq j \leq n,j\neq i}\frac{x-x_j}{x_i-x_j} 从表达式中,我们不难得到以下性质:

  1. imi \neq m,则li(xm)=0l_i(x_m)=0,因为xm==xj,\exist x_m==x_j,使得分子为0
  2. i=mi=m,则li(xm)=1l_i(x_m)=1,因为分子和分母完全相同。
  3. L(x)L(x)的阶数最大为n。0jn,jixxjxixj\prod_{0 \leq j \leq n,j\neq i}\frac{x-x_j}{x_i-x_j}一共有n个x相乘。
  4. 插值多项式是唯一的,和范德蒙矩阵求出来的是一样的。

线性回归使用一个线性函数(一次函数),调节各个参数,使它尽量满足数据误差最小。 其他笔记中写过,根据最小方差准则计算可逆矩阵回归参数的公式为: α=(XTX)1XTy\vec{\alpha}=(X^TX)^{-1}X^T\vec{y}

给定n个数据点(xi,yi),i{1,,n}(x_i,y_i),i\in\{1,\dots,n\},可以建立一个m阶多项式来预测yy的值,其推到公式如下: y=β0+β1x+β2x2++βmxm+ϵy=\beta_0+\beta_1 x+\beta_2 x^2+\dotsb+\beta_m x^m+\epsilon 对于每一个二元点有: yi=β0+β1xi+β2xi2++βmxim+ϵi,1iny_i=\beta_0+\beta_1 x_i+\beta_2 x_i^2+\dotsb+\beta_m x_i^m+\epsilon_i,1 \leq i \leq n 如果把这些点组合成矩阵形式,则有: [y0y1yn]=[1x1x1m1x2x2m1xnxnm][β0β1βm]+[ϵ0ϵ1ϵm]\begin{bmatrix}y_0\\ y_1\\ \vdots\\ y_n\\ \end{bmatrix}= \begin{bmatrix}1&x_1&\dotsb&x_1^m\\ 1&x_2&\dotsb&x_2^m\\ \vdots&\vdots&\ddots&\vdots\\ 1&x_n&\dotsb&x_n^m\\ \end{bmatrix} \begin{bmatrix}\beta_0\\ \beta_1\\ \vdots\\ \beta_m\\ \end{bmatrix}+ \begin{bmatrix}\epsilon_0\\ \epsilon_1\\ \vdots\\ \epsilon_m\\ \end{bmatrix} y=Xβ+ϵ\vec{y}=X\vec{\beta}+\vec{\epsilon},根据最小方差估计,同线性回归可得多项式回归的系数: β=(XTX)1XTy(m<n)\vec{\beta}=(X^TX)^{-1}X^T\vec{y}(m<n)

其他回归算法如脊回归,Log回归等有专门笔记记录。

内插法求解以下的问题:有一未知函数在一些特定位置下的值,求未知函数在已知数值的点之间某一点的值。

外推法类似内插法,但需要知道数值的点是在其他已知数值点的范围以外。一般而言外推法的误差会大于内插法。

曲线拟合是在已知一些数据的条件下,找到一条曲线完全符合现有的数据,数据可能是一些特定位置及其对应的值,也可能是其他资料,例如角度或曲率等。

回归分析类似曲线拟合,也是根据一些特定位置及其对应的值,要找到对应曲线。但回归分析考虑到数据可能有误差,因此所得的的曲线不需要和数据完全符合。一般会使用最小方差法来进行回归分析。[1]

[1] 数值分析,维基百科https://zh.wikipedia.org/wiki/%E6%95%B0%E5%80%BC%E5%88%86%E6%9E%90

[2] 多项式插值算法与回归算法.张戎.知乎.https://zhuanlan.zhihu.com/p/33690607