1. 从牛顿法到拟牛顿法：为什么需要这条演进路线
很多人第一次接触优化算法，都是从梯度下降开始的。梯度下降简单、直观，沿着梯度的反方向走一步，步长靠学习率控制。但用久了就会发现一个问题：它在不同方向上的收敛速度差异极大，遇到病态条件的函数（比如一个方向陡峭、另一个方向平坦的峡谷形函数），梯度下降会来回震荡，收敛慢得让人抓狂。
牛顿法就是来解决这个问题的。它的核心思想是利用二阶导数信息——海森矩阵（Hessian Matrix），来修正梯度方向。具体来说，牛顿法的迭代公式是：
x_{k+1} = 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_{k+1} = x_k + α_k * d_k - 计算梯度差
y_k = g_{k+1} - g_k，位移s_k = x_{k+1} - x_k - 根据某种更新公式，从
B_k得到B_{k+1} - 重复直到收敛
这个框架里，最关键的一步就是第 7 步——如何从B_k更新到B_{k+1}。不同的更新公式，就产生了不同的拟牛顿算法：SR1、DFP、BFGS，以及后来大名鼎鼎的 L-BFGS。它们之间的区别，本质上就是用了不同的方式去逼近海森矩阵。
这里有一个核心约束条件，叫做拟牛顿方程（Secant Equation）：
B_{k+1} * s_k = y_k这个方程的来源是：我们对梯度做泰勒展开，g_{k+1} ≈ g_k + H_k * s_k，所以H_k * s_k ≈ y_k。拟牛顿法要求近似矩阵也满足这个关系。但注意，这个方程只给了 n 个约束，而B_{k+1}有 n² 个未知量，所以解不唯一。不同的算法就是在这个解空间里，用不同的策略选出一个合适的B_{k+1}。
理解了这一点，后面看 SR1、DFP、BFGS 的公式就不会觉得是凭空冒出来的了——它们都是在满足拟牛顿方程的前提下，用不同的"最小改动"原则推导出来的。
2. SR1 更新：最直观但最不稳定的那个
SR1（Symmetric Rank-One）是三种算法里形式最简单的一个。它的名字就说明了特点：每次更新时，对矩阵的修正是一个秩一矩阵（Rank-One Matrix），而且保持对称性。
SR1 的更新公式长这样：
B_{k+1} = 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_{k+1} = B_k + v_k * v_k^T / (v_k^T * s_k)这个公式的推导逻辑其实很直接。我们想找一个对称的秩一修正v * v^T，使得B_{k+1}满足拟牛顿方程。设B_{k+1} = B_k + σ * v * v^T，代入B_{k+1} * s_k = y_k，经过推导就能得到上面的形式。
SR1 最大的优点是：它对海森矩阵的近似可以非常准确。因为它没有强制要求B_{k+1}保持正定，所以它能捕捉到海森矩阵的负曲率信息——这在某些非凸优化问题里非常重要。比如在鞍点附近，海森矩阵有负特征值，SR1 能正确反映这个信息，而 BFGS 和 DFP 因为强制正定，反而会给出错误的方向。
但 SR1 的缺点同样突出：分母可能为零或接近零。当v_k^T * s_k很小的时候，更新量会爆炸，导致数值不稳定。更糟糕的是，即使分母不为零，更新后的B_{k+1}也可能不是正定的，这意味着下降方向d_k = -B_k^{-1} * g_k可能不是下降方向，算法会直接跑飞。
我在实际使用中的体会是：SR1 很少单独使用，更多是作为其他算法的补充。比如在信赖域方法里，SR1 可以用来修正近似矩阵，因为信赖域本身限制了步长，即使方向不好也不会走太远。另外，有些实现会在 SR1 更新后检查B_{k+1}的正定性，如果不正定就跳过这次更新，或者回退到 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_{k+1} = 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_{k+1} = (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_{k+1}也正定。这个性质保证了下降方向始终是下降方向，算法不会跑飞。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_{k+1} = 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_{k+1} = (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 < 1000），BFGS 几乎总是首选。它的收敛速度通常是超线性的，比梯度下降快一到两个数量级。对于大规模问题，BFGS 的内存开销（需要存储 n×n 的矩阵）会成为瓶颈，这时候就要用 L-BFGS 了。
BFGS 还有一个变种叫 BFGS-B（Bounded BFGS），用来处理带边界约束的优化问题。它的核心思想是在更新B_k时，只考虑那些不在边界上的变量，把问题降维到自由变量空间。这个变种在工程优化里用得很多，比如结构设计、参数拟合等场景。
5. 三种算法的对比与选型建议
把 SR1、DFP、BFGS 放在一起对比，能更清楚地看到它们各自的定位。
| 特性 | SR1 | DFP | BFGS |
|---|---|---|---|
| 更新秩数 | 秩一 | 秩二 | 秩二 |
| 保持正定 | 否 | 是 | 是 |
| 自校正 | 无 | 弱 | 强 |
| 数值稳定性 | 差 | 中 | 好 |
| 对线搜索敏感度 | 高 | 高 | 低 |
| 适用场景 | 信赖域、非凸问题 | 低维问题、教学 | 通用优化、实际首选 |
从表格里能看出来，BFGS 在几乎所有维度上都优于 DFP，这也是为什么现在的优化库（比如 scipy.optimize、NLopt、Ceres Solver）默认都用 BFGS 或 L-BFGS，而 DFP 基本只出现在教科书里。
SR1 的定位比较特殊。它不适合作为主算法，但在信赖域框架里作为辅助更新很有价值。比如在 trust-region 方法里，如果 SR1 更新后的B_{k+1}能保持正定，就用 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-4，c2 = 0.9。这个组合在大多数问题上表现稳定。
第四个坑：内存布局。如果你用 Python 的 numpy 实现，注意B_k的存储方式。B_k是对称矩阵，理论上只需要存一半，但 numpy 没有原生的对称矩阵类型。我试过用scipy.linalg.blas的对称矩阵乘法来加速，效果不错，但代码复杂度会上升。对于 n < 1000 的问题，直接用完整矩阵就行，没必要优化。
第五个坑：收敛判据。很多人只用梯度范数||g_k|| < ε作为收敛条件，但这在病态问题里可能过早停止。我一般同时检查三个条件：梯度范数、步长||s_k||、以及函数值的变化|f_{k+1} - f_k|。三个条件都满足才认为收敛。阈值的选择取决于问题的尺度，我通常用1e-6作为梯度范数的阈值，1e-10作为函数值变化的阈值。
下面是一个简化的 BFGS 实现框架，用 Python 写，展示了核心逻辑：
import numpy as np def bfgs(f, grad_f, x0, max_iter=1000, tol=1e-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-BFGS（Limited-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-BFGS，m 取 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，理解它们的设计动机和适用场景。这样学下来，不仅知道怎么用，还知道为什么这么用。
