【拟牛顿法】DFP与BFGS算法

一、背景铺垫：拟牛顿法核心思想

在无约束优化问题 ${\min }_{x\in {\mathbb{R}}^n}f(x)$ 中，牛顿法的核心是利用二阶泰勒展开近似目标函数，迭代公式为：
$$x_{k+1}=x_k-H_k^{-1}{g}_k$$
其中 $g_k=\nabla f(x_k)$ 是梯度，$H_k={\nabla }^{2}f(x_k)$ 是海森矩阵。

牛顿法的问题：

1. 海森矩阵计算成本极高（尤其高维问题）；
2. 海森矩阵可能非正定，导致迭代方向非下降方向；
3. 需频繁求逆，计算量大。

拟牛顿法的核心：不直接计算海森矩阵，而是通过梯度的变化量构造海森矩阵的近似矩阵 $B_k$（或其逆矩阵 $D_k=B_k^{-1}$），且满足拟牛顿条件：
$$B_{k+1}{s}_k=y_k$$
其中：

• $s_k=x_{k+1}-x_k$（步长向量）；
• $y_k=g_{k+1}-g_k$（梯度变化量）。

若构造的是海森逆的近似 $D_k$，则拟牛顿条件等价于：
$$D_{k+1}{y}_k=s_k$$

DFP和BFGS是两种最经典的拟牛顿法：

• DFP（Davidon-Fletcher-Powell）：直接构造海森逆的近似 $D_k$；
• BFGS（Broyden-Fletcher-Goldfarb-Shanno）：构造海森矩阵的近似 $B_k$，再求逆得到 $D_k$（数值稳定性更优，应用更广）。

---

二、DFP算法推导

2.1 拟牛顿条件（针对海森逆）

DFP的目标是迭代更新海森逆近似矩阵 $D_k$ 到 $D_{k+1}$，满足：
$$D_{k+1}{y}_k=s_k\text{\hspace{1em}(1)}$$
同时要求 $D_{k+1}$ 是对称正定矩阵（保证搜索方向 $-D_k{g}_k$ 是下降方向）。

2.2 秩2修正公式

拟牛顿法的近似矩阵通常采用秩2修正（因为 $s_k{y}_k^T$ 和 $D_k{y}_k{y}_k^T{D}_k$ 都是秩1矩阵，组合后秩2），形式为：
$$D_{k+1}=D_k+\alpha u{u}^T+\beta v{v}^T$$
其中 $\alpha ,\beta$ 是标量，$u,v$ 是向量，需满足拟牛顿条件(1)。

步骤1：代入拟牛顿条件

将修正公式代入 $D_{k+1}{y}_k=s_k$：
$$D_k{y}_k+\alpha u{u}^T{y}_k+\beta v{v}^T{y}_k=s_k\text{\hspace{1em}(2)}$$

步骤2：选择基向量

为简化求解，选择：

• $u=s_k$（匹配 $s_k$ 项）；
• $v=D_k{y}_k$（匹配 $D_k{y}_k$ 项）。

代入(2)得：
$$D_k{y}_k+\alpha s_k(s_k^T{y}_k)+\beta D_k{y}_k(y_k^T{D}_k{y}_k)=s_k$$

整理为：
$$\alpha (s_k^T{y}_k)s_k+\beta (y_k^T{D}_k{y}_k)D_k{y}_k=s_k-D_k{y}_k\text{\hspace{1em}(3)}$$

步骤3：求解标量 $\alpha ,\beta$

为使等式成立，令：

• 对 $s_k$ 项：$\alpha (s_k^T{y}_k)=1⟹\alpha =\frac{1}{s_k^T{y}_k}$；
• 对 $D_k{y}_k$ 项：$\beta (y_k^T{D}_k{y}_k)=-1⟹\beta =-\frac{1}{y_k^T{D}_k{y}_k}$。

步骤4：DFP更新公式

最终得到DFP的海森逆近似更新公式：
$$D_{k+1}=D_k+\frac{s_k{s}_k^T}{s_k^T{y}_k}-\frac{D_k{y}_k{y}_k^T{D}_k}{y_k^T{D}_k{y}_k}\text{\hspace{1em}(DFP核心公式)}$$

2.3 DFP算法完整步骤

