L-BFGS算法解析：从原理到高维优化实战

做数值优化的人早晚会撞上同一堵墙参数维度一上去经典方法的底气就没了。你有目标函数有梯度甚至起点也合理可参数多到几十万乃至上百万时梯度下降一步一个脚印走到怀疑人生BFGS理论上漂亮但一个n乘n的近似矩阵就足以让内存直接爆炸。我这些年处理过的实际优化任务里最终能稳定跑起来的方案几乎都落在同一种算法上——L-BFGS也就是Limited-memory BFGS常被叫作有限内存BFGS算是在BFGS骨架上做了关键改动的省内存版本。这篇文章我想从一个手工实现开始讲起再聊到它的参数选择、常见坑和实际工程经验。适合谁看如果你正在做大规模参数估计、训练传统机器学习模型、搞物理反演或数值模拟中的参数标定或者只是被高维优化问题搞得头疼那这篇内容应该能帮上忙。我不打算堆一堆教科书公式而是尽量用“这算法到底在干什么”的视角拆给你看。1. 高维优化场景里经典算法为什么集体失灵1.1 你迟早会撞上的那一堵墙高维优化最典型的场景比如逻辑回归模型的参数训练、一个带几十万个参数的数值模拟系统校准、或者是某种反演问题里需要同时标定一大批未知系数。这时候你要做的事情本质上都一样找一个参数向量w让某个目标函数f(w)尽可能小。问题难在不同维度的尺度差异巨大、变量之间相互耦合目标函数曲面就像一条被拉伸得极长的山谷有的方向非常陡峭有的方向异常平缓。梯度下降在这种地形里的表现很糟糕。沿着陡峭方向稍微走大一点步长函数值就震荡甚至发散步长调小平缓方向的收敛又慢得让人绝望。我记得有次做一个参数标定任务优化了几千轮函数值还在缓慢下降每轮算梯度都要几分钟那感觉就是拿一个漏水的桶去浇一片菜地。理论上梯度下降的收敛速度取决于目标函数条件数。条件数越大收敛越慢。而在真实的高维模型里条件数差几个数量级太常见了。所以纯一阶方法在这种场景只能作为兜底的最后手段很少作为主力方案。1.2 牛顿法太贵BFGS又太重牛顿法之所以强大是因为它不只看梯度的方向还看曲率。它把目标函数局部近似成一个二次曲面直接一步跳到这个二次曲面的最低点附近。这个“曲率信息”藏在Hessian矩阵里也就是目标函数对参数所有二阶偏导组成的那张大矩阵。问题也出在这里。参数维度n如果是1万Hessian矩阵就是1万乘1万双精度浮点数存储需要约800MBn到10万直接是80GB要是到了百万量级比如常见的深度学习模型这个矩阵大到完全没法落地。就算能存下每一步要求解线性方程组计算量同样是灾难。BFGS算是牛顿法的一个折中它不直接计算Hessian而是利用每次迭代的梯度变化和参数变化逐步构造一个对Hessian逆矩阵的近似。好处是避免了显式求二阶导和矩阵求逆仍然保留曲率信息收敛速度远超梯度下降。但坏处是它仍然要维护一个n乘n的稠密矩阵。说白了BFGS在内存开销上和牛顿法是一路货色。所以在高维场景里你需要的是“有曲率感知能力但又不存储大矩阵”的算法。L-BFGS就是冲着这个需求来的。2. 省内存的核心机制从矩阵到向量对2.1 BFGS把问题藏在了哪里BFGS的核心迭代思路是维护一个n乘n的近似逆Hessian矩阵H_k每一步用当前梯度信息和历史信息去修正它。修正的方式是外积叠加相当于在矩阵上做低秩更新。你可以把它理解成一个不断记录“参数移动方向和梯度变化方向之间的关系”的黑板。BFGS问题在于它保存的是整个关系矩阵。一个n维问题所有维度两两之间的关系都要记录所以存储是O(n²)。这就像在一个班里BFGS要维护一张记录每两个同学之间熟悉程度的N乘N表格表格会随着班级人数成平方增长很快就把教室堆满了。那有没有可能不维护这张完整表格而是只记住最近几次活动里发生过什么当然可以。这就是有限记忆的思路。2.2 每一次迭代真正需要记住的就两件事每次迭代里算法实际上只用到两类信息参数从x_k走到x_{k1}的变化量以及梯度从g_k变到g_{k1}的变化量。这两个向量分别记为s_k和y_ks_k 参数位移y_k 梯度变化这一对向量被叫作曲率对。它物理上的含义是我沿着s方向走了一步梯度在s方向上发生了y的变化这个比例关系反映着目标函数在当前区域的弯曲程度。BFGS的近似公式本质上就是把这些曲率对按一定规则叠进那个n乘n矩阵里。L-BFGS的核心决定是既然要省内存那干脆不存那个巨大的近似矩阵只保留最近m次迭代的曲率对。当新的曲率对产生时旧的那对被丢弃就像一个滑动窗口。内存开销从O(n²)直接降到O(mn)。m一般取5到20所以实际内存开销和纯梯度下降处在同一个量级远小于BFGS。2.3 两循环递归是怎么“凭空”算出方向的省掉了矩阵之后问题变成每次迭代我还需要得到一个“近似Hessian逆矩阵乘以梯度”的方向向量没有矩阵怎么算这个东西答案就是一段经典的两循环递归算法。它的名字叫two-loop recursion是L-BFGS的灵魂。简单描述一下流程输入当前梯度先拿它从最近的曲率对开始依次往最旧的曲率对方向“倒推”每一步计算一个系数alpha并用alpha乘上梯度变化向量去修正当前向量。这个过程结束后用一个缩放系数gamma当作初始近似矩阵的倍数再让向量反过来从最旧的曲率对往最新的方向“正推”不断叠加参数位移向量的修正最终得到一个向量这个向量就是近似Hessian逆矩阵和梯度的乘积取负号就是搜索方向。代码很短逻辑确实有点绕。我一开始手写时也搞反过一次方向结果优化起来完全不对。一个比较好记的心法是“叠盘子”第一循环是从最上层盘子往下记住每层的压力第二循环是从最底层往上把这些压力再传回来。中间那个gamma相当于给整体压力乘一个比例。理解了这层直觉代码就不会写错。2.4 m参数到底在控制什么m是L-BFGS里唯一真正需要我们手动调的核心超参数含义是保留多少历史曲率对。它同时影响内存和收敛质量。m太小保留的信息太少近似矩阵质量差算法表现退化成接近梯度下降m太大内存上升更重要的是太旧的信息已经不能反映当前区域的曲率特征反而拖累收敛。经验上m取5到20之间都算合理。小型问题或者快速验证代码时取3到5一般的高维优化取10病态问题或对精度要求很高时取15到20。我通常先用10跑通再看收敛曲线效果不理想再调整。m对最终结果的影响没有终止容差那么敏感所以不必反复纠结。3. 动手实现一个轻量级L-BFGS3.1 完整实现不到一百行代码跑通为了真正看懂L-BFGS我建议你至少手写一遍核心逻辑。下面这个版本去掉了所有花哨的东西保留了两循环递归和最简单的回溯线搜索适合直接复制运行。import numpy as np def lbfgs(fun, x0, m10, max_iter1000, gtol1e-5): x x0.copy() n len(x) s_list [] y_list [] f, g fun(x) d -g for _ in range(max_iter): if np.linalg.norm(g) gtol: break # 第一循环从最新的曲率对往最旧的算 q g.copy() alphas [] for s, y in reversed(list(zip(s_list, y_list))): alpha (s q) / (y s) q q - alpha * y alphas.append(alpha) # 用最近一次曲率信息估计初始缩放因子 r q.copy() if len(s_list) 0: s, y s_list[-1], y_list[-1] gamma (s y) / (y y) r * gamma # 第二循环从最旧的曲率对往最新的算 for s, y, alpha in zip(s_list, y_list, reversed(alphas)): beta (y r) / (y s) r r (alpha - beta) * s d -r # 回溯线搜索Armijo条件 step 1.0 c1 1e-4 f_new, g_new None, None while step 1e-12: x_new x step * d f_new, g_new fun(x_new) if f_new f c1 * step * (g d): break step * 0.5 # 更新历史曲率对 s x_new - x y g_new - g if len(s_list) m: s_list.pop(0) y_list.pop(0) s_list.append(s) y_list.append(y) x x_new f, g f_new, g_new return x, f, g这里fun函数要求返回两个值函数值和梯度。比如用高维Rosenbrock函数来测梯度公式是现成的def rosenbrock(x): f np.sum(100 * (x[1:] - x[:-1] ** 2) ** 2 (1 - x[:-1]) ** 2) g np.zeros_like(x) g[:-1] -400 * x[:-1] * (x[1:] - x[:-1] ** 2) - 2 * (1 - x[:-1]) g[1:] 200 * (x[1:] - x[:-1] ** 2) return f, g x0 np.random.randn(5000) x, f, g lbfgs(rosenbrock, x0, m10, max_iter500) print(np.linalg.norm(x - 1.0), f)5000维的Rosenbrock函数L-BFGS在这种规模下完全没压力内存开销极小。我实测跑下来几百次迭代就能把梯度范数压到很低这在纯梯度下降下几乎不可能做到。3.2 用现成求解器快速落地手写版本主要用于理解原理工程上通常直接用成熟的求解器。最常用的是scipy.optimize里带的L-BFGS-B接口没错这个B指的是支持边界约束英文Bound constraints。使用方式非常简单from scipy.optimize import minimize res minimize( rosenbrock, x0, methodL-BFGS-B, jacTrue, options{ maxiter: 500, maxls: 30, ftol: 1e-12, gtol: 1e-6, }, )这里jacTrue表示我们提供的目标函数会同时返回梯度如果没有解析梯度也可以设置jacFalse让求解器用数值差分但这种情况慎用因为数值梯度本身引入的误差会影响L-BFGS的曲率估计。另外如果参数有取值范围可以传入bounds参数格式是每个维度的(min, max)列表某个维度不设边界就写None。3.3 判断收敛的三个信号用L-BFGS跑任务光看return的res对象不够还得知道什么状态才算正常收敛。我一般看三个信号第一个是梯度范数。L-BFGS终止条件里最核心的是梯度范数低于容差这在res.message里会显示成类似“CONVERGENCE: NORM_OF_PROJECTED_GRADIENTPGTOL”的提示。看到这个消息基本可以认为已经收敛到局部稳定点。第二个是函数值下降曲线。正常情况前几轮函数值应该大幅下降然后逐渐变缓。如果函数值长时间不降即使梯度范数达到容差也要怀疑是不是卡在平台或鞍点。第三个是迭代轮数。如果直接达到maxiter而没触发收敛条件说明问题比预想难或者参数需要调整。这时候我会把历史信息打出来看是线搜索频繁失败还是梯度范数始终降不下去两者对应的处理方式完全不同。4. 参数调优和实操经验4.1 历史长度m到底应该设多大这个参数在上一节提过这里展开说说实际调法。我按照问题规模和处理需求给一个经验表场景m建议快速验证代码或小规模问题3到5一般高维优化10万到100万维8到15病态问题、高精度需求15到25内存极其紧张且问题结构简单3到5m调大的直接变化是每次迭代多存几个向量内存只线性增长影响不大。真正需要关注的是历史信息是否还有参考价值。如果目标函数区域在优化过程中变化剧烈太旧的曲率对反而像过期地图会把方向带偏。所以在非凸问题里m不是越大越好。我有一次处理一个强烈非线性反演问题默认m10时优化始终在某个局部区域震荡后来把m降到4反而收敛快了。原因就是目标函数的二阶特性在全空间变化太大保留太多旧信息只会增加方向估计的偏差。这也说明参数调优不是死守默认值而是要观察收敛曲线的实际反馈。4.2 终止容差怎么定才不“过早收敛”L-BFGS的终止条件主要有两个函数值变化容差ftol和梯度范数容差gtol。默认值在不同库里有差异但经常是ftol1e-9到1e-12gtol1e-5到1e-6。实际项目里我倾向于把gtol设得比ftol更严格一些因为函数值下降变缓时梯度可能还没到最小点。特别提醒一点梯度范数接近零时并不代表你找到了全局最优。它可能是一个鞍点或者一个很平坦的局部区域。我之前做过一个参数联合标定问题目标函数在某个平台上的梯度范数确实很小但函数值本身离期望结果还差得远。后来我发现问题是目标函数设计不当。解决办法是不要只看最后一个收敛标志要把优化过程中的函数值下降曲线存下来结合业务指标判断当前结果是否合理。4.3 线搜索参数和边界约束maxls参数控制每次迭代中最大线搜索次数默认值一般是20到30。如果线搜索次数被耗尽求解器会报“线搜索无法在maxls步内充分减小函数值”。这个错误在高维问题里很常见尤其当目标函数在某些区域特别陡峭或者初始点选得很差时。遇到这个情况我的第一反应不是盲目调大maxls而是先怀疑梯度和函数值的一致性。用手写的函数很容易出现梯度符号错误或者漏项这会让算法误判方向线搜索自然一路失败。确认梯度没问题后再把maxls调到50同时减小初始步长通常就能解决。边界约束方面初始点一定要满足bounds范围。如果x0里某个分量越界求解器会直接报错。另外注意L-BFGS-B处理边界约束是靠投影和“有效集”机制它的收敛判断也会基于投影梯度因此在边界附近的收敛行为和内部点略有区别这种场景下认真读返回message里的收敛类型很有必要。5. 常见问题与排查技巧5.1 梯度里有噪声时L-BFGS最先崩L-BFGS对梯度质量的敏感度比多数一阶方法高得多。它的核心是靠曲率对推断Hessian近似结构如果梯度本身带噪声得到的y向量就有污染近似矩阵慢慢就被“教坏”了。我在机器学习任务里踩过这个坑用mini-batch估计的梯度去喂L-BFGS跑几轮就出现函数值不降反升或者梯度的范数反复震荡。排查思路很直接先用单纯的梯度下降配精确梯度跑一小段确认目标函数、梯度实现本身没问题再把L-BFGS放上去。如果纯梯度下降能稳定下降而L-BFGS不行基本可以断定是梯度噪声问题。这也就是为什么L-BFGS非常适合全批量小模型训练但很少被用来训大型深度网络——后者的梯度本身就是小批量估计出来的噪声避不开。5.2 NaN、除零和线搜索失败的真相L-BFGS在高维场景里的NaN问题通常有三个来源。一个是目标函数计算溢出中间结果出现inf梯度也跟着变成NaN。这种情况要在目标函数内部加保护把数值钳制到合理范围。另一个是曲率对的内积ys等于零这意味着两次迭代的梯度完全相同算法无法构造有效的曲率信息需要加一个极小值保护或者重置历史。还有一种可能线搜索步长踩到无界区域目标函数直接发散。排查时我习惯在目标函数里临时插入一段检查如果f或者g不是有限数值就把当前x、目前迭代轮数、以及历史曲率对的统计信息打印出来。这样基本能定位是哪个环节出了问题。大部分情况下都指向初始点离开合理范围太远让算法在早期就走上爆炸路径。5.3 高维场景内存开销实测对比很多人问L-BFGS到底省了多少内存这里给一组直观数字。以下按float64计算BFGS需要存一个n乘n近似矩阵L-BFGS按m10存20个长度为n的向量参数量nBFGS显式矩阵内存L-BFGS(m10)内存1万800 MB1.6 MB10万80 GB16 MB100万800 TB不可行160 MB这个差距是数量级的。在n为100万的问题里BFGS连矩阵都建不出来L-BFGS加上方向向量临时变量总体占用也就几百MB普通服务器完全能跑。所以“省内存”不是锦上添花而是让某些高维优化问题从“不能做”变成“可以做”的关键。5.4 非凸问题怎么办L-BFGS本质还是局部优化器它没有能力在连绵山脉里找到你想要的全局最低点。如果目标函数非凸直接跑L-BFGS很容易陷入一个局部极小值而且你浑然不觉因为它照样会输出一个“收敛成功”的标志。我的做法是把它和多起点策略配在一起用。第一步先做粗粒度采样或者网格搜索选出几个函数值明显低的初始点第二步在每个候选点上跑L-BFGS进行局部精调第三步比较各局部最优点的函数值和实际业务指标选最终结果。这个方法不能保证找到全局最优但它把L-BFGS的长处发挥得最充分而且实现成本很低。6. 多年实战沉淀下来的几个习惯6.1 小规模复现永远是最高优先级的自检接手任何新的优化问题我第一件事永远是把它缩到很小规模去复现。参数维度降到几十边界条件放宽跑几个L-BFGS迭代把每一步的x、梯度、函数值打印出来对照。这一步能揪出大量低级错误比如梯度索引错位、目标函数和梯度不匹配、初始点范围不合理等等。我有过好几次直接在大规模问题里开跑的教训结果梯度函数里一个符号写反白算了几个小时最后才发现问题在早期就用小规模复现一眼就能看出来。现在我把“小规模复现”当作和写代码同步进行的习惯而不是出了问题再回头查。6.2 梯度检查要写进流程无论梯度是手动推导还是用自动微分生成的我都建议至少做一次数值梯度检查。常用公式是中心差分g_i ≈ (f(x h e_i) - f(x - h e_i)) / (2h)h取太小会受浮点误差影响取太大会让差分近似偏差变大。一般取h sqrt(machine_epsilon) * max(1, |x_i|)然后检查解析梯度和数值梯度的相对误差。如果某些维度相对误差异常大多半是解析梯度公式有误。这个检查不用每次跑优化都做但每次新写目标函数都要过一遍。尤其当你准备用L-BFGS这种靠梯度构造曲率信息的算法时梯度的准确性直接影响收敛质量甚至决定整个优化能不能跑下去。6.3 给算法一个“干净”的起点L-BFGS对初始点不是完全无感一个离谱的初始点会白白消耗大量早期迭代。我通常会在跑之前做三件小事把输入特征的尺度标准化到相近范围把参数初始化到合理先验附近的区域而不是随机撒点先跑一小段短迭代看损失下降趋势如果前面几步函数值不降立刻停掉排查原因而不是麻木地等它跑完几千轮。最直观的体验是一个干净起点往往能让L-BFGS在十几轮内降到接近最优值附近剩下的事只是慢慢精调。而一个糟糕起点可能跑几百轮还在“迷路”两者的实际成本差别非常大。我这些年做优化最深的体会是L-BFGS不是万能的但它是一个兜底能力极强、在绝大多数高维连续光滑问题上都值得先试的工具。它有成熟实现有清晰的内存优势收敛速度又远胜梯度下降。只要你把梯度的准确性和参数的尺度处理好它能帮你省掉大量调参和等待的时间。希望这篇从原理到踩坑的梳理能帮你把它顺顺当当地用起来。