拟合简介

拟合(fitting)是数理科学中一个重要的方法,被广泛应用于数学、物理、统计学、金融等各个领域。形象的说,拟合就是把平面上一系列的点,用一条光滑的曲线连接起来。因为这条曲线有无数种可能,从而有各种拟合方法。拟合的曲线一般可以用函数表示,根据这个函数的不同有不同的拟合名字。
在科学实验或社会活动中,人们常常需要观测很多数据的规律,通过实验或者观测得到量x与y的一组数据对,其中x_i是彼此不同的。人们希望用一类与数据本质规律相适应的解析表达式来反映量x与y之间的依赖关系,其中为未知参数。拟合即为通过观测数据寻求参数b_i的最佳估计值,即寻求最佳的理论函数

线性模型与非线性模型

若理论函数的线性函数,则称为线性拟合模型,否则称为非线性拟合模型。
在实际中,许多现象之间的关系往往并不是线性的,而是呈现某种曲线关系。如服药后血药浓度与时间的关系;病毒剂量与致死率的关系;化学反应的反应物浓度与反应速度的关系。这就产生了曲线拟合,用连续曲线近似地刻画或比拟平面上离散点组所表示的坐标之间的函数关系。
拟合插值逼近是数值分析的三大基础工具,通俗意义上它们的区别在于:
  • 拟合:已知点列,从整体上靠近它们
  • 插值:已知点列并且插值函数完全经过点列
  • 逼近:已知曲线或者点列,通过逼近使得构造的函数无限靠近它们。
曲线拟合问题一般都使用最小二乘法来解决。最小二乘法通过最小化误差的平方和寻找数据的最佳函数匹配,可以用来简便地求得未知的数据,并使得这些求得的数据与实际数据之间误差的平方和为最小。

历史发展

1801年,意大利天文学家朱赛普·皮亚齐发现了第一颗小行星谷神星。经过40天的跟踪观测后,由于谷神星运行至太阳背后,使得皮亚齐失去了谷神星的位置。随后全世界的科学家利用皮亚齐的观测数据开始寻找谷神星,但是根据大多数人计算的结果来寻找谷神星都没有结果。时年24岁的高斯也计算了谷神星的轨道。奥地利天文学家海因里希·奥尔伯斯根据高斯计算出来的轨道重新发现了谷神星。
高斯使用的最小二乘法的方法发表于1809年他的著作《天体运动论》中,而法国科学家勒让德于1806年独立发明“最小二乘法”,但因不为世人所知而默默无闻。勒让德曾与高斯为谁最早创立最小二乘法原理发生争执。
1829年,高斯提供了最小二乘法的优化效果强于其他方法的证明,因此被称为高斯-马尔可夫定理。
如今对于曲线拟合问题仍然使用最小二乘法来解决,不过发展出了大量优化最小二乘问题的算法,例如梯度下降法,牛顿法,高斯-牛顿法,Levenberg-Marquardt算法等等。

最小二乘法(Least Square)

最小二乘法通过最小化误差的平方和寻找数据的最佳函数匹配。利用最小二乘法可以简便地求得未知的数据,并使得这些求得的数据与实际数据之间误差的平方和为最小。最小二乘法被广泛应用于曲线拟合等数学优化问题中。
假设现在有一系列的数据点,那么由给出的拟合函数得到的估计量就是。那么怎么评估拟合函数与实际待求解的函数的拟合程度比较高呢?这里需要引入范数的概念。

范数

对于在上的向量,有三种常用范数:
  • 范数:
  • 范数:
  • 范数:
前两种范数不利于进行微分运算,在数据量很大的情况下计算量太大,不具有可操作性。因此一般使用的是2-范数。
引入残差的概念
其中h(x)为拟合函数。拟合程度,就是指拟合函数h(x)与待求解的函数y之间的相似性。残差的范数越小,相似性就越高。

定义

由此,我们可以写出最小二乘法的定义:
对于给定的数据,在取定的假设空间中,求解,使得残差范数最小,即
从几何上讲,就是寻找与给定点距离平方和最小的曲线,其中称为拟合函数或者最小二乘解,求解拟合函数h(x)的方法称为曲线拟合的最小二乘法。 当的形式时,这样的拟合称为线性拟合,其它的被称为非线性拟合。特别地,当这样多项式的形式时,被称为多项式拟合

代码实现

由于最小二乘法与线性拟合和非线性拟合并不属于同一个级别的概念,拟合都是可以基于最小二乘法来得到拟合曲线的,所以此处的最小二乘拟合是一个定义了拟合函数模型的拟合程序,程序中以为例。
Python实现:
MATLAB实现:

