最小二乘曲线拟合原理,平面曲线与空间曲线,并附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)}$$

PS：但是这个目标公式最终是有问题的!!! 普通的直线可以拟合，但是接近与y轴的直线拟合的效果不好，这是个通病，所以这种拟合只适用于大部分情况。后文会提出解决的方法。 （详见第七部分）

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拟合效果图

6.封装函数的使用

6.1 二维点

# 1.准备数据
# x_for_fitting_list=[] #都是为python的列表格式的点
# y_for_fitting_list=[] #都是为python的列表格式的点
# pow_value=666 #带拟合的n次幂

#核心函数,求拟合点的曲线方程的函数
mat_A=least_square_fitting_2D(pow_value,x_for_fitting_list,y_for_fitting_list) #求出的mat_A 是多项式的系数矩阵

#根据拟合的系数求待画的直线,一般是给定x的起始位置和终点位置,来求y值,最终得到两个列表
# pow_value还是待拟合的次幂
# start, end, step 分别为x坐标轴的 开始 结束和 步长
x_axis_for_draw,y_axis_for_draw=solution_of_curve_points_2D(mat_A,pow_value,start, end, step) #2D画图调用
# x_axis_for_draw,y_axis_for_draw 列表里面的点可以用作画图,直线和曲线都可以.

6.2 三维点

##三维点的逻辑同二维是一样的,只是输入数据多了一维
mat_A,mat_B=least_square_fitting_3D(pow_value_xy,pow_value_yz,x.tolist(),y.tolist(),z.tolist())
x_axis_for_draw,y_axis_for_draw,z_axis_for_draw=solution_of_curve_points_3D(mat_A,mat_B,pow_value_xy,pow_value_yz,start, end, step) #2D画图调用

6.3 另一个函数

还有一个函数叫 get_least_squares_Fx()
他的作用主要是根据二维多项式的系数和维度求对应x的y值的,具体使用为: get_least_squares_Fx(mat_A,pow,in_x) . 属于单点计算。

7.补充完善算法

利用主成分分析方法（PCA）来改进以上最小二乘中，只有在单一方向上（y方向或者x方向）平方差值最小的问题。
特别备注:此方法只适用于直线的拟合，曲线拟合不太适用。可尝试 RANSAC（Random Sample Consensus，随机采样一致性） 等其他方法。

7.1 主成分分析方法的介绍

主成分分析（PCA）是一种通过线性变换将原始数据转换为一组新的、相互正交的变量的方法，这些新变量称为主成分。主成分的顺序是按照方差大小排列的，即第一个主成分解释了数据中最大的方差，第二个主成分解释了次大的方差，依此类推。
简单来说，这种方法就像在数据点云中找到最"扁"的方向，让所有点到这条直线的垂直距离最短。 就像把一把尺子放在散乱的点群中，调整尺子的方向和位置，使得所有点到尺子的垂直距离平方和最小。

7.2 主成分分析方法的推理步骤

具体步骤：

1. 找中心点：
   先把所有数据点的中心位置算出来（取x坐标平均值和y坐标平均值）。这个中心点就是最优直线必须经过的点。就像先找到点群的"重心"，把直线钉在这个点上再调整方向。

2. 计算数据分布特征：
   用数学公式计算数据在x方向的波动程度、y方向的波动程度，以及x和y之间的联动关系。这构成一个2x2的协方差矩阵，它就像一张"数据分布体检报告"。

3. 找主方向：
   通过矩阵的特征分解，找到两个互相垂直的特殊方向：

• 最长胖方向（主成分方向）： 数据在这个方向上伸展得最开，这就是最优直线的方向。
• 最瘦扁方向（垂直方向）： 数据在这个方向上最紧凑，这个方向就是直线法线方向。

4. 确定直线方程：
   用法线方向的向量结合中心点坐标，写出直线方程。例如法线方向是(a,b)，中心点是(avg_x,avg_y)，则直线方程为：

$$a(x-av{g}_x)+b(y-av{g}_y)=0$$

为什么这方法有效？

• 数据在主方向上的投影能保留最多信息（方差最大）
• 而在垂直方向上的波动最小（方差最小），自然就使垂直距离平方和最小

类比理解：
想象把数据点云压扁成一个椭圆，最优直线就是这个椭圆的长轴方向。法线方向则是短轴方向，这个方向数据最集中。

与普通最小二乘法的区别：

