[最优化方法笔记] 拟牛顿法 SR1, BFGS, DFP

1. 拟牛顿法

1.1 回顾牛顿法

牛顿法（经典牛顿法）的迭代表达式：

\[x^{k + 1} = x^k - \nabla^2 f(x^k)^{-1} \nabla f(x^k) \]

但是，牛顿法过程中 \(\text{Hessian}\) 矩阵 \(\nabla^2 f(x^k)\) 的计算和存储的代价很高，对于条件数较多的问题很难求解。因此，引入 拟牛顿法。

1.2 拟牛顿法

拟牛顿法 的核心思路在于，在牛顿法的迭代过程中，用 近似解 计算第 \(k\) 次迭代下的 \(\text{Hessian}\) 矩阵 \(\nabla^2 f(x^k)\)，近似值记为 \(B^k\)，即有 \(B^k \approx \nabla^2 f(x^k)\)，称为 拟牛顿矩阵。

用 近似值 \(B^k\) 代替牛顿法中的 \(\nabla^2 f(x^k)\)，得：

\[x^{k + 1} = x^k - (B^k)^{-1} \nabla f(x^k) \]

在近似 \(\text{Hessian}\) 矩阵时，也需要通过 某种映射关系 并 不断迭代 得到。但是依然需要求近似矩阵的逆，为了避免计算逆矩阵的开销，我们可以 直接近似 \(\text{Hessian}\) 矩阵的逆，记 \(H^k = (B^k)^{-1}\)。故我们有：

\[\begin{split} x^{k + 1} &= x^k - H^k \nabla f(x^k) \\\\ H^{k + 1} &= g(H^k) \end{split} \]

其中 \(g\) 为 近似 \(\text{Hessian}\) 矩阵的逆 的映射函数。一般有 \(H^{k + 1} = H^k + C^k\)，其中 \(C^k\) 被称为 修正矩阵。

1.3 拟牛顿法基本过程

拟牛顿法：

• 令 \(H^0 = I\)，任选初始点 \(x^0 \in \mathbb{R}^n\)，令 \(k = 0\)

• 计算 梯度 \(\nabla f(x^k)\)，如果满足终止条件 \(||\nabla f(x^k)|| < \varepsilon\)，取 \(x^* = x^k\)，并结束整个算法

• 计算 搜索方向 \(d^k = - H^k \nabla f(x^k)\)，\(H^k\) 为当前 \(x^k\) 处的 \(\text{Hessian}\) 矩阵的近似

• 迭代更新 \(x\)：\(x^{k + 1} = x^k + d^k\)

• 更新 \(H\)：\(H^{k + 1} = g(H^k)\) 根据 \(x^k\) 点的信息进行简单修正

当然，也可以用线搜索增加一个步长 \(\alpha\)，变成 拟阻尼牛顿法：

• 令 \(H^0 = I\)，任选初始点 \(x^0 \in \mathbb{R}^n\)，令 \(k = 0\)

• 计算 梯度 \(\nabla f(x^k)\)，如果满足终止条件 \(||\nabla f(x^k)|| < \varepsilon\)，取 \(x^* = x^k\)，并结束整个算法

• 计算 搜索方向 \(d^k = - H^k \nabla f(x^k)\)，\(H^k\) 为当前 \(x^k\) 处的 \(\text{Hessian}\) 矩阵的近似

• 计算 步长 \(\alpha\)，通过 线搜索 确定当前第 \(k\) 次迭代的步长（精确搜索、直接搜索、\(\text{Wolfe, Goldstein, Armijo}\) 非精确准则等）

• 迭代更新 \(x\)：\(x^{k + 1} = x^k + \alpha_k d^k\)

• 更新 \(H\)：\(H^{k + 1} = g(H^k)\) 根据 \(x^k\) 点的信息进行简单修正

---

2. 拟牛顿法 \(H^k\) 的确定

2.1 割线方程

设 \(f(x)\) 是二阶连续可微函数，对 \(\nabla f(x)\) 在点 \(x^{k + 1}\) 处进行一阶泰勒近似，得：

