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

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

数值优化方法Matlab一维线搜索

在（一）中我们提到过下降算法即是按照迭代, 其中为步长，为下降方向。在射线上寻求合适的步长，即所谓的一维线搜索，问题可表述为：

称在此最优意义下的步长为最优步长或精确线搜索步长，称这种确定步长的方法为精确线搜索。上述问题的最优性条件为, 即是(正交性条件)

由于最优步长的的精确解不一定能得出或难以计算，因此考虑求解最优步长的数值解法，主要有
试探法
  二分法
  成功-失败法
函数逼近法（用简单函数的极小点来近似原函数的极小点）
 牛顿法
 抛物线法

1. 二分法

二分法是很经典的逼近方法，其要求目标函数是一维且可导的，且在试探区间上. 算法的思路简单，如下所示：

算法 1 二分法
Step 1. 设, , , 给定终止条件, .

Step 2. 令, 若, 令并停止, 否则转下步.

Step 3. 若, 置, 否则. 转下步.

Step 4. 若, 则记并停止，否则转下步.

Step 5. 若, 令并停止，否则并转Step 2.

上述二分法是书中给出的，但是有的地方没有给出必要的说明。首先算法中对没有凸性要求，那么, 情况下二分法可能收敛到局部最优解上. 此外，如果, ，若同时是凸函数，则最小值肯定在两个端点处，而若不是凸函数，那么上述二分法找到的是极大点.

附上对应的Matlab代码，经过改进（自动找到满足条件的ab）

1. % 一维最优化问题的二分搜索 
   

   ---

2. % fun 函数名称, 或匿名函数 
   

   ---

3. % a，b区间，其中a<b 
   

   ---

4. % epsil 指定终止条件（精度） 
   

   ---

5. % maxIt 最大迭代次数 
   

   ---

6. % bifact 表示加权系数，默认0.5 
   

   ---

7. function [xk xlog] = BisectionInt(fun, a, b, epsil, maxIt, bifact) 
   

   ---

8. if nargin < 6 
   

   ---

9. bifact = 0.5 
   

   ---

10. end 
    

    ---

11. alpha = bifact; 
    

    ---

12. beta = 1 - alpha; 
    

    ---

13. diff_fun = diff(fun); 
    

    ---

14. diff_fa = double(diff_fun(a)); 
    

    ---

15. diff_fb = double(diff_fun(b)); 
    

    ---

16. % 如果端点满足导数为0的条件 
    

    ---

17. if abs(diff_fa) < epsil  
    

    ---

18. xk = a; 
    

    ---

19. xlog = a; 
    

    ---

20. elseif abs(diff_fb) < epsil  
    

    ---

21. xk = b; 
    

    ---

22. xlog = b; 
    

    ---

23. elseif diff_fa < 0 && diff_fb > 0 
    

    ---

24. % 找到一个导数为0的点，此点是局部最优解 
    

    ---

25. [xk xlog] = Bisection(a, b, epsil, maxIt, diff_fun, alpha, beta); 
    

    ---

26. else 
    

    ---

27. % 非标准二分法 
    

    ---

28. % 先尝试能否找到满足标准二分法的端点 
    

    ---

29. h = (b-1)/1000; 
    

    ---

30. ha = a; 
    

    ---

31. hb = b; 
    

    ---

32. while double(diff_fun(ha)) > 0 && ha < b 
    

    ---

33. ha = ha + h; 
    

    ---

34. end 
    

    ---

35. while double(diff_fun(hb)) < 0 && hb > a 
    

    ---

36. hb = hb - h; 
    

    ---

37. end 
    

    ---

38. if hb > ha 
    

    ---

39. %找到了这样的端点 
    

    ---

40. [xk xlog] = Bisection(ha, hb, epsil, maxIt, diff_fun, alpha, beta); 
    

    ---

41. else 
    

    ---

42. %没找到这样的端点 
    

    ---

43. fa = fun(a); 
    

    ---

44. if fun(a) < fun(b) 
    

    ---

45. xk = a; 
    

    ---

46. else 
    

    ---

47. xk = b; 
    

    ---

48. end 
    

    ---

49. xlog = xk; 
    

    ---

50. end 
    

    ---

51. end 
    

    ---

52. end 
    

    ---

53.  
    

    ---

54. function [xk xlog] = Bisection(a, b, epsil, maxIt, diff_fun, alpha, beta) 
    

    ---

55. % 标准二分法 
    

    ---

