L-BFGS算法精讲：高维优化中的省内存拟牛顿法

做优化的人，大概都经历过这种尴尬：目标函数公式不复杂，梯度也能算，可变量一上万，常规方法就顶不住了。这个时候，L-BFGS几乎是最顺手的工具。它的全称是Limited-memory BFGS，翻译过来就是“省内存的BFGS”，核心卖点只有一个：应对高维多变量优化时，用很少的内存开销，拿到接近二阶方法的收敛速度。这篇文章不打算讲教科书推导，而是从原理到能跑的代码，再到我实际调试中踩过的坑，把L-BFGS怎么用于高维问题讲清楚。

1. 高维优化为什么难：从“装不下”说起

1.1 维度一上去，直觉就失效了

当变量只有几十个的时候，很多经典方法都够用：共轭梯度、拟牛顿、甚至最老实的梯度下降，配合仔细调参都能收敛。可一旦变量个数涨到十万、百万，情况完全变了。最直观的问题就是“维度灾难”：如果用网格搜索找最优解，假设每个方向只取两个点，那一百万个变量的网格就有 (2^{1000000}) 个组合，这个数字比宇宙中基本粒子数量还大得多，根本不可能遍历。

就算不用网格搜索，梯度下降在病态问题里也极其痛苦。什么叫病态？简单说就是目标函数的等高线在某个方向收得很窄、在另一个方向拉得极长，像一个压扁的椭圆。梯度下降在这种地形上会来回震荡，明明离最优解不远，却走出一条“之”字形，半天不前进。解决办法是用二阶信息，也就是Hessian矩阵（目标函数对所有变量的二阶偏导），它告诉我们每个方向上的曲率。可惜曲率也是有代价的：一个 (n \times n) 的Hessian矩阵，当 (n=100) 时是 (100 \times 100)，没问题；当 (n=10^6) 时需要存储约 (8 \times 10^{12}) 字节，也就是8TB。这还只是“存下来”，还没算矩阵求逆的代价。普通工作站上根本不可能。

所以在一段时间里，工程界急需一种“既要二阶收敛速度、又不能背太重行李”的方法。拟牛顿法就是这个方向的答案。

1.2 BFGS是好东西，但它有个改不掉的毛病

BFGS是拟牛顿法里最经典的一种，名字来自四位提出者的姓氏首字母。它的核心思想很聪明：不直接计算二阶导，而是通过相邻两步的梯度差来“猜测”曲率信息。每次迭代时，它用一个近似矩阵 (B_k) 模拟Hessian，更新公式只需要两个向量：

[ B_{k+1} = B_k - \frac{B_k s s^T B_k}{s^T B_k s} + \frac{y y^T}{y^T s} ]

其中 (s = x_{k+1} - x_k) 是迭代点的变化，(y = g_{k+1} - g_k) 是梯度的变化。从计算角度看，这个更新非常优雅：完全避开了二阶导，只需要目标函数值和梯度。只要线搜索满足条件，这个矩阵还天然保持正定，方向质量很有保证。

但BFGS也有一个铁打的软肋：不管公式怎么优化，它始终要维护一个 (n \times n) 的近似矩阵。 (n=10^5) 时矩阵是80GB，(n=10^6) 时是8TB。也就是说，BFGS解决了一部分计算二阶导的麻烦，但没解决内存问题。很多年前我在一个变量规模超过十万的参数反演任务上试过BFGS，跑了一小会儿内存直接吃满，机器开始疯狂交换内存，最后只能换思路。这让我切身明白了“BFGS只适合中小规模问题”这件事。

2. 省内存的秘密：L-BFGS到底在记什么

2.1 不存整张矩阵，只留最近几步的历史

L-BFGS的思路，用一句俗话概括就是：与其把整本笔记背在身上，不如只带最近几页。BFGS的更新过程，本质上是在用最近发生的“位移 (s)”和“梯度变化 (y)”不断修正对曲率的估计。也就是说，真正有用的信息并不是那个巨大的矩阵本身，而是最近的若干组 ((s, y)) 向量对。

L-BFGS就只保存最近 (m) 组这样的向量对，通常 (m) 取5到20。当第 (m+1) 组出现时，最早的那组直接丢弃。这样一来，存储量从BFGS的 (O(n^2)) 降到了 (O(mn))。我们算一笔账：当变量数 (n=10^6)、(m=10) 时，需要保存20个长度为一百万的向量，每个双精度数8字节，总共约160MB。这个内存量在常规服务器上很轻松，和8TB完全不是一个量级。

