[原创]：坐标轮换法程序中匿名函数的使用心得
本帖最后由 bainhome 于 2013-3-23 22:52 编辑从前写过坐标轮换法多维无约束优化程序，当时用6，没匿名函数，加之水平眼界有限，符号函数实现，又笨又慢。最近因接触现代设计方法这本书，看到这个算法流程，来了兴致重写一遍。
言归正传，坐标轮换法是利用函数值信息直接寻优最基础的一个算法，迭代步数较多，效能取决于目标函数的形式，目标函数等值线椭圆形长短轴与坐标轴不平行时迭代计算曲折，收敛速度慢。对很多问题搜索不到合理解，算法本身没太多值得提倡的东西，不过匿名函数的用法倒还不错，另外算法实现方面，rocwoods这货今天给出重要提示，因此拿来和大家共享。
function [ outX,outY,count ] = CoodinateOpt( fx,x0,epsx,epsf )
% ================================
% 用于解决多维函数无约束优化问题
% 算法：坐标轮换法
% 用户测试代码：
% fx=@(x)100*(x(1).^2+x(2)-5).^2+(x(1)-x(2)).^2;
% [ outX,outY,count ] = CoodinateOpt( fx,,1e-2,1e-3 )
% ================================
clc
tic
lenX=length(x0);
x00=x0(:);
xStart=x0(:);
xn1=100*eye(lenX,1);
dir1=eye(lenX);
fx0=fx(x0);
count=0;
fn=1000;
while norm(xn1-x00)>epsx&&abs(fn-fx0)>epsf&&count<=5000
xIter=xStart;
x00=xStart;
fx0=fx(xStart);
xmin=0;
for i=1:lenX
f=@(a) fx(xIter+a*dir1(:,i));
xmin=Xquad(f,xmin,.01,.001);
xIter=xIter+xmin*dir1(:,i);
end
xStart=xIter;
xn1=xIter;
fn=fx(xn1);
count=count+1;
end
outX=xn1;
outY=fx(xn1);
toc
function xmin=Xquad(f,x1,h,t)
% ***********************************************
% 二次插值法计算最优步长子程序
% ***********************************************
clc
=jintui(f,x1,h);
count=0;
a1=(ax+bx)/2;
xx=;
for i=1:3
Y(i)=f(xx(i));
end
c1=(Y(3)-Y(1))/(bx-ax);
c2=((Y(2)-Y(1))/(a1-ax)-c1)/(a1-bx);
xp=.5*(ax+bx-c1/c2);
if abs(f(xp)-f(a1))<=t
xmin=min();
fmin=f(xmin);
step_opt=xmin-x1;
else
while abs(f(xp)-f(a1))>=t
=quadCompute(f,ax,a1,xp,bx);
end
xmin=min();
fmin=f(xmin);
step_opt=xmin-x1;
count=count+1;
end
count;
function =quadCompute(f,ax,a1,xp,bx)
if f(a1)<f(xp)&&xp>a1
bx=xp;
elseif f(a1)>f(xp)&&xp>a1
ax=a1;
a1=xp;
elseif f(a1)<f(xp)&&xp<a1
ax=xp;
elseif f(a1)>f(xp)&&xp<a1
bx=a1;
a1=xp;
end
xx=;
for i=1:3
Y(i)=f(xx(i));
end
c1=(Y(3)-Y(1))/(bx-ax);
c2=((Y(2)-Y(1))/(a1-ax)-c1)/(a1-bx);
xp=.5*(ax+bx-c1/c2);
function =jintui(f,x1,h)
% ================================
% 进退法确定搜索区间
% ================================
f1=f(x1);
x2=x1+h;
f2=f(x2);
countjintui=0;
if f1<f2
h=-h;
x2=x1;
f2=f1;
x3=x2+h;
f3=f(x3);
else
h=2*h;
x3=x2+h;
f3=f(x3);
end
while f2>f3
if countjintui<5000
x1=x2;
x2=x3;
x3=x2+h;
f1=f2;
f2=f3;
f3=f(x3);
countjintui=countjintui+1;
else
disp('该搜索区间不存在“高－低－高”的正确搜索形态，请修改后再试！')
return
end
end
a=min(x1,x3);
b=max(x1,x3);基本思想是将一个无约束n维优化问题转化为依次沿着相应的n个坐标轴方向的一维优化问题求解。
以二维目标函数f(X)=f(x1,x2)为例，从某个初始点X^(0)作为第一轮迭代起始点X0^(1)。
1.以x1坐标轴的单位向量$e_1=^T作为搜索方向，用一维优化方法确定其最优步长k1^(1)，获得第一轮第一个迭代点X1^(1)=X0^(1)+k1^(1)e1；
2.以X1^(1)为新起点，改以x2坐标轴的单位向量e2=^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*。
难点：例如如下二维无约束问题目标函数应如何使其被降为一维搜索问题，以确定每一次的最优步长因子。
fx=@(x) 100*(x(1).^2+x(2)-5).^2+(x(1)-x(2)).^2;解决办法是构造如下匿名函数：
fx=@(x)100*(x(1).^2+x(2)-5).^2+(x(1)-x(2)).^2;
d1=';
d2=.5*d1;
x=arrayfun(@(i) fminsearch(@(a) fx(d1(i,:)+a*d2(i,:)),0),1:length(d1))计算结果：
x =
1.5826 -0.20869 -0.80581 -1.1044 -1.2835匿名函数的确是MATLAB中比较有特色的编程方法。
用该算法计算上述二维无约束极值问题，结果如下：
Elapsed time is 0.429963 seconds.
outX =
1.7913
1.7914
outY =
9.6461e-08
count =
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:
Elapsed time is 0.251491 seconds.
fminunc采用拟牛顿法17次迭代到1e-8容差，自编程序迭代274次，再次证明梯度法收敛速度的确不是盖的。
ps1：从前倒看过很多优化问题书籍，很多优化算法仍在使用inline函数定义目标函数，这个命令的效率比较低。
ps2：很久没来，冒个泡再急速下潜...
页:
[1]