最优化学习笔记：无约束优化算法（5）

6.5 拟牛顿类算法

牛顿类算法是二阶算法，其理论完备，实践中经常被用于解决规模不大的优化问题，优势是其求解效率高。但对于大规模问题就会产生困难。
拟牛顿法能够在每一步以较小的计算代价生成近似矩阵，并且使用近似矩阵代替海瑟矩阵而产生的迭代序列仍具有 Q-超线性收敛的性质。拟牛顿方法不计算海瑟矩阵 ${\nabla }^{2}f(x)$，而是构造其近似矩阵 $B^k$ 或其逆的近似矩阵 $H^k$。

6.5.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})+\mathcal{O}(\Vert x-x^{k+1}{\Vert }^2)(\mathrm{6.5.1})$$

令 $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+\mathcal{O}(\Vert s^k{\Vert }^2)=y^k(\mathrm{6.5.2})$$

忽略高阶项 $\Vert s^k{\Vert }^2$，只希望海瑟矩阵的近似矩阵 $B^{k+1}$ 满足方程
$$y^k=B^{k+1}{s}^k(\mathrm{6.5.3})$$

或者其逆的近似矩阵 $H^{k+1}$ 满足方程
$$s^k=H^{k+1}{y}^k(\mathrm{6.5.4})$$

并称 (6.5.3) 式和 (6.5.4) 式为割线方程。

还可以从另一个角度理解割线方程 (6.5.3). 牛顿法本质上是对目标函数 $f(x)$ 在迭代点 $x^k$ 处做二阶近似然后求解. 考虑 $f(x)$ 在点 $x^{k+1}$ 处的二阶近似
$$m_{k+1}(d)=f(x^{k+1})+\nabla f(x^{k+1}{)}^{\mathrm{T}}d+\frac{1}{2}{d}^{\mathrm{T}}{B}^{k+1}d(\mathrm{6.5.5})$$

我们要求 $m_{k+1}(d)$ 在 $d=-s^k$ 和 $d=0$ 处的梯度与 $f(x)$ 在 $x=x^k$ 和 $x=x^{k+1}$ 处的梯度分别保持一致. 注意到 $\nabla m_{k+1}(0)=\nabla f(x^{k+1})$ 是自然满足的，为了使得第一个条件满足只需
$$\nabla m_{k+1}(-s^k)=\nabla f(x^{k+1})-B^{k+1}{s}^k=\nabla f(x^k)(\mathrm{6.5.6})$$

整理即可得到 (6.5.3) 式。
另外，注意到近似矩阵 $B^k$ 的正定性是一个很关键的因素，在 (6.5.3) 式两边同时左乘 $(s^k{)}^{\mathrm{T}}$ 可得 $(s^k{)}^{\mathrm{T}}{B}^{k+1}{s}^k=(s^k{)}^{\mathrm{T}}{y}^k$，因此条件

$$(s^k{)}^{\mathrm{T}}{y}^k>0(\mathrm{6.5.7})$$

为 $B^{k+1}$ 正定的一个必要条件。曲率条件: 在迭代过程中始终满足 $(s^k{)}^{\mathrm{T}}{y}^k>0$。对于一般的目标函数 $f(x)$，需要使用 Wolfe 准则线搜索来保证曲率条件成立。实际上，根据 Wolfe 准则的条件 (6.1.4b) 有 $\nabla f(x^{k+1}{)}^{\mathrm{T}}{s}^k⩾c_2\nabla f(x^k{)}^{\mathrm{T}}{s}^k$，其中 $c_2\in (0,1)$，两边同时减去 $\nabla f(x^k{)}^{\mathrm{T}}{s}^k$ 得，

$$(y^k{)}^{\mathrm{T}}{s}^k⩾(c_2-1)\nabla f(x^k{)}^{\mathrm{T}}{s}^k>0$$

这是因为 $c_2-1<0$ 且 $s^k={\alpha }_k{d}^k$ 是下降方向。
注意，仅仅使用 Armijo 准则不能保证曲率条件成立。

在通常情况下，近似矩阵 $B^{k+1}$ 或 $H^{k+1}$ 是由上一步迭代加上一个修正得到的，并且要求满足割线方程 (6.5.3)。拟牛顿方法的一般框架 (算法 6.5)，下一小节将讨论一些具体的矩阵更新方式。

# 算法 6.5 拟牛顿算法(Wolf准则和BFGS拟牛顿更新方法)
def line_search_Wolf(x, alpha_hat, d, gamma = 0.5, c1 = 0.1, c2 = 0.9):
    alpha = alpha_hat
    while (f(x + alpha * d) > f(x) + c1 * alpha * np.dot(grad_f(x), d)) and (grad_f(x + alpha * d) < c2 * np.dot(grad_f(x), d)) :
        alpha = gamma * alpha
    return alpha
def quasi_newton_algorithm(x0, B0, alpha, c, gamma, tol=1e-6):
    x = x0
    B = B0
    k = 0
    x_history = [x]
    while not np.linalg.norm(grad_f(x)) < tol: # 判断当前点梯度的范数是否小于阈值 tol
        grad = grad_f(x)
        d= -np.dot(np.linalg.inv(B), grad)  # or d = -np.dot(H, gradient)
        alpha = line_search_Wolf(x, alpha, d)
        x = x + alpha * d
        x_history.append(x)
        s = alpha * d
        grad_next = grad_f(x)
        y = grad_next - grad
        B = B + np.outer(y, y) / np.dot(y, s) - np.outer(np.dot(B, s), np.dot(B, s)) / np.dot(s, np.dot(B, s))
        # np.outer(y, y)计算了向量y的外积
        k = k + 1
    return x_history, k
# 目标函数及其梯度、海瑟矩阵的计算(同例6.3)
def f(x):
    return x[0]**2 + 10*x[1]**2
def grad_f(x):
    return np.array([2*x[0], 20*x[1]])
def hess_f(x):
    return np.array([[2, 0], [0, 20]])