56. k = 0; 
    

    ---

57. while k <= maxIt 
    

    ---

58. xk = alpha * a + beta * b; 
    

    ---

59. diff_fx = double(diff_fun(xk)); 
    

    ---

60. if diff_fx < 0 
    

    ---

61. a = xk; 
    

    ---

62. else 
    

    ---

63. b = xk; 
    

    ---

64. end 
    

    ---

65. k = k + 1; 
    

    ---

66. xlog(k) = xk;  
    

    ---

67. if abs(diff_fx) < epsil 
    

    ---

68. break; 
    

    ---

69. end 
    

    ---

70. if abs(a-b) < epsil 
    

    ---

71. break; 
    

    ---

72. end 
    

    ---

73. end 
    

    ---

74. end 
    

    ---

测试结果

1. syms f(x) 
   

   ---

2. f(x) = sin(x) + cos(x)-0.1*x; 
   

   ---

3. bifact = 0.5 
   

   ---

4. maxIt = 1000; 
   

   ---

5. epsil = 0.001; 
   

   ---

6.  
   

   ---

7. x = 0:0.01:10; 
   

   ---

8. plot(x,double(f(x))); 
   

   ---

9. hold on 
   

   ---

10. a = 2 
    

    ---

11. b = 6 
    

    ---

12. [xk xlog] = BisectionInt(f, a, b, epsil, maxIt, bifact) 
    

    ---

13. scatter(xk, double(f(xk))) 
    

    ---

14. hold on 
    

    ---

15.  
    

    ---

16. a = 0 
    

    ---

17. b = 2 
    

    ---

18. [xk xlog] = BisectionInt(f, a, b, epsil, maxIt, bifact) 
    

    ---

19. scatter(xk, double(f(xk))) 
    

    ---

20. hold on 
    

    ---

21.  
    

    ---

22. a = 2 
    

    ---

23. b = 8 
    

    ---

24. [xk xlog] = BisectionInt(f, a, b, epsil, maxIt, bifact) 
    

    ---

25. scatter(xk, double(f(xk))) 
    

    ---

26. hold off 
    

    ---

enter description here

如果函数复杂，比如分段函数，那么直接使用diff符号求导可能不合适，则需要使用diff的数值求导.
*可以看到二分法并不简单啊~~

2. 成功-失败法

二分法要求函数可导，而成功失败发不要求初始区间和目标函数可导性，其算法格式如下
算法 2 成功-失败法
Step 1. 给定初始点, 搜索步长, 计算精度, 置.
Step 2. 令, 计算. 若, 则称探索成功, . 若, 则称探索失败，, 且, 转下步.
Step 3. 检验, 若成立，则, 否则转 Step 3.

算法中标红部分是书中没有写的，若没有这一步骤算法将会周期震荡。

代码如下

1. % 一维最优化问题的二分搜索 
   

   ---

2. % fun 目标函数, 可以是符号函数或函数句柄 
   

   ---

3. % epsil 指定终止条件（精度） 
   

   ---

4. % maxIt 最大迭代次数 
   

   ---

5. function [xk xlog] = SuccFa(fun, x0, h, epsil, maxIt) 
   

   ---

6. if ~isa(fun,'function_handle') 
   

   ---

7. fun = matlabFunction(fun); 
   

   ---

8. end 
   

   ---

9. k = 0; 
   

   ---

10. xlog = x0; 
    

    ---

11. while k <= maxIt 
    

    ---

12. xk = x0 + h; 
    

    ---

13. if fun(xk) < fun(x0) 
    

    ---

14. h = 2 * h; 
    

    ---

15. x0 = xk; %如果这一行放在判断之后，则不收敛 
    

    ---

16. else 
    

    ---

17. h = - h / 4; 
    

    ---

18. end 
    

    ---

19. k = k + 1; 
    

    ---

20. xlog(k) = xk; 
    

    ---

21. if abs(h) < epsil 
    

    ---

22. break 
    

    ---

23. end 
    

    ---

24. end 
    

    ---

25. end 
    

    ---

测试结果

1. syms f(x) 
   

   ---

2. f(x) = sin(x) + cos(x); 
   

   ---

3. tf = matlabFunction(f); 
   

   ---

4. [xk xlog] = SuccFa(f, 1, 0.001, 1e-7, 1e4); 
   

   ---

5. plot(xlog) 
   

   ---

6. plot(0:0.01:10,tf(0:0.01:10)) 
   

   ---

如图所示