拟牛顿法

拟牛顿法 (Quasi-Newton Methods) 详解

1. 概念的直观理解 (The “Why” and “What”)

1.1 这是什么？我们为什么需要它？

在牛顿法优化这一篇中，我们学习了牛顿法。牛顿法之所以收敛快，是因为它使用了函数的二阶导数信息（Hessian 矩阵 $\mathbf{H}(\mathbf{x})$），能够更准确地捕捉函数曲率，从而预测到最优解的方向。

但是，牛顿法有两个主要缺点，限制了它在大型机器学习问题中的应用：

1. 计算Hessian矩阵的计算量巨大：如果模型有 $N$ 个参数，Hessian 矩阵的大小是 $N\times N$。计算这 $N^2$ 个二阶偏导数在计算上非常昂贵，尤其当 $N$ 达到数千、数万甚至数百万时。
2. 计算Hessian逆的计算量巨大：即使计算出Hessian矩阵，求其逆矩阵的计算复杂度是 $O(N^3)$，这在 $N$ 很大时是不可接受的。

拟牛顿法正是为了解决这些问题而诞生的。它的核心思想是：我们不直接计算和存储精确的Hessian矩阵 $\mathbf{H}(\mathbf{x})$ 及其逆 $\mathbf{H}(\mathbf{x}{)}^{-1}$，而是用一个近似矩阵来替代它们。这个近似矩阵会随着迭代逐步更新，以更好地模拟真实的Hessian（或其逆）。

1.2 直观类比：盲人摸象的升级版

还用我们“下山找最低点”的例子：

• 梯度下降法：你闭着眼睛，沿着脚下最陡峭的方向走一步。

• 牛顿法：你拥有超能力，能瞬间扫描整个地形，精确知道每个位置的坡度和弯曲度，然后直接跳向最低点。

• 拟牛顿法：你没有牛顿法那样瞬间扫描地形的超能力，但你比梯度下降法聪明。你虽然不能精确知道当前位置的曲率，但你可以根据你走过的路和感受到的坡度变化来“猜”地形的弯曲程度。

  • 当你走一步，记录下你的起始位置、结束位置，以及在这两点感受到的坡度（梯度）。
  • 根据这些有限的信息，你更新你对地形弯曲度的“猜测”。这个“猜测”不需要完全精确，但它会越来越好。
  • 然后，你利用这个“猜测”出的弯曲度来指引你下一步如何走，以期更快地到达最低点。

这个“猜测”Hessian（或其逆）的过程，就像是一个聪明的盲人摸象。每一次摸索（迭代）都会收集新的信息，然后根据这些信息更新他对大象形状（Hessian）的“整体印象”，而不是每次都去完整地摸一遍大象（计算完整的Hessian）。

核心思想：通过梯度变化来近似Hessian 我们知道，Hessian矩阵是梯度变化的导数。如果我们在点 ${\mathbf{x}}_k$ 和点 ${\mathbf{x}}_{k+1}$ 处计算梯度 $\nabla f({\mathbf{x}}_k)$ 和 $\nabla f({\mathbf{x}}_{k+1})$，那么梯度的变化量 ${\mathbf{y}}_k=\nabla f({\mathbf{x}}_{k+1})-\nabla f({\mathbf{x}}_k)$ 应该与我们移动的步长 ${\mathbf{s}}_k={\mathbf{x}}_{k+1}-{\mathbf{x}}_k$ 之间存在某种关系。 这个关系就是：$\mathbf{H}(\mathbf{x}){\mathbf{s}}_k\approx {\mathbf{y}}_k$。这个被称为 “割线方程 (Secant Equation)”或“拟牛顿条件 (Quasi-Newton Condition)”。 拟牛顿法就是利用这个割线方程来构建和更新Hessian的近似矩阵 ${\mathbf{B}}_k$（或其逆 ${\mathbf{H}}_k$）。

---

2. 预备数学知识与符号解释 (Mathematical Preliminaries & Notation)

除了上次提到的梯度、Hessian和泰勒展开，我们还需要了解：

2.1 牛顿法回顾

