拟牛顿法

拟牛顿法

小象实战讲义 · 人工智能数学基础

在上一节中，我们学习了牛顿法这一强大的二阶优化算法，它利用目标函数的二阶信息（Hessian矩阵）实现了快速的收敛速度。然而，牛顿法在实际应用中面临两大挑战：Hessian矩阵求逆的高计算成本，以及迭代点远离最优解时Hessian矩阵可能非正定，导致算法失效。本节介绍的拟牛顿法（Quasi-Newton Methods）正是为了解决这些问题而诞生的一类方法。它通过构造一个正定矩阵来近似Hessian矩阵或其逆矩阵，在保持较快收敛速度的同时，显著降低了每步迭代的计算复杂度。学完本节，你将掌握最主流的拟牛顿法—BFGS方法的原理、更新公式及其实现框架。

💡 核心导读

本节将带你深入理解拟牛顿法的核心思想与实现，重点包括：

1. 从牛顿法的痛点出发：回顾牛顿法在计算复杂度和矩阵正定性方面的两大缺陷，理解拟牛顿法提出的动机。
2. 拟牛顿法的基本思路：学习如何用一个易于更新的正定矩阵来近似Hessian矩阵或其逆矩阵，并理解关键的“拟牛顿条件”。
3. BFGS更新公式的推导：了解如何通过求解一个矩阵范数最小化问题，得到BFGS这一最有效的拟牛顿更新公式。
4. BFGS算法框架：掌握完整的BFGS算法流程，包括方向计算、步长搜索和矩阵更新。
5. 收敛性分析：了解BFGS方法在凸函数下的全局收敛性，以及在目标函数二阶导数满足Lipschitz连续条件下的局部超线性收敛性。

从牛顿法的痛点到拟牛顿法的思路

我们首先回顾牛顿法。在点 $x_k$ 处，牛顿法的搜索方向 $p_k$ 由下式给出： $$p_k = -[\nabla^2 f(x_k)]^{-1} \nabla f(x_k)$$ 其中 $\nabla^2 f(x_k)$ 是目标函数 $f$ 在 $x_k$ 处的 Hessian 矩阵。

牛顿法虽然具有二阶收敛速度，但存在两个显著的缺点：

1. 计算复杂度高：每一步都需要计算 Hessian 矩阵 $\nabla^2 f(x_k)$ 及其逆矩阵。矩阵求逆在最坏情况下是 $O(n^3)$ 复杂度的运算（$n$ 为变量维度），对于高维问题计算代价巨大。
2. Hessian 矩阵可能非正定：当迭代点 $x_k$ 距离最优解 $x^*$ 较远时，Hessian 矩阵 $\nabla^2 f(x_k)$ 可能不是正定矩阵。这会导致牛顿方向 $p_k$ 可能不是下降方向，破坏算法的收敛性，同时非正定矩阵求逆的数值稳定性也较差。

拟牛顿法的核心思想是：放弃直接计算复杂的 Hessian 矩阵，转而构造一个正定矩阵 $B_k$ 来近似 Hessian 矩阵 $\nabla^2 f(x_k)$，或者构造一个正定矩阵 $H_k$ 来近似 Hessian 矩阵的逆 $[\nabla^2 f(x_k)]^{-1}$。这样，搜索方向变为： $$p_k = -H_k \nabla f(x_k) \quad \text{或等价地，通过求解} \quad B_k p_k = -\nabla f(x_k)$$ 其中 $H_k \approx [\nabla^2 f(x_k)]^{-1}$， $B_k \approx \nabla^2 f(x_k)$，且 $H_k = B_k^{-1}$。

接下来的关键问题是：如何构造并更新这个近似矩阵 $B_k$ 或 $H_k$，使得它既能较好地逼近真实的二阶信息，又能在迭代中高效更新（避免 $O(n^3)$ 的求逆运算）并保持正定性？

拟牛顿条件与 BFGS 更新公式的推导