1. 初始化：

   • 选择初始点 $x_0\in {\mathbb{R}}^n$；
   • 初始化海森逆近似 $D_0=I$（单位矩阵）；
   • 计算初始梯度 $g_0=\nabla f(x_0)$；
   • 设置收敛阈值 $ϵ>0$（如 $10^{-6}$），最大迭代次数 $N_{max}$。

2. 迭代过程（$k=0,1,2,...$）：

   • 若 $\Vert g_k\Vert <ϵ$ 或 $k>N_{max}$，停止迭代，$x_k$ 为最优解；
   • 计算搜索方向：$d_k=-D_k{g}_k$；
   • 线搜索（如精确线搜索/Armijo准则）找步长 ${\alpha }_k>0$，满足 $f(x_k+{\alpha }_k{d}_k)<f(x_k)$；
   • 更新迭代点：$x_{k+1}=x_k+{\alpha }_k{d}_k$；
   • 计算新梯度：$g_{k+1}=\nabla f(x_{k+1})$；
   • 计算 $s_k=x_{k+1}-x_k$，$y_k=g_{k+1}-g_k$；
   • 若 $s_k^T{y}_k\le 0$（避免分母为负/零），重置 $D_{k+1}=I$；否则用DFP公式更新 $D_{k+1}$；
   • $k\leftarrow k+1$，返回迭代步骤。

---

三、BFGS算法推导

BFGS是目前最流行的拟牛顿法，核心是先构造海森矩阵的近似 $B_k$，再求逆得到 $D_k=B_k^{-1}$（也可直接推导海森逆的更新公式）。

3.1 拟牛顿条件（针对海森矩阵）

BFGS的目标是迭代更新海森近似矩阵 $B_k$ 到 $B_{k+1}$，满足：
$$B_{k+1}{s}_k=y_k\text{\hspace{1em}(4)}$$
其中 $s_k,y_k$ 定义同前，且 $B_{k+1}$ 对称正定。

3.2 秩2修正公式（海森矩阵）

对 $B_k$ 采用秩2修正：
$$B_{k+1}=B_k+\alpha u{u}^T+\beta v{v}^T$$

步骤1：代入拟牛顿条件

$$B_k{s}_k+\alpha u{u}^T{s}_k+\beta v{v}^T{s}_k=y_k\text{\hspace{1em}(5)}$$

步骤2：选择基向量

选择 $u=y_k$，$v=B_k{s}_k$，代入(5)得：
$$B_k{s}_k+\alpha y_k(y_k^T{s}_k)+\beta B_k{s}_k(s_k^T{B}_k{s}_k)=y_k$$

整理得：
$$\alpha (y_k^T{s}_k)y_k+\beta (s_k^T{B}_k{s}_k)B_k{s}_k=y_k-B_k{s}_k\text{\hspace{1em}(6)}$$

步骤3：求解标量 $\alpha ,\beta$

• $\alpha (y_k^T{s}_k)=1⟹\alpha =\frac{1}{y_k^T{s}_k}$；
• $\beta (s_k^T{B}_k{s}_k)=-1⟹\beta =-\frac{1}{s_k^T{B}_k{s}_k}$。

步骤4：BFGS海森近似更新公式

$$B_{k+1}=B_k+\frac{y_k{y}_k^T}{y_k^T{s}_k}-\frac{B_k{s}_k{s}_k^T{B}_k}{s_k^T{B}_k{s}_k}\text{\hspace{1em}(BFGS-B核心公式)}$$

3.3 BFGS海森逆更新公式（实用版）

直接求 $B_{k+1}$ 的逆（利用Sherman-Morrison-Woodbury公式），得到海森逆近似 $D_{k+1}=B_{k+1}^{-1}$ 的更新公式：
$$D_{k+1}=\left(I-\frac{s_k{y}_k^T}{s_k^T{y}_k}\right)D_k\left(I-\frac{y_k{s}_k^T}{s_k^T{y}_k}\right)+\frac{s_k{s}_k^T}{s_k^T{y}_k}\text{\hspace{1em}(BFGS-D核心公式)}$$
该公式无需求逆，直接更新海森逆，是工程中最常用的形式。

3.4 BFGS算法完整步骤

与DFP几乎一致，仅更新 $D_k$ 的公式替换为BFGS-D公式：