牛顿法的迭代公式是： ${\mathbf{x}}_{k+1}={\mathbf{x}}_k-\mathbf{H}({\mathbf{x}}_k{)}^{-1}\nabla f({\mathbf{x}}_k)$ 或者，我们定义搜索方向 ${\mathbf{p}}_k=-\mathbf{H}({\mathbf{x}}_k{)}^{-1}\nabla f({\mathbf{x}}_k)$，那么 ${\mathbf{x}}_{k+1}={\mathbf{x}}_k+{\mathbf{p}}_k$。 拟牛顿法就是用一个近似的逆Hessian矩阵 ${\mathbf{B}}_k^{-1}$（或者直接近似Hessian的逆矩阵 ${\mathbf{H}}_k$）来代替 $\mathbf{H}({\mathbf{x}}_k{)}^{-1}$。 所以拟牛顿法的迭代公式形式为： ${\mathbf{x}}_{k+1}={\mathbf{x}}_k-{\mathbf{H}}_k\nabla f({\mathbf{x}}_k)$ 其中 ${\mathbf{H}}_k$ 是对 $\mathbf{H}({\mathbf{x}}_k{)}^{-1}$ 的近似。

2.2 割线方程 (Secant Equation) / 拟牛顿条件

这是拟牛顿法的核心。 我们知道，对于一个足够光滑的函数 $f(\mathbf{x})$，在 ${\mathbf{x}}_k$ 附近，我们可以将梯度近似为： $\nabla f(\mathbf{x})\approx \nabla f({\mathbf{x}}_k)+\mathbf{H}({\mathbf{x}}_k)(\mathbf{x}-{\mathbf{x}}_k)$ 将 $\mathbf{x}={\mathbf{x}}_{k+1}$ 代入，我们得到： $\nabla f({\mathbf{x}}_{k+1})\approx \nabla f({\mathbf{x}}_k)+\mathbf{H}({\mathbf{x}}_k)({\mathbf{x}}_{k+1}-{\mathbf{x}}_k)$ 重新排列，定义步长向量 ${\mathbf{s}}_k={\mathbf{x}}_{k+1}-{\mathbf{x}}_k$ 和梯度变化向量 ${\mathbf{y}}_k=\nabla f({\mathbf{x}}_{k+1})-\nabla f({\mathbf{x}}_k)$： ${\mathbf{y}}_k\approx \mathbf{H}({\mathbf{x}}_k){\mathbf{s}}_k$ 这就是割线方程。拟牛顿法的目标就是找到一个近似Hessian的矩阵 ${\mathbf{B}}_{k+1}$（或近似逆Hessian的矩阵 ${\mathbf{H}}_{k+1}$），使得它满足这个方程： ${\mathbf{B}}_{k+1}{\mathbf{s}}_k={\mathbf{y}}_k$ 如果我们要近似的是逆Hessian矩阵 ${\mathbf{H}}_{k+1}$，那么对应的割线方程是： ${\mathbf{s}}_k={\mathbf{H}}_{k+1}{\mathbf{y}}_k$

2.3 低秩更新 (Low-Rank Update)

拟牛顿法在每次迭代中更新Hessian近似矩阵时，通常采用低秩更新的方式。这意味着新的近似矩阵是旧的近似矩阵加上一个（或两个）简单的矩阵（秩为1或2的矩阵）。这种更新方式计算效率高，能够保持矩阵的对称性和正定性。

例如，一个秩1更新的形式可能是 $M_{new}=M_{old}+{\mathbf{uv}}^T$，其中 $\mathbf{u},\mathbf{v}$ 是向量。一个秩2更新的形式可能是 $M_{new}=M_{old}+\mathbf{u}{\mathbf{v}}^T+\mathbf{z}{\mathbf{w}}^T$。

---

3. 更新公式

拟牛顿法中有很多不同的算法，它们的主要区别在于如何构造和更新Hessian（或其逆）的近似矩阵。最著名的包括 DFP、BFGS 和 L-BFGS。