x0 = np.array([-10, -1])
B0 = np.eye(2)
alpha = 0.085
c = 0.1
gamma = 0.5
quasi_Newton_history, k = quasi_newton_algorithm(x0, B0, alpha, c,gamma)
print("Optimal solution:", quasi_Newton_history[-1])
'''Optimal solution: [-3.84435347e-07  3.05062008e-08]'''
print("迭代次数：", k)
'''迭代次数： 188'''

# 绘制等高线图
x = np.linspace(-12, 12, 100)
y = np.linspace(-4, 4, 100)
X, Y = np.meshgrid(x, y)
Z = X**2 + 10*Y**2 
plt.figure(figsize=(10, 5))
plt.contour(X, Y, Z, levels=np.logspace(-2, 3, 20), cmap='jet')  # levels 参数控制等高线的密集程度
plt.plot([x[0] for x in M_Newton_history], [x[1] for x in M_Newton_history], markersize = 3, marker='o', markerfacecolor='none', linewidth=1, color='r', label='modified_newton')
plt.plot([x[0] for x in quasi_Newton_history], [x[1] for x in quasi_Newton_history], markersize = 3, marker='x', markerfacecolor='none', linewidth=1, color='b', label='Quasi-Newton')
plt.legend(loc='best', markerscale=1.5)  # 设置标签位置为最佳位置，调整图例的大小
plt.show()

在实际应用中基于 $H^k$ 的拟牛顿法更加实用，这是因为根据 $H^k$ 计算下降方向 $d^k$ 不需要求解线性方程组。基于 $B^k$ 的拟牛顿法有比较好的理论性质，产生的迭代序列比较稳定。如果有办法快速求解线性方程组，我们也可采用基于 $B^k$ 的拟牛顿法。此外在某些场景下，比如有些带约束的优化问题的算法设计，由于需要用到海瑟矩阵的近似，$B^k$ 的使用也很常见。

6.5.2 拟牛顿矩阵更新方式

这一小节介绍一些常见的拟牛顿矩阵的更新方式 秩一更新 (SR1)、BFGS 公式和DFP 公式。

1. 秩一更新 (SR1)
秩一更新 (SR1) 公式是结构最简单的拟牛顿矩阵更新公式．设 $B^k$ 是第 $k$ 步的近似海瑟矩阵，我们通过对 $B^k$ 进行秩一修正得到 $B^{k+1}$，使其满足割线方程 (6.5.3)．为此使用待定系数法求出修正矩阵，并设
$$B^{k+1}=B^k+au{u}^{\mathrm{T}}(\mathrm{6.5.8})$$

其中 $u\in {\mathbb{R}}^n,a\in \mathbb{R}$ 待定. 根据割线方程 (6.5.3)
$$B^{k+1}{s}^k=(B^k+au{u}^{\mathrm{T}})s^k=y^k$$

进而得到
$$au{u}^{\mathrm{T}}{s}^k=(a\cdot u^{\mathrm{T}}{s}^k)u=y^k-B^k{s}^k$$

由于 $a\cdot u^{\mathrm{T}}{s}^k$ 是一个标量，因此 $u$ 和 $y^k-B^k{s}^k$ 方向相同.
令 $u=y^k-B^k{s}^k$，代入上式可知
$$a((y^k-B^k{s}^k{)}^{\mathrm{T}}{s}^k)(y^k-B^k{s}^k)=y^k-B^k{s}^k$$

如果假设 $(y^k-B^k{s}^k{)}^{\mathrm{T}}{s}^k\ne 0$，可以得到 $a=\frac{1}{(y^k-B^k{s}^k{)}^{\mathrm{T}}{s}^k}$，最终得到更新公式为
$$B^{k+1}=B^k+\frac{(y^k-B^k{s}^k)(y^k-B^k{s}^k{)}^{\mathrm{T}}}{(y^k-B^k{s}^k{)}^{\mathrm{T}}{s}^k}(\mathrm{6.5.9})$$

我们称 (6.5.9) 式为基于 $B^k$ 的 SR1 公式. 由完全一样的过程我们可以根据割线方程 (6.5.4) 得到基于 $H^k$ 的 SR1 公式：

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

基于 $H^k$ 的 SR1 公式推导(待定系数法):
设 $H^k$ 是第 $k$ 步的海瑟矩阵逆的近似矩阵，通过对 $H^k$ 进行秩一修正得到 $H^{k+1}$，使其满足割线方程 (6.5.4)．使用待定系数法求出修正矩阵，并设
$$H^{k+1}=H^k+bv{v}^{\mathrm{T}}$$

其中 $v\in {\mathbb{R}}^n,b\in \mathbb{R}$ 待定. 根据割线方程 (6.5.4)
$$H^{k+1}{y}^k=(H^k+bv{v}^{\mathrm{T}})y^k=s^k$$

进而得到
$$bv{v}^{\mathrm{T}}{y}^k=(b\cdot v^{\mathrm{T}}{y}^k)v=s^k-H^k{y}^k$$

由于 $b\cdot v^{\mathrm{T}}{y}^k$ 是一个标量，因此 $v$ 和 $s^k-H^k{y}^k$ 方向相同.
令 $v=s^k-H^k{y}^k$，代入上式可知
$$b((s^k-H^k{y}^k{)}^{\mathrm{T}}{y}^k)(s^k-H^k{y}^k)=s^k-H^k{y}^k$$

假设 $(s^k-H^k{y}^k{)}^{\mathrm{T}}{y}^k\ne 0$，可以得到 $b=\frac{1}{(s^k-H^k{y}^k{)}^{\mathrm{T}}{y}^k}$，最终得到更新公式为
$$H^{k+1}=H^k+\frac{(s^k-H^k{y}^k)(s^k-H^k{y}^k{)}^{\mathrm{T}}}{(s^k-H^k{y}^k{)}^{\mathrm{T}}{y}^k}(\mathrm{6.5.10})$$

观察秩一公式，发现基于 $B^k$ 和 $H^k$ 的公式在形式上互为对偶，将公式 (6.5.9) 里的变量作如下替换：
$$B^k\to H^k,s^k\leftrightarrow y^k,$$

即可得到公式 (6.5.10)。实际上如果增加假设 $H^k=(B^k{)}^{-1}$，公式 (6.5.10) 也可由秩一更新的 SMW 公式 (B.1.8) 推出，反之亦然。

