最小二乘曲线拟合原理,平面曲线与空间曲线,并附python实现demo

1.最小二乘思想

1.1 基本描述

最小二乘法是一种在误差估计、不确定度、系统辨识及预测、预报等数据处理诸多学科领域得到广泛应用的数学工具。

1.2 定义

最小二乘法（又称最小平方法）是一种数学优化技术。它通过最小化误差的平方和寻找数据的最佳函数匹配。利用最小二乘法可以简便地求得未知的数据，并使得这些求得的数据与实际数据之间误差的平方和为最小。

1.3 应用范围

最小二乘法可用于曲线拟合，其他一些优化问题也可通过最小化能量或最大化熵用最小二乘法来表达。

以上引用于 <<百度百科>>词条:最小二乘

1.4 个人理解

所谓最小二乘的算法，主要思想就是“所有的值”的平方加起来取最小的值，从而求得最小值的时候带求未知数此时的值。而核心就是如何求得平方和最小哪个位置的点：通常方法为求偏导数，使其等于0那一点。

2.最小二乘拟合曲线

2.1 拟合基本思想

拟合的概念描述: 多个散点拟合成曲线，但前提要指定是多少次幂的曲线，通常为一元多次曲线方程。效果图入下所示：

## 2.2 数学模型 - **以平面曲线拟合为例**

已知： 散点数据,共有n个点.
$$[(x_1,y_1),(x_2,y_2),(x_3,y_3)....(x_n,y_n)]$$
待求： 一元多次曲线公式,假设最高阶为k.
$$y=a_0+a_{1}x+a_2{x}^2+...+a_n{x}^k=F(x)\text{\hspace{1em}(1)}$$

2.2.1目标函数

描述: 多点拟合成曲线，若选择出最优曲线，原则为（也是这类最小二乘的核心思想）：散点到曲线距离的和的绝对值最小。BUT！but! But! 点到直线的距离并不是像点到直线的距离公式那么求的，如果那样求的话，点到曲线的距离公式要相当的麻烦。 故，这里为了简化计算又不影响整体思想，采用x值和y值的差值的绝对值来进行衡量。即$|x_i-x|和|y_i-y|$的值取相对最小值，为了计算再简化，取未知点$(x,y)$为$(x_i,F(x_i))$,这样$X$方向的差值就为0，所以只用关心$Y$ 方向的差值$|y_i-F(x_i)|$即可.最终的目标是所有点的Y方向的差值的平方和最小.

最终的目标函数为:

$$min\sum _{i=1}^n(y_i-F(x_i){)}^2\text{\hspace{1em}(2)}$$

2.3 公式推导

• 设各个点的偏差和为R,n为散点个数,k为多项式的最高阶数.由公式2目标函数得:
  $$R^2=\sum _{i=1}^n[y_i-(a_0+a_1{x}_i+a_2{x}_i^2+...+a_n{x}_i^k){]}^2\text{\hspace{1em}(3)}$$
• 为了求得符合条件的$a_i$值，对等式右边求$a_i$偏导数,正常情况下导数为零点即为极值点,而依据所求得的一元多次方程的特性与实际曲线的存在情况,应该只有一个极小值符合情况.进一步推论得到每一个$a_i$的偏导数都应该为0,最终联合解算的$a_i$值才为最终多项式的系数.过程为如下::
  $$\frac{\partial R}{\partial a_0}=-2\sum _{i=1}^n[y_i-(a_0+a_1{x}_i+a_2{x}_i^2+...+a_k{x}_i^k)]\ast \frac{\partial a_0}{\partial a_0}=0$$
  $$\frac{\partial R}{\partial a_1}=-2\sum _{i=1}^n[y_i-(a_0+a_1{x}_i+a_2{x}_i^2+...+a_k{x}_i^k)]\ast \frac{\partial a_1}{\partial a_1}\ast x_i^1=0$$
  $$\frac{\partial R}{\partial a_2}=-2\sum _{i=1}^n[y_i-(a_0+a_1{x}_i+a_2{x}_i^2+...+a_k{x}_i^k)]\ast \frac{\partial a_2}{\partial a_2}\ast x_i^2=0$$
  $$...$$
  $$...$$
  $$...$$
  $$\frac{\partial R}{\partial a_n}=-2\sum _{i=1}^n[y_i-(a_0+a_1{x}_i+a_2{x}_i^2+...+a_k{x}_i^k)]\ast \frac{\partial a_k}{\partial a_k}\ast x_i^k=0$$

