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

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) 
    

    ---

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