线性拟合

介绍

当目标拟合函数为的形式时,这样的拟合称为线性拟合。根据最小二乘法的定义,线性拟合的目标是要找到这样的w,使得
此处令
为样本的平方损失函数,即我们需要优化的损失函数。
这是一个典型的极值问题,要使Q(w)最小,只需要分别对求偏导数,然后令偏导数为零,就可得到极值点,即:
求解以上方程组即可得到拟合曲线的结果。

代码实现

在Python中实现线性拟合可以直接利用numpy库的polyfit函数,代码中的error即损失函数Q(w)的值。
Python实现:
MATLAB实现:

多项式拟合

介绍

当目标拟合函数为的形式时,这样的拟合称为多项式拟合,线性拟合是情况下多项式拟合的特例。
根据最小二乘法的定义,我们要求一组使得损失函数Q(w)最小。即可以转化为这样的方程组求解:
从算法上来看,数据最小拟合的多项式方法是解一个超定方程组的最小二乘解。而多项式拟合所引出的正规方程组恰好是用超定方程组的系数矩阵的转置矩阵去左乘超定方程组左、右两端所得。
正规方程组的系数矩阵是一个病态矩阵,这类方程组被称为病态方程组。当系数矩阵或者是右端向量有微小的误差时,可能引起方程组准确解有很大的误差。
为了避免求解这样的线性方程组,在做多项式拟合时可以将多项式中的各次幂函数做正交化变换,使得所推出的正规方程的系数矩阵是对角矩阵。

矩阵的条件数

矩阵的病态程度可以用条件数来描述。设为非奇异矩阵,称
为矩阵的条件数。当A的条件数相对较大,时,线性方程组Ax=b时病态的,当系数矩阵或者是右端向量有微小的误差时,可能引起方程组准确解有很大的误差。条件数越大,越难以用一般的计算方法求得比较准确的解。

模型变种

由于非线性拟合包含线性函数以外的所有函数,因此可以取非多项式形式,常见的非多项式形式有以下几种:
  • 双曲线形式:
  • 双对数函数形式:
  • 线性-对数形式:
  • Logistic形式:
  • 指数形式:
  • 幂函数形式:
  • Gompertz形式:
此外,拟合函数可能由多个自变量组成,而非被限定在平面中,因此还有形如等拟合函数的拟合问题。当然,形如的拟合问题也被称为线性拟合问题。
在不同的下,拟合目标都是根据最小二乘法的定义,求使得损失函数最小的

求解方法

线性最小二乘拟合需要求解的方程存在闭式解,即;而非线性拟合没有闭式解,一般用迭代法来解决非线性拟合问题。迭代法是一种不断用变量的旧值递推新值的过程,在计算过程中每一步更新未知量从而来逼近最优解。迭代法包含以下几种常用方法:梯度下降法,牛顿法,高斯-牛顿法以及Levenberg-Marquardt法。

梯度下降法

在微积分里面,对多元函数的参数求偏导数,把求得的各个参数的偏导数以向量的形式写出来,就是梯度。例如对于函数,分别对求偏导数,可以得到梯度向量为
梯度的几何意义是函数在某处变化最快的方向,具体来说,对于函数在点,沿着梯度向量的方向是变化最快的地方。
进一步,沿着梯度向量的方向,更加容易找到函数的最大值。反之,沿着梯度向量相反的方向,也就是的方向,函数减少最快,即更加容易找到函数的最小值。而我们要求的最小二乘解即是后面一种情况。
我们在一座大山上的某处位置,由于我们不知道怎么下山,于是决定走一步算一步,也就是在每走到一个位置的时候,求解当前位置的梯度,沿着梯度的负方向,也就是当前最陡峭的位置向下走一步,然后继续求解当前位置梯度,向这一步所在位置沿着最陡峭最易下山的位置走一步。这样一步步的走下去,一直走到觉得我们已经到了山脚。
当然这样走下去,有可能我们不能走到山脚,而是到了某一个局部的山峰低处(即局部最小值)。
从以上例子可以看出梯度下降不一定能够找到全局的最优解,有可能是一个局部最优解。当然,如果损失函数是凸函数,梯度下降法得到的解就一定是全局最优解。
以下用数学公式的形式来说明梯度下降的算法过程:
设拟合函数为
其中包含参数,即函数拟合需要得到的最终结果,如线性拟合
损失函数为
我们先对参数及算法终止距离和步长进行初始化,然后进行以下步骤:
1.确定当前位置的损失函数的梯度
2.用步长乘以损失函数的梯度,得到当前位置下降的距离,即
对应于前面登山例子中的某一步。
3.确定是否所有w_i梯度下降的距离都小于,如果小于于则算法终止,当前所有的即为最终结果,否则进入步骤4。
4.更新,即
更新完毕后继续回到步骤1。
以上步骤详细描述了梯度下降法的计算过程,当然在实际操作中,对例如步长\alpha等各个参数的选择会影响到收敛的速度以及最后的结果,这就需要结合多种方法和情况来进行考虑了。
梯度下降法还有批量梯度下降法(BGD)、随机梯度下降法(SGD)、小批量梯度下降法(MBGD)等变种,都是对梯度下降法的补充与改进。