\[\nabla f(x) = \nabla f(x^{k + 1}) + \nabla ^2 f(x^{k + 1})(x - x^{k + 1}) + O(||x - x^{k + 1}||^2) \]

令 \(x = x^{k}\)，设 \(s^k = x^{k + 1} - x^k\) 为 点差，\(y^k = \nabla f(x^{k + 1}) - \nabla f(x^k)\) 为 梯度差，得：

\[\nabla^2 f(x^{k + 1})s^k + O(||s^k||^2) = y^k \]

忽略高阶项 \(O(||s^k||^2)\)，由此可以得到：

\[\nabla^2 f(x^{k + 1})s^k = y^k \]

所以，我们希望 近似 \(\text{Hessian}\) 矩阵 \(B^{k + 1}\) 满足方程：

\[B^{k + 1}s^k = y^k \]

因此 近似 \(\text{Hessian}\) 矩阵的逆 \(H^{k + 1}\) 满足：

\[H^{k + 1}y^k = s^k \]

上述的两个方程被称为 割线方程。

2.2 曲率条件

\(H^{k + 1}\) 满足上述 割线方程，且 保证 \(B^{k + 1}\) 对称正定，即满足了 牛顿法 中的必要条件。有时，该方程也被称为 拟牛顿方程。

保证 \(B^{k + 1}\) 对称正定，即满足条件：

\[(s^k)^T B^{k + 1} s^k > 0 \Rightarrow (s^k)^Ty^k > 0 \]

这个条件被称为 曲率条件，是 拟牛顿法迭代过程中的必要条件之一。

---

3. SR1 方法

3.1 SR1 定义

\(\text{SR1}\) 方法 （秩一更新 \(\text{Symmetric Rank-One}\)）的核心思路很简单，即 根据 \(x^k\) 处的信息得到修正量 \(\Delta H^k\) 来更新 \(H^k\)，即：

\[H^{k + 1} = H^k + \Delta H^k \]

我们希望 \(H^k \approx \nabla ^2 f(x^{k})^{-1}\) ，\(H^{k + 1} \approx \nabla ^2 f(x^{k + 1})^{-1}\)，故有：

\[\Delta H^k \approx \nabla ^2 f(x^{k + 1})^{-1} - \nabla ^2 f(x^{k})^{-1} \]

需要保证 \(H^k\) 和 \(H^{k + 1}\) 都是对称的，故显然 \(\Delta H^k\) 也是对称的。所以令 \(\beta \in \mathbb{R}, \; u \in \mathbb{R}^n\)，使得 \(\Delta H^k = \beta u u^T\)，故 \(H\) 的迭代更新表达式为：

\[H^{k + 1} = H^k + \beta u u^T \]

显然 \(\beta u u^T\) 是一个 \(n\times n\) 的 对称矩阵。\(\beta\) 是待定的标量，\(u\) 是待定的向量。

3.2 SR1 更新公式

根据 割线方程 \(H^{k + 1}y^k = s^k\)，代入 \(\text{SR1}\) 更新的结果，得到：

\[(H^k + \beta u u^T)y^k = s^k \]

整理可得：

\[\beta u u^T y^k = (\beta u^T y^k)u = s^k - H^k y^k \]

其中可以得出 \(\beta u^T y^k\) 是一个 标量，因此上式表明 向量 \(u\) 和 \(s^k - H^k y^k\) 同向。故有：

\[u = \frac{1}{\beta u^T y^k} (s^k - H^k y^k) \]

记 \(\frac{1}{\beta u^T y^k} = \gamma\)，得：

\[u = \gamma (s^k - H^k y^k) \]

将 \(u\) 回代到 \(\beta u u^T y^k = s^k - H^k y^k\)，得：

\[s^k - H^ky^k = \beta \gamma^2 (s^k - H^k y^k)(s^k - H^k y^k)^T y^k \]

由于 \(\beta \gamma^2\) 和 \((s^k - H^ky^k)^T y^k\) 都是 标量，上式可以写成：

\[s^k - H^k y^k = [\beta \gamma^2 (s^k - H^k y^k)^T y^k](s^k - H^k y^k) \]

