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; while abs(xa-xb)>e, k=k+1; xa=xb; xb=xa-subs(f,xa)/subs(diff(f),xa); end x=xb; b、雅克比迭代法程序 function [x,n]=Jacobi_Solve(A,b,x0,eps,varargin) %求解线性方程组的迭代法 %A为方程组的系数矩阵 %b方程组的右端项数字的向量组 %eps为精度要求，默认值为1e-5 %varargin为最大迭代次数，值为100 %x为方程组的解 %n为迭代次数（nargin） if nargin==3 eps=1.0e-6; M =200; elseif nargin<3 error return
高斯牛顿法的讲解和在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为例： #include<stdio.h> #include <math.h> //包含这个头文件，后面使用fabs void main() { double x=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,不是y printf("%f",x); }
如果我们的测量忽略噪声的影响，并且参数选择正学的话，那么这个向量的长度应该一直是1,无论它 指向那个方向，也就是
2 2 a2 x + ay + az = 1
(10)
或者是 ( x − mx δx )2 + ( y − my δy )2 + ( z − mx δx )2 =1 (11)
为了矫正我们的模型，我们只需要寻找参数mx , my , mz , δx , δy and δz ，这样这个等式就一直成立了。但 是
1
预备知识
1.3
常用迭代法的C语言编程
5
L=-tril(A,-1);%求A的下三角矩阵 U=-triu(A,1);%求A的上三角矩阵 G=(D-L)\U; f=(D-L)\b; x=G*x0+f; n=1;%迭代次数 while norm(x-x0)>=eps x0=x; x=G*x0+f; n=n+1; if(n>=M) disp(’迭代次数太多，可能不收敛。 ’); return; end end d、逐次超松弛迭代程序 function [x,n]=SOR_Solve(A,b,x0,w,eps,M) %求解线性方程组的迭代法 %A为方程组的细数矩阵 %b为方程组的右端项数字的向量组 %x0为迭代初始化向量 %w为松弛因子 %eps为精度要求，默认值为1e-5 %M为最大迭代次数，值为100 %x为方程组的解 %n为迭代次数 if nargin==4 eps=1.0e-6; M =200; elseif nargin<4 error return; elseif nargin==5 M=200; end if(w<=0||w>=2) error return;
其中t是单个变量，t = (t1 , t2 , t3 , . . . , tn )，今有n组实验观测值(tk , Ck ), k = 1, 2, . . . , n，在最小二乘意 义下确定m 个参数a1 , a2 , a3 , . . . , am。下面介绍一般解法步骤。 首先，给ai (i = 1, 2, . . . , m)一个初始值，记为ai ，并记初值与真值之差ϵi （未知） ，这时有 ai = ai + ϵi 若知则可求ai ，在ai 附近作Taylor级数展开并略去ϵi 的二次以上各项得 f (tk , a1 , a2 , . . . , , a(0) m ) ≈ fk0 + 其中fk0 = f (tk , a1 , a2 , . . . , am ), t = tk ∂f (t, a1 , a2 , . . . , am ) ∂fk0 = ∂ai ∂ai a1 = a1 . . .
1.2
高斯牛顿迭代法
高斯牛顿算法是一种解决非线性最小二乘法问题的方法，主要是为了解决非线性最小二乘拟合的 1
2
1
预备知识
问题。 它可以被看做是牛顿法的一种进化，为了寻找函数的最小值。 和牛顿法不同的是，高斯牛顿算 法只能被用于算平方和的最小值，但是它也有自身的优点：他不用计算实际编码中很难实现的函数二 阶导数的值。 高斯―牛顿迭代法的基本思想是使用泰勒级数展开式去近似地代替非线性回归模型，然后通过多 次迭代，多次修正回归系数，使回归系数不断逼近非线性回归模型的最佳回归系数，最后使原模型的 残差平方和达到最小。 高斯―牛顿法的一般步骤为：1、适当选择初始值（非常重要）2、泰勒级数展 开式3、估计修正因子4、精确度的检验5、重复迭代。一般非线性参数的确定，通常采用逐次逼近的方 法，即逐次“线性化”的方法，举例如下： 设某数学模型中待定参数a1 , a2 , a3 , . . . , am ，求得数学模型与时间的关系函数式为： C = f (t, a1 , a2 , a3 , . . . , am ) (1)
2.2
运用高斯牛顿迭代法算但是我们只用一个观测站从而只有一个方程。 这样我们会有无数多个解，怎 样选择一个？ • 假如我们使用6个观测站的话，我们可以得到6个未知数的唯一解，观测站尽量越不同越好这样才 不容易产生误差。但是我们知道读取数据的时候会有噪声，所以解并不是一个确定的值。 • 假如我们使用超过6个观测站的话，那么我们可以捕获到噪声的信息，但是我们得到了大于6个的 方程，但是我们只有6个未知数，所以是无解的。 这样看来我们只有选择正确的观测站数量，并且合理的选择参数，从而尽可能的使误差减小。 误差的 定义式如下 ( σi = x − mx δx )2 + ( y − my δy )2 + ( z − mx δx )2 −1 (i = 1, 2, . . . , N ) (12)
(0) (0) (0) (0) (0) (0) (0) (0) (0)
(2)
∂fk0 ∂fk0 ∂fk0 ϵ1 + ϵ2 + . . . + ϵm ∂a1 ∂a2 ∂am
(3)
(4)
am = am
k0 当ai 给定时，fk0 ， ∂f ∂ai 均可由t算得。
(0)
(0)
其次，列出正规方程（线性方程组） ，为了方便理解我们不妨设参数为4个即m = 4，则应列出以下 方程 b11 ϵ1 + b12 ϵ2 + b13 ϵ3 + b14 ϵ4 = B1 b21 ϵ1 + b22 ϵ2 + b23 ϵ3 + b24 ϵ4 = B2 b31 ϵ1 + b32 ϵ2 + b33 ϵ3 + b34 ϵ4 = B3 b41 ϵ1 + b42 ϵ2 + b43 ϵ3 + b44 ϵ4 = B4 其中bij (i, j = 1, 2, . . . , m)为方程系数，Bi 为常数。它们的表达式如下： bij = Bi =
4 elseif nargin==5 M =varargin{1}; end D=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;%迭代次数 while norm(x-x0)>=eps x0=x; x=B*x0+f; n=n+1; if(n>=M) disp(’迭代次数太多，可能不收敛。 ’); return; end end c、高斯_赛德尔迭代程序 function [x,n]=gaussseide(A,b,x0,eps,M) %求解线性方程组的迭代法 %A为方程组的细数矩阵 %b为方程组的右端项数字的向量组 %eps为精度要求，默认值为1e-5 %M为最大迭代次数，值为100 %x为方程组的解 %n为迭代次数 if nargin==3 eps=1.0e-6; M =200; elseif nargin==4 M=200; elseif nargin<3 error return; end D=diag(diag(A));%求A的对角矩阵