关键问题是：只有 ((s, y)) 历史，怎么算搜索方向？BFGS每次迭代都要用矩阵 (B_k) 乘梯度得到方向，现在矩阵没了，但我们有“历史”。由于BFGS的目标函数近似矩阵更新本身是一个递推过程，最近 (m) 步的“修补”操作其实已经包含了矩阵的主要信息，越久远的更新对当前曲率的影响越小，所以丢弃旧历史是合理的近似。这就是L-BFGS的理论根基。

2.2 双循环递归：整个算法的心脏

有了 ((s, y)) 历史之后，L-BFGS用一个非常漂亮的技巧计算搜索方向，叫做“双循环递归”。它做的事情，相当于把BFGS的矩阵更新公式“展开”到历史上，不显式构造 (B_k)，而是直接算出 (B_k^{-1} \nabla f(x))（即近似逆Hessian乘梯度）的结果。

我直接给出双循环的Python示意，这是整个L-BFGS最核心的部分：

def two_loop_recursion(grad, S, Y, RHO, H0_diag): """ 双循环递归：用历史 (s_i, y_i) 计算 H_k * grad grad: 当前梯度 S, Y: 历史向量列表，S[i]=s_i, Y[i]=y_i RHO: 每个历史对的 1/(s_i^T y_i) H0_diag: 初始对角矩阵的标量（通常每轮用缩放因子重新估计） """ q = grad.copy() alpha_list = [] # 第一个循环：从最新历史往最旧历史扫描 for i in range(len(S) - 1, -1, -1): alpha = RHO[i] * np.dot(S[i], q) alpha_list.append(alpha) q = q - alpha * Y[i] # 让 alpha 顺序从旧到新，方便第二个循环对号入座 alpha_list.reverse() # 初始对角矩阵作用：通常 H0 = gamma * I r = H0_diag * q # 第二个循环：从最旧历史往最新历史扫描 for i in range(len(S)): beta = RHO[i] * np.dot(Y[i], r) r = r + S[i] * (alpha_list[i] - beta) return r # 这是 H_k * grad，真正的方向要取负号

这个双循环看着有点绕，但拆开看其实很清晰。第一个循环从新到旧，把梯度 (q) 逐个投影到每条历史曲率方向上，得到一组系数 (\alpha)，同时用这些系数把 (q) 里的对应成分“剥掉”。第二个循环反过来从旧到新，用这些系数和重构的曲率信息把 (q) 一点点“补回来”。两个循环一正一反，刚好模拟出BFGS矩阵对梯度的作用效果。复杂度是 (O(mn))，比BFGS的 (O(n^2)) 小了一个量级，这就是L-BFGS又能快又能省的核心原因。

2.3 H0的缩放：不是玄学，是稳定性的关键

双循环里有个容易忽略但极其重要的细节：初始矩阵 (H_0)。在第一个循环剥完梯度、准备“补回”之前，需要先让 (q) 过一个初始矩阵。很多实现图省事直接取单位阵 (I)，也就是不做任何缩放。这在目标函数尺度比较正常时能用，但如果变量的量纲差异很大——比如有些变量是 (10^{-6}) 量级、有些是 (10^6) 量级，直接上单位阵会导致第一步方向长度完全错乱，线搜索被迫把步长压到极小，收敛慢得让人怀疑人生。

靠谱的做法是给 (H_0) 乘一个自动缩放的标量 (\gamma)，最常用的估计是最近一对历史向量的曲率比值：

[ \gamma = \frac{s^T y}{y^T y} ]

这个式子的直观理解是：在最近走出的一小段弧线上，用梯度的变化去拟合位移的变化，得到的比例近似表示当前地形沿搜索方向的“弧度”。在双循环递归中，将 (H_0 = \gamma I) 作为隐含的初始矩阵，能让方向长度与真实尺度匹配。我自己的习惯是，每轮迭代都重新算一次 (\gamma)，只用最新的一组 ((s, y))，效果非常稳。

3. 从零实现：能跑的双循环与线搜索

3.1 整体迭代骨架，先搭框架

完整的L-BFGS程序比双循环多两部分：线搜索和历史管理。整个迭代循环的骨架是这样的：

1. 给定初始点 (x_0)，计算目标函数值 (f(x_0)) 和梯度 (g_0)。
2. 用双循环递归计算搜索方向 (d_k = -H_k g_k)。
3. 做线搜索，找到一个满足条件的步长 (\alpha_k)，得到新点 (x_{k+1})。
4. 计算新的梯度，得到 (s = x_{k+1} - x_k) 和 (y = g_{k+1} - g_k)。
5. 更新历史队列：追加新一组 ((s, y))，如果队列长度超过 (m)，弹出最旧的一组。
6. 重复直到梯度范数小于阈值或达到最大迭代次数。