我们不会像推导牛顿法那样从零开始推导这些更新公式，因为它们涉及复杂的矩阵代数和优化理论（例如，寻找满足割线方程且与前一次近似矩阵“最接近”的矩阵），这超出了本篇文章的范围。但我们将重点解释这些更新公式的目的、组成部分，以及它们如何满足拟牛顿条件。

所有拟牛顿法的共同步骤：

1. 选择一个初始点 ${\mathbf{x}}_0$ 和一个初始的Hessian近似矩阵 ${\mathbf{B}}_0$（通常选择单位矩阵 $\mathbf{I}$）。
2. 循环迭代 ($k=0,1,2,\dots$)： a. 计算当前点的梯度 $\nabla f({\mathbf{x}}_k)$。 b. 计算搜索方向 ${\mathbf{p}}_k=-{\mathbf{B}}_k^{-1}\nabla f({\mathbf{x}}_k)$ （或者直接用近似的逆Hessian ${\mathbf{H}}_k$ 来计算 ${\mathbf{p}}_k=-{\mathbf{H}}_k\nabla f({\mathbf{x}}_k)$）。 c. 进行线搜索 (Line Search)：确定一个合适的步长 ${\alpha }_k>0$，使得 ${\mathbf{x}}_{k+1}={\mathbf{x}}_k+{\alpha }_k{\mathbf{p}}_k$ 能够显著降低函数值。线搜索是为了确保算法的稳定收敛性，并找到一个“足够好”的步长。 d. 计算新的点 ${\mathbf{x}}_{k+1}$。 e. 计算新的梯度 $\nabla f({\mathbf{x}}_{k+1})$。 f. 计算步长向量 ${\mathbf{s}}_k={\mathbf{x}}_{k+1}-{\mathbf{x}}_k$。 g. 计算梯度变化向量 ${\mathbf{y}}_k=\nabla f({\mathbf{x}}_{k+1})-\nabla f({\mathbf{x}}_k)$。 h. 更新Hessian近似矩阵 ${\mathbf{B}}_{k+1}$（或 ${\mathbf{H}}_{k+1}$），使其满足拟牛顿条件 ${\mathbf{B}}_{k+1}{\mathbf{s}}_k={\mathbf{y}}_k$（或 ${\mathbf{s}}_k={\mathbf{H}}_{k+1}{\mathbf{y}}_k$），并尽可能保持对称性和正定性。这是不同拟牛顿方法的核心区别。

3.1 DFP (Davidon–Fletcher–Powell) 更新公式

DFP 方法是最早提出的有效拟牛顿方法之一，它直接更新逆Hessian近似矩阵 ${\mathbf{H}}_k$。 其更新公式为： ${\mathbf{H}}_{k+1}={\mathbf{H}}_k+\frac{{\mathbf{s}}_k{\mathbf{s}}_k^T}{{\mathbf{s}}_k^T{\mathbf{y}}_k}-\frac{{\mathbf{H}}_k{\mathbf{y}}_k{\mathbf{y}}_k^T{\mathbf{H}}_k}{{\mathbf{y}}_k^T{\mathbf{H}}_k{\mathbf{y}}_k}$ 符号解释：

• ${\mathbf{H}}_{k+1}$: 下一次迭代的逆Hessian近似矩阵。
• ${\mathbf{H}}_k$: 当前迭代的逆Hessian近似矩阵。
• ${\mathbf{s}}_k={\mathbf{x}}_{k+1}-{\mathbf{x}}_k$: 从当前点到下一点的步长向量。
• ${\mathbf{y}}_k=\nabla f({\mathbf{x}}_{k+1})-\nabla f({\mathbf{x}}_k)$: 从当前点到下一点的梯度变化向量。
• ${\mathbf{s}}_k{\mathbf{s}}_k^T$: 向量的外积，结果是一个矩阵（秩为1）。
• ${\mathbf{s}}_k^T{\mathbf{y}}_k$: 向量的点积，结果是一个标量。
• ${\mathbf{H}}_k{\mathbf{y}}_k{\mathbf{y}}_k^T{\mathbf{H}}_k$: 矩阵与向量的乘积及外积的组合，结果是一个矩阵（秩为1）。