$$......\text{\hspace{1em}(4)}$$

• 将等式整理下，然后应该可以得到下面的等式:
  $$\sum _{i=1}^n[(a_0+a_1{x}_i+a_2{x}_i^2+...+a_k{x}_i^k)]\ast \frac{\partial a_0}{\partial a_0}\ast x_i^0=\sum _{i=1}^n{y}_i\ast x_i^0$$
  $$\sum _{i=1}^n[(a_0+a_1{x}_i+a_2{x}_i^2+...+a_k{x}_i^k)]\ast \frac{\partial a_1}{\partial a_1}\ast x_i^1=\sum _{i=1}^n{y}_i\ast x_i^1$$
  $$\sum _{i=1}^n[(a_0+a_1{x}_i+a_2{x}_i^2+...+a_k{x}_i^k)]\ast \frac{\partial a_2}{\partial a_2}\ast x_i^2=\sum _{i=1}^n{y}_i\ast x_i^2$$
  $$...$$
  $$...$$
  $$...$$
  $$\sum _{i=1}^n[(a_0+a_1{x}_i+a_2{x}_i^2+...+a_k{x}_i^k)]\ast \frac{\partial a_k}{\partial a_k}\ast x_i^k=\sum _{i=1}^n{y}_i\ast x_i^k$$

$$......\text{\hspace{1em}(5)}$$

• 将等式左边进行一下化简，然后应该可以得到下面的等式:
  $$n\ast a_0+a_1\sum _{i=1}^n{x}_i+...+a_k\sum _{i=1}^n{x}_i^k=\sum _{i=1}^n{y}_i\ast x_i^0$$
  $$a_0\sum _{i=1}^n{x}_i+a_1\sum _{i=1}^n{x}_i^2+...+a_k\sum _{i=1}^n{x}_i^{k+1}=\sum _{i=1}^n{y}_i\ast x^1$$
  $$a_0\sum _{i=1}^n{x}_i^2+a_1\sum _{i=1}^n{x}_i^3+...+a_k\sum _{i=1}^n{x}_i^{k+2}=\sum _{i=1}^n{y}_i\ast x_i^2$$
  $$...$$
  $$...$$
  $$...$$
  $$a_0\sum _{i=1}^n{x}_i^k+a_1\sum _{i=1}^n{x}_i^{k+1}+...+a_k\sum _{i=1}^n{x}_i^{2k}=\sum _{i=1}^n{y}_i\ast x_i^k$$
  $$......\text{\hspace{1em}(6)}$$

• 把这些等式表示成矩阵的形式，就可以得到下面的矩阵:

$$\left[\begin{array}{cccc}n& {\sum }_{i=1}^n{x}_i& \cdots & {\sum }_{i=1}^n{x}_i^k\\ {\sum }_{i=1}^n{x}_i& {\sum }_{i=1}^n{x}_i^2& \cdots & {\sum }_{i=1}^n{x}_i^{k+1}\\ ⋮& \ddots & ⋮& \cdots \\ {\sum }_{i=1}^n{x}_i^k& {\sum }_{i=1}^n{x}_i^{k+1}& \cdots & {\sum }_{i=1}^n{x}_i^{2k}\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ \cdots \\ a_k\end{array}\right]=\left[\begin{array}{c}{\sum }_{i=1}^n{y}_i\ast x_i^0\\ {\sum }_{i=1}^n{y}_i\ast x_i^1\\ \cdots \\ {\sum }_{i=1}^n{y}_i\ast x_i^k\end{array}\right]=\left[\begin{array}{c}{\sum }_{i=1}^n{y}_i\\ {\sum }_{i=1}^n{y}_i\ast x_i^1\\ \cdots \\ {\sum }_{i=1}^n{y}_i\ast x_i^k\end{array}\right]$$
$$......\text{\hspace{1em}(7)}$$