这里的历史队列不需要任何奇技淫巧，用列表就行，满了就把最前面的元素删掉。唯一要注意的是，更新历史前必须检查曲率条件，这点我在下一节单独说。

3.2 线搜索：光有下降还不够

很多初学者会直接用最简单的Armijo回溯线搜索，也就是步长从1开始，不断减半，直到目标函数满足充分下降条件：

[ f(x + \alpha d) \le f(x) + c_1 \alpha \nabla f(x)^T d ]

这个条件本身没错，但只满足它对于L-BFGS来说不够。L-BFGS的 ((s, y)) 历史要起作用，必须保证每一步的新点梯度满足曲率条件：(y^T s > 0)。如果违反这个条件，更新后的近似矩阵会丢失正定性，后面的方向质量会急剧恶化。Wolfe条件正是在Armijo条件基础上，额外要求新点的方向导数不能“反向得太多”：

[ |\nabla f(x + \alpha d)^T d| \le c_2 |\nabla f(x)^T d| ]

参数通常取 (c_1 = 10^{-4})、(c_2 = 0.9)。满足Wolfe条件的线搜索，才能保证每次更新都是“健康”的。完整且强壮的Wolfe线搜索实现通常用插值法，比简单回溯要复杂，但工程上很值得投资。如果只是快速验证算法，也可以用回溯法配合曲率检查来做近似，但心里要清楚这丢了理论保证。

实际调试中我发现，线搜索失败是L-BFGS跑崩的首要来源。很多情况下不是算法本身有问题，而是线搜索没写好，或者初始步长给得太大。我的建议是初始步长直接用1，交给线搜索去调整，这样符合拟牛顿法的常规假设——在接近最优点时，单位步长一般能直接接受，省掉大量函数评估。

3.3 更新历史的边界情况：曲率条件必须查

在迭代过程中，并不是每一组 ((s, y)) 都值得放进历史。如果新点根本不满足 (y^T s > 0)，说明这段路径上的曲率信息是“坏的”，强行放进队列只会污染后续方向。很多老练的实现会加一道保护，只有满足下面条件才更新历史：

curvature = dot(s, y) if curvature > 1e-12 * norm(s) * norm(y): 追加并管理队列 else: 跳过本次更新，历史保持不变

这个判断的本质是防止数值噪声带来的负曲率更新。我见过一个案例，目标函数里有一个轻微的数值污染，导致某几次迭代的 (y^T s) 是微小的负数，起初不以为意，结果十几轮迭代后方向彻底失控，步长开始指数级收缩，最后整个优化直接失败。加上曲率保护后，问题再也没出现过。

3.4 内存代价的真实对比：一张表说明问题

为了直观感受“省内存”三个字，下面这张表列出了不同变量规模下BFGS和L-BFGS（取 (m=10)）存储近似矩阵所需的粗略内存（按双精度8字节计算）：

这个对比非常清楚：在十万到百万量级的变量规模下，BFGS根本不是“省不省内存”的问题，而是物理上无法落地。L-BFGS则保持了内存始终可控，这也是它在高维多变量优化任务里成了事实标配的原因。

4. 参数调优实战：m、H0、线搜索的取值

4.1 m选多大：默认10，但别死守默认

记忆长度 (m) 是L-BFGS最直接的超参数。取小一点，内存和单轮计算量都小，但可用的曲率信息少，遇到病态问题时收敛会变慢；取大一点，历史信息充足，但每轮迭代的 (O(mn)) 成本跟着涨，而且旧历史不一定是好消息——它们记录的是老地形的曲率，可能跟当前区域没啥关系，甚至会引入干扰。

我在实际项目里的经验是：从 (m=10) 开始，先跑几十到几百轮，观察收敛曲线。如果目标函数值下降得很顺利，就保持不动；如果感觉方向质量不够、像是在盲目试探，把 (m) 提到20再看。绝大多数问题在 (m=10\sim20) 之间就解决得很好了，没有必要把 (m) 拉到50以上。真正需要 (m>50) 的问题，多半是目标函数极度病态，我会先回去怀疑变量归一化和梯度计算是不是出了问题，而不是加记忆长度硬扛。

4.2 H0缩放问题：第一轮迭代就暴露真相

如果说 (m) 是细调参数，那么 (H_0) 的缩放就是粗调参数，它的影响经常在第一轮迭代就暴露出来。如果第一轮的初始步长1在线搜索里被连续减半好几次——比如减到 (10^{-4}) 级别才被接受——大概率不是线搜索的问题，而是方向长度严重不匹配。这时候检查一下 (\gamma) 的计算：是不是分母 (y^T y) 太小了？是不是忘了每轮重新估计 (\gamma)？