公式目的： 这个公式的推导是为了找到一个最小化某种范数（通常是Frobenius范数）的更新，同时强制满足拟牛顿条件 ${\mathbf{s}}_k={\mathbf{H}}_{k+1}{\mathbf{y}}_k$，并保持矩阵的对称性。

3.2 BFGS (Broyden–Fletcher–Goldfarb–Shanno) 更新公式

BFGS 是目前最常用和最成功的拟牛顿方法。它更新Hessian近似矩阵 ${\mathbf{B}}_k$，而不是其逆。 其更新公式为： ${\mathbf{B}}_{k+1}={\mathbf{B}}_k+\frac{{\mathbf{y}}_k{\mathbf{y}}_k^T}{{\mathbf{y}}_k^T{\mathbf{s}}_k}-\frac{{\mathbf{B}}_k{\mathbf{s}}_k{\mathbf{s}}_k^T{\mathbf{B}}_k}{{\mathbf{s}}_k^T{\mathbf{B}}_k{\mathbf{s}}_k}$ 符号解释：

• ${\mathbf{B}}_{k+1}$: 下一次迭代的Hessian近似矩阵。
• ${\mathbf{B}}_k$: 当前迭代的Hessian近似矩阵。
• ${\mathbf{s}}_k,{\mathbf{y}}_k$: 同 DFP 定义。

公式目的： BFGS 可以看作是 DFP 的“对偶”形式。它同样通过一个秩2更新来满足拟牛顿条件 ${\mathbf{B}}_{k+1}{\mathbf{s}}_k={\mathbf{y}}_k$，并保持对称性和（如果 ${\mathbf{B}}_k$ 足够正定）正定性。BFGS 在实践中通常比 DFP 更鲁棒。

在实际应用中，我们通常不需要显式计算 ${\mathbf{B}}_k^{-1}$。 虽然 BFGS 更新的是 ${\mathbf{B}}_k$，但可以通过 Shermon-Morrison 公式来高效地计算 ${\mathbf{B}}_{k+1}^{-1}$，而无需每次迭代都进行矩阵求逆操作。这使得 BFGS 能够直接用于计算搜索方向 ${\mathbf{p}}_k=-{\mathbf{B}}_{k+1}^{-1}\nabla f({\mathbf{x}}_k)$。实际上，BFGS 算法通常是直接更新逆 Hessian 近似矩阵的公式，即：

${\mathbf{H}}_{k+1}=\left(\mathbf{I}-\frac{{\mathbf{s}}_k{\mathbf{y}}_k^T}{{\mathbf{y}}_k^T{\mathbf{s}}_k}\right){\mathbf{H}}_k\left(\mathbf{I}-\frac{{\mathbf{y}}_k{\mathbf{s}}_k^T}{{\mathbf{y}}_k^T{\mathbf{s}}_k}\right)+\frac{{\mathbf{s}}_k{\mathbf{s}}_k^T}{{\mathbf{y}}_k^T{\mathbf{s}}_k}$ 这个公式看起来更复杂，但它直接给出了更新后的逆Hessian近似矩阵 ${\mathbf{H}}_{k+1}$，避免了每次迭代求逆的昂贵操作。这是最常用的 BFGS 实现方式。

3.3 L-BFGS (Limited-memory BFGS)

问题： 即使 BFGS 避免了直接计算 Hessian，它仍然需要存储整个 $N\times N$ 的近似矩阵 ${\mathbf{B}}_k$ 或 ${\mathbf{H}}_k$。当 $N$ 非常大时（例如在深度学习中 $N$ 达到数百万），这个矩阵会占用巨大的内存（$N^2$ 空间），变得不可行。

L-BFGS 解决方案： L-BFGS（有限内存 BFGS）不存储完整的 $N\times N$ 矩阵。相反，它只存储最近 $m$ 次迭代的 ${\mathbf{s}}_k$ 和 ${\mathbf{y}}_k$ 向量对。

