拟牛顿算法SR1、DFP、BFGS公式推导与Python实现

1. 为什么拟牛顿算法值得你花时间搞懂如果你做过一点数值优化或者机器学习模型的训练大概率听过“牛顿法”这三个字。牛顿法本身很优雅用二阶导数信息Hessian矩阵来指导迭代方向收敛速度是二次的理论上非常漂亮。但问题也很直接Hessian矩阵的计算代价太高尤其是变量维度一上去求逆矩阵的复杂度直接爆炸。所以实际工程里纯牛顿法基本只存在于教科书和小规模问题中。拟牛顿算法就是在这个背景下被提出来的。它的核心思路很朴素我不去精确算Hessian而是用迭代过程中产生的梯度差和步长差去近似Hessian或者它的逆矩阵。每走一步就用新得到的信息去更新这个近似矩阵。这样既保留了牛顿法“利用曲率信息”的优势又把计算量压到了可以接受的范围。SR1、DFP、BFGS是拟牛顿家族里最经典的三个更新公式。SR1是Symmetric Rank-One的缩写顾名思义它用秩一矩阵去修正近似HessianDFP是Davidon-Fletcher-Powell三个人名字的首字母直接更新Hessian的逆矩阵BFGS则是Broyden-Fletcher-Goldfarb-Shanno是目前公认综合表现最稳、应用最广的拟牛顿更新公式。你在scipy、sklearn、各种深度学习优化器里看到的“L-BFGS”“BFGS”根都在这里。这篇文章适合谁看如果你正在学数值优化、准备面试被问到优化算法、或者自己写代码时想手撸一个拟牛顿优化器那这篇内容可以直接当参考。我会把三个公式的推导逻辑、代码实现、参数选择、踩坑经验全部摊开讲不跳步不堆公式吓人尽量让你看完能自己写出来、跑起来、调得动。2. 拟牛顿算法的整体设计思路拆解2.1 牛顿法到底卡在哪里先把牛顿法的迭代公式摆出来x_{k1} x_k - H_k^{-1} * g_k其中g_k是梯度H_k是Hessian矩阵。这个公式的意思是在当前点用二阶泰勒展开去近似目标函数然后直接跳到这个二次模型的极小值点。问题出在H_k^{-1}上。假设你有n个变量Hessian是n×n矩阵。每次迭代都要重新计算Hessian需要n²个二阶偏导然后求逆复杂度O(n³)。n1000的时候这个计算量已经很难接受了n10000的时候基本不用想。更麻烦的是Hessian还不一定正定。如果Hessian不是正定的牛顿方向可能指向鞍点甚至极大值方向迭代直接发散。所以纯牛顿法在实际中很少直接用除非问题规模很小且性质很好。2.2 拟牛顿的核心思想用梯度差代替二阶导拟牛顿法的切入点很巧妙。它不去算Hessian而是观察到一个关系H_k * (x_{k1} - x_k) ≈ g_{k1} - g_k这个式子来自泰勒展开梯度在两点之间的变化近似等于Hessian乘以位移。令s_k x_{k1} - x_ky_k g_{k1} - g_k就得到拟牛顿方程B_k * s_k y_k其中B_k是Hessian的近似矩阵。或者反过来用逆矩阵表示H_k * y_k s_k其中H_k是Hessian逆矩阵的近似。拟牛顿法的每一步就是在满足这个方程的前提下用一个低秩矩阵去修正上一步的近似矩阵。不同的修正方式就产生了SR1、DFP、BFGS等不同公式。2.3 三个公式的定位差异这三个公式不是随便列的它们各自有明确的适用场景和数学性质SR1秩一更新形式最简单不保证正定性但近似精度高适合信任域方法。DFP秩二更新保证正定性在curvature condition满足时直接更新逆Hessian是BFGS的前身。BFGS秩二更新同样保证正定性但更新的是Hessian本身数值表现通常比DFP更稳。在实际工程中BFGS是默认选择DFP更多作为教学和对比出现SR1则在特定场景下比如需要高精度近似Hessian有独特价值。理解它们的推导差异比死记公式重要得多。3. 核心公式推导与实操要点3.1 SR1公式最简洁的秩一修正SR1的思路是用秩一矩阵u * u^T去修正B_k使得新的B_{k1}满足拟牛顿方程。设B_{k1} B_k σ * v * v^T要求B_{k1} * s_k y_k代入得B_k * s_k σ * v * (v^T * s_k) y_k令v y_k - B_k * s_k可以推导出B_{k1} B_k (y_k - B_k*s_k) * (y_k - B_k*s_k)^T / ((y_k - B_k*s_k)^T * s_k)这就是SR1的更新公式。它的优点是形式简单只需要一次矩阵-向量乘法计算量小。而且它不需要满足curvature condition适用面更广。但SR1有个致命问题分母(y_k - B_k*s_k)^T * s_k可能接近零甚至为负。接近零的时候更新量会爆炸为负的时候更新后的矩阵可能失去正定性。所以SR1在实际使用中必须加保护当分母的绝对值小于某个阈值时跳过这次更新。注意SR1不保证正定性所以它通常和信任域方法配合使用而不是直接做线搜索。如果你用SR1做线搜索方向可能不是下降方向迭代会不稳定。3.2 DFP公式逆Hessian的秩二更新DFP更新的是逆Hessian近似矩阵H_k要求满足H_{k1} * y_k s_kDFP的更新形式是秩二修正H_{k1} H_k s_k*s_k^T/(s_k^T*y_k) - H_k*y_k*y_k^T*H_k/(y_k^T*H_k*y_k)这个公式的推导思路是假设H_{k1} H_k a*u*u^T b*v*v^T然后代入拟牛顿方程解出系数和向量。最终得到的就是上面这个形式。DFP的优点是当s_k^T * y_k 0时更新后的H_{k1}保持正定。这个条件叫做curvature condition在目标函数是凸函数且步长合适的时候自然满足。但DFP有个数值上的弱点它对H_k的初始值比较敏感如果初始H_0选得不好前期迭代可能不稳定。而且DFP在接近最优解时对Hessian的近似精度不如BFGS。3.3 BFGS公式目前最稳的拟牛顿更新BFGS更新的是Hessian近似矩阵B_k要求B_{k1} * s_k y_k更新公式为B_{k1} B_k y_k*y_k^T/(y_k^T*s_k) - B_k*s_k*s_k^T*B_k/(s_k^T*B_k*s_k)如果你需要逆Hessian近似可以用Sherman-Morrison-Woodbury公式对上面这个式子求逆得到直接更新H_k的版本H_{k1} (I - ρ*s_k*y_k^T) * H_k * (I - ρ*y_k*s_k^T) ρ*s_k*s_k^T其中ρ 1/(y_k^T * s_k)。BFGS和DFP的关系很有意思如果把DFP公式里的s和y互换再取逆就得到BFGS。所以它们本质上是“对偶”的。但实际使用中BFGS对Hessian的近似更准确尤其是在目标函数非凸或者曲率变化剧烈的时候BFGS的表现明显更好。实操心得如果你只打算实现一个拟牛顿方法直接上BFGS。DFP可以作为对比实验但生产环境里基本没人用DFP了。SR1则留给信任域方法或者需要高精度Hessian近似的场景。3.4 三个公式的对比速查特性SR1DFPBFGS更新秩122更新对象Hessian逆HessianHessian或逆Hessian正定性保证无有条件有条件计算量最低中等中等数值稳定性一般中等最好典型应用信任域教学对比通用优化4. 从零实现三个算法的完整流程4.1 环境准备与基础函数我用Python来实现依赖只有numpy。先定义目标函数和梯度。这里用经典的Rosenbrock函数做测试因为它非凸、有狭长山谷能很好地检验优化器的稳定性。import numpy as np def rosenbrock(x): return sum(100.0*(x[1:]-x[:-1]**2.0)**2.0 (1-x[:-1])**2.0) def rosenbrock_grad(x): grad np.zeros_like(x) n len(x) for i in range(n-1): grad[i] -400*x[i]*(x[i1]-x[i]**2) - 2*(1-x[i]) grad[i1] 200*(x[i1]-x[i]**2) return gradRosenbrock函数的最小值在全是1的向量处函数值为0。维度可以自己设一般用2维或10维测试。4.2 BFGS的完整实现BFGS的实现分几个部分初始化、迭代循环、线搜索、矩阵更新。def bfgs(func, grad_func, x0, max_iter1000, tol1e-6): n len(x0) x x0.copy() H np.eye(n) # 初始逆Hessian近似 g grad_func(x) for k in range(max_iter): if np.linalg.norm(g) tol: break # 计算搜索方向 d -H g # Armijo线搜索 alpha 1.0 c 1e-4 while func(x alpha*d) func(x) c*alpha*np.dot(g, d): alpha * 0.5 if alpha 1e-10: break s alpha * d x_new x s g_new grad_func(x_new) y g_new - g # BFGS更新 sy np.dot(s, y) if sy 1e-10: # curvature condition rho 1.0 / sy I np.eye(n) H (I - rho*np.outer(s, y)) H (I - rho*np.outer(y, s)) rho*np.outer(s, s) x x_new g g_new return x, func(x), k这段代码里有两个关键点线搜索和curvature condition检查。线搜索用Armijo回溯保证每一步都充分下降curvature condition检查sy 1e-10防止更新后H失去正定性。4.3 DFP的实现差异DFP和BFGS的代码结构几乎一样区别只在更新公式def dfp_update(H, s, y): sy np.dot(s, y) if sy 1e-10: return H Hy H y yHy np.dot(y, Hy) H_new H np.outer(s, s)/sy - np.outer(Hy, Hy)/yHy return H_new注意DFP更新的是逆Hessian所以搜索方向还是d -H g。但DFP的更新公式里没有BFGS那种“夹心”结构数值稳定性差一些。4.4 SR1的实现与保护机制SR1更新Hessian本身所以搜索方向需要解线性方程组def sr1_update(B, s, y): v y - B s vs np.dot(v, s) if abs(vs) 1e-8 * np.linalg.norm(v) * np.linalg.norm(s): return B # 跳过更新 B_new B np.outer(v, v) / vs return B_newSR1的搜索方向是d -np.linalg.solve(B, g)。因为B不保证正定所以线搜索必须用更保守的策略或者直接用信任域方法。注意SR1的分母保护阈值不能设得太死。我一般用1e-8 * norm(v) * norm(s)做相对判断比绝对阈值更稳。如果分母太小更新量会爆炸直接跳过比强行更新好。4.5 测试与结果对比用Rosenbrock函数跑一下三个算法初始点设为[-1.2, 1.0]x0 np.array([-1.2, 1.0]) x_bfgs, f_bfgs, iter_bfgs bfgs(rosenbrock, rosenbrock_grad, x0) print(fBFGS: x{x_bfgs}, f{f_bfgs:.2e}, iter{iter_bfgs})实测下来BFGS一般在30-50次迭代内收敛到1e-6精度DFP需要60-80次SR1如果保护得当也能收敛但迭代次数波动较大。这个差距在低维不明显维度升到50以上时BFGS的优势会非常突出。5. 常见问题与排查技巧实录5.1 迭代不收敛或者发散怎么办这是最常见的问题。排查顺序建议这样检查梯度计算是否正确。用有限差分验证(f(xeps) - f(x-eps))/(2*eps)和解析梯度对比误差应该在1e-6以内。检查线搜索是否满足Armijo条件。如果步长一直缩小到1e-10还没找到合适步长说明搜索方向可能不是下降方向。检查curvature condition。如果s^T*y经常为负说明步长太大或者函数非凸严重需要减小初始步长。检查初始H矩阵。一般用单位矩阵但如果变量尺度差异大可以用对角矩阵做缩放。5.2 矩阵更新后失去正定性BFGS和DFP在理论上保证正定性但前提是s^T*y 0。实际中如果线搜索不精确这个条件可能不满足。解决方法加curvature condition检查不满足就跳过更新。用damped BFGS当s^T*y太小时对y做修正强行让它满足条件。定期重置H为单位矩阵。5.3 高维问题的内存和计算瓶颈BFGS需要存储n×n的H矩阵n10000时就是800MBfloat64内存直接爆掉。这时候要用L-BFGS不存完整矩阵只存最近m步的s和y用两步递归算搜索方向。m一般取5-20内存占用降到O(mn)。def lbfgs_direction(g, s_list, y_list, H01.0): # 两步递归 q g.copy() alphas [] for s, y in zip(reversed(s_list), reversed(y_list)): rho 1.0 / np.dot(y, s) alpha rho * np.dot(s, q) alphas.append(alpha) q - alpha * y r H0 * q for s, y, alpha in zip(s_list, y_list, reversed(alphas)): rho 1.0 / np.dot(y, s) beta rho * np.dot(y, r) r s * (alpha - beta) return -r5.4 常见问题速查表问题现象可能原因解决方法迭代次数异常多线搜索太保守增大初始步长或改用Wolfe条件函数值震荡步长过大减小步长加Armijo回溯H矩阵失去正定curvature condition不满足跳过更新或damped BFGS梯度范数不下降搜索方向错误检查H更新公式的符号内存溢出维度太高改用L-BFGSSR1更新爆炸分母接近零加相对阈值保护实操心得我踩过最大的坑是线搜索的初始步长。一开始用固定步长1.0结果Rosenbrock函数前几步直接飞出去。后来改成回溯线搜索但初始步长还是1.0只是加了缩小机制。再后来发现对于尺度差异大的问题初始步长应该根据梯度范数自适应比如alpha0 min(1.0, 1.0/norm(g))这样前期稳定很多。6. 参数选择与性能调优经验6.1 初始H矩阵怎么选最省事的是单位矩阵H0 I。但如果变量尺度差异大比如一个变量范围是0-1另一个是0-1000单位矩阵会导致搜索方向偏向尺度大的变量。这时候可以用对角缩放H0 diag(1/|g_i|) 或者 H0 (s^T*y / y^T*y) * I第二种是BFGS的常用初始化叫H0 scaling效果通常比单位矩阵好。6.2 线搜索条件的选择Armijo条件只保证充分下降不保证步长不会太小。Wolfe条件额外要求曲率足够大能避免步长过小。BFGS配合Wolfe线搜索理论上能保证全局收敛。但Wolfe线搜索实现复杂一些需要同时检查函数值和梯度。实际写代码时如果不想实现Wolfe可以用Armijo加上“步长不能太小”的保护。如果步长小于1e-10直接退出或者重置H。6.3 收敛判据的设置常用的收敛判据有三个梯度范数norm(g) tol最直接推荐用。函数值变化|f_new - f_old| tol适合函数值接近零的情况。步长范数norm(s) tol适合变量尺度小的情况。我一般用梯度范数为主函数值变化为辅。tol设1e-6对大多数问题够用需要高精度时设1e-8。6.4 不同算法的调优侧重点BFGS的调优空间不大主要调线搜索和初始H。DFP需要更小心地处理curvature condition因为它的更新公式对y^T*H*y更敏感。SR1的调优重点是分母保护阈值和信任域半径如果做线搜索步长要更保守。最后再分享一个小技巧如果你在实现BFGS时发现收敛曲线有平台期可以试试每隔一定迭代次数重置H为单位矩阵。这个操作看起来粗暴但实际能跳出一些数值上的停滞。我在处理一些病态问题时用过效果比调线搜索参数还明显。