标签:
杂谈 |
分类: 数学 |
https://4.bp.blogspot.com/-WzBN3Jv8gS8/VrLXUB3j9_I/AAAAAAAAA0Q/MdJPTI11P9I/s1600/8.jpgLeast Square)" TITLE="非线性最小二乘学习(Nonline Least Square)" />
最小二乘问题可以描述如下:
f(x)为一个函数,输入一个n维的向量x得到一个一维的数,最小二乘问题可以简单的解释为:找到一个函数关系,使得对于给定的(x,y)拟合程度最高,或误差最小,一般当观测值多于变量个数时采用最小二乘方法进行误差分配,使得总体误差最小。
如下图所示为一个典型的非线性最小二乘问题:
http://s11/middle/002YDHC3zy6ZHianLFEba&690
得到所有离散点的最佳拟合曲线,而曲线的方程给定如下:http://s8/middle/002YDHC3zy6ZHiak7dB67&690其中x为要求的参数,t为自变量,函数为M,则根据函数关系可知,此函数有四个参数,分别为x1,x2,x3,x4,则假设有如下关系:





即使得f(x)最小,以上所描述的问题就是一个典型的非线性最小二乘的问题,现在主要问题转变为如何求解参数x(x1,x2,x3,x4)的问题。
对于以上这个具体问题,我们可以将其变得更加抽象化,对于一个给定的函数F为输入一个n维自变量得到一个一维的因变量的函数写作:

EQ(1)

如果是在整个定义域内获取最小值,则我们将其称之为求全局最优解,然而通常来说求全局最优解的的过程比较复杂,且难以实现,因此在实际使用的过程中会讨论一些比较简单的问题,即讨论函数F的局部最优解,则关系可以描述如下:
EQ(2)

即对于一个给定的小区间内找到http://s12/middle/002YDHC3zy6ZHigQXd97b&690使得在这个给定的小区间内函数F取得极小值。相比与全局的极小值问题,局部极小值通过一个极小值区间的约束将问题大大简化了。(对于全局极小值来说,不存在数值解法能够辨理整个区间进行求解,通常的迭代方法会陷入局部极小值陷阱,可能没办法得到真正意义上的全局最优解,当然如果能够证明函数在定义域上存在唯一的极值,则可以通过局部极小值的迭代算法求出其唯一极值作为全局的极小值)
2.极小值的存在性
对于代价函数F,我们认为F是可微的,将F进行二阶泰勒展开可得:EQ(3)

其中g为F的一节偏导,H为代价函数F的二阶偏导矩阵,假设http://s2/middle/002YDHC3zy6ZHiihgS5f1&690处的函数值,我们有如下极值存在定理,即:
EQ(4)
上述公式表示在若一点为极小值点,则其一阶导数一定恒等于0,当然反过来并不一定成立,此定理称之为极值存在定理,又可以称为驻点定理。当然一阶导为0只是极值存在的必要条件,而一阶导为0除了可能是极小值外还可能是极大值,可能是马鞍点,因此在这里我们引入函数的而二阶导,假设函数的一阶导为0且二阶导矩阵是正定矩阵,则说明函数的极小值存在。
3.局部极小值求法:
3.1.梯度下降法
根据极小值的存在性的条件我们可以得,得到极小值点首先需要证明其是极值点,再证明其是局部极小值点,而极值点的证明十分简单,一阶导为0的点,在实际求解过程中采用数值解法,求解到一阶导误差小于给定限差范围则说明找到了极值点,而梯度下降法就是基于这一思路,在各个方向上按照一定的步长获取梯度下降的方向认为其是极值点,最终找到布局极值点,梯度下降法伪代码描述如下:
简单的来说,左边的伪代码可以描述如下,首先获取初始位置,当然初始位置理论上是任意给定,但是要是初始位置给的比较好的话能够极大的减小运算量,减少算法的迭代次数。获取初始位置后获取初始方向,然后计算在初始方向上一阶导是否下降,如果是下降的,则迈出一步,如果不是下降的则另选方向,一直到此方向上获取不到更小值则换方向,进行判断,最后得到各个方向上都是最小的点认为是极小值点。
梯度下降法是我们能够想到的最简单的方法,其实这种方法还存在这比较大的问题:1.步长的选择,由于是给定步长,在距离极小值比较远的初始点收敛比较慢,而距离极值点比较近的位置有会错过极值点而产生震荡现象。2.需要对各个方向进行尝试,运算次数较多,运算复杂。
3.2.牛顿迭代法
不管何种迭代方法,其基本思路都是找到迭代的步长和迭代的方法,带着这个基本的思路去分析各个不同的迭代方法我们就可以知道,其方法的区别只在于迭代方向的选择和迭代步长的确定,选取较好的迭代方向和迭代步长能够极大的简化运算,减少迭代次数提高运算效率。对于一个非线性函数,进行泰勒展开并舍去极小值项:
EQ(5)
由于一阶导为0,则有如下等式:
在我们得到下一次的迭代位置为:EQ(7)