显然只有在 \(\beta \gamma^2 (s^k - H^k y^k)^T y^k = 1\) 时，等式成立。

因此，我们可以得到：

\[\beta \gamma^2 = \frac{1}{(s^k - H^k y^k)^T y^k} \]

将上式 \(\beta \gamma^2\) 回代到 迭代更新表达式 \(H^{k + 1} = H^k + \beta u u^T\):

\[\begin{split} H^{k + 1} &= H^k + \beta u u^T \\\\ &= H^k + \beta \gamma^2 (s^k - H^ky^k) (s^k - H^ky^k)^T \\\\ &= H^k + \frac{(s^k - H^ky^k) (s^k - H^ky^k)^T}{(s^k - H^k y^k)^T y^k} \end{split} \]

记 \(v = s^k - H^ky^k\)，那么上述更新表达式可以化简为：

\[H^{k + 1} = H^k + \frac{v v^T}{v^T y^k} \]

由此得到了最终 \(\text{SR1}\) 方法 的 更新公式。

3.3 SR1 的缺点

• 在迭代过程中 无法保证 \(B^k\) 正定，也就是说 搜索方向不一定下降。而且即使 \(B^k\) 正定，也 不一定保证 \(B^{k + 1}\)

• 无法保证 \(v^Ty^k\) 恒大于 0，因此也可能会导致后续的 \(B^{k + 1}\) 非正定。

由于上述缺点，\(\text{SR1}\) 方法一般很少被使用。

---

4. BFGS 方法

4.1 BFGS 定义

\(\text{BFGS}\) 方法考虑的是 对 \(B^k\) 进行秩二更新。对于拟牛顿矩阵 \(B^k \in \mathbb{R}^{n\times n}\)，设 \(u \not = 0, \; v \not = 0, \; u, v \in \mathbb{R}^n\) 以及 \(a, b \in \mathbb{R}\)，其中设定的向量和标量都是待定的，则有 秩二更新表达式：

\[B^{k + 1} = B^k + auu^T + bvv^T \]

显然 \(auu^T\) 和 \(bvv^T\) 都是对称的。

4.2 BFGS 更新公式

根据 割线方程 \(B^{k + 1}s^k = y^k\)，代入 待定参量，得：

\[B^{k + 1}s^k = (B^k + auu^T + bvv^T)s^k = y^k \]

整理可得：

\[auu^Ts^k + bvv^Ts^k = (au^Ts^k)u + (bv^Ts^k)v = y^k - B^k s^k \]

可以得出 \(au^Ts^k\) 和 \(bv^Ts^k\) 为 标量，不妨取 \((au^Ts^k)u = y^k, \; (bv^Ts^k)v = - B^ks^k\)，所以可以得到如下取值：

\[au^Ts^k = 1, \; u = y^k, \; bv^Ts^k = -1, \; v = B^ks^k \]

化简可得所有 待定参量的取值：

\[\begin{split} a &= \frac{1}{u^Ts^k} = \frac{1}{(y^k)^T s^k} \\\\ b &= - \frac{1}{v^Ts^k} = - \frac{1}{(B^ks^k)^Ts^k} = -\frac{1}{(s^k)^TB^k s^k} \end{split} \]

PS: \(B^k\) 是 对称的大小为 \(n\times n\) 的方阵，所以有 \(B^k = (B^k)^T\)

将上述取值回代到 更新表达式 \(B^{k + 1} = B^k + auu^T + bvv^T\)，得:

\[\begin{split} B^{k + 1} = B^k + \frac{y^k (y^k)^T}{(y^k)^Ts^k} - \frac{B^ks^k(s^k)^TB^k}{(s^k)^TB^ks^k} \end{split} \]

借助 \(\text{SMW}\) 公式 sherman-morrison-woodbury 公式，可以得到 \(H^{k + 1}\) 的 BFGS更新表达式：

\[H^{k + 1} = H^k + [1 + \frac{(y^k)^TH^ky^k}{(s^k)^Ty^k}] \frac{s^k(s^k)^T}{(s^k)^Ty^k} - [\frac{s^k(y^k)^TH^k + H^ky^k(s^k)^T}{(s^k)^Ty^k}] \]

