[原创]：坐标轮换法程序中匿名函数的使用心得

bainhome

从前写过坐标轮换法多维无约束优化程序，当时用6，没匿名函数，加之水平眼界有限，符号函数实现，又笨又慢。最近因接触现代设计方法这本书，看到这个算法流程，来了兴致重写一遍。

言归正传，坐标轮换法是利用函数值信息直接寻优最基础的一个算法，迭代步数较多，效能取决于目标函数的形式，目标函数等值线椭圆形长短轴与坐标轴不平行时迭代计算曲折，收敛速度慢。对很多问题搜索不到合理解，算法本身没太多值得提倡的东西，不过匿名函数的用法倒还不错，另外算法实现方面，rocwoods这货今天给出重要提示，因此拿来和大家共享。

1. 
   
2. 
   
3. 
   
4. function [ outX,outY,count ] = CoodinateOpt( fx,x0,epsx,epsf )
   
5. % ================================
   
6. % 用于解决多维函数无约束优化问题
   
7. % 算法：坐标轮换法
   
8. % 用户测试代码：
   
9. % fx=@(x)100*(x(1).^2+x(2)-5).^2+(x(1)-x(2)).^2;
   
10. % [ outX,outY,count ] = CoodinateOpt( fx,[1,1],1e-2,1e-3 )
    
11. % ================================
    
12. clc
    
13. tic
    
14. lenX=length(x0);
    
15. x00=x0(:);
    
16. xStart=x0(:);
    
17. xn1=100*eye(lenX,1);
    
18. dir1=eye(lenX);
    
19. fx0=fx(x0);
    
20. count=0;
    
21. fn=1000;
    
22. while norm(xn1-x00)>epsx&&abs(fn-fx0)>epsf&&count<=5000
    
23.     xIter=xStart;
    
24.     x00=xStart;
    
25.     fx0=fx(xStart);
    
26.     xmin=0;
    
27.     for i=1:lenX
    
28.         f=@(a) fx(xIter+a*dir1(:,i));
    
29.         xmin=Xquad(f,xmin,.01,.001);
    
30.         xIter=xIter+xmin*dir1(:,i);
    
31.     end
    
32.     xStart=xIter;
    
33.     xn1=xIter;
    
34.     fn=fx(xn1);
    
35.     count=count+1;
    
36. end
    
37. outX=xn1;
    
38. outY=fx(xn1);
    
39. toc
    
40. function xmin=Xquad(f,x1,h,t)
    
41. % ***********************************************
    
42. % 二次插值法计算最优步长子程序
    
43. % ***********************************************
    
44. clc
    
45. [ax,bx,countjintui]=jintui(f,x1,h);
    
46. count=0;
    
47. a1=(ax+bx)/2;
    
48. xx=[ax,a1,bx];
    
49. for i=1:3
    
50.     Y(i)=f(xx(i));
    
51. end
    
52. c1=(Y(3)-Y(1))/(bx-ax);
    
53. c2=((Y(2)-Y(1))/(a1-ax)-c1)/(a1-bx);
    
54. xp=.5*(ax+bx-c1/c2);
    
55. if abs(f(xp)-f(a1))<=t
    
56.     xmin=min([xp,a1]);
    
57.     fmin=f(xmin);
    
58.     step_opt=xmin-x1;
    
59. else
    
60.     while abs(f(xp)-f(a1))>=t
    
61.        [ax,a1,xp,bx]=quadCompute(f,ax,a1,xp,bx);
    
62.     end
    
63.      xmin=min([xp,a1]);
    
64.      fmin=f(xmin);
    
65.      step_opt=xmin-x1;
    
66.     count=count+1;
    
67. end
    
68. count;
    
69. function [ax,a1,xp,bx]=quadCompute(f,ax,a1,xp,bx)
    
70. if f(a1)<f(xp)&&xp>a1
    
71.     bx=xp;
    
72. elseif f(a1)>f(xp)&&xp>a1
    
73.     ax=a1;
    
74.     a1=xp;
    
75. elseif f(a1)<f(xp)&&xp<a1
    
76.     ax=xp;
    
77. elseif f(a1)>f(xp)&&xp<a1
    
78.     bx=a1;
    
79.     a1=xp;   
    
80. end
    
81. xx=[ax,a1,bx];
    
82. for i=1:3
    
83.     Y(i)=f(xx(i));
    
84. end
    
85. c1=(Y(3)-Y(1))/(bx-ax);
    