• 当需要计算搜索方向 ${\mathbf{p}}_k=-{\mathbf{H}}_k\nabla f({\mathbf{x}}_k)$ 时，L-BFGS 使用一种称为“两环递归 (two-loop recursion)”的算法，通过存储的这些 ${\mathbf{s}}_i,{\mathbf{y}}_i$ 对，来间接地计算矩阵向量乘积 ${\mathbf{H}}_k\nabla f({\mathbf{x}}_k)$。
• 它不需要显式构建和存储 ${\mathbf{H}}_k$ 矩阵。
• 通常 $m$ 的值是一个小的常数（例如 5 到 20）。这意味着 L-BFGS 的内存需求与 $N$ 成线性关系（$O(mN)$），而不是 $N^2$，这使得它能够应用于大规模优化问题。

L-BFGS 的优势：

• 内存效率高：解决了 BFGS 在高维问题上的内存瓶颈。
• 计算效率高：避免了矩阵乘法和求逆。
• 收敛性好：虽然比全量 BFGS 稍慢，但通常比梯度下降法快得多。

---

4. 数值化实例 (Concrete Numerical Example)

我们来演示一个非常简化的一步 BFGS 更新。由于完整的矩阵更新和线搜索过程复杂，我们只关注如何从 ${\mathbf{s}}_k$ 和 ${\mathbf{y}}_k$ 来更新近似矩阵。

假设我们仍在优化上一份报告中的多变量函数： $f(x_1,x_2)=x_1^2+2x_2^2-4x_1-8x_2+10$

我们已知：

• 真值梯度：$\nabla f(\mathbf{x})=\left[\begin{array}{c}2x_1-4\\ 4x_2-8\end{array}\right]$
• 真值Hessian：$\mathbf{H}(\mathbf{x})=\left[\begin{array}{cc}2& 0\\ 0& 4\end{array}\right]$

假设当前迭代 $k$：

• 当前点 ${\mathbf{x}}_k=\left[\begin{array}{c}0\\ 0\end{array}\right]$
• 初始Hessian近似矩阵 ${\mathbf{B}}_k=\mathbf{I}=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]$ (在实践中，通常从单位矩阵开始)
• 计算当前梯度 $\nabla f({\mathbf{x}}_k)=\nabla f(\left[\begin{array}{c}0\\ 0\end{array}\right])=\left[\begin{array}{c}-4\\ -8\end{array}\right]$

步骤 1：计算搜索方向 ${\mathbf{p}}_k$ ${\mathbf{p}}_k=-{\mathbf{B}}_k^{-1}\nabla f({\mathbf{x}}_k)$ 因为 ${\mathbf{B}}_k=\mathbf{I}$，所以 ${\mathbf{B}}_k^{-1}=\mathbf{I}$。 ${\mathbf{p}}_k=-\mathbf{I}\left[\begin{array}{c}-4\\ -8\end{array}\right]=\left[\begin{array}{c}4\\ 8\end{array}\right]$

步骤 2：进行线搜索（假设我们得到步长 ${\alpha }_k=0.5$） ${\mathbf{x}}_{k+1}={\mathbf{x}}_k+{\alpha }_k{\mathbf{p}}_k=\left[\begin{array}{c}0\\ 0\end{array}\right]+0.5\left[\begin{array}{c}4\\ 8\end{array}\right]=\left[\begin{array}{c}2\\ 4\end{array}\right]$

步骤 3：计算 ${\mathbf{s}}_k$ 和 ${\mathbf{y}}_k$

• 步长向量 ${\mathbf{s}}_k={\mathbf{x}}_{k+1}-{\mathbf{x}}_k$ ${\mathbf{s}}_k=\left[\begin{array}{c}2\\ 4\end{array}\right]-\left[\begin{array}{c}0\\ 0\end{array}\right]=\left[\begin{array}{c}2\\ 4\end{array}\right]$
• 梯度变化向量 ${\mathbf{y}}_k=\nabla f({\mathbf{x}}_{k+1})-\nabla f({\mathbf{x}}_k)$ 首先计算 $\nabla f({\mathbf{x}}_{k+1})$： $\nabla f(\left[\begin{array}{c}2\\ 4\end{array}\right])=\left[\begin{array}{c}2(2)-4\\ 4(4)-8\end{array}\right]=\left[\begin{array}{c}4-4\\ 16-8\end{array}\right]=\left[\begin{array}{c}0\\ 8\end{array}\right]$ 现在计算 ${\mathbf{y}}_k$: ${\mathbf{y}}_k=\left[\begin{array}{c}0\\ 8\end{array}\right]-\left[\begin{array}{c}-4\\ -8\end{array}\right]=\left[\begin{array}{c}0-(-4)\\ 8-(-8)\end{array}\right]=\left[\begin{array}{c}4\\ 16\end{array}\right]$

