拟牛顿法之BFGS算法

什么是拟牛顿法？

拟牛顿法是在牛顿法的基础上引入了Hessian矩阵的近似矩阵，避免每次迭代都计算Hessian矩阵的逆，它的收敛速度介于梯度下降法和牛顿法之间。拟牛顿法跟牛顿法一样，也是不能处理太大规模的数据，因为计算量和存储空间会开销很多。
拟牛顿法虽然每次迭代不像牛顿法那样保证是最优化的方向，但是近似矩阵始终是正定的，因此算法始终是朝着最优化的方向在搜索。具有全局收敛性和超线性收敛速度

BFGS公式推导

BFGS(Broyden，Fletcher，Goldfarb，Shanno四个人)算法是使用较多的一种拟牛顿方法，故称为BFGS校正。

将$x$写成$x=(x_1,x_2,\dots,x_n)$。对函数$f(x)$在$x=x_{k+1}$处进行泰勒展开到二阶：

$$f(x)=f(x_{k+1})+f'(x_{k+1})(x-x_{k+1})+\frac{1}{2}f''(x_{k+1})(x-x_{k+1})^2+R_n(x)\\ \approx f(x_{k+1})+f'(x_{k+1})(x-x_{k+1})+\frac{1}{2}f''(x_{k+1})(x-x_{k+1})^2$$
对上式求导并令其为0，由于 $f(x)$中的 $x$是一个向量，$f(x)$对 $x$求导意味着对$x$向量中的每个值求偏导。即， $f(x)$对 $x$的一阶导数为一个向量，对$x$的二阶导数为一个 $n*n$的矩阵
$$f'(x)=\left( \frac{\partial f(x)}{\partial x_1}， \frac{\partial f(x)}{\partial x_2},\dots, \frac{\partial f(x)}{\partial x_n}\right)\\ f''(x)=\left[\frac{\partial^2f(x)}{\partial x_i\partial x_j}\right]_{n*n}\\ $$
求导后得:
$$f'(x)=f'(x_{k+1})+f''(x_{k+1})(x-x_{k+1})$$
即：
$$\nabla f(x_k)=\nabla f(x_{k+1})+G_{k+1}(x_k-x_{k+1})$$
可以化简为:
$$\nabla f(x_{k+1})-\nabla f(x_k)=G_{k+1}(x_k-x_{k+1})$$
令 $B_{k+1}\triangleq G_{k+1}$,则可得: $B_{k+1}(x_k-x_{k+1})=\nabla f(x_{k+1})-\nabla f(x_k)$
在BFGS校正方法中，假设：
$$B_{k+1}=B_k+E_k$$

BFGS校正公式的推导

令$E_k=\alpha u_k u_k^T+\beta v_k v_k^T$，其中$u_k,v_k$均为$n *1$的向量。$y_k=\nabla f(x_{k+1})-\nabla f(x_k),s_k=x_{k+1}-x_k$.

那么$B_{k+1}(x_k-x_{k+1})=\nabla f(x_{k+1})-\nabla f(x_k)$
可以化简为：

$$B_{k+1}s_k=y_k$$
将 $B_{k+1}=B_k+E_k$代入上式得：
$$(B_k+E_k)s_k=y_k$$
将 $E_k=\alpha u_k u_k^T+\beta v_k v_k^T$代入上式得：
$$(B_k+\alpha u_k u_k^T+\beta v_k v_k^T)s_k=y_k$$

即：

$$\alpha u_k (u_k^Ts_k)+\beta v_k (v_k^Ts_k)=y_k-B_k s_k$$

$u_k^Ts_k,v_k^Ts_k$皆为实数，$y_k-B_k s_k$为$n*1$的向量，上式中，参数$\alpha$和$\beta$解的可能性有很多，我们取特殊的情况，假设$u_k=rB_ks_k,v_k=\theta y_k$。则:

$$E_k=\alpha r B_ks_k^TB_k+\beta\theta y_ky_k^T$$
代入上式：
$$\Rightarrow \alpha[(rB_ks_k)^Ts_k](rB_ks_k)+\beta[(\theta y_k)^Ts_k](\theta y_k)=y_k-B_ks_k\\ \Rightarrow [\alpha r^2(s_k^TB_ks_k)+1](B_ks_k)+[\beta \theta^2( y_k^Ts_k)-1](y_k)=0$$
令 $\Rightarrow [\alpha r^2(s_k^TB_ks_k)+1](B_ks_k)=0,\beta \theta^2( y_k^Ts_k)-1=0$,则：
$$\alpha r^2=-\frac{1}{s_k^TB_ks_k}\\ \beta\theta^2=\frac{1}{y_k^Ts_k}$$
最终的BFGS校正公式为：
$$B_{k+1}=B_k-\frac{B_ks_ks_k^TB_k}{s_k^TB_ks_k}+\frac{y_ky_k^T}{y_k^Ts_k}$$

BFGS校正的算法流程

设$B_k$对称正定，$B_{k+1}$由上述的BFGS校正公式确定，那么$B_{k+1}$对称正定的充要条件是$y_k^Ts_k\gt0$。

