最优化方法——信頼域算法

信頼域方法

信頼域方法程序

function[d,val,lam,i]=trustq(fk,gk,Bk,deltak)
n=length(gk);
beta=0.6;
sigma=0.2;
mu0=0.05;
lam0=0.05;
gamma=0.05;
d0=ones(n,1);
z0=[mu0,lam0,d0']';
zbar=[mu0,zeros(1,n+1)]';
i=0;
z=z0;
mu=mu0;
lam=lam0;
d=d0;
while(i<=150)
    H=dah(mu,lam,d,gk,Bk,deltak);
    if(norm(H)<=1.e-8)
        break;
    end
    J=JacobiH(mu,lam,d,Bk,deltak);
    b=psis(mu,lam,d,gk,Bk,deltak,gamma)*zbar-H;
    dz=J\b;
    dmu=dz(1);
    dlam=dz(2);
    dd=dz(3:n+2);
    m=0;
    mi=0;
    while(m<20)
        t1=beta^m;
        Hnew=dah(mu+t1*dmu,lam+t1*dlam,d+t1*dd,gk,Bk,deltak);
        if(norm(Hnew)<=(1-sigma*(1-gamma*mu0)*beta^m)*norm(H))
            mi=m;
            break;
        end
        m=m+1;
    end
    alpha=beta^mi;
    mu=mu+alpha*dmu;
    lam=lam+alpha*dlam;
    d=d+alpha*dd;
    i=i+1;
end
val=fk+gk'*d+0.5*d'*Bk*d;

function p=phi(mu,a,b)
p=a+b-sqrt((a-b)^2+4*mu^2);

function H=dah(mu,lam,d,gk,Bk,deltak)
n=length(d);
H=zeros(n+2,1);
H(1)=mu;
H(2)=phi(mu,lam,deltak^2-norm(d)^2);
H(3:n+2)=(Bk+lam*eye(n))*d+gk;

function J=JacobiH(mu,lam,d,Bk,deltak)
n=length(d);
J=zeros(n+2,n+2);
t2=sqrt((lam+norm(d)^2-deltak^2)^2+4*mu^2);
pmu=-4*mu/t2;
thetak=(lam+norm(d)^2-deltak^2)/t2;
J=[1,0,zeros(1,n);
pmu,1-thetak,-2*(1+thetak)*d';
zeros(n,1),d,Bk+lam*eye(n)];

function si=psis(mu,lam,d,gk,Bk,deltak,gamma)
H=dah(mu,lam,d,gk,Bk,deltak);
si=gamma*norm(H)*min(1,norm(H));

测试文件

clear;
clc;
fk=-5;
gk=[400 -200]';
Bk=[1202 -400;-400 200];
deltak=5;
[d,val,lam,i]=trustq(fk,gk,Bk,deltak)

测试文件

clear;
clc;
fk=-5;
gk=[400 -200]';
Bk=[1202 -400;-400 200];
deltak=5;
[d,val,lam,i]=trustq(fk,gk,Bk,deltak)

运行结果

作业

测试文件
clear;
clc;
fk=52;
gk=[52 24]’;
Bk=[194 -4;-4 8];
deltak=5;
[d,val,lam,i]=trustq(fk,gk,Bk,deltak)

牛顿型信頼域方法

牛顿型信頼域方法程序

function[k,x,val]=trustm(x0,epsilon)
n=length(x0);
eta1=0.1;
eta2=0.75;
tau1=0.5;
tau2=2.0;
delta=1;
dtabar=2.0;
x=x0;
Bk=Hess(x);
k=0;
while(k<50)
    fk=fun(x);
    gk=gfun(x);
    if(norm(gk)<epsilon)
        break;
    end
    [d,val,lam,i]=trustq(fk,gk,Bk,delta);
    deltaq=fk-val;
    deltaf=fun(x)-fun(x+d);
    rk=deltaf/deltaq;
    if(rk<=eta1)
        delta=tau1*delta;
    else if(rk>=eta2&&norm(d)==delta)
            delta=min(tau2*delta,dtabar);
        else
            delta=delta;
        end
    end
    if(rk>eta1)
        x=x+d;
        Bk=Hess(x);
    end
    k=k+1;
end
val=fun(x);

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))];

function He=Hess(x)
He=[1200*x(1)^2-400*x(2)+2,-400*x(1);-400*x(1),200];

测试文件

clear;
clc;
x0=[0;0];
epsilon=1e-6;
[k,x,val]=trustm(x0,epsilon)

牛顿型方法有调用前面的m文件。

输出结果