步骤 4：使用 BFGS 公式更新 ${\mathbf{B}}_{k+1}$ 回顾 BFGS 更新公式： ${\mathbf{B}}_{k+1}={\mathbf{B}}_k+\frac{{\mathbf{y}}_k{\mathbf{y}}_k^T}{{\mathbf{y}}_k^T{\mathbf{s}}_k}-\frac{{\mathbf{B}}_k{\mathbf{s}}_k{\mathbf{s}}_k^T{\mathbf{B}}_k}{{\mathbf{s}}_k^T{\mathbf{B}}_k{\mathbf{s}}_k}$

计算分母项：

• ${\mathbf{y}}_k^T{\mathbf{s}}_k$: ${\mathbf{y}}_k^T{\mathbf{s}}_k=\left[\begin{array}{cc}4& 16\end{array}\right]\left[\begin{array}{c}2\\ 4\end{array}\right]=(4)(2)+(16)(4)=8+64=72$
• ${\mathbf{s}}_k^T{\mathbf{B}}_k{\mathbf{s}}_k$: ${\mathbf{s}}_k^T{\mathbf{B}}_k{\mathbf{s}}_k=\left[\begin{array}{cc}2& 4\end{array}\right]\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]\left[\begin{array}{c}2\\ 4\end{array}\right]$ $=\left[\begin{array}{cc}2& 4\end{array}\right]\left[\begin{array}{c}(1)(2)+(0)(4)\\ (0)(2)+(1)(4)\end{array}\right]=\left[\begin{array}{cc}2& 4\end{array}\right]\left[\begin{array}{c}2\\ 4\end{array}\right]=(2)(2)+(4)(4)=4+16=20$

计算分子项和矩阵项：

• $\frac{{\mathbf{y}}_k{\mathbf{y}}_k^T}{{\mathbf{y}}_k^T{\mathbf{s}}_k}$: ${\mathbf{y}}_k{\mathbf{y}}_k^T=\left[\begin{array}{c}4\\ 16\end{array}\right]\left[\begin{array}{cc}4& 16\end{array}\right]=\left[\begin{array}{cc}(4)(4)& (4)(16)\\ (16)(4)& (16)(16)\end{array}\right]=\left[\begin{array}{cc}16& 64\\ 64& 256\end{array}\right]$ $\frac{{\mathbf{y}}_k{\mathbf{y}}_k^T}{{\mathbf{y}}_k^T{\mathbf{s}}_k}=\frac{1}{72}\left[\begin{array}{cc}16& 64\\ 64& 256\end{array}\right]=\left[\begin{array}{cc}16/72& 64/72\\ 64/72& 256/72\end{array}\right]=\left[\begin{array}{cc}2/9& 8/9\\ 8/9& 32/9\end{array}\right]\approx \left[\begin{array}{cc}0.222& 0.889\\ 0.889& 3.556\end{array}\right]$
• $\frac{{\mathbf{B}}_k{\mathbf{s}}_k{\mathbf{s}}_k^T{\mathbf{B}}_k}{{\mathbf{s}}_k^T{\mathbf{B}}_k{\mathbf{s}}_k}$: 首先计算 ${\mathbf{B}}_k{\mathbf{s}}_k$: ${\mathbf{B}}_k{\mathbf{s}}_k=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]\left[\begin{array}{c}2\\ 4\end{array}\right]=\left[\begin{array}{c}2\\ 4\end{array}\right]$ 然后计算 $({\mathbf{B}}_k{\mathbf{s}}_k)({\mathbf{B}}_k{\mathbf{s}}_k{)}^T$: $\left[\begin{array}{c}2\\ 4\end{array}\right]\left[\begin{array}{cc}2& 4\end{array}\right]=\left[\begin{array}{cc}(2)(2)& (2)(4)\\ (4)(2)& (4)(4)\end{array}\right]=\left[\begin{array}{cc}4& 8\\ 8& 16\end{array}\right]$ $\frac{{\mathbf{B}}_k{\mathbf{s}}_k{\mathbf{s}}_k^T{\mathbf{B}}_k}{{\mathbf{s}}_k^T{\mathbf{B}}_k{\mathbf{s}}_k}=\frac{1}{20}\left[\begin{array}{cc}4& 8\\ 8& 16\end{array}\right]=\left[\begin{array}{cc}4/20& 8/20\\ 8/20& 16/20\end{array}\right]=\left[\begin{array}{cc}1/5& 2/5\\ 2/5& 4/5\end{array}\right]=\left[\begin{array}{cc}0.2& 0.4\\ 0.4& 0.8\end{array}\right]$

