最优化理论——信赖域方法

算法思想

算法步骤

子问题

算法步骤

代码

function  [d,val,lam,k]=trustq(gk,Bk,dta)
% 功能: 求解信赖域子问题:  min qk(d)=gk'*d+0.5*d'*Bk*d, s.t. ||d||<=delta
%输入:  gk是xk处的梯度, Bk是第k次近似Hesse阵, dta是当前信赖域半径
%输出:  d, val分别是子问题的最优点和最优值, lam是乘子值, k是迭代次数.
n=length(gk);   gamma=0.05; 
epsilon=1.0e-6;  rho=0.6;  sigma=0.2;   
mu0=0.05;  lam0=0.05;
d0=ones(n,1);   z0=[mu0,lam0,d0']'; 
u0=[mu0,zeros(1,n+1)]';
k=0;   %k为迭代次数
z=z0;
mu=mu0; lam=lam0; d=d0;
while ( k<=150)     %Step1 of the algorithm
    dh=dah(mu,lam,d,gk,Bk,dta);
     if(norm(dh)<epsilon)
         break;
     end
     A=JacobiH(mu,lam,d,Bk,dta); b=beta(mu,lam,d,gk,Bk,dta,gamma)*u0-dh;
     B=inv(A);      dz=B*b;
     dmu=dz(1); dlam=dz(2); dd=dz(3:n+2);
     m=0;  mk=0;
     while (m<20)
         dhnew=dah(mu+rho^m*dmu,lam+rho^m*dlam,d+rho^m*dd,gk,Bk,dta);
         if(norm(dhnew)<=(1-sigma*(1-gamma*mu0)*rho^m)*dh)
             mk=m;           
             break;
         end
         m=m+1;
     end
     alpha=rho^mk;
     mu=mu+alpha*dmu;
     lam=lam+alpha*dlam;
     d=d+alpha*dd;
     k=k+1;
end
 val=gk'*d+0.5*d'*Bk*d;
%%%%%%%%%%%%%%%%%%%%%%%%%%
function p=phi(mu,a,b)
p=a+b-sqrt((a-b)^2+4*mu);
%%%%%%%%%%%%%%%%%%%%%%%%%%
function dh=dah(mu,lam,d,gk,Bk,dta)
n=length(d);
dh(1)=mu;  dh(2)=phi(mu,lam, dta^2-norm(d)^2);
mh=(Bk+lam*eye(n))*d+gk;
for(i=1:n)
    dh(2+i)=mh(i);
end
dh=dh(:);
%%%%%%%%%%%%%%%%%%%%%%%%%%
function bet=beta(mu,lam,d,gk,Bk,dta,gamma)
dh=dah(mu,lam,d,gk,Bk,dta);
bet=gamma*norm(dh)*min(1,norm(dh));
%%%%%%%%%%%%%%%%%%%%%%%%%%
function A=JacobiH(mu,lam,d,Bk,dta)
n=length(d);
A=zeros(n+2,n+2);
pmu=-4*mu/sqrt((lam+norm(d)^2-dta^2)^2+4*mu^2);
thetak=(lam+norm(d)^2-dta^2)/sqrt((lam+norm(d)^2-dta^2)^2+4*mu^2);
A=[1,0,zeros(1,n);pmu,1-thetak,-2*(1+thetak)*d';zeros(n,1),  d,Bk+lam*eye(n)];

示例

在 Matlab 命令窗口依次输入下列命令

gk=[400 -200]’;
Bk=[1202 -400; -400 200];
dta=5;
[d,val,lam,k]=trustq(gk,Bk,dta)

原问题代码

Matlab代码如下：

function [xk,val,k]=trustm(x0)
%功能: 牛顿型信赖域方法求解无约束优化问题 min f(x)
%输入: x0是初始迭代点
%输出: xk是近似极小点, val是近似极小值, k是迭代次数
n=length(x0);  x=x0; dta=1;
eta1=0.15; eta2=0.75;  dtabar=2.0;
tau1=0.5; tau2=2.0; epsilon=1e-6;
k=0;  Bk=Hess(x);  %Bk=eye(n);  
while(k<150)
    gk=gfun(x);   
    if(norm(gk)<epsilon)
        break;
    end
    [d,val,lam,ik]=trustq(gk,Bk,dta);
    deltaq=-qk(x,d);
    deltaf=fun(x)-fun(x+d);
    rk=deltaf/deltaq;
    if(rk<=eta1)
        dta=tau1*dta;
    else if (rk>=eta2&norm(d)==dta)
            dta=min(tau2*dta,dtabar);
        else
            dta=dta;
        end
    end
    if(rk>eta1)
        x0=x;     x=x+d;    
       % sk=x-x0;  yk=gfun(x)-gfun(x0);  
        %vk=sqrt(yk'*Bk*yk)*(sk/(sk'*yk)-Bk*yk/(yk'*Bk*yk));
        %Bk=Bk-Bk*yk*yk'*Bk/(yk'*Bk*yk)+sk*sk'/(sk'*yk)+vk*vk'
        %pause
        Bk=Hess(x);
    end
    k=k+1;
end
xk=x;
val=fun(xk);

示例

考虑无约束优化问题
$$minf(x)=100(x_1^2-x_2{)}^2+(x_1-1{)}^2$$
该问题有精确解x=(1,1)T，$f(x)=0$.

fun函数文件：

%目标函数
function f=fun(x)
f=100*(x(1)^2-x(2))^2+(x(1)-1)^2;

gfun函数文件：

%梯度
function gf=gfun(x)
gf=[400*x(1)*(x(1)^2-x(2))+2*(x(1)-1), -200*(x(1)^2-x(2))]’;

目标函数的Hesse阵：

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

子问题目标函数：

function qd=qk(x,d)
gk=gfun(x);  Bk=Hess(x);
qd=gk'*d+0.5*d'*Bk*d;

交互界面输入：

x0=[-1,1]’;
[x,val,k]=[xk,val,k]=trustm(x0)

结果：