526互联

2023-07-22 《数值优化方法》-庞丽萍，肖现涛-无约束最优化（七）.md

发布时间 2023-07-22 17:00:46作者: XNEF

2023-07-22 《数值优化方法》-庞丽萍，肖现涛-无约束最优化（七）

数值优化方法Matlab牛顿法

在前面我们研究了共轭方向法和共轭梯度法，两种方法都有二次终止性，那么是否可以在每次迭代的时候都用一个二次函数去近似目标函数呢？这就是牛顿法的基本思想。

我们知道函数在处的二阶泰勒展开式为


其中. 考虑的极小点，即, 我们有

得到

即得到牛顿法迭代格式

其中下降方向为.

定理 1.11 设是二阶连续可微的, 黑塞矩阵在的局部极小点附近是利普西茨连续的(常数为), 正定. 则牛顿法产生的点列满足

  1. 若初始点充分靠近, 则;
  2. 迭代点列二阶收敛于;
  3. 二阶收敛于0.

证明：由于在附近利普西茨连续且正定，因此存在, 在上是可逆的, 满足


取, 当时，有

由于（这个不等式需要一定技术）

因此


即结论2成立.

由的取值范围可知


因此对任意的, 都有.

反复带入上式，得到


即结论1成立.

注意到


我们有

由于, 所以

由于牛顿法的步长始终为1，在某些情况下可能不收敛，因此可以结合一维线搜索方法得到全局收敛性的牛顿法（阻尼牛顿法）
算法8 阻尼牛顿法
Step 1. 取初始点, 给定精度, 置.
Step 2. 检验, 若成立，则停，否则转下步.
Step 3. 计算牛顿方向


Step 4. 求一维线搜索
,
解得, 令, , 转Step 2.

修订了前面数值梯度算法中的一个错误

  1. %% 数值梯度 
  2. function [gd] = My_Gradient(f, x) 
  3. gd = x; 
  4. epsil = 1e-5; 
  5. d = [-2* epsil, -epsil 0 epsil 2*epsil]; 
  6. tz = [x x x x x]; 
  7. fx = [0,0,0,0,0]; 
  8. for i = 1:length(x) 
  9. tx = tz; 
  10. tx(i,:) = tx(i,:) + d; 
  11. for h = 1:5 
  12. fx(h) = f(tx(:,h)); 
  13. end 
  14. gf = gradient(fx) / epsil; 
  15. gd(i) = gf(3); 
  16. end 
  17. end 
  18.  
  19. %% 数值黑塞矩阵 
  20. function [mz] = My_Hessian(f, x) 
  21. % 计算以当前点为中心点的一张网格 
  22. % 根据网格数据计算最终Hessian矩阵 
  23. n = length(x); 
  24. mz = zeros(n,n); 
  25. mx = zeros(5,5); 
  26. epsil = 1e-5; 
  27. d = [-2* epsil, -epsil 0 epsil 2*epsil]; 
  28. for i = 1:n 
  29. for h = 1:n 
  30. % 针对i h方向生成网格 
  31. for j = 1:5 
  32. for k = 1:5 
  33. y = x; 
  34. y(i) = y(i) + d(j); 
  35. y(h) = y(h) + d(k); 
  36. mx(j,k) = f(y); 
  37. end 
  38. end 
  39. my = mx; 
  40. for j = 1:5 
  41. my(j,:) = gradient(mx(j,:)) / epsil; 
  42. end 
  43. myy = gradient(my(:,3)) / epsil; 
  44. mz(i,h) = myy(3); 
  45. end 
  46. end 
  47. end 
  48.  
  49. %% 主函数 
  50. % Fletcher-Reeves 共轭梯度法 
  51. function [x xlog] = Newton(f, x0, epsilon) 
  52. % f 目标函数，函数句柄 
  53. % g 梯度函数 函数句柄 
  54. % epsilon 精度要求 
  55. % method 线搜索方法 
  56. k = 0; 
  57. iter = 0; 
  58. maxIt = 1e4; 
  59. n = length(x0); 
  60. while iter <= maxIt 
  61. d1 = - inv(My_Hessian(f,x0)) * My_Gradient(f, x0); 
  62. [alpha tx] = SuccFa(f, 1, x0, d1, 1, epsilon, 1e4); 
  63. x1 = x0 + alpha * d1; 
  64. d2 = My_Gradient(f, x1); 
  65. xlog(iter+1) = norm(d2); 
  66. x = x1; 
  67. x0 = x1; 
  68. if norm(d2) < epsilon 
  69. break 
  70. end 
  71. iter = iter + 1; 
  72. end 
  73. end 
  74.  
  75. function [alphak xlog] = SuccFa(fun, alpha, x0, diff_x, h, epsil, maxIt) 
  76. k = 0; 
  77. xlog = alpha; 
  78. while k <= maxIt 
  79. alphak = alpha + h; 
  80. if fun(x0 + alphak * diff_x) < fun(x0 + alpha * diff_x) 
  81. h = 2 * h; 
  82. alpha = alphak; 
  83. else 
  84. h = - h / 4; 
  85. end 
  86. k = k + 1; 
  87. xlog(k) = alphak; 
  88. if abs(h) < epsil 
  89. break 
  90. end 
  91. end 
  92. end 

测试用例

  1. f = @(x) sin(x(1)) + cos(x(2)); 
  2. x0 = [1, 1]'; 
  3. epsilon = 1e-6; 
  4. [x xlog] = Newton(f, x0, epsilon) 
  5.  
  6. %% 书上的例子 
  7. f = @(x) 2*x(1)^4 - 4 * x(1)^2*x(2) + x(1)^2 + 2* x(2)^2-2*x(1)+5; 
  8. x0 = [0 ,0]'; 
  9. epsilon = 1e-6; 
  10. [x xlog] = Newton(f, x0, epsilon) 

牛顿法的收敛速度比前面的共轭梯度法更快，因为使用了二阶梯度信息以及线搜索方法，但是数值黑塞矩阵的计算比较复杂，如果提前知道黑塞矩阵的具体表达式则可接受.