拟牛顿法的更新策略基于对目标函数局部二次模型的近似。在 $x_k$ 处，我们用二次模型 $m_k(p)$ 近似 $f(x_k + p)$： $$f(x_k + p) \approx m_k(p) = f(x_k) + \nabla f(x_k)^T p + \frac{1}{2} p^T B_k p$$ 其中 $B_k$ 是对 Hessian 矩阵的当前估计。根据牛顿法思想，我们沿着 $p_k = -H_k \nabla f(x_k)$（$H_k = B_k^{-1}$）方向进行搜索，并通过线搜索确定步长 $\alpha_k$，得到下一个迭代点： $$x_{k+1} = x_k + \alpha_k p_k = x_k + s_k, \quad \text{其中 } s_k = \alpha_k p_k = x_{k+1} - x_k.$$

在 $x_{k+1}$ 处，我们同样可以建立一个二次模型 $m_{k+1}(p)$。一个合理的要求是，这两个相邻的模型 $m_k$ 和 $m_{k+1}$ 在它们共同的定义域内应该具有某种一致性。一个自然的选择是要求这两个模型在 $x_k$ 处具有相同的一阶导数。即： $$\nabla m_{k+1}(-s_k) = \nabla m_k(0).$$ 对 $m_{k+1}(p) = f(x_{k+1}) + \nabla f(x_{k+1})^T p + \frac{1}{2} p^T B_{k+1} p$ 求导并令 $p = -s_k$，可得： $$\nabla m_{k+1}(-s_k) = \nabla f(x_{k+1}) - B_{k+1} s_k.$$ 而 $\nabla m_k(0) = \nabla f(x_k)$。代入等式 $\nabla m_{k+1}(-s_k) = \nabla m_k(0)$，得到： $$\nabla f(x_{k+1}) - B_{k+1} s_k = \nabla f(x_k).$$ 整理后，我们得到著名的 拟牛顿条件（或割线方程）： $$B_{k+1} s_k = y_k, \quad \text{其中 } y_k = \nabla f(x_{k+1}) - \nabla f(x_k).$$ 如果使用逆 Hessian 近似 $H_{k+1} = B_{k+1}^{-1}$，则拟牛顿条件等价地写为： $$H_{k+1} y_k = s_k.$$ 这个条件要求新的近似矩阵 $B_{k+1}$（或 $H_{k+1}$）必须满足一次线性方程，它反映了梯度变化 $y_k$ 与变量变化 $s_k$ 之间应满足的二次模型关系。

现在，我们需要一个更新规则，从旧的近似 $B_k$（或 $H_k$）得到新的 $B_{k+1}$（或 $H_{k+1}$），使其满足拟牛顿条件。同时，我们希望更新是“最小”的改动，即 $B_{k+1}$ 与 $B_k$ 在某种度量下尽可能接近，并且保持对称正定性。

对于逆 Hessian 近似 $H_k$ 的更新，BFGS 方法通过求解以下优化问题来得到 $H_{k+1}$： $$\begin{aligned} \min_{H} &\quad | H - H_k |$${W} \ \text{s.t.} &\quad H = H^T, \quad H y_k = s_k. \end{aligned}Hmin​s.t.​∥H−Hk​∥W​H=HT,Hyk​=sk​.​ 这里 $| \cdot |${W} 是一种加权 Frobenius 范数。求解这个优化问题（具体推导超出本课程范围），即可得到 BFGS 更新公式： $$H_{k+1} = \left(I - \rho_k s_k y_k^T \right) H_k \left(I - \rho_k y_k s_k^T \right) + \rho_k s_k s_k^T,$$ 其中 $\rho_k = \dfrac{1}{y_k^T s_k}$。为了保证更新有意义，我们要求 $y_k^T s_k > 0$，这可以通过使用 Wolfe 条件的线搜索来保证。

