牛顿算法及其改进【阻尼牛顿法、修正牛顿法】

牛顿算法

对于优化函数\(f(x)\),\(x=(x_1;x_2;...;x_n)\),二阶连续可导

在\(x_k\)处泰勒展开，取前三项，即对于优化函数二阶拟合

\[f(x)=f(x_k)+g_k(x-x_k)+\frac{1}{2}(x-x_k)G_k(x-x_k) \]

其中\(g_k=\nabla f(x_k)\),为函数梯度；\(G_k=\nabla^2f(x_k)\)，为函数的Hesse矩阵
当\(G_k\)正定时，上式存在极小值，使得\(\nabla f(x)=0\)

\(\nabla f(x)=g_k+G_k(x-x_k)=0\)
可得牛顿法迭代公式：

\[x_{k+1}=x_k-G_k^{-1}g_k \]

1. 可见，对于牛顿法，需要计算二阶偏导数（Hesse矩阵），且Hesse矩阵必须可逆、正定
2. 并且，牛顿法对于迭代初始值\(x_0\)也有要求，当\(x_0\)距离最优解\(x*\)足够近时，算法才收敛

阻尼牛顿法

阻尼牛顿法解决了第二个问题，使得算法全局收敛
主要方法是引入线搜索技术，使得算法满足收敛性条件，关于线搜索技术
线搜索

算法：

\(step0:\)
给定迭代初始值\(x_0\),和容许误差\(\epsilon\)

\(step1:\)
计算梯度\(g_k=\nabla f(x_k)\)
if \(||g_k||<\epsilon\),break；输出当前\(x_k\)
else 解方程：

\[G_kd_k=-g_k \]

求出迭代方向$d_k$, to step 2

\(step2:\)
利用Armijo准则算法，计算迭代步长\(\alpha_k\)

\[x_{k+1}=x_k+\alpha_k d_k \]

k=k+1;to step 1

修正的牛顿算法

上面算法有前提条件1（Hesse矩阵正定）
为了扩大算法的适用范围，对算法修正，解决条件1 的强制条件，有两种方法

方法1：

对于阻尼牛顿算法求得的迭代方向\(d_k\),检查是否满足收敛条件：

\[g_k^Td_k<0 \]

if 满足
\(d_k=d_k\)
else
\(d_k=-g_k\)

就是牛顿算法和梯度下降算法的混合算法，当牛顿算法求得迭代方向不满足收敛条件，使用负梯度方向为迭代方向

关于为什么不直接使用梯度下降法

可以证明：牛顿算法是二阶收敛，梯度下降线性收敛【牛顿法收敛速度快】

方法2

引进阻尼因子\(\mu_k\),对Hesse矩阵修正,选择适当\(\mu_k\),使得

\[G_k=G_k+\mu_kI \]

正定

方法二matlab代码

function [x,val,fun_t] = revise_newton(fun,gfun,hesse,x0,iterate)
%revise_newton - Description
%
% Syntax: [x,val,fun_t] = revise_newton(fun,gfun,hesse,x0)
%
% 修正牛顿法

    n=length(x0);maxk=iterate;
    rho=0.55;Sigma=0.4;tau=0.0;
    k=0;epsilon=1e-5;
    fun_t=zeros(1,maxk);
    while k<maxk

        k=k+1;
        fun_t(1,k)=fun(x0);

        gk=gfun(x0);
        muk=norm(gk)^(1+tau);
        Gk=hesse(x0);
        Ak=Gk+muk*eye(n);
        dk=-Ak\gk;
        if norm(gk)<epsilon
            break;
        end
        m=0;mk=0;
        while m<20
            if fun(x0+rho^m*dk)<fun(x0)+Sigma*rho^m*gk'*dk
                mk=m;
                break;
            end
            m=m+1;
        end
        x0=x0+rho^mk*dk;
    end

    x=x0;
    val=fun(x);
end

main functin

clc;
close all;

fun=@(x) 100*(x(1)^2-x(2))^2+(x(1)-1)^2;
gfun=@(x) [400*(x(1)^2-x(2))*x(1)+2*(x(1)-1);-200*(x(1)^2-x(2))];
hesse=@(x) [1200*x(1)^2-400*x(2)+2,-400*x(1);-400*x(1),200];

x0=[0;0];
iterate=30;

[x,val,fun_t] = revise_newton(fun,gfun,hesse,x0,iterate);

disp(x);
disp(val);
figure(1);
plot(1:iterate,fun_t);
set(get(gca, 'XLabel'), 'String', 'number of iterations');
set(get(gca, 'YLabel'), 'String', 'function value');

result

conclusion

上述图像为函数值随迭代次数变化，可见收敛速度较快；且计算结果较精确；

reference

《最优化方法及其matlab程序设计》