牛顿法

牛顿法是用于方程求根的一种方法,其原理是利用泰勒公式,在处展开,且展开到一阶,即,求解可得,这里求得的更接近某个零点。进行迭代,进而推出
通过迭代,这个式子收敛于使得
需要注意的是,当函数有多个零点时,迭代收敛于哪个零点取决于迭代初值的选取。
除牛顿法外,数值求解方程根的方法还有简化牛顿法(平行弦法)牛顿下山法改进牛顿法等。
拟合问题可以转化为损失函数Q(\overrightarrow{w})的极小值问题,即方程求解问题Q'=0,与牛顿法求根方法相似。
假设需要求解f'(x)=0的根,把f对x泰勒展开到二阶形式
需要满足\Delta x为小量。当时,上式等价为,解得,从而推出迭代公式:
一般认为牛顿法可以利用到曲线本身的信息,比梯度下降法更容易收敛,即迭代更少次数。
以上讨论的是二维情况,高维情况的牛顿迭代公式是:
其中H为Hessian矩阵,定义为:
高维情况依然可以用牛顿迭代求解,但是问题是Hessian矩阵引入的复杂性,使得牛顿法迭代求解的难度大大增加,但是使用拟牛顿法可以用来解这个问题,不再直接计算Hessian矩阵,而是每一步的时候使用梯度向量更新Hessian矩阵的近似,拟牛顿法在这里不详细赘述。

高斯-牛顿法

高斯牛顿法实际上是牛顿法的在求解非线性最小二乘问题时的一个特例。由于要求解的损失函数为,该函数是趋于0的,因此可以直接使用高维的牛顿法来求解,其中迭代公式为:
其中g为梯度分量形成的矩阵,其分量为:
残差r为列矩阵
为Hessian矩阵,其分量
在高斯-牛顿法求解中的一个小技巧是,将二次偏导忽略,于是我们得到:
我们可以将以上两个式子改写为矩阵形式:
其中是Jacobian矩阵,代入牛顿法高维迭代方程的基本形式,得到高斯-牛顿法迭代方程:
其中

Levenberg-Marquardt法

LM算法是介于牛顿法与梯度下降法之间的一种非线性优化方法,对于过参数化问题不敏感,能有效处理冗余参数问题,使代价函数陷入局部极小值的机会大大减小,这些特性使得LM算法在计算机视觉等领域得到广泛应用。
在高斯-牛顿法中,。而在Levenberg-Marquardt法中,等式改为:
其好处就是在于可以对进行调节,即如果下降太快,使用较小的,使之更接近高斯-牛顿法;如果下降太慢,使用较大的,使之更接近梯度下降法。

代码实现

多项式拟合的拟合函数形式为
MATLAB实现:
 

参考文献

[1]  Wiki contributors. Wikipedia for Linear least squares. 2018.
[2]  李庆扬. 数值分析(第5版). 清华大学出版社. 2008.
[3]  yhao浩. 最小二乘法及其python实现. 2016.
[4]  北京交通大学. 数据拟合方法研究. 2014.
[6]  luoleicn. 牛顿法. 2011.
[8]  Jorge J Moré. The levenberg-marquardt algorithm: Implementation and theory. Lecture Notes in Mathematics, 630:105–116, 1978.
[9]  王子健. 优秀基金挑战. MCM related materials. 2017.
 

Loading...
公告
🎉大学生科技协会(2024)网站试运行🎉
-- 感谢您的支持 --- 科协获:
小平科技创新团队 (青少年科技创新领域国家级最高荣誉,全校唯一)
全国高校百强社团 中国大学生ican物联网创新创业实践教育基地 校十佳社团 “一院一品”智能空间 科技竞赛优秀组织单位
阅读科协简介了解更多吧!