证明hn为梯度下降方向,假设二阶导H为正定矩阵,则对于任意非0向量u,有http://s3/middle/002YDHC3zy6ZHilZshc02&690,因此将等式两边都乘以hn的转置,有如下公式不等式:
EQ(8)
因此hn为梯度下降方向。
可以将牛顿迭代法和梯度下降法结合起来,主要改进的地方为方向和步长的选择:
EQ(9)

由以上对于牛顿迭代法的分析我们知道,对于二阶导为正的情况,牛顿迭代法可以保证是梯度下降方向,而对于二阶导小于0的情况另h等于hsd,其中hsd为梯度下降方向。
3.3.线性搜索确定步长
对于一个给定点x和梯度下降方向h,下一次迭代的步长为x到x+h,如下图所示:
其中http://s2/middle/002YDHC3zy6ZHinYl8de1&690,等式表示如果a足够小,则对迭代解算的情况就十分清楚了,通常会给出一个a的初始值,假设a的初始值为1,使用牛顿迭代法,则可能初夏如下三种情况:
1.a很小,目标值的增加很小,因此在这种情况下需要加大步长a;
2.a太大了,导致越过了极小值的区域,为了满足梯度下降的条件应该减小a的值;
3.a接近极小值,则认为a是可以接受的
线性搜索的迭代过程如下:
1.先给定一个初始的ak:
EQ(10)
2.如果要进行精确的线性搜索需要花费大量的时间因此没有必要获取准确的phia的极小值,只要搜索方向大致相同就可以了,此时这种方法称为松弛线性搜索方法,当我们获取到一个a值不在以上1,2两个条件之中,则使用严格的梯度下降法:
EQ(11)

以上方程保证第2(a不越过极值区域)这个条件得到满足,对于情况1(步长太小的情况)有如下公式:
EQ(12)

如果初始给定的a值满足上述条件,则接受其作为as,否则需要进行迭代找到合适的线性搜索方向,具体方法可以参考3.5章节的描述。
3.4置信域阻尼法
假设有模型L,对于函数F满足如下关系:EQ(13)

其中c是一个n维向量,B是一个n*n的对称矩阵,通常来说此方法的基本思路就是使用一种形式的函数去拟合另一个函数,最典型的情况为使用F的二阶泰勒展开式对F进行模拟,如2章中公式所述,L(h)是F的近似值,但是公式满足的前提是h是一个极小值。在置信域方法中我们假设已知一个正值delta,则模型在这个delta的范围内足够精确了,假设这个delta范围的中心为x,则步长可以通过下式获取:
EQ(14)

或者通过阻尼方法获取:
EQ(15)
http://s15/middle/002YDHC3zy6ZHiqb6d88e&690
则阻尼边量为u,

算法的中心思想为:
EQ(16)
http://s15/middle/002YDHC3zy6ZHiqEnzE5e&690
EQ(16)满足公式:
EQ(17)