SR1 公式虽然结构简单，但是有一个重大缺陷：它不能保证矩阵在迭代过程中保持正定。

2. BFGS 公式
为了克服 SR1 公式的缺陷，现在考虑对 $B^k$ 的秩二更新．同样地，采用待定系数法来推导此公式．设
$$B^{k+1}=B^k+au{u}^{\mathrm{T}}+bv{v}^{\mathrm{T}}(\mathrm{6.5.11})$$

其中 $u,v\in {\mathbb{R}}^n,a,b\in \mathbb{R}$ 待定. 根据割线方程(6.5.3)，

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

整理可得

$$(a\cdot u^{\mathrm{T}}{s}^k)u+(b\cdot v^{\mathrm{T}}{s}^k)v=y^k-B^k{s}^k$$

通过选取 $u$ 和 $v$ 让以上等式成立即可.
一种最直接的取法是让上面等式左右两边的两项分别对应相等，即
$$\begin{array}{rlr}u=y^k,& & a\cdot u^{\mathrm{T}}{s}^k=1,\\ v=B^k{s}^k,& & b\cdot v^{\mathrm{T}}{s}^k=-1.\end{array}$$

将上述参量带入割线方程，得到更新方式
$$B^{k+1}=B^k+\frac{u{u}^{\mathrm{T}}}{(s^k{)}^{\mathrm{T}}u}-\frac{v{v}^{\mathrm{T}}}{(s^k{)}^{\mathrm{T}}v}=B^k+\frac{y^k(y^k{)}^{\mathrm{T}}}{(s^k{)}^{\mathrm{T}}{y}^k}-\frac{B^k{s}^k(B^k{s}^k{)}^{\mathrm{T}}}{(s^k{)}^{\mathrm{T}}{B}^k{s}^k}(\mathrm{6.5.12})$$

格式 (6.5.12) 被称为基于 $B^k$ 的 BFGS 公式，它是由 Broyden, Fletcher, Goldfarb, Shanno 四人名字的首字母组成.

根据 SMW 公式 (B.1.8) 并假设 $H^k=(B^k{)}^{-1}$, 可推出基于 $H^k$ 的 BFGS 公式：
$$H^{k+1}=(I-{\rho }_k{y}^k(s^k{)}^{\mathrm{T}}{)}^{\mathrm{T}}{H}^k(I-{\rho }_k{y}^k(s^k{)}^{\mathrm{T}})+{\rho }_k{s}^k(s^k{)}^{\mathrm{T}}(\mathrm{6.5.13})$$

其中 ${\rho }_k=\frac{1}{(s^k{)}^{\mathrm{T}}{y}^k}$.

容易看出，若要 BFGS 公式更新产生的矩阵 $H^{k+1}$ 或 $B^{k+1}$ 正定，一个充分条件是曲率条件 $(s^k{)}^{\mathrm{T}}{y}^k>0$ 成立且上一步更新矩阵 $H^k$ 或 $B^k$ 正定. 在问题求解过程中，曲率条件不一定会得到满足，可以使用 Wolfe 准则的线搜索来迫使曲率条件成立.

BFGS 格式 (6.5.13) 还有更深刻的含义，它其实满足了某种逼近的最优性. 具体来说，基于 $H^k$ 的 BFGS 公式恰好是如下优化问题的解：
$$\begin{array}{rl}\underset{H}{\min }& \Vert H-H^k{\Vert }_W\\ s.t.& H=H^{\mathrm{T}}\\ & H{y}^k=s^k\end{array}(\mathrm{6.5.14})$$

这里 $\Vert \cdot {\Vert }_W$ 是加权范数，定义为
$$\Vert H{\Vert }_W=\Vert W^{1/2}H{W}^{1/2}{\Vert }_F$$

其中 $W$ 可以是任意满足割线方程 $W{s}^k=y^k$ 的矩阵.
这个优化问题的含义是在满足割线方程 $s^k=H^{k+1}{y}^k$ 的对称矩阵中找到离 $H^k$ 最近的矩阵 $H$.
BFGS 公式是目前最有效的拟牛顿更新格式之一，它有比较好的理论性质，实现起来也并不复杂．对格式 (6.5.13) 进行改动可得到有限内存 BFGS 格式 (L-BFGS)，它是常用的处理大规模优化问题的拟牛顿类算法，在本节的末尾介绍这一算法．

3. DFP 公式
在 BFGS 公式的推导中，如果利用割线方程 $s^k=H^{k+1}{y}^k$ 对 $H^k$ 推导秩二修正的拟牛顿修正，将得到基于 $H^k$ 的拟牛顿矩阵更新
$$H^{k+1}=H^k-\frac{H^k{y}^k(H^k{y}^k{)}^{\mathrm{T}}}{(y^k{)}^{\mathrm{T}}{H}^k{y}^k}+\frac{s^k(s^k{)}^{\mathrm{T}}}{(y^k{)}^{\mathrm{T}}{s}^k}(\mathrm{6.5.15})$$

这种迭代格式首先由 Davidon 发现，此后由 Fletcher 以及 Powell 进一步发展，因此被称为 DFP 公式. 根据 SMW 公式以及 $B^k=(H^K{)}^{-1}$ 可得其关于 $B^k$ 的更新格式
$$\begin{array}{r}{B}^{k+1}=(I-{\rho }_k{s}^k(y^k{)}^{\mathrm{T}}{)}^{\mathrm{T}}{B}^k(I-{\rho }_k{s}^k(y^k{)}^{\mathrm{T}})+{\rho }_k{y}^k(y^k{)}^{\mathrm{T}}\\ \\ =(I-\frac{s^k(y^k{)}^{\mathrm{T}}}{(s^k{)}^{\mathrm{T}}{y}^k}{)}^{\mathrm{T}}{B}^k(I-\frac{s^k(y^k{)}^{\mathrm{T}}}{(s^k{)}^{\mathrm{T}}{y}^k})+\frac{y^k(y^k{)}^{\mathrm{T}}}{(s^k{)}^{\mathrm{T}}{y}^k}\end{array}(\mathrm{6.5.16})$$

