最小二乘法(又称最小平方法)是一种数学优化技术。它通过最小化误差的平方和寻找数据的最佳函数匹配。利用最小二乘法可以简便地求得未知的数据,并使得这些求得的数据与实际数据之间误差的平方和为最小。
“最小二乘法”是对过度确定系统,即其中存在比未知数更多的方程组,以回归分析求得近似解的标准方法。在这整个解决方案中,最小二乘法演算为每一方程式的结果中,将残差平方和的总和最小化。
以上来源维基百科。
从200多年前,最小二乘法一直是科技工作者整理实验数据、建立数学模型的有力工具。
最小二乘法主要应用在曲线拟合方面。
1.直线拟合
已知x-y平面上有一些点,怎样作出一条直线\(y=a+bx\)去拟合这些数据点,反映这组数据的总趋势?
从观察可以看出直线不可能通过所有点,在\(x_i\)处有误差:
$$e_i=y_i-(a+bx_i)(i=1,2…,n).$$
称\(e_i\)为直线\(y=a+bx\)在\(x_i\)处的残差,它反映了这条拟合直线在点\(x_i\)处拟合的好坏。怎样衡量总体上拟合的好坏?常用的方法有三种:
(1)残差的最大绝对值\(\max \limits_{1\le i \le n} |e_i| \);
(2) 残差的绝对值之和\(\sum _{i=1}^n|e_i|\)
(3) 残差的平方和\(\sum _{i=1}^ne^2_i\)
以上三种准则1,2两种由于含有绝对值计算比较复杂,一般都使用3准则。
那么原问题等价于求$$Q=\sum _{i=1}^n(y_i-(a+bx_i))^2$$的最小值。
由于Q与一次式有关,也就是与一次式的系数a,b有关,所以Q是a,b的二元函数。求一次式使Q最小,就是求二元函数Q的最小点。用微分法球Q的最小值点。令\(\frac{\partial Q}{\partial a}=0,\frac{\partial Q}{\partial b}=0\)
最后计算得到如下方程组。
$$
\begin{cases}
na+(\sum _{i=1}^nx_i)b=\sum _{i=1}^ny_i,\\
(\sum _{i=1}^nx_i)a+(\sum _{i=1}^nx_i^2)b=\sum _{i=1}^ny_ix_i\end{cases}
$$
直接计算求和公式的值,然后解方程组就得到a,b的值了。
哈哈,记得初高中那时候不知道偏微分,拟合直线直接记公式:
2.多项式拟合
多项式拟合的最小二乘法:
求一个m次多项式\(P_m(x)= \sum _{j=0}^ma_jx^j\)使残差平方和
$$Q=\sum _{i=1}^n(y_i- \sum _{j=0}^ma_jx^j_i)^2$$为最小,称这样作出的\(P_m(x)\)为多项式拟合的最小二乘法解。
求\(P_m(x)\)的过程同直线拟合一样(直线拟合不就是一次多项式拟合嘛),求m+1元函数的最小值点(使m+1元函数取得最小值的\(a_j,j=0,1,…m\))。令\(\frac{\partial Q}{\partial a_k}=0(k=0,1,…,m)\),计算偏导数后得:
$$\sum _{i=1}^n(y_i- \sum _{j=0}^ma_jx^j_i)x_i^k=0(k=0,1,…m)$$移项得
$$\sum _{i=1}^n(\sum _{j=0}^ma_jx^j_i)x_i^k=\sum _{i=1}^ny_ix_i^k(k=0,1,…m)$$交换顺序得$$\sum _{j=0}^m(\sum _{i=1}^nx_i^{j+k})a_j=\sum _{i=1}^ny_ix_i^k(k=0,1,…m)$$上面等式可以改写成如下矩阵相乘的形式
这样看着能舒服点,然后就是解方程组。分别计算\(\sum _{i=1}^nx_i^j(j=0,1,2,…2m)\)和\(\sum _{i=1}^nx_i^jy_i(j=0,1,2…m)\)带入方程组求解。(原式子可以写成矩阵形式XA=Y 那么A=\((X^TX)^{-1}X^TY\,如果不是方阵逆矩阵用伪逆矩阵代替,这个公式可以大胆去用),上面的方程组也可以用高斯消元去解.
事实上python numpy包中已经提供了用于解线性方程组XA=Y的方法numpy.linalg.solve(),或者matlab中直接X\Y求得。当然matlab中也提供了曲线拟合的方法lsqcurvefit()。
3.指数类型曲线拟合
已知一组数据点(采样点)(xi,yi)(i=1,2,3…n),且已知这组数据满足某一函数原型(经验公式)\(y=F(a,x)\),其中向量a待定.曲线拟合的最小二乘法就是要确定待定向量a,使残差平方和$$Q=\sum _{i=1}^n(y_i-F(a,x_i))^2$$最小。称使总误差最小的F(a,x)为曲线拟合的最小二乘解。
这里假设\(F(a,x)=a_1e^{a_2/x}\)。即\(y=a_1e^{a_2/x}\)
由于F(a,x)是指数函数,所以可以先把它化为线性模型,再按残差平方和最小规则确定向量a。具体做法如下:
对于函数y两边取对数,得:
$$lny=lna_1+\frac{a_2}{x}等价于c_0+c_1x^{-1}$$这是一个线性模型,按直线拟合的最小二乘法,求\(c_0,c_1\),使总误差$$Q=\sum _{i=1}^n(c_0+c_1x_i^{-1}-lny_i)^2$$为最小。
$$
\begin{cases}
nc_0+(\sum _{i=1}^nx_i^{-1})c_1=\sum _{i=1}^nlny_i,\\
(\sum _{i=1}^nx_i^{-1})c_0+(\sum _{i=1}^n(x_i^{-1})^2)c_1=\sum _{i=1}^nlny_ix_i^{-1}\end{cases}
$$把数据带入,求\(c_0,c_1\)然后再求\(a_1,a_2\)最后求得原函数。
4.平面拟合
平面方程的一般表达式为:
$$Ax+By+Cz+d=0(C\neq0)$$
$$z=-\frac{A}{C}-\frac{B}{C}-\frac{D}{C}$$
记:\(a0=-\frac{A}{C},a1=-\frac{B}{C},a2=-\frac{D}{C}\)则\(z=a0x+a1y+a2\)
平面方程拟合:
对于一系列的n个点\((n \ge3)\):
\((x_i,y_i,z_i),i=0,1,…n-1\)
使残差的平方和
$$S=\sum _{i=1}^n(a_0x+a_1y+a_2-z)^2$$最小,应满足:\(\frac{\partial S}{\partial a_k}=0(k=0,1,2)\)
推导后可以得到如下方程组
\(\left[ \begin{matrix} \sum x_i^2 \ \sum x_iy_i \ \sum x_i \\ \sum x_iy_i \ \sum y_i^2 \ \sum y_i \\ \sum x_i \ \sum y_i \ n
\end{matrix}\right]\)\(\left[ \begin{matrix} a_0 \\ a_1 \\
a_2 \end{matrix}\right]\)=\(\left[ \begin{matrix} \sum x_iz_i \\ \sum y_iz_i \\ \sum z_i \end{matrix}
\right]\)
可以看成X*A=Y已知XY求A,由于X存在逆矩阵那么代码就好写多了。下面是之前写的代码纯c++.