将所有项代入 BFGS 公式： ${\mathbf{B}}_{k+1}={\mathbf{B}}_k+\frac{{\mathbf{y}}_k{\mathbf{y}}_k^T}{{\mathbf{y}}_k^T{\mathbf{s}}_k}-\frac{{\mathbf{B}}_k{\mathbf{s}}_k{\mathbf{s}}_k^T{\mathbf{B}}_k}{{\mathbf{s}}_k^T{\mathbf{B}}_k{\mathbf{s}}_k}$ ${\mathbf{B}}_{k+1}=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]+\left[\begin{array}{cc}0.222& 0.889\\ 0.889& 3.556\end{array}\right]-\left[\begin{array}{cc}0.2& 0.4\\ 0.4& 0.8\end{array}\right]$ ${\mathbf{B}}_{k+1}=\left[\begin{array}{cc}1+0.222-0.2& 0+0.889-0.4\\ 0+0.889-0.4& 1+3.556-0.8\end{array}\right]$ ${\mathbf{B}}_{k+1}=\left[\begin{array}{cc}1.022& 0.489\\ 0.489& 3.756\end{array}\right]$

验证： 我们知道真实的Hessian是 $\left[\begin{array}{cc}2& 0\\ 0& 4\end{array}\right]$。可以看到，经过一次迭代，我们从单位矩阵 $\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]$ 得到了 $\left[\begin{array}{cc}1.022& 0.489\\ 0.489& 3.756\end{array}\right]$，它已经开始向真实Hessian靠近。随着迭代次数的增加，${\mathbf{B}}_k$ 会越来越精确地近似真正的Hessian。

---

5. 代码实现 (Code Implementation)

在实践中，我们很少会自己从头实现 BFGS 或 L-BFGS，因为它们涉及复杂的数值稳定性考虑和线搜索算法。Python 中的 scipy.optimize 库提供了非常成熟和高效的实现。

我们将使用 scipy.optimize.minimize 函数，它允许我们指定使用不同的优化方法，包括 ‘BFGS’ 和 ‘L-BFGS-B’。

import numpy as np
from scipy.optimize import minimize
 
# --- 1. 定义目标函数和其梯度 ---
 
def f_multivariate(x_vec):
    """
    多变量目标函数：f(x1, x2) = x1^2 + 2x2^2 - 4x1 - 8x2 + 10
    参数:
        x_vec (np.array): 包含 x1, x2 的数组
    返回:
        float: 函数值
    """
    x1, x2 = x_vec[0], x_vec[1]
    return x1**2 + 2*x2**2 - 4*x1 - 8*x2 + 10
 
def grad_multivariate(x_vec):
    """
    多变量目标函数的梯度向量：nabla f(x) = [2x1 - 4, 4x2 - 8]^T
    参数:
        x_vec (np.array): 包含 x1, x2 的数组
    返回:
        np.array: 梯度向量
    """
    x1, x2 = x_vec[0], x_vec[1]
    return np.array([2*x1 - 4, 4*x2 - 8])
 
