拟牛顿法计算实现的除零问题
(北京理工大学数学系 北京 100081)
摘要: 拟牛顿法计算修正矩阵时,为避免由计算误差引起的分母为零问题的出现,对拟牛顿法的计算实现进行了一个改进,并给出数值实验结果验证其有效性。
关键词:拟牛顿法,算法改进,分母为零问题
中图分类号:O241
1 引 言
拟牛顿法( Quasi_newton Algorithms)是求解无约束优化问题的常用算法。许多约束问题也常常通过一些方法变为无约束问题用拟牛顿法求解。
拟牛顿法的理论基础和具体算法都是在连续系统下建立的,由于计算机是按照离散的方式运算处理。算法的实现有一个从连续系统到离散系统的转变过程,这个过程常因算法的截断误差(方法误差)和计算机的有效字长的舍入误差,导致算法在具体的编程实现中出现问题。于是必须对原有算法进行误差处理方面的改进。
在拟牛顿法中,搜索方向
(1)
计算得出,公式(1)中的修正矩阵 由下列的修正公式给出
(2)
其中
(3)
在实际的编程测试时,对部分经典算例的计算中经常在(2)式出现分母为零问题以致程序无法继续执行。由于上述分析,在用拟牛顿法进行计算时,需要对分母为零的情况进行处理。
2 算法改进
拟牛顿法的分母为零问题的出现主要是由于 引起的。 是由(3)式中相邻两次迭代导数做差得到的。如果相邻两点的导数值非常接近,可能是由于步长过小,导致两迭代点过于相近,使得导数值相近;或者步长不小,两迭代点的导数值也可能非常相近,那么 可能小于某个量。在这种情况下,将 代入(2)计算,在实际的编程测试中发现(2)的规模过大且计算复杂(主要指多个矩阵之间的运算),再加上计算机字长的有效性,数字低位的略去频繁发生。这样就使得在(2)式作除法运算时,分母已经被略为0。(即使这时分母没有被略为零,也会是一个非常小的除数,这种除法计算同样会带来巨大的误差。)
出现分母为零问题的情况多种多样,它不仅和算法结构、计算机有效字长、舍入原则有关,而且和函数本身的性质相关。由于这个问题过于复杂,如果按一般解决问题的思路(先找到分母为零现象的理论依据,再以此理论依据去解决分母为零问题。)是不现实的。
避免分母为零问题出现的方法有很多,比如改变公式(2)的运算次序以此减小截断误差对运算结果的影响,或提高计算机的计算精度等。但是他们都有一定的局限性,改变公式(2)的运算次序的方法,它可能对某些函数减小了截断误差对运算结果的影响,但也可能对某些函数增大了截断误差对运算结果的影响;计算机精度在短时间内不可能无限的提高,可能在 小于某个值后,这种方法也就无效了。
最直接方法是在计算(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位双精度)编程实现的。此算法实践表明,改进算法较好地处理了由于摄入误差产生的分母为零情况,且改进的算法具有较强的适应性。
参考文献:
[1] 解可新,韩立形,林友联.最优化方法.天津:天津大学出版社,1998.
[2] 袁亚湘,孙文瑜,最优化理论与应用.北京:科学出版社,2001.
[3] 徐成贤,陈志平,李乃成.近代优化方法. 北京:科学出版社,2002.
[4] J.M.Ortega.Iterative solution of nonlinear equations in
Several variables. New York : ACADEMIC PRESS, 1970.
[5] 徐世良.计算机常用算法.北京:清华大学出版社,1999.
[6] Mathwork.com Matlab 7 Programming