BFGS 公式具有一个非常好的结构：它是 $H_k$ 的秩-2 校正。如果 $H_k$ 对称正定，且 $y_k^T s_k > 0$，则可以证明 $H_{k+1}$ 也保持对称正定性。这个公式只涉及向量内积和外积运算，计算复杂度为 $O(n^2)$，远低于直接求逆的 $O(n^3)$。

相应地，也存在对 $B_k$ 的直接 BFGS 更新公式： $$B_{k+1} = B_k - \frac{B_k s_k s_k^T B_k}{s_k^T B_k s_k} + \frac{y_k y_k^T}{y_k^T s_k}.$$ 在实际算法中，更常用的是逆 Hessian 近似 $H_k$ 的更新公式，因为它可以直接用于计算搜索方向 $p_k = -H_k \nabla f(x_k)$，无需解线性方程组。

BFGS 算法框架与实现

基于上述推导，我们可以给出完整的 BFGS 算法框架。

算法：BFGS 拟牛顿法

1. 初始化：给定初始点 $x_0$，收敛阈值 $\epsilon > 0$，初始逆 Hessian 近似 $H_0$（通常取为单位矩阵 $I$）。
2. 循环：对于 $k = 0, 1, 2, \dots$，执行以下步骤，直到 $|\nabla f(x_k)| < \epsilon$： a. 计算搜索方向： $p_k = -H_k \nabla f(x_k)$。 b. 线搜索：使用满足 Wolfe 条件的非精确线搜索，计算步长 $\alpha_k > 0$。 c. 更新迭代点： $x_{k+1} = x_k + \alpha_k p_k$。 d. 计算差分向量： $s_k = x_{k+1} - x_k$， $y_k = \nabla f(x_{k+1}) - \nabla f(x_k)$。 e. 更新逆 Hessian 近似： 计算 $\rho_k = 1 / (y_k^T s_k)$，然后按 BFGS 公式更新： $$H_{k+1} = (I - \rho_k s_k y_k^T) H_k (I - \rho_k y_k s_k^T) + \rho_k s_k s_k^T.$$
3. 输出：近似最优解 $x_k$。

初始矩阵 $H_0$ 的选择很重要，单位矩阵是一个简单而常用的选择，它意味着算法最初采用最速下降方向。 Wolfe 条件线搜索能确保 $y_k^T s_k > 0$，从而保证 $H_k$ 的正定性。

下面我们用 Python 代码演示 BFGS 更新公式的计算，并验证其满足拟牛顿条件。

import numpy as np

def bfgs_update(H_k, s_k, y_k):
    """
    根据BFGS公式更新逆Hessian近似矩阵 H_k。
    
    参数:
    H_k: 当前的逆Hessian近似矩阵 (n x n)
    s_k: 变量增量向量， x_{k+1} - x_k (n,)
    y_k: 梯度增量向量， grad_{k+1} - grad_k (n,)
    
    返回:
    H_{k+1}: 更新后的逆Hessian近似矩阵
    """
    s_k = s_k.reshape(-1, 1)  # 转换为列向量 (n, 1)
    y_k = y_k.reshape(-1, 1)  # 转换为列向量 (n, 1)
    
    rho_k = 1.0 / (y_k.T @ s_k).item()  # 标量
    
    I = np.eye(H_k.shape[0])
    term1 = I - rho_k * (s_k @ y_k.T)
    term2 = I - rho_k * (y_k @ s_k.T)
    H_next = term1 @ H_k @ term2 + rho_k * (s_k @ s_k.T)
    
    return H_next

# 示例：验证BFGS更新满足拟牛顿条件 H_{k+1} y_k = s_k
np.random.seed(42)
n = 3
# 随机生成初始 H_k (对称正定)
A = np.random.randn(n, n)
H_k = A @ A.T + np.eye(n) * 0.1  # 确保正定

# 随机生成 s_k 和 y_k (模拟一次迭代后的变化)
s_k = np.random.randn(n)
y_k = np.random.randn(n)
# 为确保 y_k^T s_k > 0 (模拟Wolfe条件)，我们进行简单处理
if np.dot(y_k, s_k) <= 0:
    y_k = y_k + 1.5 * s_k / np.linalg.norm(s_k)  # 调整y_k方向