其中 ${\rho }_k=\frac{1}{(s^k{)}^{\mathrm{T}}{y}^k}$，同 (6.5.13) 式.
可以看到，DFP 公式 (6.5.15) (6.5.16) 和 BFGS 公式 (6.5.12) (6.5.13) 分别呈对偶关系. 将 BFGS 格式 (6.5.13) 中的 $H^k$ 换成 $B^k$, $s^k$ 与 $y^k$ 对换便得到了 DFP 格式 (6.5.16). 不仅如此，在逼近性上也有这样的对偶现象. 基于 $B^k$ 的 DFP 格式是如下优化问题的解：
$$\begin{array}{rl}\underset{B}{\min }& \Vert B-B^k{\Vert }_W\\ s.t.& B=B^{\mathrm{T}}\\ & B{s}^k=y^k\end{array}(\mathrm{6.5.17})$$

其中 $\Vert \cdot {\Vert }_W$ 是加权范数，定义为 $\Vert B{\Vert }_W=\Vert W^{1/2}B{W}^{1/2}{\Vert }_F$ 且 $W$ 为任意满足 $W{y}^k=s^k$ 的矩阵. 和 BFGS 格式类似，DFP 格式要求 $B^{k+1}$ 为满足割线方程 (6.5.3) 的对称矩阵中离 $B^k$ 最近的矩阵，它也暗含某种最优性.
遗憾的是，尽管 DFP 格式在很多方面和 BFGS 格式存在对偶关系，但从实际效果来看，DFP 格式整体上不如 BFGS 格式. 因此在实际使用中人们更多使用 BFGS 格式.

拟牛顿矩阵更新方式：

method	$B^{k+1}$	$H^{k+1}$
SR1	$B^k+\frac{(y^k-B^k{s}^k)(y^k-B^k{s}^k{)}^{\mathrm{T}}}{(y^k-B^k{s}^k{)}^{\mathrm{T}}{s}^k}$	$H^k+\frac{(s^k-H^k{y}^k)(s^k-H^k{y}^k{)}^{\mathrm{T}}}{(s^k-H^k{y}^k{)}^{\mathrm{T}}{y}^k}$
BFGS	$B^k+\frac{y^k(y^k{)}^{\mathrm{T}}}{(s^k{)}^{\mathrm{T}}{y}^k}-\frac{B^k{s}^k(B^k{s}^k{)}^{\mathrm{T}}}{(s^k{)}^{\mathrm{T}}{B}^k{s}^k}$	${\left(I-\frac{y^k(s^k{)}^{\mathrm{T}}}{(s^k{)}^{\mathrm{T}}{y}^k}\right)}^{\mathrm{T}}{H}^k\left(I-\frac{y^k(s^k{)}^{\mathrm{T}}}{(s^k{)}^{\mathrm{T}}{y}^k}\right)+\frac{s^k(s^k{)}^{\mathrm{T}}}{(s^k{)}^{\mathrm{T}}{y}^k}$
DFP	${\left(I-\frac{s^k(y^k{)}^{\mathrm{T}}}{(s^k{)}^{\mathrm{T}}{y}^k}\right)}^{\mathrm{T}}{B}^k\left(I-\frac{s^k(y^k{)}^{\mathrm{T}}}{(s^k{)}^{\mathrm{T}}{y}^k}\right)+\frac{y^k(y^k{)}^{\mathrm{T}}}{(s^k{)}^{\mathrm{T}}{y}^k}$	$H^k-\frac{H^k{y}^k(H^k{y}^k{)}^{\mathrm{T}}}{(y^k{)}^{\mathrm{T}}{H}^k{y}^k}+\frac{s^k(s^k{)}^{\mathrm{T}}}{(y^k{)}^{\mathrm{T}}{s}^k}$

6.5.3 拟牛顿法的全局收敛性

本小节介绍拟牛顿法的收敛性质. 首先利用 Zoutendijk 条件得到拟牛顿法基本的收敛性，之后简要介绍收敛速度.

定理 6.8 (BFGS 全局收敛性) 假设初始矩阵 $B^0$ 是对称正定矩阵，目标函数 $f(x)$ 是二阶连续可微函数，且下水平集
L = { x ∈ R n ∣ f ( x ) ⩽ f ( x 0 ) } {\mathcal{L}}=\{x\in\mathbb{R}^{n}\mid f(x)\leqslant f(x^{0})\} L={x∈Rn∣f(x)⩽f(x0)}

是凸的，并且存在正数 $m$ 以及 $M$ 使得对于任意的 $z\in {\mathbb{R}}^n$ 以及任意的 $x\in \mathcal{L}$ 有
$$m\Vert z{\Vert }^2⩽z^{\mathrm{T}}{\nabla }^{2}f(x)z⩽M\Vert z{\Vert }^2(\mathrm{6.5.18})$$

则采用 基于 $B^k$ 的 BFGS 公式并结合 Wolfe 线搜索的拟牛顿算法全局收敛到 $f(x)$ 的极小值点 $x^{\ast }$.

证明 为了方便，首先定义 $$m_k=\frac{(y^k{)}^{\mathrm{T}}{s}^k}{(s^k{)}^{\mathrm{T}}{s}^k},M_k=\frac{(y^k{)}^{\mathrm{T}}{y}^k}{(y^k{)}^{\mathrm{T}}{s}^k}$$
由 (6.5.18) 式以及 $y^k={\int }_0^1{\nabla }^{2}f(x^k+t(x^{k+1}-x^k))s^k\text{dt}$ 可得 $$m_k⩽m,M_k⩾M$$
根据基于 $B^k$ 的 BFGS 公式(6.5.12)，两边同时求迹，得到$$\text{Tr}(B^{k+1})-\frac{\Vert B^k{s}^k{\Vert }^2}{(s^k{)}^{\text{T}}{B}^k{s}^k}+\frac{\Vert y^k{\Vert }^2}{(y^k{)}^{\text{T}}{s}^k}(\mathrm{6.5.19})$$

