第五章拟牛顿法

╰半夏微凉° 2022-10-08 14:22 338阅读 0赞

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

# 拟牛顿法及其性质 #

**基本思想**：在基本牛顿法的步2中用Hess矩阵 G k G\_k Gk的某个**近似矩阵 B k B\_k Bk(对称正定)取代 G k G\_k Gk**  
j将 d k = − G k − 1 g k d\_k=-G\_k^\{-1\}g\_k dk=−Gk−1gk转变成

d k = − B k − 1 g k d\_k=-B\_k^\{-1\}g\_k dk=−Bk−1gk（**BFGS**）或 d k = − H k g k d\_k=-H\_kg\_k dk=−Hkgk（**DFP**）

## 对称秩1算法 ##

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1MDYwNjc0_size_16_color_FFFFFF_t_70]

## 对称秩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实现 #

## 算法 ##

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1MDYwNjc0_size_16_color_FFFFFF_t_70 1]

## 程序 ##

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 G\_k Gk的近似替代及其对应的**校正公式上（ s k s\_k sk和 y k y\_k yk的位置互换了**）略有不同。

## 算法 ##

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1MDYwNjc0_size_16_color_FFFFFF_t_70_pic_center]

## 程序 ##

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

[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1MDYwNjc0_size_16_color_FFFFFF_t_70]: /images/20221005/9940477173274e8fb02c3a7ac74ec6d3.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1MDYwNjc0_size_16_color_FFFFFF_t_70 1]: /images/20221005/c9c55d9b253b4faab67764231345ba4d.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1MDYwNjc0_size_16_color_FFFFFF_t_70_pic_center]: /images/20221005/548338f54365417f88c9908ed3882aaf.png