优化算法笔记01：牛顿法的实现

牛顿法是梯度下降方法的重要成员，其算法流程如下。

对于需要最小化的目标函数

1. 给定初值

2. 代入递推公式进行迭代：

3. 判断终止条件，若满足条件则终止迭代，得到最小值点，终止条件我一般选择如下：

其中，取步长

以函数为例，给出牛顿法的MATLAB代码如下。

function [x,iters]=newton()
% Objective function, first and second Derivatives
f=@(x)3.*x.^4+4;
g=@(x)12.*x.^3;
h=@(x)36.*x.^2;
% Iterative formula
iter=@(x)x-h(x)\g(x);
% set x0
xkd=1000;
% compute x1
xk=iter(xkd);
count=0;
% iterations
while(abs(f(xkd)-f(xk))>abs((1e-8)*f(xkd)))
    count=count+1;
    xkd=xk;
    xk=iter(xkd);
end
x=xk;
iters=count;

输出：

>> [x,iters]=newton()

x =

    0.0052

iters =

    29

如果采用最速下降法，步长不能选择为常量1，否则很可能产生不收敛的状况。一般来说，我会在最速下降法中加入步长收缩的Armijo条件。

将牛顿法的收敛速度与最速下降法做比较，以同一个目标函数为例，下面给出最速下降法的MATLAB代码。

function [x,iters]=sd()
f=@(x)3.*x.^4+4;
g=@(x)12.*x.^3;
% set x0
xkd=1000;
% compute x1
xk=xkd-g(xkd);
count=0;
% iterations
while(abs(f(xkd)-f(xk))>abs((1e-8)*f(xkd)))
    count=count+1;
    alpha=1;
    while(f(xk)-f(xk-alpha*g(xk))<(0.0001*alpha*g(xk)'*g(xk)))
        alpha=alpha*0.8;       
    end
    xkd=xk;
    xk=xk-alpha*g(xk);
end
x=xk;
iters=count;

输出：

>> [x,iters]=sd()

x =

    0.0228

iters =

   106

可见最速下降法的迭代次数远多于牛顿法，得到的最小值点也不如牛顿法准确。

 

总结一下，与最速下降法相比，牛顿法的特点：

1. 迭代次数少、收敛速度快；
2. 得到的最小值点比较准确（至少在我的实验中如此）；
3. 没有选取步长的麻烦；
4. 缺点是需要计算目标函数的二阶梯度，也就是Hessian矩阵，可能计算量较大。

 

下面给出Python实现牛顿法的代码：

def Newton():
    f=lambda x:3*x**4+4
    g=lambda x:12*x**3
    h=lambda x:36*x*x
    iteration=lambda x:x-g(x)/h(x)
    xkd=1000
    xk=iteration(1000)
    c=0
    while abs(f(xkd)-f(xk))>abs(f(xkd)*1e-8):
        c+=1
        xkd=xk
        xk=iteration(xkd)
    return xk,c

x,c=Newton()
print(x,c)

 输出：

0.005215095050846564 29

现在Python代码实现的至少标量函数的牛顿法，出现向量和矩阵运算的情况我之后再研究。