则a等于1,否则a等于0,移动的步长的x+ah。然而实际上找到了x是不够的除非x正好为极小值点,通过对delta或u进行合适的修正,使用此方法的主要目的在于在下一次迭代是能够更好的接近于极小值点。
在h足够小的情况下可以认为L(h)对F(x+h)也是充分逼近,如果不能满足以上条件则说明h太大了需要减小h,另外,如果满足公式17的条件,则说明可以在下一步的时候扩大步长,从而达到减小迭代次数的目的。
关键的计算步骤是得到增益率:
EQ(18)

增益率是实际上和预测的值的下降的比值,通过构造使得分母为正,则如果此过程没有向梯度下降的方向移动则分子为负数,则整个比例为负数,这说明这一步并不是下降方向步长太大了,需要减小步长。如果比例为正数,则说明在下一步步长需要扩大,而步长如何进行扩大按照下面的规则进行处理:
EQ(19)

两个计算因子的选取很重要,一个适合的计算因子可以有效的避免delta出现震荡的情况,在以上公式中给出了默认的2和3,在实际进行处理的过程中可以根据实际问题给出适当的调节因子。
上式是采用公式14计算h得到的结果,而采用阻尼方式得到的h的结果为正好相反,当p的值比较小时需要增大阻尼因子,当p的值比较大时需要减小阻尼因子,则公式描述如下:
获取阻尼调节因子后通过步长为驻点的函数:
EQ(20)
这说明最佳步长为:
EQ(21)
http://s7/middle/002YDHC3zy6ZHitd7NA86&690
则根据L函数的定义13,可知
EQ(22)
其中B为二阶导矩阵,I为单位矩阵,h表示L的极小值,由于B是正定的,则B+uI为正定矩阵,因此可得hdm为L的极小值。
4.非线性最小二乘问题
对于一个向量函数f,
EQ(23)
由公式23可知,最小二乘的问题可以使用常规的最优化方法去解决,但是在解决的过程中会出现效率较低,迭代次数较多的问题,因此需要针对非线性最小二乘问题采用新的方法去解决。
首先函数f必须满足一定的条件:1.函数f二阶连续的函数,则函数f的二阶导是存在的,因此我们可以讲函数f按泰勒公式展开得:
EQ(24)
其中J为f的雅克比矩阵,包含了f对于各个x的一阶导数,对于F来说,由公式23可得,其关于x的一阶导为:
EQ(25)
http://s7/middle/002YDHC3zy6ZHiuuN6e36&690
则F的导数可以简化描述为:
EQ(26)
http://s1/middle/002YDHC3zy6ZHiuMrL230&690
由于在求取极小值时还要对函数的二阶导的正定性质进行判断,因此需要求出F关于x的二阶偏导,对于x的二阶偏导,http://s15/middle/002YDHC3zy6ZHiuZ9Qabe&690,则F关于x的二阶偏导可以描述为:
EQ(27)

下面举一个最简单的例子进行说明
Example1:假设f的形式为:http://s3/middle/002YDHC3zy6ZHiwKMD0c2&690,根据公式可以求出为极小值时x的值。
4.1高斯牛顿迭代法
高斯牛顿迭代方法是一个十分有效的方法,其主要思想为使用泰勒一阶展开作为一个极小值量,不断迭代逼近最优解的过程,公式描述如下:EQ(28)

将公式28代入公式23可得
EQ(29)
高斯牛顿迭代法的步长为使得L(h)最小的hgn,公式表示为:
EQ(30)
由公式26可得,L的一阶导和二阶导分别为:
EQ(31)
由公式26和公式31可得:http://s8/middle/002YDHC3zy6ZHiU44Z137&690,另外可以看到L的二阶导矩阵与步长变量h无关且如果J是列独立的矩阵,则说明二阶导恒大于0,L的极小值就是L的一阶导为0的点,则根据公式31中的第一个公式可得:
EQ(32)
梯度下降方向:
EQ(33)
因此计算hgn的步骤可以分为:
首先求解公式:

例2:

