拟牛顿法计算实现的除零问题

拟牛顿法计算实现的除零问题

（北京理工大学数学系 北京 100081）

摘要: 拟牛顿法计算修正矩阵时，为避免由计算误差引起的分母为零问题的出现，对拟牛顿法的计算实现进行了一个改进，并给出数值实验结果验证其有效性。

关键词：拟牛顿法，算法改进，分母为零问题

中图分类号：O241

1 引   言

拟牛顿法( Quasi_newton Algorithms)是求解无约束优化问题的常用算法。许多约束问题也常常通过一些方法变为无约束问题用拟牛顿法求解。

拟牛顿法的理论基础和具体算法都是在连续系统下建立的，由于计算机是按照离散的方式运算处理。算法的实现有一个从连续系统到离散系统的转变过程，这个过程常因算法的截断误差（方法误差）和计算机的有效字长的舍入误差，导致算法在具体的编程实现中出现问题。于是必须对原有算法进行误差处理方面的改进。

在拟牛顿法中，搜索方向

（1）

计算得出，公式（1）中的修正矩阵 由下列的修正公式给出

（2）

其中

（3）

在实际的编程测试时，对部分经典算例的计算中经常在（2）式出现分母为零问题以致程序无法继续执行。由于上述分析，在用拟牛顿法进行计算时，需要对分母为零的情况进行处理。

2 算法改进

拟牛顿法的分母为零问题的出现主要是由于 引起的。 是由（3）式中相邻两次迭代导数做差得到的。如果相邻两点的导数值非常接近，可能是由于步长过小，导致两迭代点过于相近，使得导数值相近；或者步长不小，两迭代点的导数值也可能非常相近，那么 可能小于某个量。在这种情况下，将 代入（２）计算，在实际的编程测试中发现（２）的规模过大且计算复杂（主要指多个矩阵之间的运算），再加上计算机字长的有效性，数字低位的略去频繁发生。这样就使得在（２）式作除法运算时，分母已经被略为０。（即使这时分母没有被略为零，也会是一个非常小的除数，这种除法计算同样会带来巨大的误差。）

出现分母为零问题的情况多种多样，它不仅和算法结构、计算机有效字长、舍入原则有关，而且和函数本身的性质相关。由于这个问题过于复杂，如果按一般解决问题的思路（先找到分母为零现象的理论依据，再以此理论依据去解决分母为零问题。）是不现实的。

避免分母为零问题出现的方法有很多，比如改变公式（２）的运算次序以此减小截断误差对运算结果的影响，或提高计算机的计算精度等。但是他们都有一定的局限性，改变公式（２）的运算次序的方法，它可能对某些函数减小了截断误差对运算结果的影响，但也可能对某些函数增大了截断误差对运算结果的影响；计算机精度在短时间内不可能无限的提高，可能在 小于某个值后，这种方法也就无效了。

最直接方法是在计算（2）式时限制 的范数不能过小。在数值实验中最有效的改进方法是在计算修正公式（2）前，先对（3）得到的 进行一次判断，如果 小于一定的值,在计算修正公式（2）时很可能发生分母为零的情况，那么就得不到修正公式 。

由于分母为零问题的出现，就计算不出下次迭代所需的修正矩阵，按拟牛顿法的收敛原理只要以一个对称正定矩阵为初始修正矩阵，算法就会收敛。那么可以以单位矩阵为下次迭代的修正矩阵，继续迭代过程。此迭代过程生成的点列仍会收敛于一个局部最优点，原迭代过程的收敛性不会改变。

这种算法的改进与其它方法相比具有很强的适应性，在编程实践中发现它不会因为具体函数的不同而改变它的有效应，而且易于编程实现。

   具体算法流程图如下：

图1  改进算法流程

3数值实验

    作为数值计算实例，我们考虑Rosen Broke函数：

有唯一极小点（1，1）,目标函数值为0，

    在实际的测试中发现，如果以（-2，2）为迭代初始点，在迭代过程中极易出现分母为零问题。

测试环境Windows XP2  CPU P-M 1.6G   RAM 512M   Matlab7.0

参数设置：控制误差0.05, 求导差分0.0000001（中心差商）， 取0.001。

下表中initial 初始点，opt-point得到的最有优点，Iteration是指牛顿代次数,Obj_Value为最优点的函数值，time为CPU运行时间。

测试A：进退法，黄金分割法，拟牛顿（DFP修正公式）

Golden-finder进退法确定（golden_a, golden_b）黄金分割法初始区间.

 表1 计算结果（1）

	DFP-Golden
Initial	Opt-point	Iteration	Obj_Value	time
(-2,2)	1.0003    1.0006	15	0.0000	0.030 s

图2 函数图形（1）

表2  迭代过程数据

No.	1	2	3	4	5	6	7	8
X1	-2.0000	1.7065	1.7078	1.5576	1.4423	1.4057	1.3085	1.2372
X2	2.0000	2.9232	2.9176	2.4067	2.0545	1.9809	1.7006	1.5141
F	409.00	0.5115	0.5011	0.3489	0.2627	0.1671	0.1085	0.0839
No.	9	10	11	12	13	14	15
X1	1.1841	1.1218	1.0742	1.0512	1.0201	0.9986	1.0003
X2	1.4071	1.2531	1.1463	1.1065	1.0390	0.9968	1.0006
F	0.0365	0.0176	0.0114	0.0028	0.0007	0.0000	0.0000

如图2和表2所示第2迭代点是可能分母为零问题出现区。

 

测试B：Wolfe不精确一维搜索，拟牛顿（BFGS修正公式）

表3 计算结果（2）

	BFGS-Wolfe
Initial	Opt-point	Iteration	Obj_Value	time
(-2,2)	0.9887    0.9775	16	0.0000	0.025 s

图3  函数图形（2）

表4  迭代过程数据

No.	1	2	3	4	5	6	7	8
X1	-2.0000	-1.2158	-0.8274	-0.8274	-0.9220	-0.7402	-0.7402	-0.6633
X2	2.0000	2.1953	0.8184	0.8184	0.7561	0.4266	0.4266	0.4740
F	409.0000	56.3322	5.1307	5.1307	4.5792	4.4982	4.4982	2.8824
No.	9	10	11	12	13	14	15	16
X1	-0.4932	-0.4932	1.0265	1.0265	0.9822	0.9883	0.9883	0.9887
X2	0.2302	0.2302	0.2302	0.9434	0.9649	0.9776	0.9776	0.9775
F	2.2468	2.2468	1.2156	1.2156	0.0003	0.0002	0.0002	0.0001

如图3和表4所以第2，6，14迭代点为分母为零问题的出现区。

 

4小结

　　计算过程是在普通计算机上，用Matlab7.0（64位双精度）编程实现的。此算法实践表明，改进算法较好地处理了由于摄入误差产生的分母为零情况，且改进的算法具有较强的适应性。

 

参考文献：

［１］  解可新,韩立形,林友联.最优化方法.天津：天津大学出版社，1998.　

［２］  袁亚湘,孙文瑜,最优化理论与应用.北京：科学出版社，2001.

［３］  徐成贤,陈志平,李乃成.近代优化方法. 北京：科学出版社，2002.

［４］  J.M.Ortega.Iterative solution of nonlinear equations in

Several variables. New York : ACADEMIC PRESS, 1970.

［５］  徐世良.计算机常用算法.北京：清华大学出版社，1999.

［６］  Mathwork.com  Matlab 7 Programming