第五章最小二乘问题

1、第五章 最小二乘问题5.1 引言 在数字处理中经常遇到寻求回归方程的问题，即根据一在数字处理中经常遇到寻求回归方程的问题，即根据一组实验数据建立两个或多个物理量（俗称因素）之间的在组实验数据建立两个或多个物理量（俗称因素）之间的在统计意义上的依赖关系式。统计意义上的依赖关系式。 这类问题的数学模型如下：这类问题的数学模型如下： 设物理量设物理量 y 与物理量与物理量 t1,t2,tl 之间的依赖关系式，设之间的依赖关系式，设其方程为：其方程为： y=F(t1,tl,x1xn) (1) 其中其中 x1xn为待定参数。我们的问题是如何通过为待定参数。我们的问题是如何通过m(n)个实验点个实验点 t

2、1(i) ,t2(i) ,tl(i), y(i)T i=1,2m 确定（确定（1）中）中n个个参数参数x1,x2xn.从而建立从而建立回归方程回归方程。 ),(1)()(2)(1)(niliiixxtttFy 对于实际一组参数对于实际一组参数x1,x2xn的值，（的值，（1）给出）给出l+1维空间中的一个超曲面。第维空间中的一个超曲面。第i个实验点个实验点（ t1(i) ,t2(i) ,tl(i) )在（在（1）中就确定超曲面上一个点即）中就确定超曲面上一个点即相应的函数值相应的函数值：这个函数值这个函数值 与测量值与测量值y(i)之差的绝对值之差的绝对值)(iy)(21)()(2)(1,in

3、iliiyxxxtttF 就是第就是第i个实验点到该曲面的一种个实验点到该曲面的一种“距离距离”。为计算方便，通常把为计算方便，通常把)2(),(),(21)(1)()(121miinilinyxxttFxxxs作为作为m个实验点到该曲面个实验点到该曲面“总距离总距离”的度量的度量。 如何选择参数如何选择参数 x1xn 使（使（2）达到极小这就是最）达到极小这就是最小二乘法问题。上述问题用向量形式记为：小二乘法问题。上述问题用向量形式记为：2-)(minYxF则mTmRxfxfxf)(),()(1( )( )( )( , ),iiif xF txy若i=1m.nxR其中i=1mTmxtFxtF

4、xtFxF),(),(),()()()2()1(=myRlTiliiiRtttt)()(2)(1)(,则上面问题可记为则上面问题可记为:min f(x)Tf(x) (3) (3) 即为最小二乘法问题一般形式。即为最小二乘法问题一般形式。 当当f(x)为线性向量值函数时，称（为线性向量值函数时，称（3）为线性最小）为线性最小二乘法问题。二乘法问题。 否则，原问题称为非线性最小二乘法问题否则，原问题称为非线性最小二乘法问题。mininmnnxxxfxxxfxxxfxxxf121221212211,min0,0,0,nmmmnnxfxfxfxfxfxfxfxfxfxfA212221212111)(

5、（3）是有是有n个变量的无约束极小化问题，一般可个变量的无约束极小化问题，一般可以用前面介绍的最优化方法求解。考虑到（以用前面介绍的最优化方法求解。考虑到（3）的特）的特殊形式，可以考虑更有效、更简单的方法求解。殊形式，可以考虑更有效、更简单的方法求解。f 的的Jacobi矩阵：矩阵：21212( )( )( )( )2( )( )( )( )2( )( )2( )( )( )( )( )( )2( ) ( )2 ( )()()()2() ()2 ()TTmTiiimiiiTkkkTkkkF xf xf xg xf xf xF xHessianG xf xf xfxf xS xfxf xG x

6、Ax A xS xNewtonf xxxf xAxA xS x 则的梯度向量而的矩阵为：若令则有先考虑无约束最优化的法：则在此处有12()()() ()()() ()TkkkTTkkkkkkkkkPf xf xAxA xS xPAxf xXXP 即 主要计算量是主要计算量是 Sk 的计算，尽管的计算，尽管Sk对称，也包含对称，也包含(1/2)mn(n+1)个二阶偏导数，但个二阶偏导数，但Hesse矩阵中第一项只含矩阵中第一项只含一阶导数的信息。因此为简化计算，我们或者忽略一阶导数的信息。因此为简化计算，我们或者忽略Sk，或者用一阶导数的信息逼近或者用一阶导数的信息逼近Sk。 由（由（4）可知，

7、当）可知，当 接近接近0或或 fi(x)接近线性从而接近线性从而 接接近于近于0，此时才可以忽略，此时才可以忽略 Sk，因此这类算法又称为，因此这类算法又称为最小余最小余量算法量算法。 而称逼近而称逼近Sk的一类算法为的一类算法为大余量算法大余量算法。2( )if x5.2 线性最小二乘法问题的解法则（则（3）为：）为：min|Ax-b|2 (mn) (6) 定理定理1 x*是（是（6）的极小点的充要条件是）的极小点的充要条件是x*满足向量组：满足向量组：ATAx*=ATb (7)当当f(x)取线性形式取线性形式 即即f(x)=Ax-b.A是是mn矩阵，矩阵，mRb证证：2*.( ) |(,)

8、2( )22*( )( *)0TTTTTTTTF xAx bAx b Ax bXA AXb AXb bF xA AXA bxF xF xA AXA b必要性对求导为：若是的极小点。则必有由此得 称形如（称形如（7）的方程组为最小二乘问题（）的方程组为最小二乘问题（6）的法方程组）的法方程组. 可见求解线性最小二乘问题等价于求解它的法方程组可见求解线性最小二乘问题等价于求解它的法方程组.）的极小点。是（可见则因为计算对任意向量）即：满足（充分性：设6- -.0-0)-(2-)(-.0)-(7*2*2222*22*2*2*XbAXbAVAZAZbAXbAXAZAZbAXbZXAbAVRZXVbAX

9、AxTTnT0,2AVAVAVRVn有故故 至少是半正定的至少是半正定的.AA又因为又因为推论推论1 当当 （8） 为（为（6）的唯一最小二乘解）的唯一最小二乘解.bAAAXnrankATT1)(*,则0detAAAATT正定推论推论3mrankAAAT正定推论推论2nrankAAA正定定理定理2 A是是mn矩阵（矩阵（mn)则则 但实际中实际中A的秩不可能事先知道的秩不可能事先知道,而求而求rankA与求解与求解线性方程组几乎等价线性方程组几乎等价,因而因而 的正定性也不能事先确的正定性也不能事先确定定,因此因此(8)仅具有理论意义仅具有理论意义,而且即使而且即使 正定正定,也不用也不用(8

10、)去求去求(6)的解的解. AATAATAAT 在在 正定性不能确定时正定性不能确定时,可用可用QR分解法求解分解法求解.AAT 在在 已知正定时已知正定时,一般可用一般可用Cholesky分解求线分解求线性方程组性方程组(7). 下面讨论下面讨论 (9)的非线性最小二乘问题的求解的非线性最小二乘问题的求解. 2)(min)()(min)(minxfxfxfxST 与与Newton法不同的是法不同的是, 至少半正定的至少半正定的,而当而当 满秩时满秩时, Pk 为一下降方向为一下降方向,但与但与Newton法一样法一样,并不能保证并不能保证 ,而而且只有当初始点充分接近于极小点且只有当初始点充

11、分接近于极小点 时时,才能使算法收敛才能使算法收敛kTkAAkAkkFF1*X在在Newton迭代公式中迭代公式中(5)忽略忽略 得得:kSkkkkTkkkTkPXXfAPAA1而而(10)称为称为法方程法方程(11)从而构成从而构成(9)的的Gauss-Newton法法:)()()()(kkTkkkTxfxAPxAxA(10) 5.3 Gauss-Newton法定理定理3：设 的局部极小点，且 为正定，若GN迭代法产生点列 收敛于X*，则当G(x)与 在X*的邻域内Lipschitz连续时，有 其中ffxFXmiCxfT)(*),1()( :2为)(kx1)()(xAxAT2(1)1()()

12、( *)( *)( *)()kTkkhA xA xS xhoh*)()(xxhkk*)(*)(xAxAT定理定理4：设 满足上面定理条件，若在极小点X*处F(x*)=0,则当初始点充分接近于X*时，GN法产生的点列 收敛于X*，且收敛速度为2阶的.miixfxF12)()()(kx22222232min( )(1)(1) ()1000010( )(1),(1)( )( )(1,21) ,( )( )1(21)( ) ( )2322TTTTF xxxxxRxFFGNF xxxxA xf xxA xA xxAx f xxxx 当时，为极小点，这是因为： （ ）， （ ）下面用算法求的极小点，因为:

13、f(x)=(例: 求解最小二乘问题: 2321223221() ()(10)() ()23221 (21)22(*)1 (21)0*0kkTkkkkkkkkkkkkkkkkA xf xPAxA xxxxxxpxxxxxxGNxx由可得：可见，若，则法只须一步迭代即找到最优解 而这个现象很清楚：因=0时,f1,f2均为线性函数,从而GN法即化为Newton法. 而当0时,若 充分小,则(*)可改写为: 因而可见其收敛速度是线性的,其原因是 而0 ,则(10)化为: (12) 当 即使奇异时,只要取值充分大,则可以保证 正定,从而(12)有解,这个解与参数有关,记为 .kTkkTkfAPIAA)(

14、kTkAAkTkAAIAAkTk)(kP 特别地=0时.(12)即为(10).并且当 非奇异时, 存在.即为GN方向:AATk)0(kPkTkkTkkfAAAP1)()0( 当 奇异时, 就不存在.这时考虑增大,使得 的每个分量与相比都充分小,以致可以忽略时,(12)即为 (13)kTkAA)0(kPkTkAAkTkkkTkkfAPfAP1)(.)(即 这表明很大, 接近于F(x)的 处的负梯度方向,因此可以想象,从零连续增大到无穷大时, 将从F(x)在 处的GN方向连续地转向负梯度方向. )(kPnx)(kPkx 若F(x)是x的二次函数,即线性最小二乘问题.则GN方法(11)相当于Newt

15、on法,因而对于任意 一次迭代即可,找到F(x)的极小点x*,即这表明F(x)为二次函数时,F(x)在 处GN方向不但指向极小点x*,而且步长 恰好等于 .kx)0()(*1kkkTkkTkkPxfAAAxxkx)0(kPkxx * 当增大时, 偏离 ,向负梯度方向 靠近. 可见起着使步长缩小的作用,称为阻尼作用.)(kPkTkkTkkfAAAP1)()0(kTkfA0)(limkP满足 (15)(kkpkTkkkkkTkfApIAA)()()(1kkkkpxx(14) 由此构成阻尼最小二乘迭代格式: 可见每次迭代怎样选择阻尼因子是很关键问题. 不同阻尼因子的选择法就导致不同阻尼最小二乘法.

16、下面介绍的一种边迭代边调整的方法(即大了就减少,小了就增大) 算法过程:mnkijkkkikijjiTmaAxxfanjmixxxfxAJacobixfxfxf)(,)(1,1. 2),10, 4 , 2(),01. 0(. 1)()()(),.()()()(0001构成计算对初点可取整数可取选定降及其已知 为避免 奇异， 必须为正，且 越大，（15）的条件越好，但步长 越小，收敛速度减慢，因此 不能选的太大 . 但 太小,由于F(x)的复杂性,不能保证(14)迭代是下降的,即可能出现IAAkkTkk)(kkpkk)()(kkkkxFpxFkmikikkkkkkkkTkkkkkkxfFpxxp

17、pgpIAArFFXXH1121111)(,)(. 68, 6,:,. 5.,*. 4计算记为求出由否则转转若终止，否则满足，则若终止后则mikikkkkkkkkTkxfFpxxpCholeskypgpIAA11211)(,)(. 3计算设解为法求）（可用解出由方程组mikiknkjkmijkikikjxfFggxxfxfgnj121)(1)()()()()(,1计算构成向量计算对.,)(:,)(,.2 , 1,. 8)(11(:)12)(1)()(1)()(,)(1为止使的值直到有某个再计算设解为并求解计算则对若klkmikikikkikikkikkTkkiikkkFFlixfFpxxpgpIAAriFF=+=+=+=+. 21,. 9)(11)(1转+=+kkxxlk