• 普通最小二乘法只考虑垂直或水平距离
• 这个方法考虑真正的垂直距离，适合x和y都有误差的情况

实际应用场景：

• 处理传感器数据（如同时存在x,y测量误差）
• 图像处理中识别物体的主方向
• 数据降维时保留主要特征

示例：
假设有三个点：(1,1), (2,2), (3,3)

1. 中心点是(2,2)
2. 协方差矩阵计算后会发现主方向正好是45度斜线
3. 最优直线就是y=x，完美符合直觉

这种方法通过数学保证了找到全局最优解，是处理二维数据拟合的经典方法。

7.2 进行直线拟合

好的，我们通过分步公式推导来揭示主成分分析（PCA）法求解最优直线系数的过程。以下用 6个关键步骤 展现推导链条：

---

步骤1：数据预处理

假设有n个数据点$(x_i,y_i)$，首先计算 中心点坐标：
$$x_o=\frac{1}{n}\sum x_i,y_o=\frac{1}{n}\sum y_i$$
将每个点平移到以中心为原点的新坐标：
$${\tilde{x}}_i=x_i-x_o,{\tilde{y}}_i=y_i-y_o$$

---

步骤2：构建协方差矩阵

计算 平移后数据 的协方差矩阵$S$（除以n-1是无偏估计，不影响方向）：
$$S=\frac{1}{n-1}\left[\begin{array}{cc}\sum {\tilde{x}}_i^2& \sum {\tilde{x}}_i{\tilde{y}}_i\\ \sum {\tilde{x}}_i{\tilde{y}}_i& \sum {\tilde{y}}_i^2\end{array}\right]=\left[\begin{array}{cc}{S}_{xx}& S_{xy}\\ S_{xy}& S_{yy}\end{array}\right]$$
这里$S_{xx}$是x方向的方差，$S_{yy}$是y方向的方差，$S_{xy}$是协方差。

---

步骤3：特征值分解

求协方差矩阵的特征值$\lambda$和特征向量$\mathbf{v}$。解特征方程：
$$\det (S-\lambda I)=0\Rightarrow {\lambda }^2-(S_{xx}+S_{yy})\lambda +(S_{xx}{S}_{yy}-S_{xy}^2)=0$$
解得两个特征值：
$${\lambda }_{\text{max}},{\lambda }_{\text{min}}=\frac{1}{2}\left[(S_{xx}+S_{yy})\pm \sqrt{(S_{xx}-S_{yy}{)}^2+4S_{xy}^2}\right]$$

---

步骤4：确定法线方向