86. c2=((Y(2)-Y(1))/(a1-ax)-c1)/(a1-bx);
    
87. xp=.5*(ax+bx-c1/c2);
    
88. function [a,b,countjintui]=jintui(f,x1,h)
    
89. % ================================
    
90. % 进退法确定搜索区间
    
91. % ================================
    
92. f1=f(x1);
    
93. x2=x1+h;
    
94. f2=f(x2);
    
95. countjintui=0;
    
96. if f1<f2
    
97.     h=-h;
    
98.     x2=x1;
    
99.     f2=f1;
    
100.     x3=x2+h;
     
101.     f3=f(x3);
     
102. else
     
103.     h=2*h;
     
104.     x3=x2+h;
     
105.     f3=f(x3);
     
106. end
     
107. while f2>f3
     
108.     if countjintui<5000
     
109.         x1=x2;
     
110.         x2=x3;
     
111.         x3=x2+h;
     
112.         f1=f2;
     
113.         f2=f3;
     
114.         f3=f(x3);
     
115.         countjintui=countjintui+1;
     
116.     else
     
117.         disp('该搜索区间不存在“高－低－高”的正确搜索形态，请修改后再试！')
     
118.         return
     
119.     end
     
120. end
     
121. a=min(x1,x3);
     
122. b=max(x1,x3);

复制代码

基本思想是将一个无约束n维优化问题转化为依次沿着相应的n个坐标轴方向的一维优化问题求解。
以二维目标函数f(X)=f(x1,x2)为例，从某个初始点X^(0)作为第一轮迭代起始点X0^(1)。

1.以x1坐标轴的单位向量$e_1=[1,0]^T作为搜索方向，用一维优化方法确定其最优步长k1^(1)，获得第一轮第一个迭代点X1^(1)=X0^(1)+k1^(1)e1；
2.以X1^(1)为新起点，改以x2坐标轴的单位向量e2=[0,1]^T为搜索方向，用一维优化方法确定其最优步长k2^(1)，获得第一轮第二个迭代点X2^(1)=X1^(1)+k2^(1)e2
3.对二维问题，经过依次沿着e1和e2的两次搜索，完成第一轮迭代，得到第一轮迭代极小点X^(1)=X2^(1)；
4.如经过第一轮迭代计算后仍然没有得到收敛到目标函数f(X)的极小点X*，则以第一轮迭代的极小点X^{1}作为第二轮迭代的起始点X2^(0)=X^(1)，依次再沿着e1和e2方向搜索，按上述方法做第3,4,...搜索，直至迭代逼近目标函数f(X)的极小点X*。

难点：例如如下二维无约束问题目标函数应如何使其被降为一维搜索问题，以确定每一次的最优步长因子。

1. 
   
2. fx=@(x) 100*(x(1).^2+x(2)-5).^2+(x(1)-x(2)).^2;

复制代码

解决办法是构造如下匿名函数：

1. fx=@(x)100*(x(1).^2+x(2)-5).^2+(x(1)-x(2)).^2;
   
2. d1=[1:5;1:5]';
   
3. d2=.5*d1;
   
4. x=arrayfun(@(i) fminsearch(@(a) fx(d1(i,:)+a*d2(i,:)),0),1:length(d1))

复制代码

计算结果：

1. x =
   
2. 1.5826 -0.20869 -0.80581 -1.1044 -1.2835

复制代码

匿名函数的确是MATLAB中比较有特色的编程方法。

用该算法计算上述二维无约束极值问题，结果如下：

1. Elapsed time is 0.429963 seconds.
   
2. outX =
   
3.        1.7913
   
4.        1.7914
   
5. outY =
   
6.    9.6461e-08
   
7. count =
   
8.    274

复制代码

fminunc命令的计算结果：
x =
1.7913 1.7913
fval =
7.3758e-12
exitflag =
1
output =
iterations: 17
funcCount: 69
stepsize: 1
firstorderopt: 5.3476e-08
algorithm: 'medium-scale: Quasi-Newton line search'
message: [1x438 char]
Elapsed time is 0.251491 seconds.

fminunc采用拟牛顿法17次迭代到1e-8容差，自编程序迭代274次，再次证明梯度法收敛速度的确不是盖的。

ps1：从前倒看过很多优化问题书籍，很多优化算法仍在使用inline函数定义目标函数，这个命令的效率比较低。
ps2：很久没来，冒个泡再急速下潜...