又因为 $$\det B^{k+1}=\det B^k\frac{(y^k{)}^{\mathrm{T}}{s}^k}{(s^k{)}^{\mathrm{T}}{B}^k{s}^k}(\mathrm{6.5.20})$$

对于向量 $u,v,x,y$ 利用分块矩阵乘法及迹的性质，有
$$\begin{array}{rl}\det (I+u{v}^{\mathrm{T}}+x{y}^{\mathrm{T}})& =\det \left[\begin{array}{ccc}I+u{v}^{\mathrm{T}}+x{y}^{\mathrm{T}}& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]\\ & =\det \left[\begin{array}{ccc}I& -u& -x\\ v^{\mathrm{T}}& 1& 0\\ y^{\mathrm{T}}& 0& 1\end{array}\right]\\ & =\det \left[\begin{array}{ccc}I& -u& -x\\ 0& 1+v^{\mathrm{T}}u& v^{\mathrm{T}}x\\ 0& y^{\mathrm{T}}u& 1+y^{\mathrm{T}}x\end{array}\right]\\ & =(1+y^{\mathrm{T}}x)(1+v^{\mathrm{T}}u)-(v^{\mathrm{T}}x)(y^{\mathrm{T}}u).\end{array}$$

因此
$$\begin{array}{cc}\det B^{k+1}& =\det \left(B^k+\frac{y^k{\left(y^k\right)}^{\mathrm{T}}}{{\left(s^k\right)}^{\mathrm{T}}{y}^k}-\frac{B^k{s}^k{\left(B^k{s}^k\right)}^{\mathrm{T}}}{{\left(s^k\right)}^{\mathrm{T}}{B}^k{s}^k}\right)\\ & =\det \left(B^k(I+\frac{(B^k{)}^{-1}{y}^k{\left(y^k\right)}^{\mathrm{T}}}{{\left(s^k\right)}^{\mathrm{T}}{y}^k}-\frac{s^k{\left(B^k{s}^k\right)}^{\mathrm{T}}}{{\left(s^k\right)}^{\mathrm{T}}{B}^k{s}^k})\right)\\ & =\det B^k\det (I+u{v}^{\mathrm{T}}+x{y}^{\mathrm{T}}).\end{array}$$

其中
$$u=\frac{(B^k{)}^{-1}{y}^k}{{\left(s^k\right)}^{\text{T}}{y}^k},v=y^k,x=\frac{s^k}{{\left(s^k\right)}^{\text{T}}{B}^k{s}^k},y=-B^k{s}^k.$$

利用上述性质，
$$\det B^{k+1}=\left((1+y^{\mathrm{T}}x)(1+v^{\mathrm{T}}u)-(v^{\mathrm{T}}x)(y^{\mathrm{T}}u)\right)\det B^k=\det B^k\frac{(y^k{)}^{\mathrm{T}}{s}^k}{(s^k{)}^{\mathrm{T}}{B}^k{s}^k}$$

接下来，定义 $$\cos {\theta }_k=\frac{(s^k{)}^{\mathrm{T}}{B}^k{s}^k}{\Vert s^k\Vert \Vert B^k{s}^k\Vert },q_k=\frac{(s^k{)}^{\mathrm{T}}{B}^k{s}^k}{(s^k{)}^{\mathrm{T}}{s}^k},$$
所以 $$\frac{\Vert B^k{s}^k{\Vert }^2}{(s^k{)}^{\text{T}}{B}^k{s}^k}=\frac{\Vert B^k{s}^k{\Vert }^2\Vert s^k{\Vert }^2}{((s^k{)}^{\text{T}}{B}^k{s}^k{)}^2}\frac{(s^k{)}^{\text{T}}{B}^k{s}^k}{\Vert s^k{\Vert }^2}=\frac{q_k}{{\cos }^2{\theta }_k}$$

$$\det B^{k+1}=\det B^k\frac{(y^k{)}^{\mathrm{T}}{s}^k}{(s^k{)}^{\mathrm{T}}{s}^k}\frac{(s^k{)}^{\mathrm{T}}{s}^k}{(s^k{)}^{\mathrm{T}}{B}^k{s}^k}=\det B^k\frac{m_k}{q_k}$$

再引进矩阵辅助函数 $$\psi (B)=\mathrm{Tr}(B)-\ln \det B,$$

那么，$$\begin{array}{rl}& \psi (B^{k+1})\\ & =\mathrm{Tr}(B^k)+M_k-\frac{q_k}{{\cos }^2{\theta }_k}-\ln \det B^k-\ln m_k+\ln q_k\\ & =\psi (B^k)+(M_k-\ln m_k-1)+\left(1-\frac{q_k}{{\cos }^2{\theta }_k}+\ln \frac{q_k}{{\cos }^2{\theta }_k}\right)+\ln {\cos }^2{\theta }_k\end{array}(\mathrm{6.5.21})$$

又因为 $\psi (B)>0,1-t+lnt⩽0,\forall t$ 可以得到
$$0<\psi (B^{k+1})⩽\psi (B^k)+(M_k-\ln m_k-1)+\ln {\cos }^2{\theta }_k⩽\psi (B^0)+(k+1)c+\sum _{j=0}^k\ln {\cos }^2{\theta }_j(\mathrm{6.5.22})$$

这里 $c=M-\ln m-1$ 且假设 $c>0$. 注意到 $s^k=-{\alpha }_k(B^k{)}^{-1}\nabla f(x^k)$ 是搜索方向，那么 $\cos {\theta }_k$ 是搜索方向与梯度方向夹角的余弦. 根据定理 6.1, 可知 $\Vert \nabla f(x^k)\Vert$ 大于某个非零常数仅当 $\cos {\theta }_k\to 0.$ 因此，为了证明 $\Vert \nabla f(x^k)\Vert \to 0$, 我们仅需证明 $\cos {\theta }_k\to 0$ 不成立，下面用反证法证明这一结论. 假设 $\cos {\theta }_k\to 0$, 那么存在 $k_1>0$, 对于任意的 $j>k_1$,
$$\begin{array}{r}\ln {\cos }^2{\theta }_j<-2c.\end{array}$$

