最优化方法——线搜索

精确线搜索

黄金分割法

待求解函数

function [y]=phi(x)
y=3*x^2-2*tan(x);

黄金分割法MATLAB实现

function [i,s,phis,ds,dphi,G]=golds(phi,a,b,epsilon,delta)  %  Golden Section
t=(sqrt(5)-1)/2;   
h=b-a;  %Define interval length
phia=feval(phi,a);  %  Write function value with m file
phib=feval(phi,b);  %  Using feval to calculate function value
p=a+(1-t)*h;  %  Calculate the initial trial point p
q=a+t*h;  %  Calculate the initial trial point q
phip=feval(phi,p);  %  Function value of trial point
phiq=feval(phi,q);
i=1;  % Number of operations initialization
G(i,:)=[a,p,q,b];  
while(abs(phib-phia)>delta)||(h>epsilon)  %Shutdown criteria
    if(phip<=phiq)
        b=q;
        phib=phiq;
        q=p;
        phiq=phip;
        h=b-a;
        p=a+(1-t)*h;
        phip=feval(phi,p);
    else
        a=p;
        phia=phip;
        p=q;
        phip=phiq;
        h=b-a;
        q=a+t*h;
        phiq=feval(phi,q);
    end
    i=i+1;G(i,:)=[a,p,q,b];
end
    if(phip<=phiq)
        s=p;
        phis=phip;
    else
        s=q;
        phis=phiq;
    end
    ds=abs(b-a);
    dphi=(phib-phia);

测试函数

clear;
clc;
phi=@(x)3*x^2-2*tan(x);
[i,s,phis,ds,dphi,G]=golds(phi,0,1,1e-4,1e-5)

while(j<maxj & abs(h)>epsilon & cond==0)
if (phi0<=phi1),
s2=s1; phi2=phi1; h=0.5*h;
s1=s0+h; phi1=feval(phi,s1);
else if (phi2<phi1),
s1=s2; phi1=phi2; h=2*h;
s2=s0+2*h; phi2=feval(phi,s2);
else

输出结果：

抛物线法

待求解函数

function [y]=phi(x)
y=3*x^2-2*tan(x);

抛物线法MATLAB实现【一段有错误的代码，需要debug】

function [i,s,phis,ds,dphi,S]=qmin(phi,a,b,epsilon,delta)
%Exact line search of parabola method
s0=a; 
i=1; 
S(i)=s0; 
maxi=30; 
maxj=20;
big=1e6; 
err=1; 
cond=0; 
h=1; 
ds=0.00001;
if (abs(s0)>1e4), 
    h=abs(s0)*(1e-4); 
end
while (i<maxi && err>delta && cond~=5)
    f1=(feval(phi,s0+ds)-feval(phi,s0-ds))/(2*ds);
    if(f1>0) 
        h=-abs(h); 
    end
    s1=s0+h; 
    s2=s0+2*h; 
    bars=s0;
    phi0=feval(phi,s0); 
    phi1=feval(phi,s1);
    phi2=feval(phi,s2); 
    barphi=phi0; 
    cond=0;
    while(j<maxj && abs(h)>epsilon && cond==0)
        if (phi0<=phi1)
            s2=s1; 
            phi2=phi1; 
            h=0.5*h;
            s1=s0+h; 
            phi1=feval(phi,s1);
        else if (phi2<phi1)
                s1=s2; 
                phi1=phi2; 
                h=2*h;
                s2=s0+2*h; 
                phi2=feval(phi,s2);
            else
                j=0;
                cond=-1;
            end
        end
        j=j+1;
        if(abs(h)>big || abs(s0)>big)
            cond=5; 
        end
    end
    if(cond==5)
        bars=s1; 
        barphi=feval(phi,s1);
    %Calculating phi value by quadratic interpolation
    else
        d=2*(2*phi1-phi0-phi2);
        if(d<0)
            barh=h*(4*phi1-3*phi0-phi2)/d;
        else
            barh=h/3; 
            cond=4;
        end
        bars=s0+barh;
        barphi=feval(phi,bars);
        h=abs(h); 
        h0=abs(barh);
        h1=abs(barh-h); 
        h2=abs(barh-2*h);
            %Determine the H value of the next iteration
            if(h0<h)
                h=h0; 
            end
            if(h1<h)
                h=h1; 
            end
            if(h2<h)
                h=h2;
            end
            if(h==0)
                h=barh; 
            end
            if(h<epsilon)
                cond=1; 
            end
            if(abs(h)>big || abs(bars)>big)
                cond=5; 
            end
            err=abs(barphi-phi1);
            s0=bars; 
            i=i+1; 
            S(i)=s0;
    end
        if(cond==2 && h<epsilon)
            cond=3; 
        end
