【优化】一篇入门之-拟牛顿法-bfgs

DFP和BFGS两种算法都是拟牛顿算法，但一般最常用的还是BFGS算法，BFGS算法的名称其实来自四个人名：Broyden、Fletcher、Goldfarb、Shanno，合起来就是BFGS，今天就来带大家看看BFGS是怎么一回事，开搞！

一、BFGS-拟牛顿算法

BFGS算法是目前最常用、最具效果的拟牛顿算法，要了解BFGS算法， 需要先了解拟牛顿算法，下面简单回顾一下拟牛顿算法的相关公式，相关原理就不再详细讲解了，请参考《拟牛顿算法》。

1.1. 拟牛顿法-回顾

首先，拟牛顿法的迭代公式如下：

其中，是目标函数在处的梯度，而对于，则需满足如下条件：

(1) 需要满足拟牛顿条件：

 其中，     ，即第次的迭代量  

                       ，即第次迭代后梯度的增量

(2) 是正定矩阵

       要求是正定矩阵，是为了保障迭代量一定是下降量

的取值方法不是固定的，只需要满足上述两点条件就可以了，对于不同的拟牛顿算法，就有不同的取值方法，不过一般都是通过一些特殊方法来构造出一个，使它满足以上两点条件。

1.2. BFGS算法

回顾完拟牛顿算法，我们就可以开始来看看BFGS算法了，它就是拟牛顿算法的一种具体实现，所以核心也就在于看它是怎么构造出拟牛顿算法中的的。

在拟牛顿法中，必须保证B满足拟牛顿条件，即：

也就是：

而BFGS算法，则假设由如下三个矩阵组成 ：

因此，拟牛顿条件就变为：

可知，将所在项分别凑出其余两项时，上式就可成立，即：

因此，BFGS算法构造了以下的公式，以满足上述两条公式，从而满足拟牛顿条件

只要简单将它们代入(1)(2) 式，消掉分子分母，就可以验证它们是符合(1)(2) 式的。

当初始值为正定矩阵时，按上式的迭代，会一直都是正定矩阵，具体证明略。

1.3. BFGS算法-迭代公式-总结

好了，说了这么多，有些凌乱了，让我们先来梳理、总结一下BFGS算法的迭代公式。

总的来说，BFGS的迭代公式为：

其中，是目标函数在处的梯度，

而的初始值为一个正定矩阵，它的迭代公式如下：

其中，

二、BFGS-算法流程

下面我们来正式地看看BFGS的具体算法流程。

如我们所知，BFGS算法是一种拟牛顿算法，它提供的迭代量h仅仅是一个下降量，如果只按照h进行迭代，会很不可控，所以在实际应用中，一般需要搭配搜索算法来确定迭代步长。

2.1. BFGS算法流程

BFGS算法的具体流程如下：

一、初始化

1. 初始化解，矩阵

    B必须初始化为一个正定矩阵，例如单位矩阵 

2. 计算初始梯度  

二、迭代

1. 计算迭代量：

2. 线性搜索迭代步长

    (1) 线性搜索迭代步长

       下面我们再来说lr_search的流程

    (2) 如果找到合适的学习率，则更新迭代量：

     (3) 如果没找到合适的学习率，则退出训练

3. 更新

4. 更新梯度

    (1) 先保存当前的G为，再计算当前梯度G

    (2) 计算G的增量：

5. 更新

6. 检查退出条件

    (1) 是否已达到最大迭代次数

    (2) 梯度是否过小

三、输出

    输出最终的解x

2.2. 线性搜索的算法流程

好了，在DFP算法中用到了线性搜索，下面我们再来说说它的算法流程。

的输入如下：

    1. 目标函数：

    2. 当前搜索方向： 

    3. 当前梯度：

的流程如下：

一、设置搜索参数 

    初始学习率：

    设置衰减系数：

    设置最大搜索次数：

    设置搜索阈值：

二、线性搜索学习率

1. 初始化 

    (1) 将学习率置为初始学习率：

    (2) 初始化搜索标记： 