由此得到了最终 \(\text{BFGS}\) 方法 的 更新公式。

使用 \(\text{SMW}\) 公式证明可以看这篇文章 Broyden类算法：BFGS算法的迭代公式推导（应用两次Sherman-Morrison公式）

4.3 BFGS 有效性

\(\text{BFGS}\) 使得拟牛顿矩阵正定的充分条件可以是：

• \(B^k\) 或 \(H^k\) 正定

• 满足 曲率条件 \((s^k)^Ty^k > 0, \; \forall k \in \mathbb{N}^+\)

基于 \(\text{BFGS}\) 方法 的 更新公式，可以得出 \(B^{k + 1}\) 及其逆 \(H^{k + 1}\) 均正定。

对于 拟阻尼牛顿法，采用 \(\text{Wolfe}\) 准则进行线搜索，即可满足 曲率条件。

综上，\(\text{BFGS}\) 方法是可以使 \(B^k\) 保持正定的，是一个有效的算法。

---

5. DFP 方法

5.1 DFP 定义

\(\text{DFP}\) 方法考虑的是 对 \(H^k\) 进行秩二更新。思路和前面都大致一致，（此处省略一些待定参量的声明），有：

\[H^{k + 1} = H^k + \beta uu^T + \gamma vv^T \]

5.2 DFP 更新公式

根据 割线方程 \(H^{k + 1}y^k = s^k\)，代入待定参量，得：

\[(H^k + \beta uu^T + \gamma vv^T)y^k = s^k \]

整理得：

\[\beta uu^T y^k + \gamma vv^T y^k = s^k - H^ky^k \]

可以得出 \(\beta u^T y^k\) 和 \(\gamma v^T y^k\) 都是 标量，所以可以将上式写成：

\[(\beta u^T y^k)u + (\gamma v^T y^k)v = s^k - H^ky^k \]

不妨取 \((\beta u^Ty^k)u = s^k, \; (\gamma v^Ty^k)v = - H^ky^k\)，所以可以得到如下取值：

\[\beta u^Ty^k = 1, \; u = s^k, \; \gamma v^Ty^k = -1, \; v = H^ky^k \]

化简可得所有 待定参量的取值：

\[\begin{split} \beta &= \frac{1}{u^Ty^k} = \frac{1}{(s^k)^T y^k} \\\\ \gamma &= - \frac{1}{v^Ty^k} = - \frac{1}{(H^ky^k)^Ty^k} = -\frac{1}{(y^k)^TH^k y^k} \end{split} \]

PS: \(H^k\) 是 对称的大小为 \(n\times n\) 的方阵，所以有 \(H^k = (H^k)^T\)

将上述取值回代到 更新表达式 \(H^{k + 1} = H^k + \beta uu^T + \gamma vv^T\)，得:

\[\begin{split} H^{k + 1} = H^k + \frac{s^k (s^k)^T}{(s^k)^Ty^k} - \frac{H^ky^k(y^k)^TH^k}{(y^k)^TH^ky^k} \end{split} \]

借助 \(\text{SMW}\) 公式 sherman-morrison-woodbury 公式，同样也可以得到 \(B^{k + 1}\) 的 DFP更新表达式：

\[B^{k + 1} = B^k + [1 + \frac{(s^k)^TB^ks^k}{(y^k)^Ts^k}] \frac{y^k(y^k)^T}{(y^k)^Ts^k} - [\frac{y^k(s^k)^TB^k + B^ks^k(y^k)^T}{(y^k)^Ts^k}] \]

由此得到了最终 \(\text{DFP}\) 方法 的 更新公式。

5.3 DFP 方法的劣势

尽管 \(\text{DFP}\) 格式与 \(\text{BFGS}\) 对偶, 但从实际效果而言, \(\text{DFP}\) 格式的 求解效率 整体 不如 \(\text{BFGS}\) 格式。

---

参考

刘浩洋，户将， 李勇锋，文再文，《最优化：建模、算法与理论》

最优化方法复习笔记（四）拟牛顿法与SR1,DFP,BFGS三种拟牛顿算法的推导与代码实现

sherman-morrison-woodbury 公式