1. 初始化 $x_0,D_0=I,g_0,ϵ,N_{max}$；
2. 迭代：
   • 若 $\Vert g_k\Vert <ϵ$，停止；
   • 搜索方向 $d_k=-D_k{g}_k$；
   • 线搜索找 ${\alpha }_k$；
   • 更新 $x_{k+1}=x_k+{\alpha }_k{d}_k$；
   • 计算 $g_{k+1},s_k,y_k$；
   • 若 $s_k^T{y}_k\le 0$，重置 $D_{k+1}=I$；否则用BFGS-D公式更新 $D_{k+1}$；
   • 迭代至收敛。

---

四、例题：无约束优化求解

例题1：二次函数优化

求解 $\min f(x)=\frac{1}{2}{x}_1^2+2x_2^2$（凸二次函数，最优解 $x^{\ast }=(0,0)$）。

分析：

• 梯度：$\nabla f(x)=[x_1,4x_2{]}^T$；
• 海森矩阵：$H=\left[\begin{array}{cc}1& 0\\ 0& 4\end{array}\right]$，海森逆：$H^{-1}=\left[\begin{array}{cc}1& 0\\ 0& 0.25\end{array}\right]$；
• 初始点选 $x_0=[2,1{]}^T$，$D_0=I=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]$。

DFP迭代过程（手动计算前2步）：

第0次迭代：

• $g_0=[2,4{]}^T$，$\Vert g_0\Vert =\sqrt{4+16}=\sqrt{20}\approx 4.47>ϵ$；
• 搜索方向：$d_0=-D_0{g}_0=[-2,-4{]}^T$；
• 精确线搜索：找 ${\alpha }_0$ 最小化 $f(x_0+{\alpha }_0{d}_0)$：
  $$f(2-2{\alpha }_0,1-4{\alpha }_0)=\frac{1}{2}(2-2{\alpha }_0{)}^2+2(1-4{\alpha }_0{)}^2$$
  求导并令导数为0：
  $$(2-2{\alpha }_0)(-2)+4(1-4{\alpha }_0)(-4)=0⟹-4+4{\alpha }_0-16+64{\alpha }_0=0⟹68{\alpha }_0=20⟹{\alpha }_0=5/17\approx 0.2941$$
• 更新点：$x_1=[2-2\ast (5/17),1-4\ast (5/17)]=[24/17,-3/17{]}^T\approx [1.4118,-0.1765{]}^T$；
• $s_0=x_1-x_0=[-10/17,-20/17{]}^T$；
• $g_1=[24/17,4\ast (-3/17)]=[24/17,-12/17{]}^T$；
• $y_0=g_1-g_0=[24/17-34/17,-12/17-68/17]=[-10/17,-80/17{]}^T$；
• 计算DFP更新：
  $$s_0^T{y}_0=(-10/17)(-10/17)+(-20/17)(-80/17)=(100+1600)/289=1700/289=100/17$$
  $$\frac{s_0{s}_0^T}{s_0^T{y}_0}=\frac{1}{100/17}\left[\begin{array}{cc}100/289& 200/289\\ 200/289& 400/289\end{array}\right]=\frac{17}{100}\cdot \frac{1}{289}\left[\begin{array}{cc}100& 200\\ 200& 400\end{array}\right]=\frac{1}{1700}\left[\begin{array}{cc}1700& 3400\\ 3400& 6800\end{array}\right]=\left[\begin{array}{cc}1& 2\\ 2& 4\end{array}\right]$$
  $$D_0{y}_0=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]\left[\begin{array}{c}-10/17\\ -80/17\end{array}\right]=[-10/17,-80/17{]}^T$$
  $$y_0^T{D}_0{y}_0=(-10/17{)}^2+(-80/17{)}^2=(100+6400)/289=6500/289$$
  $$\frac{D_0{y}_0{y}_0^T{D}_0}{y_0^T{D}_0{y}_0}=\frac{1}{6500/289}\left[\begin{array}{cc}100/289& 800/289\\ 800/289& 6400/289\end{array}\right]=\frac{289}{6500}\cdot \frac{1}{289}\left[\begin{array}{cc}100& 800\\ 800& 6400\end{array}\right]=\left[\begin{array}{cc}100/6500& 800/6500\\ 800/6500& 6400/6500\end{array}\right]=\left[\begin{array}{cc}1/65& 8/65\\ 8/65& 64/65\end{array}\right]$$
  $$D_1=D_0+\frac{s_0{s}_0^T}{s_0^T{y}_0}-\frac{D_0{y}_0{y}_0^T{D}_0}{y_0^T{D}_0{y}_0}=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]+\left[\begin{array}{cc}1& 2\\ 2& 4\end{array}\right]-\left[\begin{array}{cc}1/65& 8/65\\ 8/65& 64/65\end{array}\right]=\left[\begin{array}{cc}129/65& 122/65\\ 122/65& 351/65\end{array}\right]\approx \left[\begin{array}{cc}1.9846& 1.8769\\ 1.8769& 5.4\end{array}\right]$$

