修正牛顿法及其matlab实现

由上篇文章知道，牛顿法具有二阶收敛速度，收敛速度比较快，但是用这个算法的一个必要条件是，目标函数的Hess阵必须是正定的，否则难以保证他的下降方向是正确的，所以我们可以采用牛顿法与最速下降法结合，即修正牛顿法。

修正牛顿法的基本原理是：当Hess阵正定时，采用牛顿方向作为搜素方向；若不正定，则采用最速下降法，即负梯度方向作为搜索方向。算法的步骤为：

步0：确定终止误差e=（0~1）,给定初始值，rho=(0~1),sigma=(0~0.5)，给定初始值x0，令k=0

步1：计算gk=f(xk)的值，若||gk||<=e停止运算，输出xk作为近似最优解，否则转步2

步2：计算Gk=^2f(xk)，并解方程Gk*dk+gk=0，求出dk，并且要求gk'*dk<0(就是满足正定)，若满足，则转步3,。否则令dk=-gk，转步3

步3：用线搜索技术来确定步长因子k

步4：令Xk+1=xk+k*dk，k=k+1，转步1

具体代码与前两篇文章类似，这里不再赘述，下面给出另一种修正牛顿法，即引进一个阻尼因子uk>=0,在每一迭代步选取适当的参数uk，使得Ak=G（xk）+uk*I正定，I是单位矩阵，具体的算法步骤为：

步0：确定终止误差e=（0~1）,给定初始值，=(0~1),=(0~0.5)，=[0,1],给定初始值x0，令k=0

步1：计算gk=f(xk)，uk=||gk||^(1+),若||gk||<=e,停止运算，输出xk为近似最优解，否则转步2

步2： 计算Hess阵，Gk=^2f(xk)，并解线性方程（Gk+uk*I）*dk=-gk,解得dk

步3：利用Armijo搜索技术确定步长，步长k=^mk，m的值从0开始

若满足不等式f(xk+ ^m*dk)<=f(xk)+*^m*gk'dk

则 mk=m，步长 k=^mk，若不满足上式，则m=m+1，直到满足上述不等式为止

步4：得到步长后令Xk+1=xk+ k*dk，k=k+1，转步1

 具体代码如下：

1.修正牛顿法函数

function [x,val,k] = revisenm(fun,gfun,Hess,x0)
%功能：用修正牛顿法来求解无约束问题：minif（x）
%revisenm是修正的意思
%输入：fun,gfun,Hess分别是目标函数，梯度，二阶导数，x0是初始点
%输出：x，val分别是近似极小点和近似最优值，k是迭代次数
maxk=100;
rho=0.55;
sigma=0.4;
tau=0.0;%这是引进的阻尼因子
n=length(x0);%这是为了后面生成单位矩阵做铺垫
k=0;
e=1e-5;%精度要求
while(k<maxk)
    gk=feval(gfun,x0);
    muk=norm(gk)^(1+tau);
    if(norm(gk)<=e),break;end
    Gk=feval(Hess,x0);
    Ak=Gk+muk*eye(n);
    dk=-Ak\gk;
    m=0;
    mk=0;
    while(m<20)
        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+dk*rho^m;
     k=k+1;      
end
x=x0;
val=feval(fun,x0);
end

2.目标函数

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

 3.目标函数梯度

function  g=gfun(x)
%目标函数的梯度
g=[400*x(1)*(x(1)^2-x(2))+2*(x(1)-1),-200*(x(1)^2-x(2))]';
end

4.目标函数Hess阵

function f = Hess(x)
%n=length(x);
%f=zero(n,n);
f=[1200*x(1)^2-400*x(2)+2,-400*x(1);
    -400*x(1),             200];
end

5.运行结果

 

 与上一篇文章比较，修正牛顿法的收敛速度比不上阻尼牛顿法，这是因为（不告诉你，噗）。

下一篇文章将要介绍共轭梯度法，共轭梯度法具有超线性的收敛速度，且算法比较简单，避免了Hess阵的计算，所以共轭梯度法是一种求解无约束问题的比较实用的算法，敬请后续更新...........