我在复现别人代码时见过一个很典型的错误：有人实现了L-BFGS的双循环，但把 (H_0) 写死成单位阵，还在目标函数值下降到一定程度之前都不缩放。结果就是前几十轮迭代的步长一直在 (10^{-5}) 附近徘徊，优化慢得像在爬。加上动态 (\gamma) 之后，同一个问题二十轮就收敛了。所以遇到收敛慢，先别怀疑算法，先把第一轮迭代的“方向范数”和“接受步长”打印出来看一眼，这个习惯能省下很多调试时间。

4.3 线搜索参数取值的常见区间

Wolfe条件里的两个参数 (c_1)、(c_2) 在准牛顿框架下的取值选择相对固定。(c_1) 控制Armijo条件的“充分下降”程度，取太小基本没有约束力，取太大会要求每步都大幅下降，反而限制步长，一般用 (10^{-4}) 比较稳。(c_2) 控制曲率条件的严格程度，L-BFGS这类准牛顿方法一般取 (0.9)，意思是允许方向导数有较大的变化幅度；如果取 (0.1)，线搜索会非常严格，每步都要找精确的最优步长，函数评估次数会明显增多，但迭代步数可能会少一点。在目标是“总运行时间最短”的工程场景里，宽松一点（(c_2=0.9)）通常更合适。

4.4 梯度精度会直接毁掉L-BFGS

这一点我要单独强调：L-BFGS比梯度下降更依赖梯度的精确性。因为算法从“相邻两次梯度的差值 (y)”中提取曲率信息，如果梯度本身有比较明显的噪声，那么提取出来的曲率基本就是噪声，后面的方向计算全部失去意义。

实际操作中，如果目标函数的梯度是用有限差分近似出来的，差分步长要仔细选；如果目标函数内部本身有随机性（比如深度学习里的小批量采样），L-BFGS直接硬上会非常痛苦，因为每次评估的梯度都在跳变，历史里的 (y) 根本无法反映真实的曲率。这种场景下应该考虑随机拟牛顿方法，或者干脆用适合随机优化的算法。记住一个判断标准：如果两次相邻迭代的梯度之差和当前梯度本身在同一个量级，那L-BFGS大概率会失控。

5. 踩坑实录：收敛失败与NaN排查

5.1 一迭代就NaN，先检查线搜索初始步长

L-BFGS跑出NaN，最常见的元凶不是算法本身，而是第一步线搜索的初始步长太大，直接跨出了目标函数的合法定义域。比如目标函数里有 (\log)、开根号、除法这类操作，步长过大就会让中间变量变成负数或零，一算就炸。

排查思路是这样的：先关掉所有历史更新，只做一次梯度下降式的迭代，看看步长1是否会触发NaN。如果会，把初始步长改成更小的值（比如 (10^{-3})），或者检查目标函数内部有没有需要在数学域上进行约束的地方。另一个原因是目标函数或者梯度的实现本身有bug，这种情况下任何优化算法都没辙，只能先在几个随机的初始点上做梯度检查（数值梯度与分析梯度对比），确认没问题再上L-BFGS。

5.2 线搜索步长不断收缩：方向出问题了

另一种常见情况是线搜索启动了，但每次的初始步长1都被反复减半，减到 (10^{-8}) 才勉强接受，导致每轮迭代只在原地挪一小步。这种问题几乎可以断定是搜索方向质量太差。什么事情会让方向质量差？第一，梯度计算错误，方向里带着大量噪声；第二，历史队列被“坏曲率”污染了，方向越来越离谱；第三，目标函数存在不可微或非光滑区域，某些方向上的下降完全是假的。

我调试的时候会打印两个量：一是方向与负梯度方向夹角的余弦值，二是线搜索最终接受步长。如果余弦值长期低于0.5，说明方向和梯度下降几乎垂直，这种“优化”基本是在兜圈子。这时候先怀疑梯度和历史队列，别去调线搜索参数。

5.3 收敛慢得令人怀疑：可能只是没归一化

有个印象很深的例子：某个参数反演任务里，变量横跨了好几个数量级，有的接近 (10^{-6})，有的接近 (10^6)。L-BFGS跑起来之后，目标函数值虽然一直在降，但降得非常慢，几十轮迭代都看不到数量级改善。后来我把变量做了一个简单的归一化变换，让每个变量都落在差不多的尺度上，同一个问题十几轮就收敛了。