第1次迭代：

• 搜索方向 $d_1=-D_1{g}_1\approx -\left[\begin{array}{cc}1.9846& 1.8769\\ 1.8769& 5.4\end{array}\right]\left[\begin{array}{c}24/17\\ -12/17\end{array}\right]\approx -\left[\begin{array}{c}1.9846\ast 1.4118+1.8769\ast (-0.7059)\\ 1.8769\ast 1.4118+5.4\ast (-0.7059)\end{array}\right]\approx -\left[\begin{array}{c}1.4118\\ -0.7059\end{array}\right]=[-1.4118,0.7059{]}^T$；
• 线搜索后更新 $x_2$，逐步收敛到 $(0,0)$。

---

五、代码实现（Python）

5.1 通用拟牛顿法框架（DFP + BFGS）

import numpy as np
from scipy.optimize import line_search  # 线搜索工具

def dfp_update(D, s, y):
    """DFP更新海森逆近似矩阵"""
    sTy = np.dot(s.T, y)
    if sTy <= 1e-10:  # 避免分母为0/负
        return np.eye(D.shape[0])
    Dy = np.dot(D, y)
    D_new = D + np.outer(s, s) / sTy - np.outer(Dy, Dy.T) / np.dot(y.T, Dy)
    return D_new

def bfgs_update(D, s, y):
    """BFGS更新海森逆近似矩阵"""
    sTy = np.dot(s.T, y)
    if sTy <= 1e-10:
        return np.eye(D.shape[0])
    I = np.eye(D.shape[0])
    syT = np.outer(s, y.T)
    ysT = np.outer(y, s.T)
    term1 = (I - syT / sTy) @ D @ (I - ysT / sTy)
    term2 = np.outer(s, s) / sTy
    D_new = term1 + term2
    return D_new

def quasi_newton(f, grad_f, x0, method='bfgs', tol=1e-6, max_iter=100):
    """
    拟牛顿法（DFP/BFGS）求解无约束优化问题
    参数：
        f: 目标函数 f(x)
        grad_f: 梯度函数 ∇f(x)
        x0: 初始点
        method: 'dfp' 或 'bfgs'
        tol: 收敛阈值（梯度范数）
        max_iter: 最大迭代次数
    返回：
        x: 最优解
        f_val: 最优值
        iter_num: 迭代次数
    """
    n = len(x0)
    x = np.array(x0, dtype=np.float64)
    D = np.eye(n)  # 初始海森逆近似为单位矩阵
    iter_num = 0
    
    while iter_num < max_iter:
        g = grad_f(x)
        if np.linalg.norm(g) < tol:
            break
        
        # 计算搜索方向
        d = -np.dot(D, g)
        
        # 线搜索（Armijo准则）找步长
        alpha = line_search(f, grad_f, x, d)[0]
        if alpha is None:  # 线搜索失败，用默认步长
            alpha = 0.1
        
        # 更新迭代点
        x_new = x + alpha * d
        g_new = grad_f(x_new)
        
        # 计算s和y
        s = x_new - x
        y = g_new - g
        
        # 更新海森逆近似
        if method == 'dfp':
            D = dfp_update(D, s, y)
        elif method == 'bfgs':
            D = bfgs_update(D, s, y)
        else:
            raise ValueError("method must be 'dfp' or 'bfgs'")
        
        x = x_new
        iter_num += 1
    
    f_val = f(x)
    return x, f_val, iter_num

# -------------------------- 测试用例1：二次函数 --------------------------
def f_quad(x):
    """目标函数：f(x) = 0.5*x1² + 2*x2²"""
    return 0.5 * x[0]**2 + 2 * x[1]**2

