GaussNewton

高斯牛顿法的讲解和在APM上的应用李智迅1预备知识1.1牛顿迭代法牛顿法是一种求解方程f(x)=0的方法。多数方程不存在求根公式，从而求精确根非常困难，甚至不可能，从而寻找方程的近似根就显得特别重要。迭代法是用某个固定公式反复地计算，用以校正方程所得根的近似值，使之逐步精确化，最后得到的精度要求的结果。牛顿迭代法是用切线来逼近零点的，收敛的速度很快，但要求的条件也高。首先要有一个区间，在这个区间的端点函数值反向，其次第一次迭代点不能随意取，否则第一次迭代后的点可能跑出原先的区间，收敛性就不一定能保证了（即有的情况也能有收敛性，有的情况就没有收敛性了，全看函数的性能）。比如：设r是f(x)=0的根，选取x0作为r初始近似值，过点(x0,f(x0))做曲线y=f(x)的切线L，L的方程为y=f(x0)+f′(x0)(x−x0),求出L与x轴交点的横坐标x1=x0−f(x0)/f′(x0),称x1为r的一次近似值，过点(x1,f(x1))做曲线y=f(x)的切线，并求该切线与x轴的横坐标x2=x1−f(x1)/f′(x1)称x2为r的二次近似值，重复以上过程，得r的近似值序列Xn,其中Xn+1=Xn−f(Xn)/f′(Xn),称为r的n+1次近似值。上式称为牛顿迭代公式。以f(x)=2x3−4x2+3x−6为例：#includestdio.h#includemath.h//包含这个头文件，后面使用fabsvoidmain(){doublex=1.5,y,y1;do{y=2*x*x*x-4*x*x+3*x-6;y1=6*x*x-8*x+3;x=x-y/y1;}while(fabs(y/y1)1e-6);//是y/y1,不是yprintf(%f,x);}1.2高斯牛顿迭代法高斯牛顿算法是一种解决非线性最小二乘法问题的方法，主要是为了解决非线性最小二乘拟合的121预备知识问题。它可以被看做是牛顿法的一种进化，为了寻找函数的最小值。和牛顿法不同的是，高斯牛顿算法只能被用于算平方和的最小值，但是它也有自身的优点：他不用计算实际编码中很难实现的函数二阶导数的值。高斯―牛顿迭代法的基本思想是使用泰勒级数展开式去近似地代替非线性回归模型，然后通过多次迭代，多次修正回归系数，使回归系数不断逼近非线性回归模型的最佳回归系数，最后使原模型的残差平方和达到最小。高斯―牛顿法的一般步骤为：1、适当选择初始值（非常重要）2、泰勒级数展开式3、估计修正因子4、精确度的检验5、重复迭代。一般非线性参数的确定，通常采用逐次逼近的方法，即逐次“线性化”的方法，举例如下：设某数学模型中待定参数a1,a2,a3,...,am，求得数学模型与时间的关系函数式为：C=f(t,a1,a2,a3,...,am)(1)其中t是单个变量，t=(t1,t2,t3,...,tn)，今有n组实验观测值(tk,Ck),k=1,2,...,n，在最小二乘意义下确定m个参数a1,a2,a3,...,am。下面介绍一般解法步骤。首先，给ai(i=1,2,...,m)一个初始值，记为a(0)i，并记初值与真值之差ϵi（未知），这时有ai=a(0)i+ϵi(2)若知则可求ai，在a(0)i附近作Taylor级数展开并略去ϵi的二次以上各项得f(tk,a(0)1,a(0)2,...,,a(0)m)≈fk0+∂fk0∂a1ϵ1+∂fk0∂a2ϵ2+...+∂fk0∂amϵm(3)其中fk0=f(tk,a(0)1,a(0)2,...,a(0)m),∂fk0∂ai=∂f(t,a1,a2,...,am)∂ait=tka1=a(0)1...am=a(0)m(4)当a(0)i给定时，fk0，@fk0@ai均可由t算得。其次，列出正规方程（线性方程组），为了方便理解我们不妨设参数为4个即m=4，则应列出以下方程b11ϵ1+b12ϵ2+b13ϵ3+b14ϵ4=B1(5)b21ϵ1+b22ϵ2+b23ϵ3+b24ϵ4=B2b31ϵ1+b32ϵ2+b33ϵ3+b34ϵ4=B3b41ϵ1+b42ϵ2+b43ϵ3+b44ϵ4=B4其中bij(i,j=1,2,...,m)为方程系数，Bi为常数。它们的表达式如下：bij=n∑k=1∂fk0∂ai·∂fk0∂aj(6)Bi=n∑k=1∂fk0∂ai(Ck−fk0)(7)1.3常用迭代法的C语言编程3上述方程组采用高斯消元法求解ϵi，继而算出ai的值。最后，将算得的ai作为初值，重复上述步骤，反复迭代和修正ϵi的值，直至∥ϵi∥小于允许误差或符合残差平方的要求。残差平方和是检查计算是否达到要求的重要指标，其表达式如下所示ξ=n∑k=1[f(tk,a1,a2,...,am)−Ck]2(8)上述高斯-牛顿迭代法对于初始值依赖高，常先用残数法求出初始值。若初始值偏离真值太远，往往迭代数次后又会发散，导致线性回归失败。1.3常用迭代法的C语言编程a、牛顿迭代法程序：function[x,k]=newton(f,p0,e)%求f(x)=0在给定p0的根。%f为所求的函数f(x)，p0为迭代初始值，e为迭代精度。k为迭代次数%diff(f)为对函数求导，subs是赋值函数，用数值替代符号变量替换函数例xa=p0;xb=xa-subs(f,xa)/subs(diff(f),xa);k=1;whileabs(xa-xb)e,k=k+1;xa=xb;xb=xa-subs(f,xa)/subs(diff(f),xa);endx=xb;b、雅克比迭代法程序function[x,n]=Jacobi_Solve(A,b,x0,eps,varargin)%求解线性方程组的迭代法%A为方程组的系数矩阵%b方程组的右端项数字的向量组%eps为精度要求，默认值为1e-5%varargin为最大迭代次数，值为100%x为方程组的解%n为迭代次数（nargin）ifnargin==3eps=1.0e-6;M=200;elseifnargin3errorreturn41预备知识elseifnargin==5M=varargin{1};endD=diag(diag(A));%求A的对角矩阵L=-tril(A,-1);%求A的下三角矩阵U=-triu(A,1);%求A的上三角矩阵B=D\(L+U);f=D\b;x=B*x0+f;n=1;%迭代次数whilenorm(x-x0)=epsx0=x;x=B*x0+f;n=n+1;if(n=M)disp('迭代次数太多，可能不收敛。');return;endendc、高斯_赛德尔迭代程序function[x,n]=gaussseide(A,b,x0,eps,M)%求解线性方程组的迭代法%A为方程组的细数矩阵%b为方程组的右端项数字的向量组%eps为精度要求，默认值为1e-5%M为最大迭代次数，值为100%x为方程组的解%n为迭代次数ifnargin==3eps=1.0e-6;M=200;elseifnargin==4M=200;elseifnargin3errorreturn;endD=diag(diag(A));%求A的对角矩阵1.3常用迭代法的C语言编程5L=-tril(A,-1);%求A的下三角矩阵U=-triu(A,1);%求A的上三角矩阵G=(D-L)\U;f=(D-L)\b;x=G*x0+f;n=1;%迭代次数whilenorm(x-x0)=epsx0=x;x=G*x0+f;n=n+1;if(n=M)disp('迭代次数太多，可能不收敛。');return;endendd、逐次超松弛迭代程序function[x,n]=SOR_Solve(A,b,x0,w,eps,M)%求解线性方程组的迭代法%A为方程组的细数矩阵%b为方程组的右端项数字的向量组%x0为迭代初始化向量%w为松弛因子%eps为精度要求，默认值为1e-5%M为最大迭代次数，值为100%x为方程组的解%n为迭代次数ifnargin==4eps=1.0e-6;M=200;elseifnargin4errorreturn;elseifnargin==5M=200;endif(w=0||w=2)errorreturn;62高斯牛顿法在APM中的应用endD=diag(diag(A));%求A的对角矩阵L=-tril(A,-1);%求A的下三角矩阵U=-triu(A,1);%求A的上三角矩阵B=inv(D-L*w)*((1-w)*D+w*U);%inv为矩阵求逆f=w*inv((D-L*w))*b;x=B*x0+f;n=1;%迭代次数whilenorm(x-x0)=epsx0=x;x=B*x0+f;n=n+1;if(n=M)disp('迭代次数太多，可能不收敛。');return;endend2高斯牛顿法在APM中的应用2.1问题的提出在APM加速度计校准的环节中，单点校准法并不适用，因为这种方法我们只在一个方向的一个轴进行测量。我们需要分别寻找x,y,z三个轴上的0G点标记mx,my,mz和相应的灵敏度δx,δy,δz。其中G为重力加速度。所以我们需要有6个需要测量的元素。从单点测量法的测量数据上来看，测量精度是非常差的，所以我们需要更精确的数学方法，在无法直接算出确切结果的前提下，使理论值和真实值的残差平方和最小。无论我们怎样在读取数据前摆放我们的传感器，我们假设它是静止的并且在某些方向测量值是1G.如果我们假设x为X方向的读取值，y为Y方向的读取值，z为Z方向的读取值，那么加速度向量a=(ax,ay,az)的表达式如下所示：ax=x−mxδxay=y−myδyaz=z−mxδz(9)如果我们的测量忽略噪声的影响，并且参数选择正学的话，那么这个向量的长度应该一直是1,无论它指向那个方向，也就是a2x+a2y+a2z=1(10)或者是(x−mxδx)2+(y−myδy)2+(z−mxδx)2=1(11)为了矫正我们的模型，我们只需要寻找参数mx,my,mz,δx,δyandδz，这样这个等式就一直成立了。但是2.2运用高斯牛顿迭代法算出最小残差7•我们有6个未知数，但是我们只用一个观测站从而只有一个方程。这样我们会有无数多个解，怎样选择一个？•假如我们使用6个观测站的话，我们可以得到6个未知数的唯一解，观测站尽量越不同越好这样才不容易产生误差。但是我们知道读取数据的时候会有噪声，所以解并不是一个确定的值。•假如我们使用超过6个观测站的话，那么我们可以捕获到噪声的信息，但是我们得到了大于6个的方程，但是我们只有6个未知数，所以是无解的。这样看来我们只有选择正确的观测站数量，并且合理的选择参数，从而尽可能的使误差减小。误差的定义式如下σi=(x−mxδx)2+(y−myδy)2+(z−mxδx)2−1(i=1,2,...,N)(12)那么我们要找到合适的参数，使得所有的σi都很小。怎样评判综合误差最小呢？我们需要用到残差平方和，它的表达形式如下N∑i=1σ2i=[(x−mxδx)2+(y−myδy)2+(z−mxδx)2−1]2(i=1,2,...,N)(13)2.2运用高斯牛顿迭代法算出最小残差在下面代码中，我们使用向量beta来代替要计算的六个未知数mx=beta(1)(14)my=beta(2)mz=beta(3)1/δx=beta(4)1/δy=beta(5)1/δz=beta(6)constintNUM_READINGS=6;//theseconstantsdescribethepins.Theywon'tchange:constintxpin=A1;//x-axisoftheaccelerometerx轴加速度计constintypin=A2;//y-axisy轴constintzpin=A3;//z-axis(onlyon3-axismodels)z轴//datacollectionstructures采集数据unsignedint*data;//dynamicallyallocatedarrayfordatastorage.Don'tl