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

1. 拟牛顿法

1.1 回顾牛顿法

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

$$x^{k+1}=x^k-{\mathrm{\nabla }}^{2}f(x^k{)}^{-1}\mathrm{\nabla }f(x^k)$$

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

1.2 拟牛顿法

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

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

$$x^{k+1}=x^k-(B^k{)}^{-1}\mathrm{\nabla }f(x^k)$$

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

$$\begin{array}{rl}{x}^{k+1}& =x^k-H^k\mathrm{\nabla }f(x^k)\\ \\ H^{k+1}& =g(H^k)\end{array}$$

其中 $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$

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

• 计算 搜索方向 $d^k=-H^k\mathrm{\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$

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

• 计算 搜索方向 $d^k=-H^k\mathrm{\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)$ 是二阶连续可微函数，对 $\mathrm{\nabla }f(x)$ 在点 $x^{k+1}$ 处进行一阶泰勒近似，得：

$$\mathrm{\nabla }f(x)=\mathrm{\nabla }f(x^{k+1})+{\mathrm{\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=\mathrm{\nabla }f(x^{k+1})-\mathrm{\nabla }f(x^k)$ 为 梯度差，得：

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

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

$${\mathrm{\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{)}^T{y}^k>0$$

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

---

3. SR1 方法

3.1 SR1 定义

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

$$H^{k+1}=H^k+\mathrm{\Delta }{H}^k$$

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

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

需要保证 $H^k$ 和 $H^{k+1}$ 都是对称的，故显然 $\mathrm{\Delta }{H}^k$ 也是对称的。所以令 $\beta \in \mathbb{R},u\in {\mathbb{R}}^n$，使得 $\mathrm{\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^k{y}^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^k{y}^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{array}{rl}{H}^{k+1}& =H^k+\beta u{u}^T\\ \\ & =H^k+\beta {\gamma }^2(s^k-H^k{y}^k)(s^k-H^k{y}^k{)}^T\\ \\ & =H^k+\frac{(s^k-H^k{y}^k)(s^k-H^k{y}^k{)}^T}{(s^k-H^k{y}^k{)}^T{y}^k}\end{array}$$

记 $v=s^k-H^k{y}^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^T{y}^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\ne 0,v\ne 0,u,v\in {\mathbb{R}}^n$ 以及 $a,b\in \mathbb{R}$，其中设定的向量和标量都是待定的，则有 秩二更新表达式：

$$B^{k+1}=B^k+au{u}^T+bv{v}^T$$

显然 $au{u}^T$ 和 $bv{v}^T$ 都是对称的。

4.2 BFGS 更新公式

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

$$B^{k+1}{s}^k=(B^k+au{u}^T+bv{v}^T)s^k=y^k$$

整理可得：

$$au{u}^T{s}^k+bv{v}^T{s}^k=(a{u}^T{s}^k)u+(b{v}^T{s}^k)v=y^k-B^k{s}^k$$

可以得出 $a{u}^T{s}^k$ 和 $b{v}^T{s}^k$ 为 标量，不妨取 $(a{u}^T{s}^k)u=y^k,(b{v}^T{s}^k)v=-B^k{s}^k$，所以可以得到如下取值：

$$a{u}^T{s}^k=1,u=y^k,b{v}^T{s}^k=-1,v=B^k{s}^k$$

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

$$\begin{array}{rl}a& =\frac{1}{u^T{s}^k}=\frac{1}{(y^k{)}^T{s}^k}\\ \\ b& =-\frac{1}{v^T{s}^k}=-\frac{1}{(B^k{s}^k{)}^T{s}^k}=-\frac{1}{(s^k{)}^T{B}^k{s}^k}\end{array}$$

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

将上述取值回代到 更新表达式 $B^{k+1}=B^k+au{u}^T+bv{v}^T$，得:

$$\begin{array}{r}{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}\end{array}$$

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

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

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

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

4.3 BFGS 有效性

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

• $B^k$ 或 $H^k$ 正定

• 满足 曲率条件 $(s^k{)}^T{y}^k>0,\mathrm{\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 u{u}^T+\gamma v{v}^T$$

5.2 DFP 更新公式

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

$$(H^k+\beta u{u}^T+\gamma v{v}^T)y^k=s^k$$

整理得：

$$\beta u{u}^T{y}^k+\gamma v{v}^T{y}^k=s^k-H^k{y}^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^k{y}^k$$

不妨取 $(\beta u^T{y}^k)u=s^k,(\gamma v^T{y}^k)v=-H^k{y}^k$，所以可以得到如下取值：

$$\beta u^T{y}^k=1,u=s^k,\gamma v^T{y}^k=-1,v=H^k{y}^k$$

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

$$\begin{array}{rl}\beta & =\frac{1}{u^T{y}^k}=\frac{1}{(s^k{)}^T{y}^k}\\ \\ \gamma & =-\frac{1}{v^T{y}^k}=-\frac{1}{(H^k{y}^k{)}^T{y}^k}=-\frac{1}{(y^k{)}^T{H}^k{y}^k}\end{array}$$

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

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

$$\begin{array}{r}{H}^{k+1}=H^k+\frac{s^k(s^k{)}^T}{(s^k{)}^T{y}^k}-\frac{H^k{y}^k(y^k{)}^T{H}^k}{(y^k{)}^T{H}^k{y}^k}\end{array}$$

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

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

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

5.3 DFP 方法的劣势

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

---

参考

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

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

sherman-morrison-woodbury 公式