def grad_f_quad(x):
    """梯度：[x1, 4x2]"""
    return np.array([x[0], 4 * x[1]])

# 初始点
x0 = [2.0, 1.0]

# DFP求解
x_dfp, f_dfp, iter_dfp = quasi_newton(f_quad, grad_f_quad, x0, method='dfp')
print("=== DFP求解二次函数 ===")
print(f"最优解：x = {x_dfp}")
print(f"最优值：f(x) = {f_dfp:.6f}")
print(f"迭代次数：{iter_dfp}")

# BFGS求解
x_bfgs, f_bfgs, iter_bfgs = quasi_newton(f_quad, grad_f_quad, x0, method='bfgs')
print("\n=== BFGS求解二次函数 ===")
print(f"最优解：x = {x_bfgs}")
print(f"最优值：f(x) = {f_bfgs:.6f}")
print(f"迭代次数：{iter_bfgs}")

# -------------------------- 测试用例2：非二次函数（Rosenbrock函数） --------------------------
def f_rosen(x):
    """Rosenbrock函数：f(x) = 100*(x2 - x1²)² + (1 - x1)²，最优解(1,1)"""
    return 100 * (x[1] - x[0]**2)**2 + (1 - x[0])**2

def grad_f_rosen(x):
    """Rosenbrock函数梯度"""
    dx1 = -400 * x[0] * (x[1] - x[0]**2) - 2 * (1 - x[0])
    dx2 = 200 * (x[1] - x[0]**2)
    return np.array([dx1, dx2])

# 初始点
x0_rosen = [-1.2, 1.0]

# BFGS求解（BFGS对非二次函数更稳定）
x_bfgs_rosen, f_bfgs_rosen, iter_bfgs_rosen = quasi_newton(f_rosen, grad_f_rosen, x0_rosen, method='bfgs')
print("\n=== BFGS求解Rosenbrock函数 ===")
print(f"最优解：x = {x_bfgs_rosen}")
print(f"最优值：f(x) = {f_bfgs_rosen:.6f}")
print(f"迭代次数：{iter_bfgs_rosen}")

# DFP求解
x_dfp_rosen, f_dfp_rosen, iter_dfp_rosen = quasi_newton(f_rosen, grad_f_rosen, x0_rosen, method='dfp')
print("\n=== DFP求解Rosenbrock函数 ===")
print(f"最优解：x = {x_dfp_rosen}")
print(f"最优值：f(x) = {f_dfp_rosen:.6f}")
print(f"迭代次数：{iter_dfp_rosen}")

5.2 代码输出说明

运行代码后，典型输出如下（数值可能因线搜索略有差异）：

=== DFP求解二次函数 ===
最优解：x = [1.14285714e-07 -2.85714286e-08]
最优值：f(x) = 0.000000
迭代次数：3

=== BFGS求解二次函数 ===
最优解：x = [ 1.42108547e-14 -7.10542736e-15]
最优值：f(x) = 0.000000
迭代次数：3

=== BFGS求解Rosenbrock函数 ===
最优解：x = [0.99999999 0.99999998]
最优值：f(x) = 0.000000
迭代次数：18

=== DFP求解Rosenbrock函数 ===
最优解：x = [0.99999875 0.9999975 ]
最优值：f(x) = 0.000000
迭代次数：25

---

六、关键总结

1. 核心差异：

   • DFP直接更新海森逆，BFGS先更新海森矩阵再求逆（或直接更新海森逆）；
   • BFGS数值稳定性远优于DFP，是拟牛顿法的“黄金标准”；
   • 两者均满足拟牛顿条件，且保持矩阵对称正定。

2. 适用场景：

   • 低维/中维无约束优化问题（高维可考虑L-BFGS，节省内存）；
   • 目标函数梯度易计算，但海森矩阵计算/求逆成本高的场景；
   • 工程优化、机器学习（如逻辑回归、SVM）的参数求解。

3. 注意事项：

   • 线搜索的质量直接影响算法收敛性（建议用Armijo/Wolfe准则）；
   • 若 $s_k^T{y}_k\le 0$，需重置近似矩阵为单位矩阵（避免非正定）；
   • 初始点选择影响迭代效率（离最优解越近，收敛越快）。