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

在拟牛顿算法中，最常见的就是DFP和BFGS两种算法，它们的思路是类似的，但又略微不同，其中，DFP算法先由Davidon于1959年提出，后经Fletcher和Powell于1963年改进，最后就成了DFP算法(Davidon、Fletcher、Powell)，好了，下面我们就具体的来看看DFP算法的原理和算法流程。

一、DFP算法原理

我们都知道，DFP算法就是拟牛顿算法的一种具体实现，要了解DFP算法， 需要先了解拟牛顿算法，下面我们不妨先来简单回顾一下拟牛顿算法，再来讲DFP算法。

1.1. 拟牛顿算法-回顾

拟牛顿算法的原理我们就不详细说了，只晒一下拟牛顿算法的关键公式。

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

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

而对于，则需满足如下条件：

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

 其中，     ，即第次的迭代量  

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

(2) 是正定矩阵

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

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

1.2. DFP算法

好了，终于可以开始说我们的DFP算法了，它就是拟牛顿算法的一种具体实现，所以我们主要就是看看它是怎么构造拟牛顿算法中的的，好了，表演开始！

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

DFP算法假设由如下三个矩阵组成 ：

则拟牛顿条件变为：

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

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

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

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

1.3. DFP迭代公式-总结

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

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

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

其中，

二、DFP-算法流程

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

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

2.1. DFP-算法流程

DFP算法的具体流程如下：

一、初始化

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. 输出找到的学习率：

三、DFP算法-代码实现    

下面使用DFP算法来求函数的最小值。

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

  ，

DFP算法的具体实现代码如下：

"""
本代码展示DFP拟牛顿法求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                                           # 返回查找结果
														       
# DFP主流程	                                                   
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 = -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  = h@h.T/(h.T@dG)                                         # 计算P
    Q  = -B@dG@dG.T@B/(dG.T@B@dG)                               # 计算Q
    B  = B+P+Q                                                  # 更新二阶矩阵B

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

代码运行结果如下：

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

总结

DFP算法就是拟牛顿算法的一种实现，而拟牛顿算法最核心就是怎么构造二次项矩阵B，在DFP中，将B拆成多项，然后再让它们分别凑出符合拟牛顿条件的B，有了符合条件的B，再按拟牛顿算法的流程再实现就可以了。