非精确的一维搜索（线搜索）准则：Armijo搜索准则，搜索准则的目的是为了帮助我们确定学习率，还有其他的一些准则，如Wolfe准则以及精确线搜索等。在利用Armijo搜索准则时并不是都满足上述的充要条件，此时可以对BFGS校正公式做些许改变：

$$B_{k+1}= \begin{cases} B_k, & if \quad y_k^Ts_k\le0 \\ B_k-\frac{B_ks_ks_k^TB_k}{s_k^TB_ks_k}+\frac{y_ky_k^T}{y_k^Ts_k} ,&if \quad y_k^Ts_k\gt0\end{cases}$$

注：在李航写的那本《统计学习方法》中说是正定的，但是并没有说上述情况下会怎么样

算法

1. 给定参数$\delta \in(0,1),\sigma \in(0,0.5)$，初始化点$x_0 \in R^n$，终止误差$0 \le\epsilon \ll1$,初始化对称正定阵$B_0$，通常取为$G(x_o)$或单位阵$I_n$;令$k=0$。
2. 计算$g_k=\nabla f(x_k)$，若$\left \| g_k \right \| \ll \epsilon$,终止，输出$x_k$作为近似极小点。
3. 解线性方程组得解$d_k$:$B_kd=-g_k$.
4. 设$m_k$是满足下列不等式的最小非负整数m:
   $$f(x_k+\delta^m d_k)\le f(x_k)+\sigma \delta^m g_k^Td_k$$
   令$\alpha_k=\delta^{m_k},x_{k+1}=x_k+\alpha_kd_k$.
5. 由BFGS校正公式确定$B_{k+1}$
6. 令$k=k+1$，转向步骤“2”

求解具体优化问题

求解无约束优化问题:

$$\min f(s)=100(x_1^2-x_2)^2+(x_1-1)^2,x=(x_1,x_2)^T\in R^2$$

#coding:UTF-8  
'''
Created on 2017年4月20日

@author: zhangdapeng
'''
from numpy import *    
import matplotlib.pyplot as plt
from numpy.matrixlib.defmatrix import mat
#fun  原始函数
def fun(x):  
    return 100 * (x[0,0] ** 2 - x[1,0]) ** 2 + (x[0,0] - 1) ** 2  

#对x1，x2求导后的函数  
def gfun(x):  
    result = zeros((2, 1))  
#     对x1求导
    result[0, 0] = 400 * x[0,0] * (x[0,0] ** 2 - x[1,0]) + 2 * (x[0,0] - 1)  
    result[1, 0] = -200 * (x[0,0] ** 2 - x[1,0])  #对x2求导
    return result  
def bfgs(fun, gfun, x0):  
    result = []  
    maxk = 500  
    delta = 0.55  
    sigma = 0.4  
    m = shape(x0)[0]  
    Bk = eye(m)  
    k = 0
    epsilon=1e-10
    while (k < maxk):  
        gk = mat(gfun(x0))#计算梯度 ，mat函数将数组转化为矩阵。
#         print(gk) 
#         print(linalg.norm(gk,1))
        #axis=0,沿着纵轴方向
        if linalg.norm(gk,1)<epsilon:
            break
        dk = mat(-linalg.solve(Bk, gk))  #解矩阵方程Bk*x=gk得到x
        m = 0  
        mk = 0  
        while (m < 20):  
            newf = fun(x0 + delta ** m * dk)  
            oldf = fun(x0)  
            if (newf < oldf + sigma * (delta ** m) * (gk.T * dk)[0,0]):  
                mk = m  
                break  
            m = m + 1  

        #BFGS校正  
        x = x0 + delta ** mk * dk  
        sk = x - x0  
        yk = gfun(x) - gk  
#         print(math.isnan(yk.T * sk))
        if (yk.T * sk > 0): 
            Bk = Bk - (Bk * sk * sk.T * Bk) / (sk.T * Bk * sk) + (yk * yk.T) / (yk.T * sk)  

        k = k + 1  
        x0 = x  
        result.append(fun(x0))  

    return result  

#初始化x0  
x0 = mat([[-1.2], [1]])  
result = bfgs(fun, gfun, x0)  
print("result:",result[-1])
n = len(result)  
ax = plt.figure().add_subplot(111)  
x = arange(0, n, 1)  
y = result  
ax.plot(x,y)  

plt.show()  

输出：

result: 2.68262011582e-28

输出图片：

http://blog.csdn.net/google19890102/article/details/45867789

http://blog.csdn.net/acdreamers/article/details/44664941
http://www.codelast.com/%E5%8E%9F%E5%88%9B%E7%94%A8%E4%BA%BA%E8%AF%9D%E8%A7%A3%E9%87%8A%E4%B8%8D%E7%B2%BE%E7%A1%AE%E7%BA%BF%E6%90%9C%E7%B4%A2%E4%B8%AD%E7%9A%84armijo-goldstein%E5%87%86%E5%88%99%E5%8F%8Awo/