结合(6.5.22)式，当 $k>k_1$ 时，
$$\begin{array}{rl}& \text{0}<\psi (B^{k+1})⩽\psi (B^0)+(k+1)c+\sum _{j=0}^{k_1}\ln {\cos }^2{\theta }_j+\sum _{j=k_1}^k(-2c)(\mathrm{6.5.23})\\ & =\psi (B^0)+\sum _{j=0}^{k_1}\ln {\cos }^2{\theta }_j+2c{k}_1+c-ck(\mathrm{6.5.24})\end{array}$$

上式的右边对于充分大的 $k$ 是负的，而不等式左边是 0, 这矛盾，因此假设不成立，即存在一个子列 { j k } k = 1 , 2 , ⋯ \{j_k\}_{k=1,2,\cdots} {jk​}k=1,2,⋯​ 使得 $\cos {\theta }_{j_k}⩾\delta >0$. 根据 Zoutendijk 条件 (定理 6.1), 可以得到 $\underset{k\to \infty }{\lim }\inf \Vert \nabla f(x^k)\Vert \to 0.$ 又因为问题对$x\in \mathcal{L}$ 是强凸的，所以这可以导出 $x^k\to x^{\ast }.$ 得证（若函数的性质稍差, 即函数并非强凸的, 则全局收敛性可能将退化为局部的收敛性）

定理 6.8 叙述了 BFGS 格式的全局收敛性，但没有说明以什么速度收敛．下面这个定理介绍了在一定条件下 BFGS 格式会达到 Q-超线性收敛速度．

定理 6.9 (BFGS 收敛速度) 设 $f(x)$ 二阶连续可微，在最优点 $x^{\ast }$ 的一个邻域内海瑟矩阵利普希茨连续，且使用 BFGS 迭代格式收敛到 $f$ 的最优值点 $x^{\ast }$. 若迭代点列 { x k } \{x^k\} {xk} 满足
$$\sum _{k=1}^{\infty }\Vert x^k-x^{\star }\Vert <+\infty (\mathrm{6.5.25})$$

则 { x k } \{x^k\} {xk} 以 Q-超线性收敛到 $x^{\ast }$.

由于仅仅使用了海瑟矩阵的近似矩阵，拟牛顿法只能达到 Q-超线性收敛速度，这个速度和牛顿法相比较慢. 但由于拟牛顿法不需要每一步都计算海瑟矩阵，它在整体计算开销方面可能远远小于牛顿法，因此在实际问题中较为实用.

6.5.4 有限内存 BFGS 方法

拟牛顿法虽然克服了计算海瑟矩阵的困难，但是它仍然无法应用在大规模优化问题上. 一般来说，拟牛顿矩阵 $B^k$ 或 $H^k$ 是稠密矩阵（稠密矩阵指的是矩阵中大多数元素都是非零元素的矩阵），而存储稠密矩阵要消耗 $\mathcal{O}(n^2)$ 的内存，这对于大规模问题显然是不可能实现的. 在本小节介绍的有限内存 BFGS 方法（L-BFGS）解决了这一存储问题，从而使得人们在大规模问题上也可应用拟牛顿类方法加速迭代的收敛．
L-BFGS 方法是根据 BFGS 公式 (6.5.12) (6.5.13) 变形而来的．为了推导方便，我们以 $H^k$ 的更新公式 (6.5.13) 为基础来推导相应的 L-BFGS 公式．
首先引入新的记号重写 (6.5.13) 式：

$$H^{k+1}={\left(V^k\right)}^{\mathrm{T}}{H}^k{V}^k+{\rho }_k{s}^k(s^k{)}^{\mathrm{T}}(\mathrm{6.5.26})$$

其中
$${\rho }_k=\frac{1}{(y^k{)}^{\mathrm{T}}{s}^k},V^k=I-{\rho }_k{y}^k(s^k{)}^{\mathrm{T}}(\mathrm{6.5.27})$$

观察到 (6.5.26) 式有类似递推的性质，为此可将 (6.5.26) 式递归地展开 $m$ 次，其中 $m$ 是一个给定的整数：

$$H^k=(V^{k-m}\cdots V^{k-1}{)}^{\mathrm{T}}{H}^{k-m}(V^{k-m}\cdots V^{k-1})+{\rho }_{k-m}(V^{k-m+1}\cdots V^{k-1}{)}^{\mathrm{T}}{s}^{k-m}(s^{k-m}{)}^{\mathrm{T}}(V^{k-m+1}\cdots V^{k-1})+{\rho }_{k-m+1}(V^{k-m+2}\cdots V^{k-1}{)}^{\mathrm{T}}{s}^{k-m+1}(s^{k-m+1}{)}^{\mathrm{T}}(V^{k-m+2}\cdots V^{k-1})+\cdots +{\rho }_{k-1}{s}^{k-1}(s^{k-1}{)}^{\mathrm{T}}(\mathrm{6.5.28})$$

为了达到节省内存的目的，(6.5.26) 式不能无限展开下去，但这会发现 $H^{k-m}$ 还是无法显式求出. 一个很自然的想法就是用 $H^{k-m}$ 的近似矩阵来代替 $H^{k-m}$ 进行计算.
在拟牛顿迭代中，实际上并不需要计算 $H^k$ 的显式形式，只需要利用 $H^k\nabla f(x^k)$ 来计算迭代方向 $d^k$. 为此先直接给出一个利用展开式 (6.5.28) 直接求解 $H^k\nabla f(x^k)$ 的算法，见算法 6.6. 它充分利用了 (6.5.28) 式的结构来尽量节省计算 $H^k\nabla f(x^k)$ 的开销. 由于其主体结构包含了方向相反的两个循环，因此它也被称为双循环递归算法.

我们现在给出双循环递归算法的执行过程. 在 (6.5.28) 式中，等式左右两边同时右乘 $\nabla f(x^k)$，则等式左侧为 $-d^k$. 若只观察等式右侧，则需要计算
$$V^{k-1}\nabla f(x^k),V^{k-2}{V}^{k-1}\nabla f(x^k),\cdots ,V^{k-m}\cdots V^{k-2}{V}^{k-1}\nabla f(x^k).$$