对应 最小特征值${\lambda }_{\text{min}}$的特征向量$\mathbf{n}=(n_x,n_y)$是 直线的法线方向。具体求解：
$$(S-{\lambda }_{\text{min}}I)\mathbf{n}=0\Rightarrow \{\begin{array}{l}(S_{xx}-{\lambda }_{\text{min}})n_x+S_{xy}{n}_y=0\\ S_{xy}{n}_x+(S_{yy}-{\lambda }_{\text{min}})n_y=0\end{array}$$
特征向量方向为：
$$\mathbf{n}=\left(S_{xx}-S_{yy}-\sqrt{(S_{xx}-S_{yy}{)}^2+4S_{xy}^2},2S_{xy}\right)$$
（需归一化为单位向量）

---

步骤5：构造直线方程

已知法线方向$\mathbf{n}=(n_x,n_y)$，直线经过中心点$(x_o,y_o)$，方程为：
$$n_x(x-x_o)+n_y(y-y_o)=0$$
展开为标准形式：
$$\boxed{n_{x}x+n_{y}y-(n_x{x}_o+n_y{y}_o)=0}$$

---

步骤6：转换为斜截式（可选）

如果需要显式的斜率和截距：
$$y=-\frac{n_x}{n_y}x+\frac{n_x{x}_o+n_y{y}_o}{n_y}$$
对应斜率为$k=-n_x/n_y$，截距为$b=(n_x{x}_o+n_y{y}_o)/n_y$。

---

关键结论

• 法线向量$\mathbf{n}=(n_x,n_y)$直接给出直线方程的系数
• 特征值${\lambda }_{\text{min}}$对应垂直距离的方差最小值
• 特征向量的正交性保证了直线方向的最优性

---

几何解释图示

数据点云 → 中心平移 → 协方差矩阵 → 特征分解 → 法线方向 → 直线方程

整个过程通过矩阵的谱分解，将数据分布的主轴方向转化为直线方向，数学上严格保证了最小化垂直距离平方和。

7.3 特征值分解的解释

好的，我来详细解释这一步的原理和推导过程：

---

特征方程的来源

特征值分解的核心思想是寻找矩阵$S$的特殊方向：在这些方向上，矩阵的变换作用相当于简单的缩放。具体到我们的协方差矩阵，就是要找到两个正交的方向（特征向量），使得数据在这两个方向上的方差（特征值）达到极值。

---

详细推导步骤

1. 定义特征方程
   对于任意方阵$S$，若存在非零向量$\mathbf{v}$和标量$\lambda$，使得：
   $$S\mathbf{v}=\lambda \mathbf{v}$$
   则称$\lambda$是矩阵$S$的特征值，$\mathbf{v}$是对应的特征向量。

2. 改写为齐次方程
   将等式变形为：
   $$(S-\lambda I)\mathbf{v}=0$$
   其中$I$是单位矩阵。为了保证非零解存在，系数矩阵$S-\lambda I$必须是奇异的（不可逆），因此其行列式必须为零：
   $$\det (S-\lambda I)=0$$

3. 代入协方差矩阵形式
   对于二维协方差矩阵：
   $$S=\left[\begin{array}{cc}{S}_{xx}& S_{xy}\\ S_{xy}& S_{yy}\end{array}\right]$$
   构造$S-\lambda I$：
   $$S-\lambda I=\left[\begin{array}{cc}{S}_{xx}-\lambda & S_{xy}\\ S_{xy}& S_{yy}-\lambda \end{array}\right]$$

4. 计算行列式
   展开行列式计算：
   $$\det (S-\lambda I)=(S_{xx}-\lambda )(S_{yy}-\lambda )-S_{xy}^2$$
   展开括号：
   $$=S_{xx}{S}_{yy}-S_{xx}\lambda -S_{yy}\lambda +{\lambda }^2-S_{xy}^2$$
   整理为二次方程：
   $${\lambda }^2-(S_{xx}+S_{yy})\lambda +(S_{xx}{S}_{yy}-S_{xy}^2)=0$$

---

几何意义

• 特征值$\lambda$:
  代表数据在对应特征向量方向上的方差。
  -${\lambda }_{\text{max}}$：数据在主方向（最分散方向）的方差
  -${\lambda }_{\text{min}}$：数据在垂直方向（最集中方向）的方差

• 特征向量:
  -${\mathbf{v}}_{\text{max}}$：指向数据伸展最开的方向（最优直线的方向）
  -${\mathbf{v}}_{\text{min}}$：指向数据最压缩的方向（直线的法线方向）

---

为什么这样能求解最优直线？

1. 协方差矩阵的本质$S$描述了数据的分布形状：

• 对角线元素$S_{xx},S_{yy}$表示x和y方向的离散程度
• 非对角线元素$S_{xy}$表示x和y的协同变化

2. 特征分解的物理意义 通过找到$S$的特征方向，我们实际上在寻找：

• 最分散方向（${\mathbf{v}}_{\text{max}}$）：保留数据最多信息的方向
• 最集中方向（${\mathbf{v}}_{\text{min}}$）：数据波动最小的方向，这正是点到直线垂直距离最小的方向

---

图解辅助理解

假设数据分布如下图：

        ↗ 数据主要沿此方向分散（对应λ_max）
       /  
      /    
     · ——→ 数据在此方向最集中（对应λ_min）
  (中心点)

• 椭圆的长轴方向 →${\mathbf{v}}_{\text{max}}$
• 椭圆的短轴方向 →${\mathbf{v}}_{\text{min}}$

---

关键公式总结

1. 特征方程：
   $${\lambda }^2-(S_{xx}+S_{yy})\lambda +(S_{xx}{S}_{yy}-S_{xy}^2)=0$$
2. 特征值解：
   $$\lambda =\frac{1}{2}\left[(S_{xx}+S_{yy})\pm \sqrt{(S_{xx}-S_{yy}{)}^2+4S_{xy}^2}\right]$$
3. 最小特征值对应法线方向：
   $$\mathbf{n}=\left(S_{xx}-S_{yy}-\sqrt{(S_{xx}-S_{yy}{)}^2+4S_{xy}^2},2S_{xy}\right)$$

---

通过这种分解，我们实际上将数据分布的几何特征转化为代数计算，最终得到最优拟合直线的数学表达式。

希望这个解释更清楚地阐明了特征值分解在求解最优直线中的作用。