这是 BFGS 算法系列手把手教程的第四篇。前面三篇我们把拟牛顿法的动机、BFGS 公式的数学推导以及一个最基础的可运行骨架都讲完了代码也能在简单问题上跑通。但到这里你会发现骨架离“真正能用”还差好几块关键拼图怎么保证每一步都下降怎么判断算法该停矩阵更新遇到数值异常怎么办所以这一篇我不打算再绕原理直接把上一版骨架升级成一个完整、可复现的 BFGS 优化器补齐非精确线搜索、收敛判据、日志输出再用标准测试函数做数值实验最后把我实践中踩过的坑整理成排错清单。适合三类读者跟着系列走到这里的人、学过最优化理论但没手写过完整实现的人以及正在调试自己优化代码的人。代码基于 Python 3.8 和 NumPy环境只要一个 numpy 就够了。1. 本篇定位从公式骨架到可用优化器还差多少1.1 前三篇解决了什么这一篇要补什么前几篇的逻辑链大致是这样的第一篇从目标函数和梯度讲起说明为什么优化问题本质上是在找梯度为零的点第二篇引入拟牛顿条件解释为什么可以用逆 Hessian 的近似矩阵 (H_k) 来构造搜索方向而不是直接解牛顿方程第三篇给出了最简框架让你能看到“计算方向 - 更新 x - 更新 H”这条主循环长什么样。但那个框架距离一个真正能投入使用的优化器还有不少距离。首要问题是步长第三篇的代码通常在每步设一个固定步长这对凸二次型勉强能用一旦碰到曲率变化剧烈的函数固定步长要么让迭代震荡要么为了保证不震荡被压得极小导致几十步都挪不动。其次是停止条件如果只靠“函数值不再下降”来停机你永远分不清自己是收敛了还是卡在了数值噪声里。再一个是更新矩阵的健壮性(H_k) 更新公式本身要求 (s_k^T y_k 0)但实际数值计算中这个条件可能被破坏不做保护的话程序会直接跑飞。这一篇的目标就是逐一解决这些问题最后给出一份能在标准测试函数上稳定工作的完整代码。1.2 整体模块划分与数据流设计先把最终代码的模块拆开你会看到整个优化器其实就由四块组成目标函数模块提供 (f(x)) 和梯度 (\nabla f(x))在测试实验里我会分别实现二次型、Rosenbrock 和 Himmelblau 函数线搜索模块给定当前点、搜索方向返回一个合适的步长 (\alpha)本篇用 Armijo 回溯作为主实现BFGS 更新模块根据相邻两点的位移 (s_k) 和梯度差 (y_k)更新逆 Hessian 近似矩阵 (H_k)主循环模块负责计算搜索方向、调用线搜索、判断收敛、输出日志。数据流是固定的从 (x_k) 出发先算梯度 (g_k)得到搜索方向 (p_k -H_k g_k)然后线搜索返回步长 (\alpha_k)更新点位 (x_{k1} x_k \alpha_k p_k)接着用新梯度 (g_{k1}) 计算 (s_k x_{k1} - x_k) 和 (y_k g_{k1} - g_k)最后更新 (H_{k1})。整个循环里最容易搞反的就是 (s_k) 和 (y_k) 的顺序很多初学者会把 (y_k) 算成 (g_k - g_{k1})行列式一翻更新公式整个方向就错了。一句话总结设计思路让每个模块只干一件事主循环只做调度这样定位问题时非常快。如果你的梯度函数写错了你只需要查目标函数模块而不用在主循环里翻来翻去。1.3 环境准备与依赖这篇代码只用 NumPy绘图和 SciPy 对照放在可选项里。安装命令很简单pip install numpy matplotlib scipy如果你机器上还没有 Python 环境去官网下载 3.8 以上的安装包安装时勾选“Add Python to PATH”然后在命令行里执行上面的 pip 命令即可。Windows 用户如果遇到 wheel 安装失败多半是 Python 版本太旧升级到 3.9 基本能解决。Mac 用户如果用的是系统自带 Python我建议装一个 Homebrew 的 Python 再继续避免系统权限问题。2. BFGS更新公式的工程化实现别把NumPy公式抄错2.1 从数学公式到NumPy代码的关键映射先温习一下 BFGS 算法中最核心的逆 Hessian 近似更新公式。我们用 (H_k) 表示 Hessian 逆的近似矩阵(s_k x_{k1} - x_k)(y_k \nabla f(x_{k1}) - \nabla f(x_k))定义 (\rho_k 1 / (y_k^T s_k))则 BFGS 更新可以写成紧凑形式[ H_{k1} (I - \rho_k s_k y_k^T)^T H_k (I - \rho_k s_k y_k^T) \rho_k s_k s_k^T ]这个形式比展开成四项加减更稳定因为它在数学上等价于一个秩二修正的乘积表达。对应到 NumPy 就是几行import numpy as np def bfgs_update(H, s, y): rho 1.0 / np.dot(y, s) V np.eye(len(s)) - rho * np.outer(s, y) H_new V.T H V rho * np.outer(s, s) return H_new为什么这样写从数值角度看直接套用展开公式会涉及 (H_k s_k y_k^T y_k s_k^T H_k) 这种两个大矩阵相加再除以一个标量的操作浮点误差是逐项累积的。而紧凑形式本质上是先做两个矩阵积再做一次秩一修正偏差更小。我在实际测试里对比过两种写法在病态问题上的差异等 (\rho) 达到 (10^6) 量级时展开式算出的矩阵甚至可能出现特征值为负的情况而紧凑形式明显更稳。这里还要提醒一个细节代码中V用了np.outer(s, y)而有的教科书把 (V) 定义为 (I - \rho y s^T)两者只差一个转置。如果你以前记的是另一种写法套代码时务必用 (H_{new} V^T H V \rho s s^T) 验证一下不要 V 和转置搞混。2.2 初始H_0怎么选才合理BFGS 的迭代公式本身保证不了初始矩阵的质量所以 (H_0) 的选择直接影响前面几步的收敛速度。最简单的做法是 (H_0 I)这在大多数教学例子里够用但遇到条件数很大的问题时第一步的搜索方向和真正需要的方向可能差一个很大的比例因子导致线搜索要回溯很多次。实践中我更推荐一种启发式初始化在拿到第一步 (s_0 x_1 - x_0) 和 (y_0 g_1 - g_0) 之后把 (H_1) 的初值设为 ((y_0^T y_0) / (y_0^T s_0) \cdot I)。背后的直觉是我们希望 (H_0) 在某种意义上“贴合”目标函数的真实曲率量级而 (y_0^T y_0 / y_0^T s_0) 在二次型情形正好就等于系数矩阵的某个缩放估计。def initial_hessian_inv(x0, y0, s0): sy np.dot(s0, y0) yy np.dot(y0, y0) if sy 0 and yy 0: return (yy / sy) * np.eye(len(x0)) return np.eye(len(x0))实测下来在 Rosenbrock 函数上做对比固定 (H_0I) 大约需要 35~45 次迭代而用这种缩放初始化可以把迭代次数压到 25~30 次左右。注意这个技巧要在第一步之后才能用也就是说主循环里首次搜索用的是 (H_0I)第一次迭代结束后再重新赋值。2.3 数值稳定性处理对称性、曲率条件与溢出保护BFGS 更新公式在精确算术下能保持 (H_k) 的对称正定性但浮点计算会积累微小误差迭代几百步后(H) 矩阵的对角线附近可能出现 (10^{-8}) 量级的不对称。最省事的修复是在每次更新后做一次对称化H_new 0.5 * (H_new H_new.T)第二个必须处理的是曲率条件 (s_k^T y_k 0)。理论上只要线搜索满足 Wolfe 条件这个内积一定大于零但实际中如果目标函数非光滑或者线搜索实现得不够精细这个条件会被突破。一旦 (s_k^T y_k \le 0)(\rho) 会变成负数后续更新出来的矩阵几乎肯定不正定接着搜索方向就不下降算法直接崩溃。我的习惯是在更新入口加一个保护def safe_bfgs_update(H, s, y): sy np.dot(s, y) if sy 1e-12 * np.linalg.norm(s) * np.linalg.norm(y): return H # 跳过本次更新保持原有矩阵 return bfgs_update(H, s, y)阈值用相对值而不是绝对 0是因为不同问题的变量量级差很多绝对值判断在缩放后不靠谱。还有一个溢出场景是 (s) 或 (y) 的模长达到 (10^8) 以上直接算 (\rho) 可能溢出成 inf进一步导致矩阵出现 NaN。为了保险可以在bfgs_update开头检查np.isfinite(rho)不是有限数就直接返回原矩阵。2.4 矩阵损坏了怎么补救即使加了保护偶尔还是会遇到更新后矩阵不正定的情况。检测方法很简单算一下最小特征值eigs np.linalg.eigvalsh(H) if eigs[0] 0: H H (np.abs(eigs[0]) 1e-6) * np.eye(n)给对角加一个小的正数等价于对矩阵做正则化能很快把特征值推回正半轴。另一个更“粗暴”但实用的方案是重置 (HI)代价是丢失前面累积的曲率信息一般只用于算法已经明显出问题的情况。从工程角度说如果频繁走到这个分支你先别急着调正则项回去查线搜索是不是没满足曲率条件——那才是根因。3. 非精确线搜索手写实现Armijo回溯与Wolfe条件的落地3.1 为什么BFGS必须配线搜索不看推导只谈现象的话BFGS 恰恰是为“非精确线搜索”设计的。精确线搜索需要反复求解一维最小化问题每次都要算很多次函数值成本不低而 BFGS 只需要一个粗糙但保证下降的步长就能持续构造出足够好的曲率信息。固定步长在这个框架里尤其危险。假设某一步 (H_k) 给的搜索方向很棒但步长太大冲过了头函数值不降反升步长太小又会让 (s_k) 变得很小数值上 (y_k) 和 (s_k) 的关系失真更新出来的矩阵质量会急剧恶化。所以线搜索的作用有两个一是让目标函数充分下降二是保证新产生的 (s_k) 和 (y_k) 满足曲率条件从而让 (H_k) 能持续积累真实曲率信息。3.2 Armijo回溯的代码实现与参数选取Armijo 条件的本质是要求函数值下降量至少是线性下降量的一定比例[ f(x_k \alpha p_k) \le f(x_k) c_1 \alpha \nabla f(x_k)^T p_k ]其中 (c_1) 通常取 (10^{-4})这保证下降是“充分”的但又不至于苛刻到找不到步长。回溯策略是从 (\alpha1) 开始不满足条件就把 (\alpha) 乘以一个收缩因子 (\rho)常用 (\rho0.5) 或 (0.8)。def armijo_backtracking(f, x, p, g, alpha_start1.0, c11e-4, rho0.5, max_iter50): f0 f(x) g0 np.dot(g, p) if g0 0: raise RuntimeError(搜索方向不是下降方向请检查梯度或H矩阵) alpha alpha_start f_new f0 for _ in range(max_iter): f_new f(x alpha * p) if f_new f0 c1 * alpha * g0: return alpha, f_new alpha * rho return alpha, f_new注意开头对 (g0) 的判断。向量 (p) 与梯度内积为负是搜索方向下降的必要条件如果这一步都不满足后面的回溯纯属浪费时间而且大概率会掩盖更深层的问题梯度写错、(H) 不正定等。把这个检查放在线搜索入口能帮你第一时间定位错误。3.3 从Armijo升级到强Wolfe条件如果你只要求代码能跑出结果Armijo 回溯通常就够了。但如果你想保证 BFGS 的 (H_k) 永远正定最好让步长同时满足强 Wolfe 条件。强 Wolfe 除了 Armijo 下降条件还要求[ |\nabla f(x \alpha p)^T p| \le c_2 |\nabla f(x)^T p| ]这个条件的含义是走到新点时沿 (p) 方向的方向导数已经不再显著说明步长没有停在梯度的陡坡上。(c_2) 通常取 (0.9)。下面给出一个能工作的实践版搜索逻辑先尝试 (\alpha1)若违背 Armijo 条件就继续回溯若满足下降条件但不满足曲率条件则放大或微调 (\alpha)在有限步内找到候选步长。def wolfe_line_search(f, grad, x, p, g, alpha_start1.0, c11e-4, c20.9, max_iter100): f0 f(x) g0 np.dot(g, p) alpha alpha_start for _ in range(max_iter): f_new f(x alpha * p) if f_new f0 c1 * alpha * g0: alpha * 0.5 continue g_new grad(x alpha * p) if abs(np.dot(g_new, p)) c2 * abs(g0): return alpha, f_new alpha min(alpha * 2.0, alpha_start * 4.0) return alpha, f_new严格说正规产品级 Wolfe 搜索需要两阶段先括定区间再在区间内缩我上面给的是“在多数问题中能跑通”的版本。实际使用中你会发现对于光滑凸函数这个简单版本已经比纯 Armijo 稳定不少尤其是碰到曲率变化大的区域时步长能自动调整到合适区间。3.4 步长策略的实测对比我在二维二次型函数上做过一个直观实验同一目标函数、同一初始点分别用固定步长 0.1、纯 Armijo、简单 Wolfe 三种策略跑 BFGS。固定步长在前几轮下降看起来挺快但到接近最优值时开始来回震荡迭代 500 次梯度范数还在 (10^{-1}) 量级Armijo 花了大约 30 次迭代收敛到 (10^{-8})Wolfe 版本由于步长选择更精确几乎每一步都在做有效更新20 次以内就收敛了。这个实验结果并不意外。对 BFGS 来说精确的步长并不是必须的但步长必须能让函数值充分下降同时保留一定的“前进幅度”。Armijo 回溯在初始点离最优点很远时通常会多回溯几次而 Wolfe 条件的曲率约束恰好能避免步长卡在太小的区域。所以我的个人建议是教学示范用 Armijo工程项目尽量上 Wolfe 或者直接调用 SciPy 的快速实现。4. 收敛判据和日志系统让迭代过程一眼看清4.1 收敛判据怎么定才靠谱很多初学者只知道“函数值不再下降就停”但这在实际数值环境中非常危险。函数值在浮点精度附近会出现锯齿状抖动你可能在局部极小附近来回震荡却误以为已经收敛。更可靠的判据是以梯度范数为主最优点必须满足 (\nabla f(x)0)因此我们可以定义[ |\nabla f(x_k)|_{\infty} \le \varepsilon ]当 (\varepsilon10^{-5}) 时绝大多数平滑问题的结果已经足够精确。不过梯度范数并不总是能压到特别小尤其是问题本身有噪声或条件数极大时你可能只达到 (10^{-3}) 就卡住了。这时需要辅助判据连续两次迭代的 (x) 距离很小或者函数值相对变化很小。一套实用的停止规则是def converged(f_old, f_new, x_old, x_new, grad_new, tol_grad1e-5, tol_x1e-8, tol_f1e-12): if np.linalg.norm(grad_new, ordnp.inf) tol_grad: return True if np.linalg.norm(x_new - x_old, ordnp.inf) tol_x: return True if abs(f_new - f_old) tol_f * (1.0 abs(f_old)): return True return False注意最后一个判据用的是相对变化不是绝对变化否则不同量级的目标函数会互相影响。永远不要忘记设置最大迭代次数——这是我踩过最多次的坑没有兜底限制的循环在数值异常时会无限跑下去。4.2 日志输出怎么设计日志是排查优化器问题的第一工具。我习惯每步输出五个指标迭代数、函数值、梯度无穷范数、步长、当前点的坐标。这样从输出里一眼就能看出问题出在哪个阶段。def print_log(iteration, x, fval, grad_norm, alpha): print(fIter {iteration:4d} | f {fval: .10f} | ||g|| {grad_norm: .3e} | alpha {alpha: .3e} | x {np.round(x, 6)})实际运行时会看到类似这样的输出Iter 0 | f 24.2000000000 | ||g|| 5.318e02 | alpha 1.000e00 | x [-1.2 1. ] Iter 1 | f 4.2300000000 | ||g|| 2.928e01 | alpha 1.000e00 | x [-0.96 0.92] ... Iter 30 | f 0.0000009982 | ||g|| 9.812e-06 | alpha 1.000e00 | x [ 0.999998 0.999996]从日志判断问题的方法是如果 (f) 下降但梯度范数一直不降说明步长太小线搜索在反复回溯这时要检查初始 (H_0) 的缩放如果梯度范数在下降但迭代次数异常多说明搜索方向可能不是最优的问题多半出在更新公式或线搜索参数上。4.3 完整BFGS主循环参考实现把上面的模块拼起来就是一份可以直接复制使用的优化器class BFGSOptimizer: def __init__(self, f, grad, tol_grad1e-5, max_iter1000, verboseTrue): self.f f self.grad grad self.tol_grad tol_grad self.max_iter max_iter self.verbose verbose def minimize(self, x0): x np.array(x0, dtypefloat) H np.eye(len(x)) g self.grad(x) f self.f(x) if self.verbose: print_log(0, x, f, np.linalg.norm(g, ordnp.inf), 0.0) for k in range(1, self.max_iter 1): p -H g alpha, f_new armijo_backtracking(self.f, x, p, g) x_new x alpha * p g_new self.grad(x_new) s x_new - x y g_new - g H safe_bfgs_update(H, s, y) f, g, x f_new, g_new, x_new if self.verbose: print_log(k, x, f, np.linalg.norm(g, ordnp.inf), alpha) if np.linalg.norm(g, ordnp.inf) self.tol_grad: return x, f, g, k return x, f, g, self.max_iter这段代码不追求华丽但结构是完整的。换个目标函数你只需要改f和grad两个函数引用主循环完全不用动。5. 在三个标准测试函数上的完整数值实验5.1 正定二次型最快验证实现是否正确的试金石写优化器第一件事不是跑复杂函数而是先把一个二维正定二次型跑通。我用A np.array([[2.0, 0.5], [0.5, 1.0]]) b np.array([1.0, 0.5]) def f_quad(x): return 0.5 * x A x b x def grad_quad(x): return A x b初始点选x0 np.array([-3.0, 2.0])。理论上 BFGS 在精确线搜索下最多 (n) 步收敛对二维问题就是 2 步左右即使使用非精确 Armijo通常也在 5 次迭代内达到 (10^{-8}) 的梯度范数。如果你在这个测试函数上跑了十几步还收敛不了基本可以确定是线搜索或更新公式出了问题别急着怀疑算法理论。5.2 Rosenbrock函数经典验证场Rosenbrock 函数是优化算法验证的“标配”[ f(x, y) 100 (y - x^2)^2 (1 - x)^2 ]它的特点是有一个极狭窄的弯曲“山谷”从初始点 ((-1.2, 1.0)) 出发梯度下降法会被谷壁反弹需要成千上万次迭代才能慢慢挪向最优解 ((1,1))。而 BFGS 能利用曲率信息快速调整搜索方向通常在 30 次左右收敛。解析梯度如下def f_rosen(x): return 100.0 * (x[1] - x[0]**2)**2 (1.0 - x[0])**2 def grad_rosen(x): dx -400.0 * x[0] * (x[1] - x[0]**2) - 2.0 * (1.0 - x[0]) dy 200.0 * (x[1] - x[0]**2) return np.array([dx, dy])我用Armijo回溯版本跑出来的典型结果是前 10 步函数值从 (24.2) 降到 (10^{-1})中间进入谷道后步长会被适度压缩之后函数值稳定下降到 (10^{-10}) 以下迭代次数大约 26 到 32 之间。相比梯度法这个优势非常明显。5.3 Himmelblau函数多个局部极小时的行为观察Himmelblau 函数有四个局部极小点适合用来提醒新手“BFGS 是局部优化器”[ f(x, y) (x^2 y - 11)^2 (x y^2 - 7)^2 ]不同初始点会收敛到不同的极小所以看到收敛结果时不要默认它是全局最优。几个已知极小点大致是初始点收敛点常见(0, 0)约 (3.0, 2.0)(4, 0)约 (3.584, -1.848)(-4, 0)约 (-3.779, -3.283)(0, -4)约 (-2.805, 3.131)跑这个函数的价值在于验证你的优化器能稳定地从一个随机初始点收敛到某个局部解并且函数值都落在 0 附近。如果从某个初始点出发不降反升那大概率不是公式问题而是线搜索没有正确生效。5.4 与梯度下降的对比实验为了让你对 BFGS 的优势有体感我把梯度下降配合同样的 Armijo 线搜索和 BFGS 放到同一个 Rosenbrock 问题上对比方法迭代次数最终梯度范数备注梯度下降 Armijo大约 3000 次以上(10^{-5}) 左右受限于狭窄谷道BFGS Armijo大约 30 次(10^{-8}) 以下曲率信息持续累积梯度下降在谷道里每步都被迫走得很短因为步长一旦偏大就会冲出谷壁BFGS 则通过 (H_k) 积累“该往哪个方向拐”的信息自然能选择一条更顺滑的路径。这个差距在更高维和非二次问题上只会更夸张也是工程里更愿意用拟牛顿系列的原因。6. 手写BFGS排错实录最常踩的几个坑6.1 函数值不降反升直接发散遇到这种情况第一反应不是怀疑算法而是检查方向是否真的是下降方向。在armijo_backtracking里我们已经加了g0 0的报错但如果你的梯度是手算的很可能符号出了问题。另一个常见原因是初始步长太大而回溯收缩因子设成了接近 1 的值比如 0.95导致一次回溯需要上百次函数调用才能救回来。把 (\rho) 设成 0.5一般情况下几轮就能找到合适步长。6.2 梯度范数卡在某个值不再下降这通常有两种可能。第一是目标函数本身在最优解附近梯度很小但并非最优比如噪声函数这时要降低收敛容差或者改用其他判据。第二是你的 (H_k) 累积的曲率信息“失忆”了更新矩阵一直在循环但方向没有实质改善。我在实验里发现遇到这种情况把 (H_0) 重新缩放一次或者干脆每 50 步重置一次 (H)往往能打破僵局。6.3 解析梯度对不对用中心差分先验证手写梯度最容易错尤其是 Rosenbrock 这种链式法则不熟会写错的函数。强烈建议每个新目标函数都先用中心差分做一次验证def check_grad(f, grad, x, eps1e-6): g_num np.zeros_like(x) for i in range(len(x)): e np.zeros_like(x) e[i] eps g_num[i] (f(x e) - f(x - e)) / (2.0 * eps) return np.linalg.norm(grad(x) - g_num, ordnp.inf)如果返回的误差超过 (10^{-5})基本可以断定梯度公式写错了。别小看这一步它能帮你省掉好几小时的调试时间。6.4 (s_k^T y_k) 为负或接近 0这个坑我在第 2.3 节已经讲过但值得单独拿出来强调一旦出现这个情况说明线搜索没能产生满足曲率条件的迭代点。原因通常有两个——线搜索只是纯 Armijo而 Armijo 本身不保证曲率条件或者目标函数在某点附近不可导。处理方案就是safe_bfgs_update里的跳过逻辑不要强行更新。如果你希望根治那就把线搜索换成第 3.3 节的 Wolfe 版本。6.5 高维问题内存爆炸BFGS 需要维护一个 (n \times n) 的矩阵(n1) 万时矩阵就要占约 800MB 内存(n10) 万就是直接不可用。遇到高维问题正确的思路是转向 L-BFGS不显式存储 (H_k)只保留最近的 (m) 组 ((s_i, y_i)) 对然后用双循环递归计算搜索方向。这是自然的下一步扩展也是工程界更常用的方案。6.6 用SciPy结果做对照自写代码不管看起来多正确最好还是和 SciPy 的实现对照一次from scipy.optimize import minimize res minimize(f_rosen, x0, jacgrad_rosen, methodBFGS, options{gtol: 1e-5, disp: True})用我上面的测试函数跑一遍你会看到自写代码和 SciPy 最终函数值基本一致。注意因为线搜索和容差细节不同迭代轨迹不必完全对齐但最终精度应该在同一个量级。如果差得离谱回到梯度检查环节。最后说一个我自己的调试习惯。每次写完一个优化器我第一件事不是跑到 Rosenbrock 上看最终结果而是先拿二维正定二次型做冒烟测试。如果二次型都不能在十次以内迭代收敛那问题多半出在线搜索或更新顺序上而不是算法本身。这个习惯帮我排掉了至少一半的 bug。下一篇如果要继续往下走就是 L-BFGS、阻尼 BFGS 以及带约束问题的扩展这些在工程中的出场率比裸 BFGS 更高。