end
    s=s0; 
    phis=feval(phi,s);
    ds=abs(s-s1); 
    dphi=err;

function [i,s,phis,ds,dphi,S]=qmin(phi,a,b,epsilon,delta)
s0=a; 
i=1; 
S(i)=s0; 
maxi=30; 
maxj=20;
big=1e6; 
err=1; 
cond=0; 
h=1; 
ds=0.00001;
if (abs(s0)>1e4), h=abs(s0)*(1e-4); end
while (i<maxi && err>delta && cond~=5)
    f1=(feval(phi,s0+ds)-feval(phi,s0-ds))/(2*ds);
    if(f1>0), h=-abs(h); end
    s1=s0+h; 
    s2=s0+2*h; 
    bars=s0;
    phi0=feval(phi,s0); 
    phi1=feval(phi,s1);
    phi2=feval(phi,s2); 
    barphi=phi0; cond=0;
    j=0; 
    while(j<maxj && abs(h)>epsilon && cond==0)
        if (phi0<=phi1),s2=s1; phi2=phi1; h=0.5*h;
            s1=s0+h; phi1=feval(phi,s1);
        else if (phi2<phi1),s1=s2; phi1=phi2; 
                h=2*h;s2=s0+2*h; phi2=feval(phi,s2);
            else
                cond=-1;
            end
        end
        j=j+1;
        if(abs(h)>big || abs(s0)>big), cond=5; end
    end
    if(cond==5)
        bars=s1; 
        barphi=feval(phi,s1);
    else
        d=2*(2*phi1-phi0-phi2);
        if(d<0),barh=h*(4*phi1-3*phi0-phi2)/d;
        else
            barh=h/3; cond=4;
        end
        bars=s0+barh; barphi=feval(phi,bars);
        h=abs(h); h0=abs(barh);
        h1=abs(barh-h); h2=abs(barh-2*h);
        if(h0<h), h=h0; end
        if(h1<h), h=h1; end
        if(h2<h), h=h2; end
        if(h==0), h=barh; end
        if(h<epsilon), cond=1; end
        if(abs(h)>big || abs(bars)>big), cond=5; end
        err=abs(barphi-phi1);
        s0=bars; i=i+1; S(i)=s0;
    end
    if(cond==2 && h<epsilon), cond=3; end
end
s=s0; phis=feval(phi,s);
ds=abs(s-s1); dphi=err;

测试代码

clear;
clc;
phi=@(x)3*x^2-2*tan(x);
[i,s,phis,ds,dphi,S]=qmin(phi,0,1,1e-4,1e-5)

输出结果：

非精确线搜索

amijo算法

目标函数

function f=fun(x)
f=100*(x(1)^2-x(2))^2+(x(1)-1)^2;

梯度

function gf=gfun(x)
gf=[400*x(1)*(x(1)^2-x(2))+2*(x(1)-1);-200*(x(1)^2-x(2))];

amijo算法的MATLAB实现

function mk=armijo(xk,dk)
beta=0.5; 
sigma=0.2;
m=0; 
maxm=20;
while (m<=maxm)
    if(fun(xk+beta^m*dk)<=fun(xk)+sigma*beta^m*gfun(xk)'*dk)
        mk=m; 
        break;
    end
    m=m+1;
end
alpha=beta^mk
newxk=xk+alpha*dk
fk=fun(xk)
newfk=fun(newxk)

测试文件

clear;
clc;
xk=[-1,1]'; 
dk=[1,-2]'; 
mk=armijo(xk,dk)

输出结果：