print("初始 H_k:")
print(H_k)
print(f"\ns_k = {s_k}")
print(f"y_k = {y_k}")
print(f"检查 y_k^T s_k = {np.dot(y_k, s_k):.6f} (应 > 0)")

# 执行BFGS更新
H_next = bfgs_update(H_k, s_k, y_k)

print("\n更新后的 H_{k+1}:")
print(H_next)

# 验证拟牛顿条件: H_{k+1} y_k 应等于 s_k
lhs = H_next @ y_k
print(f"\n验证拟牛顿条件:")
print(f"H_{{k+1}} y_k = {lhs}")
print(f"s_k         = {s_k}")
print(f"差值范数 ||H_{{k+1}} y_k - s_k|| = {np.linalg.norm(lhs - s_k):.10e}")
# 由于浮点数计算，差值应非常接近0
assert np.allclose(lhs, s_k), "BFGS更新不满足拟牛顿条件！"
print("验证通过：BFGS更新满足 H_{k+1} y_k = s_k。")

# 验证对称性
print(f"\n验证对称性: ||H_{{k+1}} - H_{{k+1}}^T|| = {np.linalg.norm(H_next - H_next.T):.10e}")
assert np.allclose(H_next, H_next.T), "更新后的矩阵不对称！"

BFGS 方法的收敛性

了解一个优化算法的收敛性质至关重要。对于 BFGS 方法，其收敛性理论如下：

1. 全局收敛性：对于一致凸的函数（即 Hessian 矩阵的特征值有正的上、下界），如果采用回溯法（Backtracking）线搜索，BFGS 方法具有全局收敛性，即从任意初始点出发都能收敛到全局极小点。在实际中，使用 Wolfe 条件的线搜索通常也能保证良好的全局收敛行为。

2. 局部收敛性与收敛速度：如果目标函数 $f$ 在最优解 $x^$x∗ 的某个邻域内二阶连续可微，且其 Hessian 矩阵 $\nabla^2 f(x^$) 是正定的，并且二阶导数满足 Lipschitz 连续条件，那么当初始点 $x_0$ 足够接近 $x^$x∗ 时，BFGS 方法产生的序列 ${x_k}$ 将超线性收敛到 $x^$。

   收敛速度比较：

   • 牛顿法：在类似条件下，具有 二阶收敛速度（Quadratic Convergence）。这意味着误差大致按平方速率减少：$|x_{k+1} - x^$| \approx C |x_k - x^|^2。
   • BFGS 拟牛顿法：具有 超线性收敛速度（Superlinear Convergence）。这意味着 $\lim_{k \to \infty} \frac{|x_{k+1} - x^$|}{|x_k - x^|} = 0。虽然超线性收敛比二阶收敛慢，但它仍然远快于一阶方法（如梯度下降法）的线性收敛。

   尽管 BFGS 的理论收敛速度慢于牛顿法，但其每步迭代的计算成本（$O(n^2)$）远低于牛顿法（$O(n^3)$）。因此，对于中高维度的优化问题，BFGS 的总计算时间往往远少于牛顿法，是实践中更受欢迎的选择。

📝 动手练一练

1. BFGS 更新计算：考虑目标函数 $f(x) = x_1^2 + 2x_2^2$。假设在迭代中，$x_0 = (1, 1)^T$, $x_1 = (4/9, -1/9)^T$。初始逆 Hessian 近似为 $H_0 = I$（2x2单位矩阵）。请计算对应的 $s_0$, $y_0$，并手动推导（或编写代码验证）BFGS 更新后的矩阵 $H_1$。

2. 拟牛顿条件验证：对于上题计算得到的 $H_1$，请验证其是否满足拟牛顿条件 $H_1 y_0 = s_0$。

参考答案：