这些结果可以递推地进行，无需重复计算.
同时，在计算 $V^{k-l}\cdots V^{k-1}\nabla f(x^k)$ 的过程中恰好可以计算上一步的 ${\rho }_{k-l}(s^{k-l}{)}^{\mathrm{T}}[V^{k-l+1}$ $\cdots V^{k-1}\nabla f(x^k)]$，这是一个标量，可以记 $$\begin{array}{rl}q& =V^{k-m}\cdots V^{k-1}\nabla f(x^k),\\ & \\ {\alpha }_i& ={\rho }_{k-l}(s^{k-l}{)}^{\mathrm{T}}[V^{k-l+1}\cdots V^{k-1}\nabla f(x^k)],\end{array}$$

递归公式可化为如下的形式：
$$\begin{array}{rl}{H}^k\nabla f(x^k)=& (V^{k-m}\cdots V^{k-1}{)}^{\mathrm{T}}{H}^{k-m}q+(V^{k-m+1}\cdots V^{k-1}{)}^{\mathrm{T}}{s}^{k-m}{\alpha }_{k-m}+\\ & (V^{k-m+2}\cdots V^{k-1}{)}^{\mathrm{T}}{s}^{k-m+1}{\alpha }_{k-m+1}+\cdots +s^{k-1}{\alpha }_{k-1}.\end{array}(\mathrm{6.5.29})$$

公式 (6.5.29) 已经简化了不少，接下来的第二个循环就是自上而下合并每一项. 以合并前两项为例，它们有公共的因子 $(V^{k-m+1}\cdots V^{k-1}{)}^{\mathrm{T}}$, 提取出来之后前两项的和可以写为

$$(V^{k-m+1}\cdots V^{k-1}{)}^{\mathrm{T}}((V^{k-m}{)}^{\mathrm{T}}r+{\alpha }_{k-m}{s}^{k-m})=(V^{k-m+1}\cdots V^{k-1}{)}^{\mathrm{T}}(r+({\alpha }_{k-m}-\beta )s^{k-m})$$

这正是第二个循环的迭代格式. 注意合并后 (6.5.29) 式的结构仍不变，因此可递归地计算下去，最后变量 $r$ 就是我们期望的结果 $H^k\nabla f(x^k)$.

L-BFGS 双循环递归算法约需要 $4mn$ 次乘法运算与 $2mn$ 次加法运算，若近似矩阵 ${\hat{H}}^{k-m}$ 是对角矩阵，则额外需要 $n$ 次乘法运算. 由于 $m$ 不会很大，因此该算法的复杂度是 $\mathcal{O}(mn)$. 算法所需要的额外存储为临时变量 ${\alpha }_i$, 它的大小是 $\mathcal{O}(m)$. 综上所述，L-BFGS 双循环算法是非常高效的.

# 算法 6.6 L-BFGS 双循环递归算法
import numpy as np
def lbfgs_double_loop(x, k, m, s_hist, y_hist, rho_hist, H_hat):
    k = len(s_hist)
    q = grad_f(x)
    if k >=1:
        H_hat = np.dot(s_hist[-1], y_hist[-1]) / np.dot(y_hist[-1], y_hist[-1])
    else:
        H_hat = 1
    alpha = np.zeros(m)
    for i in range(k-1,-1,-1):
        alpha[i] = rho_hist[i] * np.dot(s_hist[i], q)
        q -= alpha[i] * y_hist[i]
    r = H_hat * q
    for i in range(k):
        beta = rho_hist[i] * np.dot(y_hist[i], r) 
        r += s_hist[i] * (float(alpha[i] - beta))  
    return r

近似矩阵 ${\hat{H}}^{k-m}$ 的取法可以是对角矩阵 ${\hat{H}}^{k-m}={\gamma }_{k}I$, 其中
$${\gamma }_k=\frac{(s^{k-1}{)}^{\mathrm{T}}{y}^{k-1}}{(y^{k-1}{)}^{\mathrm{T}}{y}^{k-1}}(\mathrm{6.5.30})$$

这恰好是 BB 方法的第一个步长.

# 算法 6.7 L-BFGS 方法
import numpy as np
import matplotlib.pyplot as plt
def line_search_Wolf(x, alpha_hat, d, gamma = 0.5, c1 = 0.1, c2 = 0.9):
    alpha = alpha_hat
    while (f(x + alpha * d) > f(x) + c1 * alpha * np.dot(grad_f(x), d)) and (grad_f(x + alpha * d) < c2 * np.dot(grad_f(x), d)) :
        alpha = gamma * alpha
    return alpha

def lbfgs_method(x0, m, alpha_hat, max_iter = 1000, tol = 1e-6):
    n = len(x0)
    x = x0
    k = 0
    s = []
    y = []
    rho = []
    alpha = []
    x_history = [x] 
    while k < max_iter:
        if k == 0:
            H_hat = np.eye(n)
        else:
            gamma_k = np.dot(s[-1], y[-1]) / np.dot(y[-1], y[-1])
            H_hat = gamma_k * np.eye(n)

        d = - lbfgs_double_loop(x, k, m, s, y, rho, H_hat)
        alpha_k = line_search_Wolf(x, alpha_hat, d)

        x_next = x + alpha_k * d
        
        if k >= m:
            s.pop(0)
            y.pop(0)
            rho.pop(0)
            alpha.pop(0)
        
        s_k = x_next - x
        y_k = grad_f(x_next) - grad_f(x)
        rho_k = 1 / np.dot(s_k, y_k)
        
        s.append(s_k)
        y.append(y_k)
        rho.append(rho_k)
        alpha.append(alpha_k)
        x = x_next
        k += 1
        x_history.append(x)
        if np.linalg.norm(grad_f(x)) < tol: # 判断当前点梯度的范数是否小于阈值 tol
            break       
    return x_history, k

# 目标函数及其梯度、海瑟矩阵的计算(同例6.3)
def f(x):
    return x[0]**2 + 10 * x[1]**2
def grad_f(x):
    return np.array([2 * x[0], 20 * x[1]])
def hess_f(x):
    return np.array([[2, 0], [0, 20]])

# 初始点
x0 = np.array([-10, -1])
alpha_hat = 0.085
m = 4 # Memory size

