线性及非线性最小二乘问题

第5章线性与非线性最小二乘问题§5.1前言1801年,意大利天文学家朱赛普·皮亚齐发现了第一颗小行星谷神星.经过40天的跟踪观测后,由于谷神星运行至太阳背后，使得皮亚齐失去了谷神星的位置。随后全世界的科学家利用皮亚齐的观测数据开始寻找谷神星,但是根据大多数人计算的结果来寻找谷神星都没有结果.时年24岁的高斯也计算了谷神星的轨道.奥地利天文学家海因里希·奥尔伯斯根据高斯计算出来的轨道重新发现了谷神星.高斯使用的最小二乘法的方法发表于1809年他的著作《天体运动论》中.而法国科学家勒让德于1806年独立发现“最小二乘法”,但因不为时人所知而默默无闻.两人曾为谁最早创立最小二乘法原理发生争执.1829年,高斯提供了最小二乘法的优化效果强于其他方法的证明,称为高斯-马尔可夫定理.无约束最优化问题最小二乘的形式22111min()||()||[()],.22nmixRifxrxrxmn12()((),(),,())Tmrxrxrxrx()(1,,)irxim称为残量函数()irx是x的线性函数22111min()||()||[()],.22nmixRifxrxrxmn线性最小二乘问题()Tiiirxaxb12(,,,)Tiiiinaaaan为维向量1ibR为纯量211min[]2mTiiiaxb21min||||2Axb12(,,,)TmAaaa12(,,,)Tmbbbb§5.2线性最小二乘问题的解法解线性最小二乘问题设为nm阶矩阵,为m维mRb向量,线性最小二乘问题为求的最优解.,21min2bAxxq其中为向量的范数.将按范数2bAx展开得bbAxbAxAxxqTTTT2121minijAa对称矩阵AAT至少半正定任一最优解都是全局最优解bbAxbAxAxxqTTTT2121min凸优化问题,0nTTxRxAAx()0TTTxAAxAxAx当向量b属于矩阵A的像空间则存在nRxAxyyA,使nRx否则问题的最优02bAx值不为零.设为问题的最优解xnm0bAx不一定为零.bAxmn()rankAm一阶必要条件bAAxAxqTT是问题的最优解TTAAxAb法方程组xA列满秩正规方程组正定AATbAbAAAxTT1方程组的解唯一且可表示为TTAAAA1矩阵A的广义逆bbAxbAxAxxqTTTT2121min作为一个二次函数的无约束优化对方程组的求解，一般采用矩阵AAT的Cholesky分解，其中L下三角矩阵TLL如果矩阵正定对于线性最小二乘问题，我们可以利用其问题的特殊结构，设计更为有效的求解方法。对于线线性最小二乘问题形成的矩阵不宜采用AAT先计算矩阵AAT再对它进行分解的方法AAT的条件数是矩阵A的条件数的平方因为矩阵对称矩阵的三角分解定理设A为n阶对称矩阵,且A的所有顺序主子式均不为零,则A可唯一分解为TALDL其中L为单位下三角矩阵,D为对角矩阵.对称正定矩阵的三角分解或Cholesky分解设A为n阶对称正定矩阵,.TALL当限定L的对角元素为正时,这种分解是唯一的.则存在一个实的非奇异下三角矩阵L使得矩阵的QR分解,()),,,,mnmnnnnnnnARmnrAnAQRQRRRRmnQRRR是列满秩矩阵(存在分解式其中列正交矩阵,非奇异上三角阵.若限定阵对角元符号,则分解式是唯一的.当时正交阵非奇异上三角阵.0:......::...QRAnnnnnmnmQRA0:......::...nmmmnmQRA0...00:......::...为避免计算乘积矩阵再进行分解，可以采用AAT对增广矩阵作QR正交分解的方法。bA设A列满秩112110RAQRQQQRbRQbA其中Q为m×m阶正交矩阵1mnQR为n×n上三角矩阵()2mmnQR1R,TbQbbAAxATT,11111nTTTTbRbQRxRR112110RAQRQQQR1mnQR()2mmnQR1,TnbQb()TTTTAAxQRQRxRQQRxTRRx()TTTTAbQRbRQbTTTRRxRQb因此最优解x*可由三角方程组nbxR1经回代确定求线性最小二乘问题最优解正交分解算法如下:步1.对增广矩阵[Ab]作QR正交分解得和向量1R步2.取b为向量nbb的前n个分量形成的向量步3.用回代解方程组得解nbxR1*x在上述QR正交分解算法的分析过程中我们假定矩阵A的列线性无关.当矩阵A的列线性相关,即A不是列满秩时,矩阵A的QR正交分解是一个r×r非奇异的上三角矩阵112110[0]00RAQRQQQR1Rr=rank(A)＜n表示矩阵A的秩数表示由rb向量b的前r个分量组成的向量这时确定最优解的方程组1rRxbmnARTAA的特征值为1,2,,Aiiin称为矩阵的奇异值.利用矩阵A的奇异值分解1210rrn设A为m×n(m＞n)阶矩阵则存在m×m阶正交矩阵U和n×n阶正交矩阵V使得TTVUUSVA000其中S为m×n阶的块对角矩阵12r为r×r阶对角矩阵120r为矩阵A的奇异值A的秩为r＜n最小二乘法解不唯一1*TriiiiubxvTVSUb12(,,,)mUuuu12(,,,)nVvvvTTVUUSVA000极小范数最小二乘解规范最小二乘解范数最小2l其中,min**Xxxx是最小二乘问题所有解的集合.*X即在所有最小二乘解中当r=n时,最小二乘解唯一11()TTTTAVSUVUAAA1*Tniiiiubxv1TVUbAb由奇异值分解所确定的矩阵§5.3非线性最小二乘的Gauss-Newton法考虑非线性最小二乘问题miiTRxnmxrxrxrxfn12,,)]([21)()(21)(min其中Tm(x)，rx)),rxrxr(),(()(21为nRx的非线性函数，称为残量函数.如果函数r(x)二阶连续可微，则f(x)的一阶和二阶导数(Hesse矩阵)分别为miiiTxrxrxrxAxf1)()()()()(221()()()()()()()mTiiifxAxAxrxrxMxSxmiiiTxrxrxSxAxAxM12)()()(),()()(12()[(),(),,()]TmAxrxrxrx牛顿法11kkkkxxGg拟牛顿法kkkdHg221()()()()()()()mTiiifxAxAxrxrxMxSx()()1min()()()2kkTTkkqfxfxB()kxx()1min()()2kTTTkkkkkqfxrrAATTkkkkAAAr(1)()()()1()kkkkTTkkkkxxxAAAr算法(Gauss-Newton法)步1.给定解的初始估计置k=1;步2.如果满足精度要求,停止迭代;步3.解方程组()TTkkkKkAAAr得；步4.置k:=k+1后转步2;(1)()(),kkkxx()kx(1)x221()()()()()()()mTiiifxAxAxrxrxMxSx()()SxMx算法(阻尼Gauss-Newton法)步1.给定解的初始估计，置k=1;步2.如果满足精度要求,停止迭代;步3.解方程组并置(1)x()kx()TTkkkKkAAAr得；()();kks步4.沿方向进行线性搜索,确定步长ka置)()()1(kkkksaxx,k:=k+1后转步2.()ks§5.4信赖域方法信赖域方法是求解最优化问题的另一类有效方法，其最初的设计思想可追溯至Levenberg和Marquardt对Gauss-Newton法的修正。离最优解较远时)(kx和)(kkTkkAAM确定的点)()()1(kkkxx满足)1(kx(1)()()()kkfxfx)(k比较大,超出了有效的邻域不能保证而且对于越大的正数,解向量的长度越短,)(k因此,必可找到适当的值,使)()(kkx在)(kx的一个较小的邻域内,从而有).()()()()(kkkxfxf()TTkkkkAAIAr()1min12kTTkkkqsfgssBsksts.其中,,kkkxfgxxs是Hesse阵kxf2的近似0k为信赖域半径．根据模型函数()kqs对目标函数xf的拟合程度来调整信赖域半径.k对于问题(1)的解,ks定义比值：()()0kkkkkkkkkfxfxsAredrPredqqs它衡量模型函数()kqs与目标函数xf的一致性程度．实际下降量预测下降量注：(1)kr越接近于１，表明模型函数与目标函数的一致性程度越好，可以增大k以扩大信赖域．(2)0kr不接近于１，可以保持k不变．(3)kr接近于零或取负值，表明模型函数与目标函数的一致性程度不好，可以减小k以缩小信赖域．步1.给定控制迭代的参数值0,10,10max2121,初始点)1(x,以及初始信赖域半径max1,置k=1;步2.如果)(kx满足终止条件,停止迭代,取作为最优解的近似;)(kx步3.对解信赖域子问题得解)(k步4.计算)(Pr)(ked)()(kAred和kr步5.如果1kr,置kk1:,转步3;步6.置)()()1(kkkxx2max1min,,kkk2krkk)(转步2.1:kk