学最优化绕不开拟牛顿法。我刚开始学那会儿看着教材里BFGS的更新公式脑子里全是问号那一大串 ( V^T H V \rho s s^T ) 到底是怎么凑出来的教材里要么直接丢结论要么搬出一堆加权Frobenius范数的最小化问题看得人昏昏欲睡。后来自己把推导一步步抄下来才发现核心逻辑其实特别朴素——拟牛顿法从头到尾只坚持了一件事用梯度信息去逼近Hessian不同方法只是逼近的“补全策略”不同而已。这篇笔记就把这条线完整捋一遍从牛顿法为什么贵讲起到割线方程、秩1修正、秩2修正再到BFGS公式的完整推导最后给一份能直接跑通的最小Python实现适合正在学最优化理论、或者想搞明白BFGS为什么长这样的同学参考。1. 先说清楚牛顿法贵在哪拟牛顿法在讨价还价什么考虑无约束光滑优化问题[ \min_{x \in \mathbb{R}^n} f(x) ]经典牛顿法每一步在点 (x_k) 处构造二次模型[ m_k(d) f_k g_k^T d \frac{1}{2} d^T \nabla^2 f(x_k) d ]其中 (g_k \nabla f(x_k))(d) 是搜索方向。对 (m_k(d)) 求极小令梯度为零得到线性方程组[ \nabla^2 f(x_k) d -g_k ]这就是牛顿方向。理论上是很好的方法——如果每一步都用精确Hessian在极小点附近能达到二次收敛迭代次数少得惊人。但工程上牛顿法有三个“贵”求Hessian贵。二阶偏导要算 (O(n^2)) 个分量对复杂的目标函数来说解析推导麻烦自动微分也烧计算量。解线性方程组贵。直接法求解 (n) 元方程组是 (O(n^3)) 的复杂度(n) 稍微大一点就扛不住。即便用迭代法也需要额外设计预条件麻烦。Hessian不保证正定。在非凸区域(\nabla^2 f(x_k)) 可能是不定矩阵牛顿方向不一定是下降方向。为了保证收敛还得加修正策略复杂度进一步提升。所以实际问题里很少有人直接用纯牛顿法去硬刚大尺度问题。那能不能只用一阶信息——也就是梯度——来间接估计曲率这是拟牛顿法出现的根本动机。拟牛顿法的思路是维护一个矩阵 (B_k) 作为Hessian (\nabla^2 f(x)) 的近似迭代过程中只靠相邻两步的梯度变化量 (y_k) 和位移量 (s_k) 去更新 (B_k)。每一步求方向只需要做一次矩阵向量乘法[ d_k -B_k^{-1} g_k ]注意这里维护的是 (B_k)但实际用的时候要的是 (B_k^{-1})所以后面大部分实现直接维护逆Hessian近似 (H_k)一步方向计算就是[ d_k -H_k g_k ]这样每步复杂度是 (O(n^2)) 的矩阵向量乘没有二阶导数没有线性方程组求解。在 (n10^6) 这种规模下牛顿法基本无望而拟牛顿法配合有限内存版本L-BFGS依然能跑。一句话总结拟牛顿法做的事情就是“用廉价的一阶信息凑出一个足够好的曲率估计”。接下来的问题就是——怎么凑。2. 割线方程拟牛顿法唯一称得上“必须满足”的约束假设我们已经从 (x_k) 走到了 (x_{k1})记[ s_k x_{k1} - x_k,\qquad y_k g_{k1} - g_k ]把梯度函数做一阶Taylor展开或者更严谨地用积分中值定理[ y_k g_{k1} - g_k \int_0^1 \nabla^2 f(x_k t s_k) s_k , dt ]可以看到梯度差 (y_k) 等于Hessian沿 (s_k) 方向在区间上的某种平均作用。如果拟牛顿矩阵 (B_{k1}) 要近似Hessian它至少应该复现这条信息[ B_{k1} s_k y_k ]这就是割线方程也叫拟牛顿方程是整个拟牛顿法唯一称得上“必须满足”的硬约束。这个约束其实很弱。一个对称矩阵在 (\mathbb{R}^n) 上有 (n(n1)/2) 个独立自由度割线方程只给定了 (n) 个线性方程。也就是说满足割线方程的对称矩阵有无穷多个。举个例子类比就像我说“有一个过点 ((1,2)) 的光滑函数”你根本没法确定它是直线、抛物线还是三角函数。割线方程也只约束了一个方向上的“平均曲率”剩下海量的自由度要靠其他偏好去填充。拟牛顿法的不同流派本质上就是给“如何补全自由度”加上了不同的偏好对称性(B_{k1}) 必须对称这是结构约束正定性希望 (B_{k1}) 正定保证搜索方向 (d_k -B_{k1}^{-1}g_k) 是下降方向这是品质约束最小变动希望 (B_{k1}) 尽量靠近 (B_k)不要因为某一步的噪声把之前积累的曲率信息全冲掉这是风格约束。SR1、DFP、BFGS都是这些约束下的不同解。下面从最简单的秩1修正开始看。3. 从秩1修正到秩2修正BFGS公式是怎么一步一步长出来的3.1 SR1为什么容易翻车最简单的补全方式是让当前矩阵加一个秩1外积修正[ B_{k1} B_k c, v v^T ]代入割线方程[ B_k s_k c, v (v^T s_k) y_k ]移项[ c, (v^T s_k), v y_k - B_k s_k ]记残差 (r_k y_k - B_k s_k)可以看到 (v) 必须与 (r_k) 共线。取 (v r_k)两边点积后解得[ c \frac{1}{r_k^T s_k} ]于是得到对称秩1更新SR1Symmetric Rank-One[ B_{k1} B_k \frac{(y_k - B_k s_k)(y_k - B_k s_k)^T}{(y_k - B_k s_k)^T s_k} ]这个公式推导起来很干净但实际用起来心里要提心吊胆。问题出在分母(r_k^T s_k) 可能等于零公式直接爆掉(r_k^T s_k 0) 时秩1项是负定方向极大概率把 (B_{k1}) 弄得不正定即便分母为正SR1也不能保证每一步都保持正定。所以SR1虽然在某些问题上收敛速度出奇地快但稳定性很差属于“时灵时不灵”的类型。工程上基本不会当作默认算法用。它最大的意义在于告诉我们秩1修正的自由度太少了控制不了正定性。3.2 秩2修正BFGS的B更新公式既然一个秩1不够自然想到秩2修正。设[ B_{k1} B_k a, u u^T b, v v^T ]其中 (u,v) 是两个待定的方向(a,b) 是待定系数。把它代入割线方程[ B_k s_k a, u (u^T s_k) b, v (v^T s_k) y_k ]整理[ a, (u^T s_k), u b, (v^T s_k), v y_k - B_k s_k ]关键的一步来了。看右边的结构(y_k - B_k s_k)这是两个向量的差。如果我们让第一项去“生成” (y_k)让第二项去“消灭” (B_k s_k)那方向的选择就很自然了[ u y_k,\qquad v B_k s_k ]这样只要令[ a (y_k^T s_k) 1,\qquad b ((B_k s_k)^T s_k) -1 ]割线方程就能直接满足。于是[ a \frac{1}{y_k^T s_k},\qquad b -\frac{1}{s_k^T B_k s_k} ]代回秩2修正表达式得到BFGS对Hessian近似 (B) 的更新公式[ B_{k1} B_k - \frac{B_k s_k s_k^T B_k}{s_k^T B_k s_k} \frac{y_k y_k^T}{y_k^T s_k} ]这个公式长得很对称第一项减去的是沿 (B_k s_k) 方向的旧曲率信息第二项加上的是沿 (y_k) 方向的新曲率观测。你可能会觉得“选 (u y_k)(v B_k s_k)”有点拍脑袋。更严格的推导可以写成变分问题在所有满足割线方程和对称性的矩阵中找一个在某种加权Frobenius范数意义下离 (B_k) 最近的矩阵解出来的结果正是这个公式。换句话说BFGS的更新是“表面像拍脑袋实际是最小变动原理的必然结果”。这个变分版本我就不展开了常见教材里都能找到。顺带一提如果在逆Hessian近似 (H_k) 上直接做同样形式的秩2构造也就是设[ H_{k1} H_k a u u^T b v v^T ]并要求 (H_{k1} y_k s_k)会得到另一个著名公式——DFP更新[ H_{k1} H_k - \frac{H_k y_k y_k^T H_k}{y_k^T H_k y_k} \frac{s_k s_k^T}{y_k^T s_k} ]BFGS和DFP之间存在对偶关系互换 (s \leftrightarrow y)、(B \leftrightarrow H)两个公式可以互相转化。BFGS在实践中表现更稳所以成了默认选择。3.3 为什么说BFGS是“最小变动”我学这个公式时最大的体会是它不是在“猜”而是在“补全”。每一步迭代我们手里已经有一个积累了大量历史信息的 (B_k)。割线方程要求新矩阵在 (s_k) 方向上的作用等于观测到的 (y_k)。满足这个条件的矩阵很多但绝大多数都会把 (B_k) 辛苦积累的信息毁掉。BFGS做的是在保持对称性和割线约束的前提下尽可能少地改变 (B_k)同时利用“减去秩1”和“加上秩1”的组合让新矩阵有机会保持正定。用大白话说旧信息别全扔新信息也别不用两边各让一步。这个朴素的分配原则就是BFGS能成为无约束光滑优化默认算法的底层原因。4. 从近似Hessian到近似逆HessianSMW公式才是幕后功臣BFGS的B更新公式理解起来直观但实际实现时我们更想要 (H_{k1} B_{k1}^{-1})因为迭代方向的计算变成[ d_k -H_k g_k ]一步矩阵向量乘完事不用解线性方程组。那问题就来了已知 (B_{k1}) 的更新怎么求逆答案是用Sherman-Morrison-Woodbury公式。先回忆一下秩1情况下的核心结论[ (A u v^T)^{-1} A^{-1} - \frac{A^{-1} u v^T A^{-1}}{1 v^T A^{-1} u} ]我们的 (B_{k1}) 更新可以看成是“加一个秩1再减一个秩1”[ B_{k1} B_k u_1 u_1^T - u_2 u_2^T ]其中[ u_1 \frac{y_k}{\sqrt{y_k^T s_k}},\qquad u_2 \frac{B_k s_k}{\sqrt{s_k^T B_k s_k}} ]对上述表达式连续使用两次SMW公式经过一顿代数整理就能得到著名的BFGS逆Hessian更新公式[ H_{k1} \left(I - \rho_k s_k y_k^T\right) H_k \left(I - \rho_k y_k s_k^T\right) \rho_k s_k s_k^T ]其中[ \rho_k \frac{1}{y_k^T s_k} ]这个形式比B更新公式更常用。为了记忆方便可以写成[ V_k I - \rho_k s_k y_k^T ][ H_{k1} V_k H_k V_k^T \rho_k s_k s_k^T ]这样结构就清晰了中间是“用 (V_k) 把旧 (H_k) 夹一下”再补一个秩1项。如果把它展开还能写成另一种等价形式方便程序里逐项验证[ H_{k1} H_k \frac{(s_k - H_k y_k)s_k^T s_k(s_k - H_k y_k)^T}{y_k^T s_k} - \frac{(s_k - H_k y_k)^T y_k}{(y_k^T s_k)^2} s_k s_k^T ]刚看到这个展开式的人容易懵其实它就是简洁形式展开后的结果。我个人推荐实现时用带 (V_k) 的简洁形式因为它既好写也天然保持了 (H_{k1}) 的对称性。有人会问怎么确认更新后的 (H_{k1}) 确实满足拟牛顿方程的逆形式 (H_{k1} y_k s_k)直接验证一下[ H_{k1} y_k V_k H_k (y_k - \rho_k y_k (s_k^T y_k)) s_k ]因为 (\rho_k (s_k^T y_k) 1)括号里是 (y_k - y_k 0)所以第一项消失剩下的正是[ H_{k1} y_k s_k ]验证通过。这说明从 (B_{k1}) 出发用SMW公式求逆得到的 (H_{k1}) 确实是对应逆Hessian的合理近似不是随便凑的。5. 最小实现验证用一份简单Python代码把BFGS跑在Rosenbrock上理论推完动手验证一下。下面这份代码刻意保持最小只依赖numpy和scipy的线搜索函数核心就是BFGS的 (V_k H_k V_k^T \rho s s^T) 更新。import numpy as np from scipy.optimize import line_search def rosen(x): return 100.0 * (x[1] - x[0]**2)**2 (1.0 - x[0])**2 def rosen_grad(x): g np.zeros_like(x) g[0] -400.0 * x[0] * (x[1] - x[0]**2) - 2.0 * (1.0 - x[0]) g[1] 200.0 * (x[1] - x[0]**2) return g def bfgs(f, grad, x0, max_iter100, tol1e-6): n len(x0) H np.eye(n) x x0.copy() g grad(x) for it in range(max_iter): if np.linalg.norm(g, np.inf) tol: break d -H g alpha, _, _, _, _, _ line_search(f, grad, x, d, g) if alpha is None: # 线搜索失败时退回到一个小步长避免程序中断 alpha 1e-3 s alpha * d x_new x s g_new grad(x_new) y g_new - g rho 1.0 / (y s) V np.eye(n) - rho * np.outer(s, y) H V H V.T rho * np.outer(s, s) x, g x_new, g_new return x, it, g if __name__ __main__: x0 np.array([-1.2, 1.0]) x_opt, n_iter, g_end bfgs(rosen, rosen_grad, x0) print(最优解:, x_opt) print(迭代次数:, n_iter) print(终止梯度范数:, np.linalg.norm(g_end, np.inf))从 ((-1.2, 1.0)) 出发用上面代码跑Rosenbrock函数通常三四十次迭代内就能让梯度无穷范数降到 (10^{-6}) 量级。作为对比纯最速下降法在这个问题上可能需要几千甚至上万次差距非常明显。这也印证了BFGS具备超线性收敛特性——虽然不如牛顿法的二次收敛快但每步代价低得多综合效益极高。代码里有几个细节值得说明我直接调用了scipy.optimize.line_search它内部实现了强Wolfe条件不是简单回溯。这是有意为之因为Wolfe条件的第二个不等式——曲率条件——正好保证 (y^T s 0)也就是 (\rho) 的分母为正。没有这个保证BFGS公式更新出来的 (H) 很容易失去正定性。更新顺序是先算 (s) 再算 (y)顺序别反。很多初学实现出错是因为拿旧梯度去算新梯度差结果 (s) 和 (y) 不是同一步的对应量。这里用的是 (H) 更新形式所以代码里看不到任何线性方程组求解。这就是拟牛顿法工程上的核心优势。6. 实现避坑数值梯度、Wolfe条件与H0初始化6.1 数值梯度的步长怎么选如果目标函数太复杂无法手推梯度很多人会选择有限差分[ \frac{\partial f}{\partial x_i} \approx \frac{f(x h e_i) - f(x - h e_i)}{2h} ]中心差分精度是 (O(h^2))但 (h) 不是越小越好。(h) 太小浮点舍入误差会淹没真实差分值(h) 太大截断误差变大。以Rosenbrock这类二次程度较高的函数为例我实测过 (h 10^{-8}) 时数值梯度经常抖动收敛稳定性明显变差取 (h 10^{-6}) 或 (10^{-5}) 反而更稳。如果你的问题解析梯度可行就优先用解析梯度。数值梯度只适合做单元测试和验证不适合长时间跑优化。6.2 Wolfe条件与 (y^T s) 的正定性BFGS推导里最容易被忽略的是前提条件(y_k^T s_k 0) 必须成立。否则 (\rho) 为负更新方向就乱了。Wolfe条件第二条曲率条件正是干这个的[ g_{k1}^T d_k \ge c_2 , g_k^T d_k,\qquad c_2 \in (0,1) ]因为 (g_k^T d_k 0)上式保证了迭代点走出足够远使得沿 (d_k) 方向的导数有显著改变。把它代入[ y_k^T s_k (g_{k1} - g_k)^T (\alpha_k d_k) ]可以得到 (y_k^T s_k 0)。这就是为什么实现BFGS时不能随便用一个“能下降就行”的简单线搜索。如果线搜索不满足曲率条件BFGS更新很快会失去正定性后面全乱套。我自己的经验第一次写BFGS时图省事用了Armijo回溯前几步凑合到后来经常出现非正定更新只能靠重置 (HI) 硬撑。后来老老实实换成带曲率条件的线搜索问题就消失了。所以别省这一步。6.3 H0的缩放和矩阵对称性磨损初始矩阵 (H_0) 的选择直接影响前几次迭代效率。最省事的是取单位阵 (I)但遇到曲率尺度差异很大的问题时收敛很慢。工程上更推荐用Barzilai-Borwein式缩放[ H_0 \gamma I,\qquad \gamma \frac{s_0^T y_0}{y_0^T y_0} ]也就是先走一步用第一步的 (s_0, y_0) 估一个整体缩放因子再继续迭代。很多算法库默认这么干收敛速度提升明显。另外一个隐蔽问题是数值对称性磨损。理论上 (H_{k1}) 是对称的但浮点运算是近似过程迭代几十步后 (H) 可能慢慢不那么对称最终影响搜索方向。保险起见我一般每隔一段时间强制对称化一次H 0.5 * (H H.T)如果发现更新后的 (y^T s \le 0)更稳妥的做法是直接跳过本轮更新保留旧 (H)相当于这一步退回最速下降方向。这比强行继续更新更安全。6.4 什么时候该重置H大规模L-BFGS里通常不轻易重置近似的逆Hessian因为历史信息很宝贵。但小规模标准BFGS里如果目标函数有强非凸性或者问题本身在迭代中切换了“曲率模式”旧矩阵反而成为包袱。此时一个简单策略是当梯度的下降方向与负梯度方向夹角过大时把 (H) 重置为单位阵或 (\gamma I)重新积累曲率信息。实测下来这个技巧对很多病态问题能救急。最后再分享一点个人体会。学拟牛顿法最关键的不是背BFGS的公式而是理解割线方程只是一个“方向上的插值”。Hessian有 (n(n1)/2) 个自由度梯度差只给了n个约束剩下所有问题的本质都是同一个在满足这条插值关系的前提下如何尽量沿用已经积累的曲率信息。SR1、DFP、BFGS、Broyden族本质全是这个问题的不同答案。抓住这一点再看那些公式会轻松很多。下一份笔记我打算写写L-BFGS的two-loop recursion——在大规模问题上真正发力的那个版本。