牛顿法

关于基本牛顿法、阻尼牛顿法、修正牛顿法的介绍见于我的百度文库（http://wenku.baidu.com/view/d272965af01dc281e53af0b2.html）。

基于Armijo准则的阻尼牛顿法

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

修正牛顿法

function [x,val,k]=revisenm(fun,gfun,Hess,x0)
%功能：用修正牛顿法求解无约束问题
%输入：
%
%输出：
%
n=length(x0);
maxk=150;
beta=0.55;
sigma=0.4;
tau=0.0;                  %
epsilon=1e-5;
k=0;
while k<maxk
    gk=feval(gfun,x0);    %计算梯度
    muk=norm(gk)^(1+tau); %阻尼因子
    Gk=feval(Hess,x0);    %计算Hess矩阵
    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+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,x);
    

上述两个算法用的fun、gfun和Hess需要自己编写相应M文件实现，下面是我自己用于检验的上述算法的M文件。

function f=fun(x)
f=100*(x(1)^2-x(2))^2+(x(1)-1)^2;

function g=gfun(x)
g=[400*x(1)*(x(1)^2-x(2))+2*(x(1)-1),-200*(x(1)^2-x(2))]';

function H=Hess(x);
n=length(x);
H=zeros(n,n);
H=[1200*x(1)^2-400*x(2)+2, -400*x(1);
    -400*x(1),200];