这件事的教训是：L-BFGS虽然比梯度下降更抗病态问题，但它并不是万能的。变量尺度差异过大时，搜索方向会因为不同方向的曲率差异太大而失衡。给变量做归一化或者预处理，往往比调任何参数都管用。如果你的目标函数可以分解成“线性变换+简单形式”，先做变量变换再优化，几乎不会有坏效果。

5.4 结果对初始点异常敏感：这是局部算法的宿命

L-BFGS本质上是局部优化算法，它只知道从当前点开始，根据地形的曲率估计找出附近最陡的方向。目标函数如果有多个局部极小点，不同初始点当然会收敛到不同位置。这不是算法缺陷，是问题本身的非凸性。

应对方法没有捷径，就是多起点跑。我一般会随机生成几十个初始点，分别用L-BFGS运行一遍，然后看最终目标函数值的分布。如果所有结果都集中在相近的值，说明找到的可能是全局最优或至少是满意的局部最优；如果结果差异很大，那就需要调整初始点生成策略或者换一种全局优化方法了。

5.5 怎么判断“真的收敛了”

很多人只看目标函数值不再下降就宣布收敛，但这是不可靠的。更稳妥的做法是同时看三个指标：梯度无穷范数是否降到阈值以下、目标函数的相对变化是否足够小、迭代步长是否也趋于零。如果目标函数曲线已经很平，但梯度范数还是很大，说明当前位置大概率不是极值点，只是线搜索卡住了或者地形太崎岖。这种情况下建议重新检查梯度和线搜索实现，不要强行终止。

6. 场景边界：L-BFGS能用在哪，别硬用在哪

6.1 高维参数估计与反问题里最能体现优势

L-BFGS最舒服的应用场景，是目标函数光滑、梯度可算、变量规模巨大的确定性优化问题。比如图像复原、全变分去噪、地震反演、气象同化这类“反问题”，一次要解的变量数动辄几十万甚至上百万，L-BFGS几乎成了标准配置。我用过一个模拟的图像复原任务，把512x512图像的每个像素当成变量，变量数超过26万，用 (m=10) 的L-BFGS跑起来，内存开销只有约40MB，几十轮迭代就能看到清晰的复原效果。如果换成BFGS尝试维护 (26万 \times 26万) 的矩阵，直接没法算。

另一个典型场景是大规模机器学习中的线性模型训练，比如逻辑回归、线性支持向量机。这类问题的目标函数通常是凸的，吃满L-BFGS的所有优点：全局最优、快速收敛、内存可控。需要注意的是，如果数据量极大且训练方式切换成了随机小批量，L-BFGS就不太合适了，因为梯度噪声会破坏 ((s, y)) 的曲率估计。

6.2 哪些情况别碰L-BFGS

L-BFGS不是万能钥匙。第一，目标函数不可微或者有大量离散变量时，别用。比如决策树里的分裂点选择、组合优化问题，L-BFGS缺乏梯度信息，根本无法启动。第二，函数评估里有随机噪声时，别用。前面反复强调过，L-BFGS的整套机制建立在对梯度的精确估计上，噪声一旦盖过曲率信号，算法会像一个喝了酒的人走路，方向越来越乱。第三，目标函数有大量局部极小、你又需要全局最优解的时候，别指望L-BFGS。它可以作为“局部精修”的一环，但全局探索应该交给其他全局优化方法。

6.3 值得留意的衍生版本

如果你实际用手头的开源数值库时发现L-BFGS遇到“变量边界限制”的问题——比如某些参数必须在0到1之间——会看到L-BFGS-B这个变体。它在L-BFGS框架上加了边界约束处理，把无约束问题变成带上下限的问题，内存代价基本不变。如果目标是稀疏模型训练里常见的L1正则，还有一个叫OWL-QN的变体，专门处理不可导的L1惩罚项。这些变体的核心思想仍然是“省内存的拟牛顿”，只是针对特殊目标函数做了相应的修正。

我在实际项目里的体会是：L-BFGS很难在纸面上显得惊艳，但它在高维多变量优化里的稳定性和通用性极其出色。调试时我有一个习惯，先在20维的小问题上把方向、梯度和线搜索全部验证一遍，再放到上百万维的实际任务里；第一轮迭代结束后就打印本轮接受步长和方向夹角，一旦异常，先怀疑梯度计算，而不是怀疑算法本身。这套流程帮我省下过太多在参数调优上瞎折腾的时间。说到底，L-BFGS的设计哲学用一个词概括就是“放下”——放下完整的Hessian，放下整张近似矩阵，只留下最近几步的曲率记忆。这个取舍，让它在高维多变量优化里始终是我最信任的默认选择。