拟牛顿法学习

拟牛顿法

拟牛顿法的基本思想就是将海森矩阵${\mathbf{G}}_k={\nabla }^{2}f\left({\mathbf{x}}_k\right)$用一个近似矩阵${\mathbf{B}}_k$来代替

通常${\mathbf{B}}_k$应该有3个特点

1. ${\mathbf{B}}_k\approx {\mathbf{G}}_k$使相应的算法产生的方向近似于牛顿方向，以确保算法具有较快的收敛速度
2. 对于所有的$k$,${\mathbf{B}}_k$对称正定，从而使得算法产生的方向是函数$f$在${\mathbf{x}}_k$处的下降方向
3. ${\mathbf{B}}_k$更新规则相对简单，即通常采用一个秩1矩阵或者秩2矩阵来矫正

设$f:{\mathbb{R}}^n\to \mathbb{R}$在开集$D\subset {\mathbb{R}}^n$上二次连续可微，则
$$f\left(\mathbf{x}\right)\approx \left({\mathbf{x}}_{k+1}\right)+{\mathbf{g}}_{k+1}^T\left(\mathbf{x}-{\mathbf{x}}_{k+1}\right)+\frac{1}{2}{\left(\mathbf{x}-{\mathbf{x}}_{k+1}\right)}^T{\mathbf{G}}_{k+1}\left(\mathbf{x}-{\mathbf{x}}_{k+1}\right)$$
求导
$$\mathbf{g}\left(\mathbf{x}\right)={\mathbf{g}}_{k+1}+{\mathbf{G}}_{k+1}\left(\mathbf{x}-{\mathbf{x}}_{k+1}\right)$$

令$\mathbf{x}={\mathbf{x}}_k,{\mathbf{s}}_k={\mathbf{x}}_{k+1}-{\mathbf{x}}_k,{\mathbf{y}}_k={\mathbf{g}}_{k+1}-{\mathbf{g}}_k$,则有
$${\mathbf{G}}_{k+1}{\mathbf{s}}_k\approx {\mathbf{y}}_k$$

现在要构造近似的矩阵${\mathbf{B}}_k$,即
$${\mathbf{B}}_{k+1}{\mathbf{s}}_k={\mathbf{y}}_k$$
令${\mathbf{H}}_{k+1}={\mathbf{B}}_{k+1}^{-1}$,则有
$${\mathbf{H}}_{k+1}{\mathbf{y}}_k={\mathbf{s}}_k$$

根据第三个特点,可以令
$${\mathbf{B}}_{k+1}={\mathbf{B}}_k+{\mathbf{E}}_k,{\mathbf{H}}_{k+1}={\mathbf{H}}_k+{\mathbf{D}}_k$$
其中${\mathbf{E}}_k,{\mathbf{D}}_k$是秩1矩阵或者秩2矩阵

秩1矫正

设${\mathbf{E}}_k=\alpha {\mathbf{u}}_k{\mathbf{u}}_k^T$,其中$\alpha \in \mathbb{R},{\mathbf{u}}_k\in {\mathbb{R}}^n$
于是
$$\left({\mathbf{B}}_k+\alpha {\mathbf{u}}_k{\mathbf{u}}_k^T\right){\mathbf{s}}_k={\mathbf{y}}_k$$
即
$$\alpha \left({\mathbf{u}}_k^T{\mathbf{s}}_k\right){\mathbf{u}}_k={\mathbf{y}}_k-{\mathbf{B}}_k{\mathbf{s}}_k$$
也就是说${\mathbf{u}}_k$平行于${\mathbf{y}}_k-{\mathbf{B}}_k{\mathbf{s}}_k$,即存在$\beta$,使得${\mathbf{u}}_k=\beta \left({\mathbf{y}}_k-{\mathbf{B}}_k{\mathbf{s}}_k\right)$

