多彩编程 多彩编程MZPH · CODE BLOG
ARTICLE DETAIL

文章详情

深耕前端与后端开发技术的一线实战笔记与踩坑复盘。

拟牛顿法推导详解:从割线方程到BFGS更新公式

拟牛顿法推导详解:从割线方程到BFGS更新公式 最优化这门课学到拟牛顿法的时候很多人会有一种“咦好像懂了但让自己推一遍就卡壳”的感觉。尤其是教材里DFP、BFGS那一串公式看起来像变魔术一样不知道当初是怎么凑出来的。这篇文章就是我整理的一份拟牛顿法推导笔记把我自己踩过的坑、绕过的弯都写清楚争取做到“推导过程不跳步、每个选择都有原因”。不管你是正在准备考试的学生还是工作中需要用到优化算法的工程师这篇笔记都能帮你把这块硬骨头啃下来。1. 从牛顿法说起拟牛顿法到底在解决什么问题1.1 最速下降法和牛顿法的“两个极端”先简单回顾一下无约束优化问题min f(x)其中f是二次连续可微的。最基本的迭代格式就是x_{k1} x_k α_k * d_k区别在于怎么选搜索方向d_k。最速下降法选 d_k -∇f(x_k)也就是负梯度方向。这个方向的好处是计算简单只用到一阶导数坏处是收敛速度慢尤其在病态问题Hessian矩阵条件数很大上会走出明显的“之”字形路径往往需要很多轮迭代才能收敛。牛顿法则用二阶信息直接取 d_k -[∇²f(x_k)]⁻¹ ∇f(x_k)。这个方向相当于用二次模型去近似原函数然后直接跳到这个二次模型的极小点。如果f本身是二次函数牛顿法一步就能收敛对一般函数只要初始点不太离谱收敛速度是二阶的非常快。但牛顿法的代价也很明显每轮都要计算Hessian矩阵∇²f(x_k)然后还要对它求逆或者解一个线性方程组。当变量维数n很大时Hessian矩阵是n×n的存储就要O(n²)内存求逆更是O(n³)的计算量。在实际问题里n动辄几千几万甚至更高这就让纯牛顿法变得不现实。更麻烦的是Hessian矩阵可能根本没法解析求出来只能靠有限差分近似那噪声和计算量就更大了。1.2 拟牛顿法的核心思想用“近似”换“效率”拟牛顿法走的是中间路线不直接计算Hessian矩阵而是通过迭代过程中收集到的梯度差信息逐步构造一个对Hessian矩阵或其逆的近似。这样既保留了二阶方法收敛快的优点又避开了二阶导数计算和大规模矩阵求逆的麻烦。打个比方牛顿法像是一个每步都精确测量地形、再规划最优路径的登山者可靠但费时最速下降法像是一个只看脚下、哪陡往哪走的登山者省力但容易绕路拟牛顿法则像是边走边积累地形经验、越走越有把握的登山者兼顾了效率和稳健性。具体来说拟牛顿法维护一个矩阵B_k用来逼近∇²f(x_k)或者维护H_k用来逼近[∇²f(x_k)]⁻¹。每次迭代用这个近似矩阵构造方向d_k -H_k ∇f(x_k)然后做线搜索得到步长α_k更新x_{k1}再用(x_k, x_{k1})之间的位移和梯度变化量去修正H_k得到H_{k1}。整个过程完全不需要二阶导数。理论上如果H_k最终能收敛到真正的Hessian逆那拟牛顿法的收敛速度就会趋近于牛顿法。实际中BFGS方法在满足一定条件时具有超线性收敛性虽然赶不上牛顿法的二阶收敛但比最速下降法快得多更重要的是每轮开销小得多。2. 拟牛顿条件整个推导的地基2.1 割线方程的来历拟牛顿法构造近似矩阵时核心依据是所谓的“割线方程”也叫拟牛顿方程。推导过程很简单但理解它背后的思想是关键。设 f(x) 是二次连续可微的在点 x_{k1} 处对梯度∇f做一阶泰勒展开∇f(x_{k1}) ≈ ∇f(x_k) ∇²f(x_k) (x_{k1} - x_k)把移项整理一下。定义两个常用记号s_k x_{k1} - x_ky_k ∇f(x_{k1}) - ∇f(x_k)于是有近似关系y_k ≈ ∇²f(x_k) s_k如果把∇²f(x_k)换成它的近似矩阵B_{k1}注意这里用k1时刻的近似就得到B_{k1} s_k y_k这就是割线方程。它说的是近似矩阵B_{k1}在沿着位移方向s_k作用时应等于梯度变化y_k。这是拟牛顿法向后传播信息的基本约束。如果用的是逆矩阵近似H_{k1} ≈ [∇²f(x_{k1})]⁻¹则割线方程等价地写成H_{k1} y_k s_k注意很多教材直接把割线方程“啪”地写出来然后开始推公式。我一开始就很困惑为什么一定是这个形式其实关键在于泰勒展开的近似——它要求相邻两个迭代点足够近梯度变化才能近似为线性关系。所以线搜索步长不宜取得太大否则割线方程的“信息质量”就会变差。2.2 为什么要求对称正定Hessian矩阵本身是对称的对二次连续可微函数由克莱罗定理保证所以近似矩阵B和H也应当保持对称性。更重要的是我们要求B_k或H_k是正定的。理由是搜索方向要用到H_k乘以负梯度。如果H_k正定那么∇f(x_k)ᵀ d_k -∇f(x_k)ᵀ H_k ∇f(x_k) 0也就是说d_k一定是下降方向。这是保证算法能收敛的基础。如果H_k不正定d_k可能根本不是下降方向迭代就会出问题。所以拟牛顿更新公式的设计必须满足两个硬约束满足割线方程对B_{k1}或H_{k1}从正定的B_k出发更新后B_{k1}仍然正定。你可能已经感觉到这两个约束加在一起可选的更新公式就不是很多了。这正是接下来推导的主线。3. 主流更新公式推导从秩一校正到BFGS3.1 SR1从秩一校正看起最简单的思路是在B_k的基础上加上一个秩为1的修正项让结果满足割线方程。设B_{k1} B_k c u uᵀ其中u是一个待定向量c是待定系数。这个修正项是秩一的外积u uᵀ的秩为1。把它代入割线方程B_{k1} s_k y_kB_k s_k c u (uᵀ s_k) y_k所以 c (uᵀ s_k) u y_k - B_k s_k这说明u必须与(y_k - B_k s_k)平行。取u y_k - B_k s_k则c (uᵀ s_k) 1所以B_{k1} B_k ((y_k - B_k s_k)(y_k - B_k s_k)ᵀ) / ((y_k - B_k s_k)ᵀ s_k)这就是SR1对称秩一校正公式。形式简洁计算量也小。但SR1有个臭名昭著的问题分母 (y_k - B_k s_k)ᵀ s_k 可能接近零甚至等于零导致校正项爆炸或者无定义。而且SR1不保证保持正定性。因此尽管SR1在某些特定问题上有奇效它并不是最“稳”的选择。3.2 DFP对逆矩阵的更新既然维护逆矩阵H更方便直接用H乘以负梯度就得到方向那就在H上做文章。设H_{k1} H_k 校正项要求满足 H_{k1} y_k s_k。DFP公式的核心思路是加两项秩一校正一项用来“补”y_k方向的信息一项用来“校正”H_k y_k方向的信息。推导如下构造H_{k1} H_k a u uᵀ b v vᵀ让它满足 H_{k1} y_k s_k。于是H_k y_k a u (uᵀ y_k) b v (vᵀ y_k) s_k很自然地取第一项 u s_k第二项 v H_k y_k。为了抵消H_k y_k对结果的影响令b (vᵀ y_k) -1即 b -1 / (y_kᵀ H_k y_k)a (uᵀ y_k) 1即 a 1 / (y_kᵀ s_k)代入整理得到DFP更新公式H_{k1} H_k - (H_k y_k y_kᵀ H_k) / (y_kᵀ H_k y_k) (s_k s_kᵀ) / (y_kᵀ s_k)这个公式在历史上比BFGS更早出现。它的优点是只要y_kᵀ s_k 0线搜索满足Wolfe条件就能保证且H_k正定那么H_{k1}也正定。这在理论上非常漂亮。但实际计算中DFP对“精确”线搜索的依赖更强处理非精确线搜索时的数值稳定性不如BFGS。这也是后来BFGS逐渐成为主流的原因之一。3.3 BFGS最常用的那个BFGSBroyden-Fletcher-Goldfarb-Shanno可以看作是DFP的“对偶”DFP直接更新逆矩阵H而BFGS直接更新Hessian近似B然后通过Sherman-Morrison-Woodbury公式反推出逆矩阵的更新式。先推导B的更新。思路和DFP完全对称设B_{k1} B_k 校正项要求 B_{k1} s_k y_k。同样用两项秩一校正取一项为y_k方向用来“补”y_k一项为B_k s_k方向用来抵消B_k s_k。设B_{k1} B_k a y_k y_kᵀ b (B_k s_k)(B_k s_k)ᵀ代入割线方程B_k s_k a y_k (y_kᵀ s_k) b (B_k s_k)(s_kᵀ B_k s_k) y_k令a (y_kᵀ s_k) 1 a 1 / (y_kᵀ s_k)b (s_kᵀ B_k s_k) -1 b -1 / (s_kᵀ B_k s_k)于是B_{k1} B_k (y_k y_kᵀ) / (y_kᵀ s_k) - (B_k s_k s_kᵀ B_k) / (s_kᵀ B_k s_k)这是B的更新式。现在需要H B⁻¹的更新式。这里直接套Sherman-Morrison-Woodbury公式即矩阵求逆引理可以推出来但推导比较繁琐。我这里直接给出结果并附一个记忆技巧H_{k1} (I - ρ_k s_k y_kᵀ) H_k (I - ρ_k y_k s_kᵀ) ρ_k s_k s_kᵀ其中 ρ_k 1 / (y_kᵀ s_k)。这个形式叫“乘积形式”product form在数值实现中非常常用因为它的结构是左右夹着H_k能最大限度地保持正定性和数值稳定性。记忆小技巧把 (I - ρ s yᵀ) 看成是对H_k做了一次“投影修正”s_k s_kᵀ是补充的新信息。我在自己推导的时候是先记住H_{k1} (I - ρ s yᵀ) H_k (I - ρ y sᵀ) ρ s sᵀ再展开验证它满足H_{k1} y_k s_k。验证割线方程也很容易把H_{k1}作用到y_k上左边第一项中因为(I - ρ y_k s_kᵀ) y_k y_k - y_k(s_kᵀ y_k)/(y_kᵀ s_k) 0注意分母相同约掉了所以只剩ρ s_k (s_kᵀ y_k) s_k。干净利落。3.4 Broyden族与统一视角看完了DFP和BFGS你可能会想这两个公式是不是某种更大框架的特例没错。定义H_{k1} (1 - φ_k) H_{k1}^{DFP} φ_k H_{k1}^{BFGS}其中φ_k是任意实数就得到Broyden族。当φ_k0时退化为DFPφ_k1时退化为BFGS。Broyden族的所有成员都满足割线方程且在φ_k取值合适时也能保持正定性。在实际工程中BFGS几乎是一边倒的选择。原因很简单在非精确线搜索尤其是Wolfe条件配合下BFGS的数值表现最稳定对步长选择的敏感度低收敛也最可靠。DFP虽然在理论上同样优雅但实践中容易在Hessian近似上积累误差导致收敛变慢。L-BFGSLimited-memory BFGS更是把BFGS推广到了大规模问题不显式存储H_k只保留最近m轮的(s_k, y_k)对用它们隐式地逼近H_k作用在梯度上的结果。如果你的问题是高维的比如超过1万维请直接考虑L-BFGS。4. 实操要点写代码前必须想清楚的事4.1 完整算法流程这里我给出一份可直接实现的BFGS逆矩阵版本算法流程给定初始点x_0初始近似逆矩阵H_0 I允许误差ε。对k 0,1,2,...重复以下步骤直到停机计算梯度 g_k ∇f(x_k)。若||g_k|| ε停止。取搜索方向 d_k -H_k g_k。沿d_k做线搜索求α_k满足Wolfe条件或Armijo条件令x_{k1} x_k α_k d_k。计算 s_k x_{k1} - x_ky_k g_{k1} - g_k。计算ρ_k 1 / (y_kᵀ s_k)。更新 H_{k1} (I - ρ_k s_k y_kᵀ) H_k (I - ρ_k y_k s_kᵀ) ρ_k s_k s_kᵀ。输出x_{k1}。第4步是很多初学实现容易忽略的关键点线搜索必须满足Wolfe条件至少满足Armijo条件否则y_kᵀ s_k可能非正导致ρ_k为负H的更新就会破坏正定性算法直接“跑飞”。4.2 初始矩阵H_0怎么选H_0通常取单位矩阵I。这样第一轮迭代实际上就是最速下降法d_0 -g_0后续迭代逐步积累曲率信息逐渐“升级”为拟牛顿方向。更精细的做法是取H_0 (s_0ᵀ y_0) / (y_0ᵀ y_0) * I这个缩放因子来自对特征值的启发式估计可以让第一轮迭代的步长更合理。实测下来在目标函数各向异性比较强时这种缩放能让收敛快不少。4.3 线搜索配合Wolfe条件为什么重要我之前在实现时偷懒用了固定步长α_k 0.001结果BFGS的表现甚至不如最速下降法。为什么因为固定小步长导致s_k非常短y_kᵀ s_k趋近于零H_k的更新信息基本是噪声近似矩阵的质量自然上不去。改用Wolfe条件线搜索后情况立刻不同。Wolfe条件包含两条Armijo条件充分下降条件f(x_k α d_k) ≤ f(x_k) c₁ α ∇f(x_k)ᵀ d_k其中c₁通常取1e-4。曲率条件∇f(x_k α d_k)ᵀ d_k ≥ c₂ ∇f(x_k)ᵀ d_k其中c₂对拟牛顿法通常取0.9对共轭梯度法取0.1。曲率条件保证步长不会太小从而确保y_kᵀ s_k 0。这一条对拟牛顿法至关重要因为它是H_k保持正定性的前提。我在实践中用的比较多的是Nocedal和Wright书里的线搜索实现或者直接用SciPy的line_search函数。实际编码时有个小坑如果直接用解析梯度线搜索时要用同一个方向d_k和同一个初始步长不能每轮都重新算一个方向。所有主流优化库都是这个逻辑。5. 常见问题与调试经验速查5.1 更新后矩阵不正定怎么办这是初学者最容易碰到的问题。原因几乎只有三个步长不满足Wolfe条件的曲率条件、数值精度导致y_kᵀ s_k被算成负数实际应该是正的、目标函数本身不可导或不光滑。排查方法打印每一轮的y_kᵀ s_k值看是否始终为正。如果出现负值先用更严格的线搜索如果线搜索没问题检查梯度是否计算正确用中心差分验证一次。5.2 计算中y_kᵀ s_k过小除零风险在极值点附近梯度的变化量y_k会趋于零y_kᵀ s_k也会变得非常小。这本身不是问题——因为此时算法已经接近收敛可以在更新前加一个判断如果y_kᵀ s_k 某个阈值比如1e-12直接跳过本轮更新保持H_k不变。这样既避免了除零又不会影响收敛精度。5.3 BFGS比最速下降法还慢如果出现这种情况十有八九是初始点离最优解太远或者目标函数有严重的非光滑区域。拟牛顿法是局部优化方法对初始点敏感。一个实用的技巧是先用少量迭代的梯度下降把点“拉”到相对正常的区域再切换到BFGS俗称“热启动”。另一个可能是你的步长选的太保守被迫走了大量小碎步这就要靠线搜索的曲率条件来解决了。5.4 收敛精度遇到瓶颈BFGS在极小点附近的收敛有时会受限于近似矩阵的精度尤其是当Hessian矩阵本身病态条件数极大时BFGS的数值误差会被放大。这时可以改用L-BFGS并配合更强精度的线搜索或者做一次二阶修正——在临近收敛时计算一次真实的Hessian逆作为H的重新初始化。下面整理成一张速查表问题可能原因检查方法解决方案H更新后不正定线搜索不满足曲率条件检查y_kᵀ s_k是否为正改用满足Wolfe条件的线搜索除零或数值爆炸y_kᵀ s_k接近零打印每轮y_kᵀ s_k加阈值跳过更新收敛极慢初始点差 / 步长太小看第一轮方向是否像负梯度热启动或换初始H_0缩放目标函数不可导梯度计算错误中心差分验证梯度修正梯度实现5.5 一个容易忽略的应用细节约束优化如果你遇到的是带约束的最优化问题不能直接套BFGS。解决方案是结合惩罚函数法、增广拉格朗日法或者用序列二次规划SQP框架把BFGS用在SQP的子问题上。很多号称“BFGS求解器”的库内部也都是在处理KKT系统的近似块。这个思路与陈宝林老师《最优化理论与算法》里讲的“罚函数法 无约束优化”一脉相承——先把约束问题转化为一系列无约束问题再用拟牛顿法去解。从实际应用来看拟牛顿法尤其是L-BFGS是机器学习、数值优化、工程建模里应用最广的二阶近似方法之一。很多深度学习训练器里也内置了L-BFGS作为小批量确定性优化的一种选择。理解了BFGS的推导你再去看L-BFGS的实现就会豁然开朗——它只是把H_k隐式化不再存储n×n矩阵而是存储最近的s和y对核心更新逻辑完全一致。这段笔记写到这里我自己的体会是拟牛顿法的关键不在于背下那几个更新公式而在于理解“割线方程 对称正定约束”这两条主线。顺着这条主线DFP、BFGS、Broyden族都是自然而然的选择而不是从天而降的魔法。如果你正在学最优化或者准备复试面试强烈建议自己动手推导一遍BFGS的更新式再配合一个小维度的二次函数做数值实验观察H_k的演化。这个过程走完拟牛顿法相关的考题基本就难不住你了。
返回列表