• 将这个范德蒙得矩阵化简后(此公式为一种理想模式,符号跟之前的没有关系,只是说明可以化简到下面这样的模式,方便讲解运算)可得到:
  $${\mathbf{J}}_F=\left[\begin{array}{cccc}1& x_1& \cdots & x_1^k\\ 1& x_2& \cdots & x_2^k\\ ⋮& \ddots & ⋮& \cdots \\ 1& x_n& \cdots & x_n^k\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ \cdots \\ a_k\end{array}\right]=\left[\begin{array}{c}{y}_1\\ y_2\\ \cdots \\ y_n\end{array}\right]$$
  $$......\text{\hspace{1em}(8)}$$

• 也就是说最终简化为矩阵运算$X\ast A=Y$求$A$的问题，那么$A=Y\ast X^{-1}$，便得到了系数矩阵$A$，同时，我们也就得到了拟合曲线.

注意: 这里$X$是已知数,$A$是未知数,其中最复杂的就是求矩阵$A$的系数$X$的过程,一定要细心.

3.拓展

3.1 空间曲线最小二乘拟合
已知: 空间内的点 $(x_1,y_1,z_1),(x_2,y_2,z_2),(x_3,y_3,z_3)...(x_n,y_n,z_n)$
所要拟合的多项式直线的最高次幂: $k$
**待求:**多项式表达式
$$z=a_0+a_{1}x+a_2{x}^2+...+a_n{x}^k+b_0+b_{1}y+b_2{y}^2+...+b_n{y}^k=F(x,y)$$
根据经典最小二乘算法可以拓展到加权最小二乘算法等,具体应用场景可以查询相关资料.
求解思路: 投影法.也就是把三维问题化成二维问题来解决.也就是将空间曲线F(x,y,z)投影到x,y平面和y,z平面即可.只用投影到两个平面即可.因为投影后空间曲线就变成空间曲面了,而两次投影生成两个曲面,而两个曲面就可以确定一个,且唯一一个空间曲线.故连个空间曲面的方程联立后就是所求的空间曲线.
**公式推导部分:**基本同上二维部分.
假设在x,y平面和y,z平面求空间曲面
注意空间点$p(x_n,y_n,z_n)$在使用时,要每次带入对应的投影面的值,比如在xy平面,已知值为$p_{xy}(x_n,y_n)$,在yz平面时,已知点值为$p_{yz}(y_n,z_n)$.
然后根据上面平面最小二乘拟合公式分别在每个投影面求得平面方程为:
$$y=a_0+a_{1}x+a_2{x}^2+...+a_n{x}^k\text{\hspace{1em}(9)}$$
$$z=b_0+b_{1}y+b_2{y}^2+...+b_n{y}^k\text{\hspace{1em}(10)}$$
最后公式9与公式10联立,既可以求得看见曲线方程,形式大致为:
$$z=a_0^{\prime }+a_1^{\prime }x+a_2^{\prime }{x}^2+...+a_n^{\prime }{x}^k+b_0^{\prime }+b_1^{\prime }y+b_2^{\prime }{y}^2+...+b_n^{\prime }{y}^k\text{\hspace{1em}(11)}$$

4.参考文献

https://blog.csdn.net/neuldp/article/details/52049367
https://www.cnblogs.com/paiandlu/p/7843236.html
https://www.datalearner.com/blog/1051539222770649

5.示例代码

点击跳转下载链接

demo里含平面曲线(2D)和空间曲线(3D)的实现源码,并通过画图软件实时的二维或者三维展示

执行代码后效果如下图所示：

2D拟合效果图

3D拟合效果图