则可以求得F的一阶导和二阶导分别为:
http://s16/middle/002YDHC3zy6ZHiVCXXFbf&690
则对于lamda小于1,F的二阶导大于0,则0为F的一个极值点,实际上0为F的一个极小值点,对于f求雅克比行列式得:
对于传统的高斯牛顿迭代方法根据公式32有:
由于lamda不等于0且xk趋近于0,因此可以将上述公式简化为:
由上式可得如果lamda的绝对值小于1,则高斯牛顿迭代方法是线性收敛的,如果lamda小于-1,则传统的高斯牛顿迭代法不能找到极小值点,如果lamda为0,则迭代一次就可以找到极小值点,因为f函数是一个线性函数。
例3:对于本文最开始提出的线性拟合的问题,我们可以得到第i行的雅克比矩阵为:

同理我们可以获取其步长,得到求解结果。
4.2LM方法
LM方法结合了高斯牛顿方法和阻尼方法,给定步长hm的定义如下:EQ(34)

在这里引入了一个阻尼参数u,u的引入可以达到如下几个效果:
1.对于对于任意u大于0,则稀疏矩阵一定是正定的,这样可以保证hlm是一个梯度下降方向;
2.如果u的值很大,则有:

3.如果u值很小,则阻尼高斯方法和高斯方法近似,可以近似的得到一个二阶收敛的迭代过程
因此阻尼变量的选取既影响了方向又对步长有影响,对于u值的选择与元素值的大小有关系
EQ(35)
变量t是用户给定的,在迭代过程中u的大小是自适应的变化的,,对于u的大小的变化是根据一个控制变量来进行处理得到的:
EQ(36)
分母是一个增益的预测值,而增益的预测值求解的公式如下:
EQ(37)
由公式37可得,L0-Lhm必定是一个正值,如果控制变量的值比较大,则说明L,是一个比较好的拟合,我们可以减少u所以下一个LM算法更加接近于高斯牛顿过程,如果控制变量比较小甚至可能是负值,这说明L是一个比较差的你和,在此情况下我们可以增加u的值,使得迭代的步骤更加接近最大梯度下降方向。而迭代的终止条件为:
如果F的一阶导为0,或者小于一个给定的变量时停止迭代,或达到最大迭代次数时停止迭代。
LM算法的具体过程为:
对于非线性最小二乘的问题以及其求解方法就介绍到这里,其实除了上述迭代方法外还有很多迭代方法,和迭代处理的技巧,其实总结以下,不管是何种迭代方法其目的只有一个,求出极值点,而各种不同的迭代方法主要区别在于迭代的梯度下降方向的选取和步长的选择,一般来说在距离最优解距离比较远的位置需要采用比较大的步长更快的接近最优解,减小迭代次数,在比较接近最优解时需要减小步长以避免出现震荡的现象。把握好了这些原则,在看那些论文对迭代算法进行改进时就比较清晰,知道算法到底是从哪些角度进行了改进,以及改进的方法和效果如何,如果对于迭代算法有兴趣,可以进一步参看以下论文:
G. Golub and V. Pereyra, (2003): Separable Nonlinear Least Squares: The Variable Projection Method and its Applications, Inverse Problems, 19, ppR1–R26.
J.E. Dennis, Jr. and R.B. Schnabel (1983): Numerical Methods for Unconstrained Optimization and Nonlinear Equations, Prentice Hall.
K. Madsen (1988): A Combined Gauss-Newton and Quasi-Newton Method for Non-Linear Least Squares. Institute for Numerical Analysis (now part ofIMM), DTU. Report NI-88-10.
H.B. Nielsen (1999): Damping Parameter in Marquardt’s Method. IMM, DTU. Report IMM-REP-1999-05. Available as http://www.imm.dtu.dk/»hbn/publ/TR9905.ps
M.J.D. Powell (1970): A Hybrid Method for Non-Linear Equations. In P. Rabinowitz(ed): Numerical Methods for Non-Linear Algebraic Equations, Gordon and Breach. pp 87ff.
P.L. Toint (1987): On Large Scale Nonlinear Least Squares Calculations. SIAM J.S.S.C. 8, pp 416–435.