最优化方法——拟牛顿法

对称秩1算法

目标函数

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

梯度

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

对称秩1算法程序

function[k,x,val]=sr1(fun,gfun,x0,epsilon,N)
if nargin<5
    N=1000;
end
if nargin<4
    epsilon=1.e-5;
end
beta=0.55;
sigma=0.4;
n=length(x0);
Hk=eye(n);
k=0;
while(k<N)
    gk=feval(gfun,x0);
    dk=-Hk*gk;
    if(norm(gk)<epsilon)
        break;
    end
    m=0;
    mk=0;
    while(m<20)
        if(feval(fun,x0+beta^m*dk)<=feval(fun,x0)+sigma*beta^m*gk'*dk)
            mk=m;
            break;
        end
        m=m+1;
    end
    x=x0+beta^mk*dk;
    sk=x-x0;
    yk=feval(gfun,x)-gk;
    Hk=Hk+(sk-Hk*yk)*(sk-Hk*yk)'/((sk-Hk*yk)'*yk);
    k=k+1;
    x0=x;
end
val=feval(fun,x0);

测试文件

clear;
clc;
x0=[1.0;1.0];
[k,x,val]=sr1('fun','gfun',x0)

输出结果

BFGS算法

BFGS算法程序

function[k,x,val]=bfgs(fun,gfun,x0,varargin)
N=1000;
epsilon=1.e-5;
beta=0.55;
sigma=0.4;
n=length(x0);
Bk=eye(n);
k=0;
while(k<N)
    gk=feval(gfun,x0,varargin{:});
    if(norm(gk)<epsilon),
        break;
    end
    dk=-Bk\gk;
    m=0;
    mk=0;
    while(m<20)
        newf=feval(fun,x0+beta^m*dk,varargin{:});
        oldf=feval(fun,x0,varargin{:});
        if(newf<=oldf+sigma*beta^m*gk'*dk)
            mk=m;
            break;
        end
        m=m+1;
    end
    x=x0+beta^mk*dk;
    sk=x-x0;
    yk=feval(gfun,x,varargin{:})-gk;
    if(yk'*sk>0)
        Bk=Bk-(Bk*(sk*sk')*Bk)/(sk'*Bk*sk)+(yk*yk')/(yk'*sk);
    end
    k=k+1;
    x0=x;
end
val=feval(fun,x0,varargin{:});

测试文件

clear;
clc;
x0=[1.0;1.0];
[k,x,val]=bfgs('fun','gfun',x0)

fun、gfun已经在最开始部分给出！

输出结果

DPF算法

DPF算法程序

function[k,x,val]=dfp(fun,gfun,x0,epsilon,N)
if nargin<5
    N=1000;
end
if nargin<4
    epsilon=1.e-5;
end
beta=0.55;
sigma=0.4;
n=length(x0);
Hk=eye(n);
k=0;
while(k<N)
    gk=feval(gfun,x0);
    if(norm(gk)<epsilon),
        break;
    end
    dk=-Hk*gk;
    m=0;
    mk=0;
    while(m<20)
        if(feval(fun,x0+beta^m*dk)<=feval(fun,x0)+sigma*beta^m*gk'*dk)
            mk=m;
            break;
        end
        m=m+1;
    end
    x=x0+beta^mk*dk;
    sk=x-x0;
    yk=feval(gfun,x)-gk;
    if(sk'*yk>0)
        Hk=Hk-(Hk*(yk*yk')*Hk)/(yk'*Hk*yk)+(sk*sk')/(sk'*yk);
    end
    k=k+1;
    x0=x;
end
val=feval(fun,x0);

测试文件

clear;
clc;
x0=[1.0;1.0];
[k,x,val]=dfp('fun','gfun',x0)

fun、gfun已经在最开始部分给出！

输出结果

Broyden族算法

Broyden族算法程序

function[k,x,val]=broyden(fun,gfun,x0,epsilon,N)
phi=0.5;
if nargin<5
    N=1000;
end
if nargin<4
    epsilon=1.e-5;
end
beta=0.55;
sigma=0.4;
n=length(x0);
Hk=eye(n);
k=0;
while(k<N)
    gk=feval(gfun,x0);
    if(norm(gk)<epsilon),
        break;
    end
    dk=-Hk*gk;
    m=0;
    mk=0;
    while(m<20)
        if(feval(fun,x0+beta^m*dk)<=feval(fun,x0)+sigma*beta^m*gk'*dk)
            mk=m;
            break;
        end
        m=m+1;
    end
    x=x0+beta^mk*dk;
    sk=x-x0;
    yk=feval(gfun,x)-gk;
    Hy=Hk*yk;
    sy=sk'*yk;
    yHy=yk'*Hk*yk;
    if(sy<0.2*yHy)
        theta=0.8*yHy/(yHy-sy);
        sk=theta*sk+(1-theta)*Hy;
        sy=0.2*yHy;
    end
    vk=sqrt(yHy)*(sk/sy-Hy/yHy);
    Hk=Hk-(Hy*Hy')/yHy+(sk*sk')/sy+phi*vk*vk';
    x0=x;
    k=k+1;
end
val=feval(fun,x0);

文件测试

clear;
clc;
x0=[1.0;1.0];
[k,x,val]=broyden('fun','gfun',x0)

输出结果