2. 逐步搜索m步：

    计算搜索判断条件：

    如果满足搜索判断条件：

        将搜索标记置为，并退出搜索

    否则：

        下降学习率：

三、输出

    1. 输出是否找到合适的学习率：

    2. 输出找到的学习率：

三、BFGS算法-代码实现

最后的最后，我们来看看BFGS算法的代码如何实现。

我们就以为目标函数好了，然后用BFGS算法来求它的最小值。

由于BFGS算法需要使用目标函数的梯度，所以需要先算出梯度，如下：

  ，

算出了梯度，就可以按算法流程来实现BFGS算法了，BFGS算法的具体实现代码如下：

"""
本代码展示BFGS拟牛顿法求y= (x1-2)^2+(x2-3)^2的最小解
本代码来自《老饼讲解-机器学习》www.bbblearn.com 
"""
import numpy as np
# 目标函数
def f(x):                                                       # 目标函数
    return (x[0,0]-2)**2+(x[1,0]-3)**2                          # 返回目标函数值
														       
# 线性搜索函数                                                 
def lr_search(f,h,G):                                           # 线性搜索函数
    m = 100                                                     # 线性搜索的最大次数
    lr_init = 1                                                 # 初始学习率
    lr_desc = 0.9                                               # 学习率衰减系数
    gamma   = 0.5                                               # 线性搜索的阈值系数
    # 线性搜索                                                 
    is_find = 0                                                 # 学习率是否有效
    lr = lr_init                                                # 初始学习率
    for t in range(m):                                          # 逐步搜索
        is_find = f(x+lr*h)<= f(x)+gamma*lr*G.T@h               # 当前学习率是否有效
        if(is_find):                                            # 如果有效
            break                                               # 退出搜索
        lr = lr * lr_desc                                       # 对学习率误差
    return is_find,lr                                           # 返回查找结果
														       
# BFGS主流程	                                                   
x = np.array([[0],[0]])                                         # 初始化x
B = np.eye(len(x))                                              # 初始化二阶矩阵B
G = np.array([[2*x[0,0]-4], [2*x[1,0]-6]])                      # 计算梯度G 	
for i in range(100):                                            # 最大迭代100次
    # -----计算迭代量并更新x-----
    h = -np.linalg.inv(B)@G                                     # 计算迭代量
    [is_find,lr] = lr_search(f,h,G)                             # 线性搜索学习率         
    if(is_find==0):                                             # 如果没有搜索到有效的学习率
        break                                                   # 退出训练
    h = lr*h                                                    # 更新迭代量
    x = x + h                                                   # 更新x
    # ----------更新B-----------                                         
    Gp = G                                                      # 备份当前的梯度
    G  = np.array([[2*x[0,0]-4], [2*x[1,0]-6]])                 # 计算当前梯度G  
    dG = G-Gp                                                   # 计算梯度的增量
    P  = dG@dG.T/(dG.T@h)                                       # 计算P
    Q  = -B@h@h.T@B/(h.T@B@h)                                   # 计算Q
    B  = B+P+Q                                                  # 更新二阶矩阵B

    print("第",i+1,"轮迭代:x=[",x[0,0],",",x[1,0],"],y=",f(x))  # 打印当前结果
    if((max(abs(G))< 0.001) ):break                             # 如果梯度过小，则退出迭代

代码运行结果如下：

可以看到，经过2轮迭代后，已经达到目标的最小解[2,3]了~

结束语

BFGS算法与DFP算法很像，差异在于DFP算法替代的是牛顿法中的，而BFGS则只是替代，由于BFGS替代的是，所以还需要进一步求逆才能得到，这使得BFGS算法的计算要比DFP更加复杂。

但是，大家公认的是，在实践中BFGS就是比DFP的效果要更加好，因此更常用的是BFGS，而不是DFP，为什么BFGS比DFP更好呢？目前这一点在理论上似乎还没有严谨的证明，老饼没有进一步去深究，如果有好的材料，大家可以推荐给老饼。