1. 首先计算梯度和差分向量：

   • $\nabla f(x) = (2x_1, 4x_2)^T$。
   • $s_0 = x_1 - x_0 = (4/9 - 1, -1/9 - 1)^T = (-5/9, -10/9)^T$。
   • $\nabla f(x_0) = (2, 4)^T$, $\nabla f(x_1) = (8/9, -4/9)^T$。
   • $y_0 = \nabla f(x_1) - \nabla f(x_0) = (8/9-2, -4/9-4)^T = (-10/9, -40/9)^T$。
   • $\rho_0 = 1 / (y_0^T s_0) = 1 / ((-10/9)$(-5/9) + (-40/9)(-10/9)) = 1 / (50/81 + 400/81) = 81/450 = 9/50。

   代入 BFGS 公式 $H_1 = (I - \rho_0 s_0 y_0^T) H_0 (I - \rho_0 y_0 s_0^T) + \rho_0 s_0 s_0^T$，其中 $H_0 = I$。 计算过程略，最终结果为： $$H_1 = \frac{1}{162}\begin{bmatrix} 169 & -22 \ -22 & 46 \end{bmatrix} \approx \begin{bmatrix} 1.0432 & -0.1358 \ -0.1358 & 0.2840 \end{bmatrix}.$$

2. 验证拟牛顿条件：

   • 计算 $H_1 y_0 = \frac{1}{162}\begin{bmatrix} 169 & -22 \ -22 & 46 \end{bmatrix} \begin{bmatrix} -10/9 \ -40/9 \end{bmatrix} = \frac{1}{162} \begin{bmatrix} -1690/9 + 880/9 \ 220/9 - 1840/9 \end{bmatrix} = \frac{1}{162} \begin{bmatrix} -810/9 \ -1620/9 \end{bmatrix} = \begin{bmatrix} -5/9 \ -10/9 \end{bmatrix}$。
   • 这正是 $s_0 = (-5/9, -10/9)^T$。因此，$H_1 y_0 = s_0$ 严格成立，验证了 BFGS 更新公式天然满足拟牛顿条件。

本章小结

本节我们深入探讨了拟牛顿法，特别是其中应用最广泛的 BFGS 方法。我们了解到拟牛顿法通过构造并迭代更新一个正定矩阵来近似 Hessian 矩阵的逆，巧妙地规避了牛顿法在计算复杂度和矩阵正定性方面的两大核心难题。

要点回顾：

• 动机：解决牛顿法 Hessian 矩阵求逆计算成本高 ($O(n^3)$) 和可能非正定的问题。
• 核心：用矩阵 $H_k \approx [\nabla^2 f(x_k)]^{-1}$ 估计逆 Hessian，搜索方向为 $p_k = -H_k \nabla f(x_k)$。
• 拟牛顿条件：更新矩阵必须满足 $H_{k+1} y_k = s_k$，这是连接一阶梯度信息和二阶曲率信息的桥梁。
• BFGS 更新：通过求解一个矩阵最小二乘问题，得到秩-2校正公式，能以 $O(n^2)$ 成本更新 $H_k$ 并保持其对称正定性。
• 算法框架：结合 Wolfe 条件线搜索，形成稳定高效的 BFGS 算法流程。
• 收敛性：对凸函数具有全局收敛性；在最优解附近具有超线性收敛速度，虽慢于牛顿法的二阶收敛，但综合计算效率更高。

行动清单：

1. 推导一遍：亲手推导一次从拟牛顿条件到 BFGS 更新公式的核心步骤，加深对“最小改动”更新思想的理解。
2. 代码实现：根据本节给出的算法框架，尝试用 Python 实现一个完整的 BFGS 求解器，并用于优化一个简单的凸函数（如 Rosenbrock 函数），观察其收敛过程。
3. 对比实验：将你实现的 BFGS 算法与梯度下降法、共轭梯度法在相同问题上进行对比，从迭代次数和计算时间两个维度直观感受不同优化算法的效率差异。

— 小象教研组