梯度法-最优化算法设计与分析之一

引言：无约束最优化问题的一般形式如下 ：。例如无约束优化问.此二维空间的最优化问题该如何求解。从图形上反应的图形为如图一所示：
缩小图形尺寸，得到的微缩图形如图二所示。从图像上可以看出，最优解为x*=(1,1)，最优值为f(x*)=0。
梯度法 ： 梯度法是求解无约束优化问题最简单和最古老的方法之一。设f(x)在附近连续可微，
为搜索方向向量，由泰勒展开式，得，那么目标函数f(x)在处沿方向下降的变化率为，式中为与的夹角。显然，对于不同的方向，函数变化率取决于它与夹角的余弦值。要使变化率最小，只有，即=时才能达到，亦即应该取。即负梯度方向是目标函数f(x)在当前点的最速下降方向，因此梯度法也称为最速下降法。最速下降法的算法流程如下：
算法一：最速下降法
step one: 选取初始点，容许误差。令
step two:计算。若，停止计算，输出作为近似极小点。
step three： 取方向。
step four:由线索方法确定步长因子
step five:令转步骤2.
其中的确定可以使用Armijo准则。
具体算法如下：
算法二：Armijo准则
给定参数对，若不等式成立，则置

牛顿法：牛顿法也是求解无约束优化问题最早使用的经典算法之一，其基本思想是：用迭代点处的一阶导数(梯度)和二阶导数(Hesse阵)对目标函数进行二次函数近似，然后把二次函数的极小点作为新的迭带点，并不断重复这一过程，直至求得满足精度的近似极小点。

牛顿法的具体算法步骤如下：
算法三：基本牛顿法

由于实际问题的精确极小点一般是不知道的，因此，初始点的选取给算法的实际操作带来了很大的困难。为了克服这一困难，可引入线搜索方法以得到大范围收敛的算法，即所谓的阻尼牛顿法。具体步骤如下
算法四：阻尼牛顿法

实验结果：
实验分别选取五组初始点 为，终止准则为
梯度法的数值结果：
阻尼牛顿法的数值结果
实验结果可知：阻尼牛顿法的实验效果更好，同样的初始点，迭代次数均小于梯度下降法，且目标函数值更接近于真实的最优值0.
附录：
梯度法程序：

function[k,x,val]=grad(fun,gfun,x0,epsilon)
%功能： 梯度法求解无约束优化问题：min f(x)
%输入： fun,gfun分别是目标函数及其梯度，x0是初始点，
% epsilon为容许误差
%输出：k是迭代次数，x,val 分别是近似最优点和最优值
maxk=5000;%最大的迭代次数
beta=0.5;sigma=0.4;
k=0;
while(k<maxk)
     gk=feval(gfun,x0);%计算梯度
     dk=-gk;%计算搜索方向
     if(norm(gk)<epsilon),break;end %检验终止准则
     m=0 ； mk=0;
     while(m<20) %用Armijo搜索求步长
        if（feval(fun,x0+beta^m*dk)<=feval(fun,x0)+sigma*beta^m*gk'*dk）
            mk=m ; break;
        end
        m=m+1;
     end
     x0=x0+beta^mk*dk;
     k=k+1;
   end
   x=x0;
   val = feval(fun,x0);    

阻尼牛顿法程序：

function[k,x,val]=dampnm(fun,gfun,Hess,x0,epsilon)
%功能：阻尼牛顿法求解无约束优化问题： min f(x)
%输入： fun,gfun,Hess 分别是目标函数及其梯度和Hesson阵，
%      x0是初始点， epsilon为容许误差
%  输出：k是迭代次数，x,val 分别是近似最优解和最优值
maxk=5000;%最大迭代次数
beta=0.5;sigma=0.4;k=0;
while(k<maxk)
    gk=feval(gfun,x0);%计算梯度
    Gk=feval(Hess,x0);%计算Hesson阵
    dk=-Gk\gk;%解方程组 Gk*dk=-gk，计算搜索方向
    if(norm(gk)<epsilon),break;end %检验终止准则
    m=0;mk=0;
    while(m<20)%用Armijo搜索求步长
      if(feval(fun,x0+beta^m*dk)<=feval(fun,x0)+sigma*beta^m*gk'*dk)
           mk=m;break;
      end
      m=m+!;
    end
    x0=x0+beta^m*dk;k=k+1;
end
x=x0;
val=feval(fun,x);           

其中：

%目标函数
function f = fun(x)
f = 4*(x(1)^2 - x(2))^2 + 3*(x(1)-1)^2;
%梯度
function gf = gfun(x)
gf = [16*x(1)*(x(1)^2 - x(2)) + 6*(x(1)-1) ; -8*(x(1)^2-x(2))];
%Hesson阵
function He = Hess(x)
He=[48*x(1)^2-16*x(2)+6,   -16*x(1);
    -16*x(1),                   8]

参考文献： 陈宝林 .最优化理论与算法设计 第二版
马昌凤.柯艺芬.谢亚君 最优化计算方法及其MATLAB程序实现