因此
$${\mathbf{E}}_k=\alpha {\beta }^2\left({\mathbf{y}}_k-{\mathbf{B}}_k{\mathbf{s}}_k\right){\left({\mathbf{y}}_k-{\mathbf{B}}_k{\mathbf{s}}_k\right)}^T$$
以及
$$\alpha {\beta }^2\left[{\left({\mathbf{y}}_k-{\mathbf{B}}_k{\mathbf{s}}_k\right)}^T{\mathbf{s}}_k\right]\left({\mathbf{y}}_k-{\mathbf{B}}_k{\mathbf{s}}_k\right)=\left({\mathbf{y}}_k-{\mathbf{B}}_k{\mathbf{s}}_k\right)$$
如果${\left({\mathbf{y}}_k-{\mathbf{B}}_k{\mathbf{s}}_k\right)}^T{\mathbf{s}}_k\ne 0$
则
$$\alpha {\beta }^2\left[{\left({\mathbf{y}}_k-{\mathbf{B}}_k{\mathbf{s}}_k\right)}^T{\mathbf{s}}_k\right]=1$$
于是
$$\alpha {\beta }^2=\frac{1}{{\left({\mathbf{y}}_k-{\mathbf{B}}_k{\mathbf{s}}_k\right)}^T{\mathbf{s}}_k}$$
然后
$${\mathbf{E}}_k=\frac{\left({\mathbf{y}}_k-{\mathbf{B}}_k{\mathbf{s}}_k\right){\left({\mathbf{y}}_k-{\mathbf{B}}_k{\mathbf{s}}_k\right)}^T}{{\left({\mathbf{y}}_k-{\mathbf{B}}_k{\mathbf{s}}_k\right)}^T{\mathbf{s}}_k}$$
所以
$${\mathbf{B}}_{k+1}={\mathbf{B}}_k+\frac{\left({\mathbf{y}}_k-{\mathbf{B}}_k{\mathbf{s}}_k\right){\left({\mathbf{y}}_k-{\mathbf{B}}_k{\mathbf{s}}_k\right)}^T}{{\left({\mathbf{y}}_k-{\mathbf{B}}_k{\mathbf{s}}_k\right)}^T{\mathbf{s}}_k}$$
类似地
$${\mathbf{H}}_{k+1}={\mathbf{H}}_k+\frac{\left({\mathbf{s}}_k-{\mathbf{H}}_k{\mathbf{y}}_k\right){\left({\mathbf{s}}_k-{\mathbf{H}}_k{\mathbf{y}}_k\right)}^T}{{\left({\mathbf{s}}_k-{\mathbf{H}}_k{\mathbf{y}}_k\right)}^T{\mathbf{y}}_k}$$

所以秩1矫正步骤
初始化：选定${\mathbf{x}}_0$，${\mathbf{H}}_0$(通常选择${\mathbf{I}}_n$)
步骤：
1.若$\Vert {\mathbf{g}}_k\Vert \le ϵ$,停止，输出${\mathbf{x}}_k$
2.计算${\mathbf{d}}_k=-{\mathbf{H}}_k{\mathbf{g}}_k$
3.求步长${\alpha }_k$
4.${\mathbf{x}}_{k+1}={\mathbf{x}}_k+{\alpha }_k{\mathbf{d}}_k$,矫正${\mathbf{H}}_k$得到${\mathbf{H}}_{k+1}$
5.k=k+1,转1

代码：

function [x,fun_val]=sr1(f,g,x0,s,alpha,...
    beta,epsilon)
% sr1 with backtracking stepsize rule
%
% INPUT
%=======================================
% f ......... objective function
% g ......... gradient of the objective function
% x0......... initial point
% s ......... initial choice of stepsize
% alpha ..... tolerance parameter for the stepsize selection
% beta ...... the constant in which the stepsize is multiplied
%             at each backtracking step (0<beta<1)
% epsilon ... tolerance parameter for stopping rule
% OUTPUT
%=======================================
% x ......... optimal solution (up to a tolerance)
%             of min f(x)
% fun_val ... optimal function value
grad=g(x0);
fun_val=f(x0);
iter=0;
H_k=eye(length(x0));
while (norm(grad)>epsilon)
    iter=iter+1;
    dk=H_k*grad;
    t=s;
    while (fun_val-f(x0-t*dk)<alpha*t*grad'*dk)
        t=beta*t;
    end
    x=x0-t*dk;
    s_k=x-x0;
    y_k=g(x)-grad;
    H_k=H_k+(s_k-H_k*y_k)*(s_k-H_k*y_k)'/((s_k-H_k*y_k)'*y_k);
    fun_val=f(x);
    grad=g(x);
    x0=x;
    fprintf('iter_number = %3d norm_grad = %2.6f fun_val = %2.6f \n',...
        iter,norm(grad),fun_val);
end