# --- 2. 设置初始点 ---
initial_x_vec = np.array([0.0, 0.0])
 
# --- 3. 使用 BFGS 方法进行优化 ---
print("--- 使用 BFGS 方法优化 ---")
# `minimize` 函数会自动处理 Hessian 矩阵的近似和更新
# 这里我们只需要提供目标函数 `fun` 和其梯度 `jac` (Jacobian，对于标量函数就是梯度)
result_bfgs = minimize(fun=f_multivariate, x0=initial_x_vec, jac=grad_multivariate, method='BFGS')
 
print("\nBFGS 优化结果:")
print(f"收敛状态: {result_bfgs.success}")
print(f"收敛消息: {result_bfgs.message}")
print(f"找到的最优参数 (x1, x2): {result_bfgs.x}")
print(f"在该点的函数值: {result_bfgs.fun}")
print(f"迭代次数: {result_bfgs.nit}")
 
# --- 4. 使用 L-BFGS-B 方法进行优化 ---
# L-BFGS-B 是 L-BFGS 的一个变种，支持参数边界约束 (Bounds)，
# 但即使没有边界，它也是一个非常常用的 L-BFGS 实现。
print("\n--- 使用 L-BFGS-B 方法优化 ---")
result_lbfgsb = minimize(fun=f_multivariate, x0=initial_x_vec, jac=grad_multivariate, method='L-BFGS-B')
 
print("\nL-BFGS-B 优化结果:")
print(f"收敛状态: {result_lbfgsb.success}")
print(f"收敛消息: {result_lbfgsb.message}")
print(f"找到的最优参数 (x1, x2): {result_lbfgsb.x}")
print(f"在该点的函数值: {result_lbfgsb.fun}")
print(f"迭代次数: {result_lbfgsb.nit}")
 
# --- 对比：牛顿法 (需要提供Hessian) ---
# 为了完整性，这里也展示牛顿法（需要提供Hessian函数）
def hessian_multivariate(x_vec):
    """Hessian矩阵：H(x) = [[2, 0], [0, 4]]"""
    return np.array([[2, 0], [0, 4]])
 
print("\n--- 使用 Newton-CG 方法优化 (需要提供Hessian) ---")
# Newton-CG 是一个利用共轭梯度法计算牛顿方向的牛顿方法变种，
# 适合处理稀疏Hessian或 Hessian-vector 乘积可用情况
result_newton_cg = minimize(fun=f_multivariate, x0=initial_x_vec, jac=grad_multivariate, hess=hessian_multivariate, method='Newton-CG')
 
print("\nNewton-CG 优化结果:")
print(f"收敛状态: {result_newton_cg.success}")
print(f"收敛消息: {result_newton_cg.message}")
print(f"找到的最优参数 (x1, x2): {result_newton_cg.x}")
print(f"在该点的函数值: {result_newton_cg.fun}")
print(f"迭代次数: {result_newton_cg.nit}")

代码与数学公式的关联：

• f_multivariate(x_vec) 对应我们的目标函数 $f(\mathbf{x})$。
• grad_multivariate(x_vec) 对应梯度 $\nabla f(\mathbf{x})$。scipy.optimize.minimize 会内部调用这个函数来获取梯度信息。
• method='BFGS' 和 method='L-BFGS-B' 告诉 scipy 使用相应的拟牛顿算法。它们内部会自动处理 BFGS 或 L-BFGS 的近似矩阵更新逻辑（如上面数值例子中的公式）。你不需要手动编写这些复杂的矩阵更新。
• x0 是初始点 ${\mathbf{x}}_0$。
• result.x 是算法找到的最优参数 ${\mathbf{x}}^{\ast }$。
• result.fun 是在最优参数处的函数值 $f({\mathbf{x}}^{\ast })$。

对于我们这个简单的二次函数例子，所有方法都能在很少的迭代次数内收敛到精确解 $\left[\begin{array}{c}2\\ 2\end{array}\right]$。在更复杂的、非二次的函数中，拟牛顿法通常比梯度下降法快，而比牛顿法慢（但更实用）。

---