有限内存BFGS算法

有限内存 BFGS 算法（limited-memory BFGS，L-BFGS）是拟牛顿法的一种实现：只保存最近若干次移动及梯度差，用它们计算逆曲率矩阵与梯度的乘积，而不显式存储完整矩阵。适用对象通常是能计算梯度的大规模光滑无约束优化问题；变量边界需另用相应变体处理。

保存的是曲率对

设 $f:{ℝ}^n\to ℝ$ 可微，第 $j$ 步留下 $$s_j=x_{j+1}-x_j,y_j=\mathrm{\nabla }f(x_{j+1})-\mathrm{\nabla }f(x_j),{\rho }_j=(s_j^{𝖳}{y}_j{)}^{-1}.$$ 以下要求每一对满足 $s_j^{𝖳}{y}_j>0$。内存参数 $m\ge 1$ 是最多保存的对数，不是向量维数。第 $m+1$ 对到来时丢弃最旧的一对；每次从一个新的正定基矩阵出发，将保留的更新按旧到新依次施加。这与“先形成全历史矩阵再删除几行”不同。

新对进入时只移除最旧的一对；有限内存按曲率对管理。

常取基矩阵 $H^{(0)}=\gamma I$，$\gamma >0$。有最近曲率对时，一种常用尺度是 $\gamma =s^{𝖳}y/(y^{𝖳}y)$，使近似尺度反映最近曲率。下面为便于手算使用 $\gamma =1$。不同尺度产生不同方向，报告算例时应说清。

双循环怎样代替矩阵乘法

目标是求 $r=Hg$，方向为 $-r$。设保留的索引按旧到新为 $0,\dots ,\ell -1$。先置 $q=g$，从新到旧计算 $${\alpha }_j={\rho }_j{s}_j^{𝖳}q,q\leftarrow q-{\alpha }_j{y}_j.$$ 接着取 $r=\gamma q$，再从旧到新计算 $${\beta }_j={\rho }_j{y}_j^{𝖳}r,r\leftarrow r+s_j({\alpha }_j-{\beta }_j).$$ 循环次序不能交换，因为不同曲率更新一般不交换。

正确性的关键可先看一对。将逆 BFGS 公式右乘 $g$，令 $\alpha =\rho s^{𝖳}g$，可将结果写成 $$H^{+}g=(I-\rho s{y}^{𝖳})H(g-\alpha y)+\alpha s.$$ 若已经算得 $r_0=H(g-\alpha y)$，只须再算 $\beta =\rho y^{𝖳}{r}_0$，便得到 $H^{+}g=r_0+s(\alpha -\beta )$。多对情形递归展开同一恒等式：第一循环依次准备各层的右侧向量；基矩阵完成最内层乘法；第二循环依次补回各层左侧修正。这证明双循环与所定义的有限次 BFGS 乘积完全相同，并非另一种经验近似。

一对曲率信息的数值轨迹

取 $s=(1,1)$、$y=(1,4)$、$g=(2,1)$、$\gamma =1$。于是 $\rho =1/5$，第一循环给出 $$\alpha =3/5,q=(2,1)-\frac{3}{5}(1,4)=(7/5,-7/5).$$ 中间乘法仍为 $r_0=q$，随后 $$\beta =\frac{1}{5}(1,4{)}^{𝖳}(7/5,-7/5)=-21/25,r=(71/25,1/25).$$ 若显式建立矩阵，结果是 $H=\frac{1}{25}\left(\begin{array}{cc}37& -3\\ -3& 7\end{array}\right)$，乘 $g$ 恰为同一向量。因此方向为 $p=(-2.84,-0.04)$，而 $g^{𝖳}p=-5.72<0$。

先从 g 消去 y 方向，再补回 s 方向，得到 r=Hg。

变式：保持这对 $s,y$ 不变，把输入梯度改为 $g=y$。第一循环的 $\alpha =1$，故 $q=0$、$\beta =0$，最后 $r=s$，正好重建割线条件。若将 $\gamma$ 改为 $5/17$，这个变式的结果仍是 $s$，但原来的 $g=(2,1)$ 会得到另一乘积。

方向正确与迭代收敛的区别

保留对均满足正曲率、基矩阵正定时，逐次 BFGS 更新保持正定，所以非零梯度产生下降方向。随后仍需线搜索选择步长。仅保存有限历史一般不会保持所有旧割线条件，也不继承全 BFGS 的所有收敛速率结论。经典全局收敛分析对目标的曲率、初始水平集及线搜索设有条件；“全局收敛”在这类文献中也不应被误读为任意非凸函数都能到达全局最优点。

不合格对不能硬除：若 $s=(1,0),y=(0,1)$，内积为零；若内积很小为正，$\rho$ 很大也可能导致数值不稳。实现可按经过说明的阈值跳过或阻尼，而不能把修改隐藏在公式中。内存过小可能反复遗忘有用方向；内存更大也不保证每个问题都更快。

不形成矩阵的实现

函数输入梯度、按旧到新排列的曲率对以及正尺度，返回 $Hg$，不改写输入。只用标准库；向量维数为 $n$、保留 $\ell \le m$ 对时，时间和额外工作空间分别为 $O(mn)$ 与 $O(n+m)$，曲率历史占 $O(mn)$。下例测试的是方向子程序；一次完整优化还要计算新梯度、选择步长并更新队列。

Python 3

def lbfgs_product(gradient, history, gamma=1.0):
    g = [float(v) for v in gradient]
    pairs = list(history)
    if gamma <= 0 or not g:
        raise ValueError("positive scale and nonempty gradient required")
    dot = lambda a, b: sum(x*y for x, y in zip(a, b))
    rho = []
    for s, y in pairs:
        if len(s) != len(g) or len(y) != len(g) or dot(s, y) <= 0:
            raise ValueError("invalid curvature pair")
        rho.append(1.0/dot(s, y))
    q, alpha = g[:], [0.0]*len(pairs)
    for j in range(len(pairs)-1, -1, -1):
        s, y = pairs[j]
        alpha[j] = rho[j]*dot(s, q)
        q = [v-alpha[j]*w for v, w in zip(q, y)]
    r = [gamma*v for v in q]
    for j, (s, y) in enumerate(pairs):
        beta = rho[j]*dot(y, r)
        r = [v+(alpha[j]-beta)*w for v, w in zip(r, s)]
    return tuple(r)

h = [((1, 1), (1, 4))]
print(tuple(round(v, 6) for v in lbfgs_product((2, 1), h)))
print(lbfgs_product((1, 4), h))
assert lbfgs_product((2, 1), [], 0.5) == (1.0, 0.5)

输出为 (2.84, 0.04) 和 (1.0, 1.0)；空历史返回 $\gamma g$，符合基矩阵的定义。