1. 从牛顿法到拟牛顿法为什么需要这条演进路线很多人第一次接触优化算法都是从梯度下降开始的。梯度下降简单、直观沿着梯度的反方向走一步步长靠学习率控制。但用久了就会发现一个问题它在不同方向上的收敛速度差异极大遇到病态条件的函数比如一个方向陡峭、另一个方向平坦的峡谷形函数梯度下降会来回震荡收敛慢得让人抓狂。牛顿法就是来解决这个问题的。它的核心思想是利用二阶导数信息——海森矩阵Hessian Matrix来修正梯度方向。具体来说牛顿法的迭代公式是x_{k1} x_k - H_k^{-1} * g_k其中H_k是当前点的海森矩阵g_k是梯度。这个公式的物理含义是用二次函数去近似当前点的局部曲面然后直接跳到这个二次函数的最小值点。对于正定的二次函数牛顿法一步就能收敛到最优解这个性质非常漂亮。但牛顿法的问题也很致命。第一计算海森矩阵本身代价很高对于 n 维问题海森矩阵有 n² 个元素每次迭代都要计算这些二阶偏导数。第二更麻烦的是求逆——海森矩阵求逆的时间复杂度是 O(n³)当 n 达到几千甚至上万时这个计算量完全不可接受。第三牛顿法要求海森矩阵正定否则下降方向可能不是下降方向算法会发散。所以拟牛顿法Quasi-Newton Methods的思路就很自然了我不直接计算海森矩阵而是用一个近似矩阵 B 来替代它并且这个近似矩阵在每次迭代中通过梯度信息来更新。这样既保留了牛顿法利用曲率信息的优势又避免了计算海森矩阵和求逆的高昂代价。拟牛顿法的通用框架是这样的初始化一个近似海森矩阵B_0通常取单位矩阵 I计算当前梯度g_k求解方向d_k -B_k^{-1} * g_k沿方向d_k做线搜索找到合适的步长α_k更新x_{k1} x_k α_k * d_k计算梯度差y_k g_{k1} - g_k位移s_k x_{k1} - x_k根据某种更新公式从B_k得到B_{k1}重复直到收敛这个框架里最关键的一步就是第 7 步——如何从B_k更新到B_{k1}。不同的更新公式就产生了不同的拟牛顿算法SR1、DFP、BFGS以及后来大名鼎鼎的 L-BFGS。它们之间的区别本质上就是用了不同的方式去逼近海森矩阵。这里有一个核心约束条件叫做拟牛顿方程Secant EquationB_{k1} * s_k y_k这个方程的来源是我们对梯度做泰勒展开g_{k1} ≈ g_k H_k * s_k所以H_k * s_k ≈ y_k。拟牛顿法要求近似矩阵也满足这个关系。但注意这个方程只给了 n 个约束而B_{k1}有 n² 个未知量所以解不唯一。不同的算法就是在这个解空间里用不同的策略选出一个合适的B_{k1}。理解了这一点后面看 SR1、DFP、BFGS 的公式就不会觉得是凭空冒出来的了——它们都是在满足拟牛顿方程的前提下用不同的最小改动原则推导出来的。2. SR1 更新最直观但最不稳定的那个SR1Symmetric Rank-One是三种算法里形式最简单的一个。它的名字就说明了特点每次更新时对矩阵的修正是一个秩一矩阵Rank-One Matrix而且保持对称性。SR1 的更新公式长这样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)如果用v_k y_k - B_k * s_k来简写公式就是B_{k1} B_k v_k * v_k^T / (v_k^T * s_k)这个公式的推导逻辑其实很直接。我们想找一个对称的秩一修正v * v^T使得B_{k1}满足拟牛顿方程。设B_{k1} B_k σ * v * v^T代入B_{k1} * s_k y_k经过推导就能得到上面的形式。SR1 最大的优点是它对海森矩阵的近似可以非常准确。因为它没有强制要求B_{k1}保持正定所以它能捕捉到海森矩阵的负曲率信息——这在某些非凸优化问题里非常重要。比如在鞍点附近海森矩阵有负特征值SR1 能正确反映这个信息而 BFGS 和 DFP 因为强制正定反而会给出错误的方向。但 SR1 的缺点同样突出分母可能为零或接近零。当v_k^T * s_k很小的时候更新量会爆炸导致数值不稳定。更糟糕的是即使分母不为零更新后的B_{k1}也可能不是正定的这意味着下降方向d_k -B_k^{-1} * g_k可能不是下降方向算法会直接跑飞。我在实际使用中的体会是SR1 很少单独使用更多是作为其他算法的补充。比如在信赖域方法里SR1 可以用来修正近似矩阵因为信赖域本身限制了步长即使方向不好也不会走太远。另外有些实现会在 SR1 更新后检查B_{k1}的正定性如果不正定就跳过这次更新或者回退到 BFGS 更新。还有一个细节值得注意SR1 的更新公式里如果v_k^T * s_k的绝对值小于某个阈值比如1e-8 * ||v_k|| * ||s_k||通常就直接跳过更新。这个阈值的选择需要根据问题的尺度来调整太小了起不到保护作用太大了会丢失有用的曲率信息。3. DFP 更新第一个实用的拟牛顿算法DFP 算法以 Davidon、Fletcher、Powell 三人的名字命名是历史上第一个被广泛认可的拟牛顿算法。它的更新公式比 SR1 复杂一些但保证了正定性。DFP 的更新公式针对海森矩阵的逆矩阵H_k B_k^{-1}是H_{k1} H_k - (H_k * y_k * y_k^T * H_k) / (y_k^T * H_k * y_k) (s_k * s_k^T) / (y_k^T * s_k)如果直接对B_k更新公式是B_{k1} (I - (y_k * s_k^T) / (y_k^T * s_k)) * B_k * (I - (s_k * y_k^T) / (y_k^T * s_k)) (y_k * y_k^T) / (y_k^T * s_k)这个公式看起来复杂但结构其实很有规律。它由两部分组成第一部分是对B_k做相似变换相当于在保持特征值不变的情况下旋转矩阵第二部分是加上一个秩一矩阵用来修正曲率信息。DFP 的关键性质是如果B_k正定且y_k^T * s_k 0那么B_{k1}也正定。这个性质保证了下降方向始终是下降方向算法不会跑飞。y_k^T * s_k 0这个条件在函数是凸函数且线搜索满足 Wolfe 条件时是自动满足的。但 DFP 有一个在实际中很要命的问题它对海森矩阵的近似容易变得病态。特别是在高维问题里H_k的条件数会越来越大导致数值不稳定。我试过在一个 1000 维的二次规划问题上用 DFP迭代到后面H_k的条件数超过了 1e12方向计算完全失去了精度。DFP 的另一个问题是它没有 BFGS 那样的自校正机制。BFGS 在数值误差导致B_k偏离真实海森矩阵时能够自动修正回来而 DFP 一旦偏离就很难恢复。这也是为什么现在几乎没有人单独用 DFP 了——它更多是作为理解 BFGS 的铺垫以及在某些特定场景下作为 BFGS 的补充。不过DFP 也不是完全没有价值。在一些低维问题比如 n 50上DFP 的表现和 BFGS 差不多而且计算量略小。另外DFP 和 BFGS 可以组合使用比如在迭代初期用 DFP后期切换到 BFGS利用两者的互补性。4. BFGS 更新拟牛顿法的实际标准BFGS 以 Broyden、Fletcher、Goldfarb、Shanno 四人的名字命名是目前公认最有效的拟牛顿算法。它的更新公式和 DFP 很像但把s_k和y_k的角色对调了。BFGS 对B_k的更新公式是B_{k1} B_k - (B_k * s_k * s_k^T * B_k) / (s_k^T * B_k * s_k) (y_k * y_k^T) / (y_k^T * s_k)对逆矩阵H_k的更新公式是H_{k1} (I - (s_k * y_k^T) / (y_k^T * s_k)) * H_k * (I - (y_k * s_k^T) / (y_k^T * s_k)) (s_k * s_k^T) / (y_k^T * s_k)这个公式和 DFP 的逆更新公式结构完全一样只是s_k和y_k互换了位置。但就是这个互换带来了本质的区别。BFGS 的核心优势在于它的自校正性质。当数值误差导致B_k偏离真实海森矩阵时BFGS 的更新会倾向于把B_k拉回来。这个性质在理论上被称为BFGS 的收敛性保证在实际中表现为即使初始B_0选得很差比如单位矩阵BFGS 也能在若干次迭代后逼近真实的海森矩阵。另一个关键点是 BFGS 对线搜索的鲁棒性。DFP 对线搜索的精度要求很高如果步长选得不好y_k^T * s_k可能接近零甚至为负导致更新失败。BFGS 对这个问题不那么敏感即使线搜索不够精确它也能保持较好的性能。这也是为什么在实际实现中BFGS 通常搭配 Armijo 线搜索或 Wolfe 线搜索而不是精确线搜索。我在实际项目里用 BFGS 的经验是对于中小规模问题n 1000BFGS 几乎总是首选。它的收敛速度通常是超线性的比梯度下降快一到两个数量级。对于大规模问题BFGS 的内存开销需要存储 n×n 的矩阵会成为瓶颈这时候就要用 L-BFGS 了。BFGS 还有一个变种叫 BFGS-BBounded BFGS用来处理带边界约束的优化问题。它的核心思想是在更新B_k时只考虑那些不在边界上的变量把问题降维到自由变量空间。这个变种在工程优化里用得很多比如结构设计、参数拟合等场景。5. 三种算法的对比与选型建议把 SR1、DFP、BFGS 放在一起对比能更清楚地看到它们各自的定位。特性SR1DFPBFGS更新秩数秩一秩二秩二保持正定否是是自校正无弱强数值稳定性差中好对线搜索敏感度高高低适用场景信赖域、非凸问题低维问题、教学通用优化、实际首选从表格里能看出来BFGS 在几乎所有维度上都优于 DFP这也是为什么现在的优化库比如 scipy.optimize、NLopt、Ceres Solver默认都用 BFGS 或 L-BFGS而 DFP 基本只出现在教科书里。SR1 的定位比较特殊。它不适合作为主算法但在信赖域框架里作为辅助更新很有价值。比如在 trust-region 方法里如果 SR1 更新后的B_{k1}能保持正定就用 SR1否则回退到 BFGS。这种混合策略在一些高级优化器里有实现。选型的时候我一般按这个逻辑走如果问题是凸的、维度不高n 500直接用 BFGS搭配 Wolfe 线搜索。如果问题是非凸的或者有鞍点考虑用 SR1 作为补充或者用信赖域方法。如果维度很高n 10000用 L-BFGS只存储最近 m 步的s_k和y_k。如果有边界约束用 L-BFGS-B 或 BFGS-B。DFP 基本不用除非是在教学场景或者需要和 BFGS 做对比实验。还有一个实际中容易忽略的点初始矩阵B_0的选择。大多数实现默认用单位矩阵但这不一定最优。如果知道问题的尺度信息可以用一个对角矩阵来缩放比如B_0 (y_0^T * s_0) / (y_0^T * y_0) * I。这个缩放能显著改善条件数减少迭代次数。我在一个参数拟合问题里试过用缩放后的B_0比单位矩阵少了将近 30% 的迭代。6. 手写实现中的关键细节与踩坑记录如果你打算自己实现一遍这三种算法有几个细节是文档里不会写、但实际会坑死人的。第一个坑y_k^T * s_k的符号检查。在 BFGS 和 DFP 里如果y_k^T * s_k 0更新公式的分母会出问题而且正定性也无法保证。这个情况在非凸问题里很常见。我的处理方式是如果y_k^T * s_k 1e-10就跳过这次更新保持B_k不变。虽然这会损失一些曲率信息但比让算法跑飞要好。第二个坑矩阵求逆的数值精度。虽然 BFGS 可以直接更新逆矩阵H_k避免了显式求逆但H_k在多次更新后可能失去正定性由于浮点误差累积。我的做法是每隔一定迭代次数比如 50 次用当前的B_k重新计算H_k B_k^{-1}或者直接用 Cholesky 分解来保证正定性。第三个坑线搜索的精度。拟牛顿法对线搜索的精度要求比梯度下降高。如果线搜索太粗糙y_k的精度不够更新公式就会引入很大的误差。我一般用 Wolfe 条件Armijo 条件 曲率条件参数取c1 1e-4c2 0.9。这个组合在大多数问题上表现稳定。第四个坑内存布局。如果你用 Python 的 numpy 实现注意B_k的存储方式。B_k是对称矩阵理论上只需要存一半但 numpy 没有原生的对称矩阵类型。我试过用scipy.linalg.blas的对称矩阵乘法来加速效果不错但代码复杂度会上升。对于 n 1000 的问题直接用完整矩阵就行没必要优化。第五个坑收敛判据。很多人只用梯度范数||g_k|| ε作为收敛条件但这在病态问题里可能过早停止。我一般同时检查三个条件梯度范数、步长||s_k||、以及函数值的变化|f_{k1} - f_k|。三个条件都满足才认为收敛。阈值的选择取决于问题的尺度我通常用1e-6作为梯度范数的阈值1e-10作为函数值变化的阈值。下面是一个简化的 BFGS 实现框架用 Python 写展示了核心逻辑import numpy as np def bfgs(f, grad_f, x0, max_iter1000, tol1e-6): n len(x0) x x0.copy() B np.eye(n) # 初始近似海森矩阵 g grad_f(x) for k in range(max_iter): if np.linalg.norm(g) tol: break # 计算方向 d -np.linalg.solve(B, g) # Wolfe 线搜索 alpha line_search_wolfe(f, grad_f, x, d, g) # 更新 s alpha * d x_new x s g_new grad_f(x_new) y g_new - g # 检查曲率条件 ys y s if ys 1e-10: # BFGS 更新 Bs B s B B - np.outer(Bs, Bs) / (s Bs) np.outer(y, y) / ys x, g x_new, g_new return x, f(x)这个实现里line_search_wolfe需要自己实现核心是满足 Armijo 条件和曲率条件。np.linalg.solve(B, g)比直接求逆np.linalg.inv(B) g更稳定也更高效。如果你要实现 SR1把更新部分换成v y - B s vs v s if abs(vs) 1e-8 * np.linalg.norm(v) * np.linalg.norm(s): B B np.outer(v, v) / vs注意 SR1 不需要检查ys 0但需要检查vs的大小避免除以接近零的数。DFP 的更新则是By B y yBy y By ys y s if ys 1e-10: B B - np.outer(By, By) / yBy np.outer(y, y) / ys注意 DFP 的公式里第一项的分母是y^T * B * y不是s^T * B * s。这个细节很容易写错写错之后算法可能还能跑但收敛速度会差很多。7. 从 BFGS 到 L-BFGS大规模问题的出路BFGS 虽然好用但它的内存开销是 O(n²)对于 n 100000 的问题光存储B_k就需要 80GB 内存双精度完全不现实。L-BFGSLimited-memory BFGS就是来解决这个问题的。L-BFGS 的核心思想是不存储完整的B_k而是存储最近 m 步的s_k和y_k通常 m 取 5 到 20然后用这些向量来隐式地表示B_k。计算方向d_k -B_k^{-1} * g_k时通过一个两循环递归Two-Loop Recursion来完成不需要显式构造矩阵。两循环递归的过程是这样的q g_k for i k-1 down to k-m: rho_i 1 / (y_i^T * s_i) alpha_i rho_i * s_i^T * q q q - alpha_i * y_i r H_0 * q # H_0 通常取 (s_{k-1}^T * y_{k-1}) / (y_{k-1}^T * y_{k-1}) * I for i k-m to k-1: beta rho_i * y_i^T * r r r s_i * (alpha_i - beta) d_k -r这个递归的计算量是 O(mn)内存开销是 O(mn)对于大规模问题非常友好。我试过在一个 50000 维的逻辑回归问题上用 L-BFGSm 取 10内存占用不到 10MB收敛速度比随机梯度下降快得多。L-BFGS 的另一个优势是它天然适合分布式计算。因为两循环递归只涉及向量运算可以很容易地并行化。Spark 的 MLlib 和 TensorFlow 的优化器里都有 L-BFGS 的实现就是看中了这一点。不过 L-BFGS 也有它的局限。它丢失了完整 BFGS 的自校正性质在病态问题上的表现可能不如完整 BFGS。另外m 的选择需要权衡m 太小曲率信息不足收敛慢m 太大内存和计算开销增加。我的经验是 m 取 10 到 20 之间比较合适具体取决于问题的维度和条件数。还有一个实际中容易忽略的点L-BFGS 的初始矩阵H_0的选择。大多数实现用H_0 γ * I其中γ (s_{k-1}^T * y_{k-1}) / (y_{k-1}^T * y_{k-1})。这个缩放能显著改善条件数特别是在问题的尺度差异很大的时候。我在一个特征尺度差异达到 1e6 的问题上试过用缩放后的H_0比单位矩阵少了将近一半的迭代次数。8. 实际应用中的性能调优与经验总结在实际项目里用拟牛顿法光知道公式是不够的还需要根据问题的特点做调优。我总结了几条经验都是踩过坑之后才明白的。第一条预处理比算法选择更重要。如果问题的变量尺度差异很大比如一个变量在 1e-3 量级另一个在 1e3 量级那么无论用 BFGS 还是 L-BFGS条件数都会很差。这时候应该先做变量缩放把所有变量归一化到相近的尺度。我一般用(x - x_mean) / x_std来做标准化或者根据问题的物理意义手动缩放。这个步骤看起来简单但效果往往比换算法更明显。第二条线搜索的参数需要调。Wolfe 条件的c1和c2不是固定的。对于大多数问题c1 1e-4、c2 0.9是安全的。但如果函数值变化很剧烈可以把c1调小到1e-6避免步长过大。如果函数很平滑可以把c2调大到0.95让线搜索更精确。我一般会先跑一遍默认参数如果收敛慢再调。第三条注意函数的计算精度。拟牛顿法依赖梯度信息如果梯度是用有限差分算的精度损失会很大。我试过在一个问题上用有限差分梯度BFGS 迭代了 500 次还没收敛换成解析梯度后30 次就收敛了。如果必须用有限差分步长要选得合适一般取sqrt(eps) * max(1, |x_i|)其中eps是机器精度。第四条监控B_k的条件数。如果条件数超过 1e12说明B_k已经严重病态继续迭代可能没有意义。这时候可以考虑重启把B_k重置为单位矩阵或者用当前的s_k和y_k重新初始化。我在一个问题上遇到过这种情况重启之后算法又恢复了正常收敛。第五条不要忽视问题的结构。如果问题有特殊结构比如稀疏性、低秩性、或者可分性应该利用这些结构来加速。比如对于稀疏问题可以用稀疏矩阵存储B_k或者用 L-BFGS 的变种来利用稀疏性。对于可分问题可以用坐标下降或者分块更新。拟牛顿法是通用方法但通用方法不一定是最优的。最后说一个我自己的体会拟牛顿法的理论很漂亮但实际用起来80% 的时间花在调试线搜索、调整参数、处理数值问题上只有 20% 的时间在享受超线性收敛的快感。但就是这 20% 的快感让拟牛顿法成为了我工具箱里最常用的优化算法之一。如果你刚开始学建议先从 BFGS 入手把线搜索和更新公式搞明白然后再去看 SR1 和 DFP理解它们的设计动机和适用场景。这样学下来不仅知道怎么用还知道为什么这么用。