lbfgs_history, k = lbfgs_method(x0, m, alpha_hat)
print("Optimal solution:", lbfgs_history[-1])
'''Optimal solution: [-4.68605798e-07  4.33419093e-09]'''
print("迭代次数：", k)
'''迭代次数： 177'''
# 绘制等高线图
x = np.linspace(-12, 12, 100)
y = np.linspace(-4, 4, 100)
X, Y = np.meshgrid(x, y)
Z = X**2 + 10*Y**2 
plt.figure(figsize=(10, 5))
plt.contour(X, Y, Z, levels=np.logspace(-2, 3, 20), cmap='jet')  
plt.plot([x[0] for x in lbfgs_history], [x[1] for x in lbfgs_history], markersize = 3, marker='o', markerfacecolor='none', linewidth=1, color='r', label='lbfgs Method')
plt.plot([x[0] for x in quasi_Newton_history], [x[1] for x in quasi_Newton_history], markersize = 3, marker='x', markerfacecolor='none', linewidth=1, color='green', label='Quasi-Newton')
plt.legend(loc='best', markerscale=1.5)  # 设置标签位置为最佳位置，调整图例的大小
plt.show()

至此基本介绍了 L-BFGS 方法的迭代格式 (见算法 6.7)，下面从另一个角度出发来重新理解这个格式，进而可以直接给出 L-BFGS 格式下拟牛顿矩阵的形式. 为了讨论问题方便，我们引入新的记号
$$S^k=[s^0,s^1,\cdots ,s^{k-1}],Y^k=[y^0,y^1,\cdots ,y^{k-1}](\mathrm{6.5.31})$$

定理 6.10 设 $H^0$ 为 BFGS 格式的初始矩阵，且是对称正定的，又设 $k$ 个向量对 { s i , y i } i = 0 k − 1 \{s^i,y^i\}_{i=0}^{k-1} {si,yi}i=0k−1​ 满足 $(s^i{)}^{\mathrm{T}}{y}^i>0,H^k$ 是 BFGS 格式 (6.5.13) 产生的拟牛顿矩阵，则有
$$H^k=H^0+\left[\begin{array}{cc}{S}^k& H^0{Y}^k\end{array}\right]\left[\begin{array}{cc}{W}^k& -(R^k{)}^{-\mathrm{T}}\\ -(R^k{)}^{-1}& 0\end{array}\right]\left[\begin{array}{c}(S^k{)}^{\mathrm{T}}\\ (Y^k{)}^{\mathrm{T}}{H}^0\end{array}\right](\mathrm{6.5.32})$$

其中
$$W^k={\left((R^k{)}^{-1}\right)}^{\mathrm{T}}(D^k+(Y^k{)}^{\mathrm{T}}{H}^0{Y}^k)(R^k{)}^{-1}$$

矩阵 $R^k$ 是 $k\times k$ 上三角矩阵，其元素为
$$(R^k{)}_{ij}=\{\begin{array}{ll}(s^{i-1}{)}^{\mathrm{T}}{y}^{j-1}& i⩽j,\\ 0,& i>j\end{array}$$

$D^k=\text{Diag}((s^0{)}^{\mathrm{T}}{y}^0,(s^1{)}^{\mathrm{T}}{y}^1,\cdots ,(s^{k-1}{)}^{\mathrm{T}}{y}^{k-1})$ 为对角矩阵.

该定理的重要意义在于可在给定 $H^0,S^k,Y^k$ 的条件下直接计算出 BFGS 迭代矩阵 $H^k$. 如果 $H^0$ 是一个近似矩阵，那么 (6.5.32) 式将给出 L-BFGS 格式迭代矩阵的显式格式. 格式 (6.5.32) 也可直接用于 L-BFGS 的计算，但总计算量要比 L-BFGS 双循环递归算法多一个常数量级. 然而 (6.5.32) 涉及矩阵操作，在现代计算机实现下可以使用 BLAS2 (BLAS3) 操作更高效地执行，实际效果很可能赶超双循环算法.

利用 SMW 公式 (B.1.8) 可以求出 BFGS 基于 $B^k$ 的块迭代格式.

定理 6.11 设 $B^0$ 为 BFGS 格式的初始矩阵，且是对称正定的，又设 $k$ 个向量对 { s i , y i } i = 0 k − 1 \{s^i,y^i\}_{i=0}^{k-1} {si,yi}i=0k−1​ 满足 $(s^i{)}^{\mathrm{T}}{y}^i>0,B^k$ 是 BFGS 格式 (6.5.12) 产生的拟牛顿矩阵，则有
$$B^k=B^0-\left[\begin{array}{cc}{B}^0{S}^k& Y^k\end{array}\right]{\left[\begin{array}{cc}(S^k{)}^{\mathrm{T}}{B}^0{S}^k& L^k\\ (L^k{)}^{\mathrm{T}}& -D^k\end{array}\right]}^{-1}\left[\begin{array}{c}(S^k{)}^{\mathrm{T}}{B}^0\\ (Y^k{)}^{\mathrm{T}}\end{array}\right](\mathrm{6.5.33})$$

其中 $L^k$ 是 $k\times k$ 严格下三角矩阵，其元素为
$$(L^k{)}_{ij}=\{\begin{array}{ll}(s^{i-1}{)}^{\mathrm{T}}{y}^{j-1},& i>j,\\ 0,& i⩽j.\end{array}$$

$D^k=\mathrm{Diag}((s^0{)}^{\mathrm{T}}{y}^0,\cdots ,(s^{k-1}{)}^{\mathrm{T}}{y}^{k-1})$ 为对角矩阵.

正因为 L-BFGS 方法的出现，人们可以使用拟牛顿类算法求解优化问题．虽然有关 L-BFGS 方法的收敛性质依然很有限，但在实际应用中 L-BFGS 方法很快成为了应用最广泛的拟牛顿类算法．尽管 DFP 公式和 BFGS 公式呈对偶关系，但极少有人研究有限内存的 DFP 格式，这可能是因为 DFP 方法本身的数值表现不如BFGS方法所致.

参考教材：《最优化：建模、算法与理论》