牛頓算法及其改進【阻尼牛頓法、修正牛頓法】

牛頓算法

對於優化函數\(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程序設計》