第三章最速下降法和牛顿法

刺骨的言语ヽ痛彻心扉 2022-10-08 12:26 265阅读 0赞

> 内容来自马昌凤编著的《最优化方法及其Matlab程序设计》，文章仅为个人的学习笔记，感兴趣的朋友详见原书。  
> 本章讨论无约束优化问题 m i n f ( x ) minf(x) minf(x)的最速下降法和牛顿法及其改进算法，前者简单而古老，虽不再具有实用性，却是研究其他无约束优化算法的基础；后者也是一种经典的无约束优化算法，具有收敛速度快和自适应性等优点。

# 最速下降方法及其Matlab实现 #

## 算法 ##

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

## 程序 ##

function [x,val,k]=grad(fun,gfun,x0)
    % 功能: 用最速下降法求解无约束问题:  min f(x)
    %输入:  x0是初始点, fun, gfun分别是目标函数和梯度
    %输出:  x, val分别是近似最优点和最优值,  k是迭代次数.
    maxk=5000;   %最大迭代次数
    rho=0.5;sigma=0.4;
    k=0;  epsilon=1e-5;
    while(k<maxk)
        g=feval(gfun,x0);  %计算梯度
        d=-g;    %计算搜索方向（精确线搜索）
        if(norm(d)<epsilon), break; end
        m=0; mk=0;
        while(m<20)   %Armijo搜索
            if(feval(fun,x0+rho^m*d)<feval(fun,x0)+sigma*rho^m*g'*d)
                mk=m; break;
            end
            m=m+1;
        end
        x0=x0+rho^mk*d;
        k=k+1;
    end
    x=x0;
    val=feval(fun,x0);

注：这里的 f u n fun fun与 g f u n gfun gfun分指目标函数与梯度，需在grad运行前准备好。（Armijo搜索部分可详见[第二章][Link 1]相关内容）

## 示例1 ##

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

## 示例2 ##

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

# 牛顿法及其Matlab实现 #

**基本思想**：**用迭代点 x k x\_k xk处的一阶导数和二阶导数对目标函数进行二次函数近似（泰勒多项式）**，然后把二次模型的极小点作为新的迭代点，并不断重复这一过程，直至求得满足精度的近似极小点。

## 算法 ##

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1MDYwNjc0_size_16_color_FFFFFF_t_70 3]  
牛顿法最突出的优点是收敛速度快，具有局部二阶收敛性。但是对初始点要求足够“靠近”极小点，否则，有可能导致算法不收敛。由于实际问题的精度极小点一般是未知的，为克服这一困难，故引入线搜索技术以得到大范围收敛的算法，即下面的阻尼牛顿法。（此处基于Armijo准则）

# 阻尼牛顿法 #

## 算法 ##

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

## 程序 ##

function [x,val,k]=dampnm(fun,gfun, Hess,x0)
    %功能: 用阻尼牛顿法求解无约束问题:  min f(x)
    %输入: x0是初始点, fun, gfun, Hess 分别是求
    %         目标函数,梯度,Hesse 阵的函数
    %输出:  x, val分别是近似最优点和最优值,  k是迭代次数.
    maxk=100;   %给出最大迭代次数
    rho=0.55;sigma=0.4;
    k=0;  epsilon=1e-5;
    while(k<maxk)
        gk=feval(gfun,x0); %计算梯度
        Gk=feval(Hess,x0);  %计算Hesse阵
        dk=-Gk\gk; %解方程组Gk*dk=-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
        x0=x0+rho^mk*dk;
        k=k+1;
    end
    x=x0;
    val=feval(fun,x); 
    %gval=norm(gfun(x));

此处除需提前建立目标函数和梯度的M文件外，还需建立Hess矩阵的求解文件

function He=Hess(x)
    n=length(x)
    He=zeros(n,n);
    He=[此处填函数的二阶求导结果]

# 修正牛顿法及其Matlab实现 #

牛顿法具有不低于二阶的收敛速度，但该算法要求目标函数的Hess矩阵在每个迭代点 x k x\_k xk处是正定的，否则，难以保证牛顿方向 d k = − G k − 1 g k d\_k=-G\_k^\{-1\}g\_k dk=−Gk−1gk是 f f f在 x k x\_k xk处的下降方向。为克服这一缺陷，故对其进行修正：将牛顿法与最速下降法结合起来，构造\*\*“牛顿-最速下降算法”\*\*。其基本思想为：**当Hess正定时，采用牛顿方向作为搜索方向；否则，采用负梯度方向作为搜索方向。**

## 算法 ##

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

## 程序 ##

function [x,val,k]=revisenm(fun,gfun,Hess,x0)
    %功能: 用修正牛顿法求解无约束问题:  min f(x)
    %输入: x0是初始点, fun, gfun, Hess 分别是求
    %         目标函数,梯度,Hesse 阵的函数
    %输出:  x, val分别是近似最优点和最优值,  k是迭代次数.
    n=length(x0); maxk=150;
    rho=0.55;sigma=0.4; tau=0.0;
    k=0;  epsilon=1e-5;
    while(k<maxk)
        gk=feval(gfun,x0); % 计算梯度
        muk=norm(gk)^(1+tau);
        Gk=feval(Hess,x0);  % 计算Hesse阵
        Ak=Gk+muk*eye(n);
        dk=-Ak\gk; %解方程组Gk*dk=-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
        x0=x0+rho^mk*dk;
        k=k+1;
    end
    x=x0;
    val=feval(fun,x); 
    %gval=norm(gfun(x));

[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1MDYwNjc0_size_16_color_FFFFFF_t_70]: /images/20221005/ab4c8de2c6914772862873e81a41c576.png
[Link 1]: https://blog.csdn.net/qq_45060674/article/details/118076054?spm=1001.2014.3001.5501
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1MDYwNjc0_size_16_color_FFFFFF_t_70 1]: /images/20221005/f2c00bb7221c4039bd6a803779b84e11.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1MDYwNjc0_size_16_color_FFFFFF_t_70 2]: /images/20221005/438ee4fdd0d34785b969f5b3586b9640.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1MDYwNjc0_size_16_color_FFFFFF_t_70 3]: /images/20221005/d10536842c4a437caab9b41c26b2be72.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1MDYwNjc0_size_16_color_FFFFFF_t_70 4]: /images/20221005/530b91399c2e4f099a47ccdd25e72368.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1MDYwNjc0_size_16_color_FFFFFF_t_70 5]: /images/20221005/6fd190af8c994f9cbb38561eee679f63.png