第五章 拟牛顿法

内容来自马昌凤编著的《最优化方法及其Matlab程序设计》，文章仅为个人的学习笔记，感兴趣的朋友详见原书。

拟牛顿法及其性质

基本思想：在基本牛顿法的步2中用Hess矩阵$G_k$的某个近似矩阵$B_k$(对称正定)取代$G_k$
j将$d_k=-G_k^{-1}{g}_k$转变成

$d_k=-B_k^{-1}{g}_k$（BFGS）或$d_k=-H_k{g}_k$（DFP）

对称秩1算法

对称秩1算法程序

function [x,val,k]=sr1(fun,gfun, x0)
%功能: 用对称秩1算法求解无约束问题:  min f(x)
%输入: x0是初始点, fun, gfun分别是目标函数及其梯度
%输出:  x, val分别是近似最优点和最优值,  k是迭代次数.
maxk=500;   %给出最大迭代次数
rho=0.55;sigma=0.4; epsilon=1e-5; 
k=0;   n=length(x0); Hk=eye(n); 
while(k<maxk)
    gk=feval(gfun,x0); %计算梯度
    dk=-Hk*gk; %计算搜索方向
    if(norm(gk)<epsilon), break; end  %检验终止准则
    m=0; mk=0;
    while(m<20)   % 用Armijo搜索求步长 
        if(feval(fun,x0+rho^m*dk)<feval(fun,x0)+sigma*rho^m*gk'*dk)
            mk=m; break;
        end
        m=m+1;
    end
    x=x0+rho^mk*dk;  
    sk=x-x0;  yk=feval(gfun,x)-gk;
    Hk=Hk+(sk-Hk*yk)*(sk-Hk*yk)'/((sk-Hk*yk)'*yk); %秩1校正
    k=k+1;     x0=x;
end
val=feval(fun,x0); 

BFGS算法及其Matlab实现

算法

程序

function [x,val,k]=bfgs(fun,gfun,x0,varargin)
%功能: 用BFGS算法求解无约束问题:  min f(x)
%输入: x0是初始点, fun, gfun分别是目标函数及其梯度;
% varargin是输入的可变参数变量, 简单调用bfgs时可以忽略它,
% 但若其它程序循环调用该程序时将发挥重要的作用
%输出:  x, val分别是近似最优点和最优值,  k是迭代次数.
maxk=500;   %给出最大迭代次数
rho=0.55; sigma1=0.4; epsilon1=1e-5; 
k=0;   n=length(x0); 
Bk=eye(n);   %Bk=feval('Hess',x0);  两种取法，任选一个即可
while(k<maxk)
    gk=feval(gfun,x0,varargin{:}); %计算梯度
    if(norm(gk)<epsilon1), break; end  %检验终止准则
    dk=-Bk\gk;  %解方程组, 计算搜索方向 注意这里是左除
    m=0; mk=0;
    while(m<20)   % 用Armijo搜索求步长 
        newf=feval(fun,x0+rho^m*dk,varargin{:});
        oldf=feval(fun,x0,varargin{:});
        if(newf<oldf+sigma1*rho^m*gk'*dk)
            mk=m; break;
        end
        m=m+1;
    end
    %BFGS校正
    x=x0+rho^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{:}); 

DFP算法及其Matlab实现

DFP与BFGS算法类似，只是在$G_k$的近似替代及其对应的校正公式上（$s_k$和$y_k$的位置互换了）略有不同。

算法

程序

function [x,val,k]=dfp(fun,gfun,x0)
%功能: 用DFP算法求解无约束问题:  min f(x)
%输入: x0是初始点, fun, gfun分别是目标函数及其梯度
%输出:  x, val分别是近似最优点和最优值,  k是迭代次数.
maxk=1e5;   %给出最大迭代次数
rho=0.55;sigma=0.4; epsilon=1e-5; 
k=0;   n=length(x0); 
Hk=inv(feval('Hess',x0));   %Hk=eye(n);
while(k<maxk)
    gk=feval(gfun,x0); %计算梯度
    if(norm(gk)<epsilon), break; end  %检验终止准则
    dk=-Hk*gk;  %解方程组, 计算搜索方向
    m=0; mk=0;
    while(m<20)   % 用Armijo搜索求步长 
        if(feval(fun,x0+rho^m*dk)<feval(fun,x0)+sigma*rho^m*gk'*dk)
            mk=m; break;
        end
        m=m+1;
    end
    %DFP校正
    x=x0+rho^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);