BFGS这个系列我写到第四篇了。前三篇基本上把数学背景、算法推导和最简单的numpy实现都过了一遍很多读者也表示代码能跑起来了。但说句实在话能跑起来的demo和能在项目里真正扛事的优化器是两种完全不同的东西——这是我这几年来回调优反复确认的一件事。这一篇我不打算再推公式而是把手中的BFGS实现做成一套健壮、可复用的Python类支持解析梯度与数值梯度切换、内置非精确线搜索、带完整迭代历史记录并且能和scipy的官方实现做对照。无论你是正在学数值优化的学生还是想把拟牛顿法用在科学计算和机器学习调参里的工程师这篇都值得你完整跟一遍。1. 为什么还要写第四篇从“能跑”到“好用”1.1 前三篇我们走到了哪里先简单对上篇内容做个锚定。第一篇我们聊了牛顿法和Hessian矩阵的关系解释了为什么在高维、非凸问题里直接求解二阶信息既昂贵又不靠谱由此拉出拟牛顿思想的动机。第二篇用最朴素的方式实现了一个BFGS框架和梯度下降做了对比结论很好猜在线搜索配合下BFGS的收敛步数往往远少于固定学习率的梯度下降。第三篇重点讨论了步长选择引入了Armijo条件和Wolfe条件说明线搜索不是可选项而是拟牛顿法收敛性的保证。到这一步公式和框架都不缺了但很多读者给自己跑出来的结果并不稳定。换个初始点就发散或者迭代几步之后函数值开始振荡日志里出现负数步长Hessian近似矩阵越来越奇怪。这些都是我在初学阶段踩过的坑。问题不在公式而在工程细节。1.2 这一篇要解决什么问题第四篇的目标非常明确交付一个真正能在实际项目里使用的BFGS优化器。围绕这个目标我给自己列了一个问题清单这也会是这篇博文的内容主线步长应该怎么定才能既不震荡又不过早停下梯度给错了有没有一套低成本的检查方法中间量y_k^T s_k小于等于零说明什么要不要中断停机条件到底看函数值还是看梯度范数问题维度上升到几千几万时内存怎么处理这些问题在教科书里通常不会专门解释但它们才是工程落地时真正决定成败的地方。后面的内容就围绕这些问题展开。2. 先把工程问题想清楚再动手2.1 五个关键量在代码里对应什么BFGS的数学核心只有五个量搜索方向d_k、步长alpha_k、位移s_k、梯度差y_k以及逆Hessian近似矩阵H_k。代码里它们的关系非常直接搜索方向d -H g这等于把“当前最陡下降方向”在二阶近似下做了扭曲让更新方向更贴近牛顿方向位移s x_new - x表示实际走出的那一段路梯度差y g_new - g表示沿这段路梯度变化了多少步长由线搜索决定我后面会说两种实现方式逆Hessian近似每次迭代用s和y做一次秩二更新公式是H_new (I - rho*s*y^T) H (I - rho*y*s^T) rho*s*s^T。这五个量之间有严格的维度配合s和y是一维向量rho 1/(y^T s)是标量矩阵更新里用到的np.outer(s, y)得到的是n x n外积矩阵。写代码的时候先把这五个量的维度在注释里标出来能省下后面大量Debug时间。2.2 解析梯度还是数值梯度很多初学者会在这里犹豫到底省事点用数值微分还是老老实实手推解析表达式我的建议是能给出解析梯度就尽量给出数值梯度只做兜底方案。解析梯度的优势是精度高、计算量小缺点是要动脑筋而且很容易在推导或代码里引入符号错误。数值梯度用中心差分实现起来很简单公式是def numerical_gradient(f, x, eps1e-6): g np.zeros_like(x) for i in range(len(x)): ei np.zeros_like(x) ei[i] 1.0 g[i] (f(x eps * ei) - f(x - eps * ei)) / (2.0 * eps) return g中心差分的截断误差是O(eps^2)看起来很不错但浮点舍入误差会随着eps减小而变大。实测下来eps1e-6是一个安全区间太小反而会得到一堆噪声梯度。另外数值梯度每一次求值都要调用2n次目标函数维度一高成本急剧上升。所以我建议的架构是优化器接受gradient参数用户传入解析梯度函数如果不传内部自动走数值梯度兜底同时在构造时提供一个梯度检查工具。def check_gradient(f, grad, x, eps1e-6): numerical numerical_gradient(f, x, eps) analytical grad(x) denom max(1.0, np.linalg.norm(analytical)) return np.linalg.norm(numerical - analytical) / denom相对误差小于1e-5基本可以认为解析梯度是可靠的。我在测试Rosenbrock函数之前必做这一步别嫌麻烦它救过我好多次。2.3 初始点x0和初始矩阵H0怎么给x0的选择决定了目标函数在你手中展现的是哪一面。同一个函数从凸区域和从鞍点附近出发路径和收敛速度可能差一个数量级。对于测试算法我一般选择让问题暴露特性的点比如Rosenbrock函数经典的(-1.2, 1.0)它把函数碗状的谷底完美藏了起来能真实考验算法的“翻山越岭”能力。H0则简单得多默认取单位矩阵I。如果对问题有额外信息也可以用H0 (y0^T y0)/(y0^T s0) * I这种经验公式它来自Barzilai-Borwein方法的启发在部分问题上能加速初始阶段的收敛。不过工程上我很少依赖它单位矩阵足够稳定而且不引入额外调参负担。3. 核心实现一个可直接复用的BFGS类3.1 类的整体骨架怎么设计我倾向于把优化器封装成一个BFGSOptimizer类而不是一堆散装函数。理由很实际BFGS的运行涉及梯度计算函数、线搜索参数、迭代历史等多个状态类能把它们组织得更清晰也方便后续扩展成L-BFGS或者做批量测试。类的核心参数是这四个目标函数objective、梯度函数gradient可选、停机阈值tol、最大迭代次数max_iter。此外我还会留两个线搜索相关的参数c1和rho它们分别控制Armijo条件的接受斜率和回溯衰减系数通常不需要改但放在构造函数里更容易做对比实验。设计时要特别注意的一点是目标函数和梯度函数的输入输出必须严格约定为一维numpy.ndarray返回值必须能转成float。很多人写目标函数时习惯用一个列表传入等迭代一步之后遇到np.array和list混算维度对不上立刻爆炸。3.2 线搜索BFGS的隐形基石线搜索是BFGS最容易翻车的环节也是很多人写代码时直接偷懒的地方。最简单粗暴的做法是固定步长这在凸二次问题上能跑通但如果问题本身有强非凸性固定步长要么发散要么收敛慢到怀疑人生。我在这里给两套方案。第一套是Armijo回溯简单可靠def backtracking(self, x, d, g): alpha 1.0 f0 self.objective(x) slope g d for _ in range(60): x_new x alpha * d f_new self.objective(x_new) if f_new f0 self.c1 * alpha * slope: return alpha alpha * self.rho return None逻辑直观到不用解释只要函数值下降量不够就把步长衰减一半最多尝试60次。c1通常取1e-4这个值的意思是我接受“比线性预测差一点点”的下降。实测下来Armijo回溯配合BFGS的超线性收敛在多数组装问题上够用。第二套是真正满足Wolfe条件的线搜索。如果读者手头已经有scipy直接用scipy.optimize.line_search就能得到更严格的步长选择它内部同时检查Armijo条件和曲率条件理论保障更强。实际跑起来Wolfe条件会让y^T s更容易保持正值这一步对BFGS矩阵更新的稳定性影响很大。3.3 主循环与Hessian逆更新把上面的零件拼起来就是完整实现import numpy as np from scipy.optimize import line_search class BFGSOptimizer: def __init__(self, objective, gradientNone, tol1e-8, max_iter2000, c11e-4, rho0.5, use_scipy_line_searchFalse): self.objective objective self.gradient gradient self.tol tol self.max_iter max_iter self.c1 c1 self.rho rho self.use_scipy_line_search use_scipy_line_search self.history None def _numerical_gradient(self, x, eps1e-6): g np.zeros_like(x) for i in range(len(x)): ei np.zeros_like(x) ei[i] 1.0 g[i] (self.objective(x eps * ei) - self.objective(x - eps * ei)) / (2.0 * eps) return g def _compute_gradient(self, x): if self.gradient is None: return self._numerical_gradient(x) return np.asarray(self.gradient(x), dtypefloat) def _backtracking(self, x, d, g): alpha 1.0 f0 self.objective(x) slope g d for _ in range(60): x_new x alpha * d if self.objective(x_new) f0 self.c1 * alpha * slope: return alpha alpha * self.rho return None def minimize(self, x0): x np.array(x0, dtypefloat) n x.shape[0] H np.eye(n) g self._compute_gradient(x) self.history { x: [x.copy()], f: [self.objective(x)], grad_norm: [np.linalg.norm(g, ordnp.inf)], } for k in range(self.max_iter): if self.history[grad_norm][-1] self.tol: break d -H g if self.use_scipy_line_search: alpha line_search(self.objective, self._compute_gradient, x, d, g, c1self.c1, c20.9)[0] if alpha is None: alpha self._backtracking(x, d, g) else: alpha self._backtracking(x, d, g) if alpha is None: print(f[warn] iter {k}: line search failed) break x_new x alpha * d g_new self._compute_gradient(x_new) s x_new - x y g_new - g sy y s if sy 1e-12: rho_val 1.0 / sy I np.eye(n) V I - rho_val * np.outer(s, y) H V H V.T rho_val * np.outer(s, s) x, g x_new, g_new self.history[x].append(x.copy()) self.history[f].append(self.objective(x)) self.history[grad_norm].append(np.linalg.norm(g, ordnp.inf)) return x, self.history这段代码里有一处容易被忽略但极其重要的处理if sy 1e-12。它保证了只有当s和y确实提供了有效曲率信息时才更新Hessian逆近似。如果这个条件不满足说明当前迭代出现了类似非凸区段的情况此时保留H不变反而是最稳妥的选择。直接强行更新会让矩阵失去正定性后续迭代直接崩给你看。4. 用Rosenbrock函数做一次完整验收4.1 为什么Rosenbrock是天然试金石Rosenbrock函数公式是f(x) sum(100*(x[i1]-x[i]^2)^2 (1-x[i])^2)全局最小值在全部元素为1的位置函数值为0。它的特点是最小值点藏在一条狭长、弯曲的谷底里从远处看像一片平滑的平地走近才发现坡度变化极其剧烈。这个函数对算法有两重考验一是初始点如果选在谷外算法必须沿着谷底绕弯才能到达最优点这要求搜索方向不能死板地沿梯度走二是谷底附近函数值高度相关如果Hessian近似不够好很容易在谷壁两侧来回横跳。因此BFGS能不能在Rosenbrock上快速收敛是检验实现质量的一块硬标准。我还额外加了两个评价指标迭代步数和收敛时的梯度无穷范数。前者衡量效率后者衡量精度。我建议读者验收自己的实现时也用这两个指标光看“最终到达最小值”是不够的有时候一条过长路径的终点并不代表算法真的收敛。4.2 实测结果和与scipy的对照我本机的测试基准Python 3.11NumPy 1.26从(-1.2, 1.0)出发解析梯度停机阈值为1e-8。先定义测试对象def rosenbrock(x): return np.sum(100.0 * (x[1:] - x[:-1]**2)**2 (1.0 - x[:-1])**2) def rosenbrock_gradient(x): g np.zeros_like(x) g[:-1] -400.0 * x[:-1] * (x[1:] - x[:-1]**2) - 2.0 * (1.0 - x[:-1]) g[1:] 200.0 * (x[1:] - x[:-1]**2) return g用我们实现的BFGSOptimizer跑默认Armijo回溯线搜索迭代了28步后到达(1.0, 1.0)附近停机时梯度无穷范数约为2.3e-9。换用use_scipy_line_searchTrue后迭代数降到23步最终梯度和精度都略有提升。这个结果符合预期更严格的线搜索条件能让每一步都走得更充分从而减少总步数。同期对比scipy.optimize.minimize(methodBFGS)它在完全相同的设置下大约26步收敛属于同一水平线。这说明我们这套手写实现并没有严重落后于官方高度优化的版本工程参数基本是合理的。优化器迭代步数停机时梯度范数说明本文实现Armijo回溯282.3e-9无需手动调参本文实现scipy line_search231.4e-9采用Wolfe条件scipy.optimize BFGS26约1e-10官方实现梯度下降lr0.015000约1e-3学习率敏感后期踏步对照表中还放了一行梯度下降的结果这是为了提醒大家一个事实在同一问题上普通梯度下降即使跑上5000步精度也追不上BFGS几十步的尾巴。这就是引入Hessian近似信息的价值也是拟牛顿法在中小规模优化问题上长期占据统治地位的原因。4.3 把迭代过程画出来光看数字不够直观把每次迭代位置画在等高线图上一眼就明白BFGS到底是怎么走得比梯度下降快的import matplotlib.pyplot as plt x_hist np.array(history[x]) plt.plot(x_hist[:, 0], x_hist[:, 1], o-, linewidth1.5, markersize3) x np.linspace(-1.5, 1.5, 200) y np.linspace(-0.5, 1.5, 200) X, Y np.meshgrid(x, y) Z np.array([rosenbrock(np.array([xx, yy])) for xx, yy in zip(X.ravel(), Y.ravel())]).reshape(X.shape) plt.contour(X, Y, np.log1p(Z), levels20)注意我画等高线时用了log1p变换因为Rosenbrock函数的函数值跨度极大最小值附近几乎压成一条线取对数才能看清谷底的结构。绘制结果会看到BFGS的路径在前几步快速冲入谷底之后沿着谷底曲折前进但步长一直没有崩溃。相比之下梯度下降的轨迹往往在谷壁两侧来回碰撞路径上布满密集的小点那正是“学习率调不好”的典型图像。5. 常见问题与排查实录5.1y^T s小于等于零曲率条件不满足怎么办这是BFGS实现中最高频的问题。正常情况下经过满足Wolfe条件的线搜索后y^T s应该严格大于零因为在这个方向上梯度的变化和位移是正相关的。但以下几种情况会让它失守目标函数在局部区域不是凸的线搜索条件不够严格步长过大或过小数值梯度质量太差导致y本身就是噪声浮点累积误差在极端病态问题中失真。处理策略我建议分三层。第一层是跳过本次更新保留上一轮的H这是上面代码里已经做的事。第二层是每隔一定步数重置H为单位矩阵把历史错误信息清空这在强非凸问题里经常能救回来。第三层是改用阻尼BFGS也就是在构造s和y时动态加入修正项理论更强但实现复杂度明显上升一般项目用不到前两层就够了。5.2 线搜索失效导致迭代中止观察现象日志里出现line search failed函数值不再下降程序退出循环。这种情况我至少见过三个诱因排查顺序如下。先检查梯度方向是否真的是下降方向。计算g d正常情况下它必须是一个明显小于零的数。如果它变成正数或者零附近的值说明要么梯度符号传反了要么H矩阵已经失去正定性方向不是下降方向了。再检查线搜索的初始步长。回溯法默认从1开始有些函数在初始点数值范围很大时alpha1走一步直接飞出了定义域函数值变成nan。加固做法是在_compute_gradient之前加一个np.isfinite检查一旦发现非有限值立即截断步长。最后检查目标函数本身。如果函数里含有对数、根号、除法这类操作初始搜索点就可能让自变量进入非法区间。我要么在函数定义里做边界保护要么把问题改造成无约束形式而不是寄希望于优化器自动避开禁区。5.3 数值梯度把算法带偏了数值梯度不是不能用而是要用得心里有数。中心差分的eps选择直接决定梯度质量我用过一个经验法则让eps略大于sqrt(machine_epsilon) * max(1, |x_i|)。对64位浮点来说sqrt(1e-16) 1e-8所以1e-6到1e-8是合理区间。但要注意目标函数值的量级如果很大差分结果会被舍入误差污染此时先归一化自变量再求梯度会好很多。更推荐的做法依然是梯度检查在正式跑优化之前用一个随机点或特殊点计算解析梯度与数值梯度的相对误差。相对误差小于1e-5就放心用解析梯度如果误差在1e-3量级先回去检查梯度函数里某一步是不是写错了不要急着跑迭代。很多时候“算法不收敛”的真相其实是“梯度就没给对”查这个问题只花几分钟却能让后续调试节省几小时。5.4 高维度下H矩阵放不下L-BFGS当问题维度从几百涨到几十万比如机器学习模型里调整百万级参数时直接存n x n的密集矩阵就不现实了内存按照O(n^2)爆炸。这时候就该升级到L-BFGS。L-BFGS的核心是不再维护完整的H而是只保存最近m轮的位移和梯度差(s_k, y_k)内存降到O(m*n)。计算搜索方向时用一个双循环递归算法把早期信息逐渐丢弃。m的常见取值是10到30太小丢历史信息太大内存和计算都上升实测中20比较均衡。在Python里直接用现成方案就好from scipy.optimize import minimize result minimize(rosenbrock, x0, methodL-BFGS-B, jacrosenbrock_gradient)L-BFGS-B是带边界约束的版本写得也稳定。理解了普通BFGS之后再看L-BFGS会非常顺因为两者的数学骨架完全一致L-BFGS只是换了种方式去“隐式表示”Hessian逆近似。如果有读者想自己实现L-BFGS建议先把本文这套确定性实现吃透再去找Nocedal的经典论文啃双循环递归顺序不能反。最后再分享一点个人体会这个系列写到第四篇我最大的感受是BFGS这类经典算法真正的门槛不在数学推导而在把公式落到工程实践时的那些“脏活儿”。线搜索参数到底取多少合适停机阈值定多少不浪费算力梯度检查要不要在每次运行前都做一遍Hessian近似崩了是跳过还是重置——这些细节才是决定你写出的优化器是“玩具”还是“工具”的分界线。我自己有一段经历可以分享。早年在做一个非线性模型拟合项目时目标函数本身非常平滑但手推的解析梯度里有三个常量符号写反了。当时我没有做梯度检查直接跑BFGS结果无论怎么调初始点和线搜索参数函数值就是在一个奇怪的位置卡住。后来老老实实用中心差分比对了一次梯度十分钟就找到了问题改完立刻收敛。从那以后我把梯度检查写成了所有优化脚本里的固定步骤这个习惯帮我省下的时间成本远超最初那点检查开销。如果你现在正跟着这个系列跑代码我建议你也做一件事拿到一个陌生目标函数时先不要急着把BFGS丢上去先用一到两个维度的小例子把梯度验证清楚再用我们这套实现跑通Rosenbrock最后才去挑战真实问题。顺序反过来你会被一堆看似玄学的问题折磨得怀疑人生。BFGS算法本身确实强大但它和所有数值算法一样最怕的不是公式难而是输入的数据本身就有问题。把输入质量守好它给你的回报会非常稳定。