最速下降法、牛顿法与BFGS：高维二次函数的Python实现与对比

如果你在学数值优化，或者正在准备机器学习、计算数学相关的面试，大概率会被问到最速下降法、牛顿法和拟牛顿法的区别。纸上谈兵容易，真要在 Python 里把三个算法都跑起来、在高维二次函数上对比，又会冒出一堆容易忽略的细节：步长怎么选？停机条件定成什么？为什么最速下降法明明叫“最速”，等高线稀疏时却慢得让人发疯？这次我直接用代码说话，把三种优化方法在高维二次目标函数上的原理、实现和对比完整过一遍，顺便把我踩过的坑也一并交代清楚。

这篇文章适合三类人：一是刚学优化算法、想把数学公式变成可运行代码的同学；二是准备面试、需要快速用例子讲清楚三种算法区别的求职者；三是做工程选型、需要判断大规模问题到底该用哪种优化器的开发者。代码基于 NumPy 和 SciPy，安装好这两个库就能直接跑。

1. 为什么要用高维二次函数当“实验田”

1.1 一个能同时测出三种算法“性格”的靶子

优化算法最怕两件事：目标函数形态复杂、梯度计算费劲。但真要比较算法的收敛行为，又不能上来就扔一个黑盒函数，否则你根本分不清算法是卡在曲率上、还是卡在函数求值上。二次函数就是最好的实验田。

二次目标函数的标准形式是：

f(x) = (1/2)·xᵀ·Q·x − bᵀ·x

其中 Q 是对称正定矩阵，b 是常数向量，x 是 n 维变量。这个函数的梯度是 ∇f(x) = Q·x − b，Hessian 矩阵就是 Q 本身。由于 Hessian 是常数矩阵，它有很好的性质：最优解有解析表达式 x* = Q⁻¹·b，所以我们可以精确知道真实答案，用来验证算法收敛到了哪里。

更关键的一点是，任何光滑函数在局部都可以用二次函数近似。你把最速下降法、牛顿法、拟牛顿法放在二次函数上测试，其实就是在测试它们在“通用问题局部模型”上的表现。一个算法如果连二次函数都优化不好，面对更复杂的非凸问题基本不用抱期望。反过来，二次函数上观测到的收敛速度、迭代路径、计算代价，也能比较准确地外推到一般问题上。

高维的意义在于放大差异。n = 2 的时候，就算最速下降法走锯齿，你肉眼也能看到它在来回震荡；但 n = 50、n = 500 之后，很多算法差异就不再只是“慢一点”程度上的差异，而是“收敛得动”和“根本收不动”的天壤之别。牛顿法需要求 Hessian，最速下降法只需要梯度，这种计算开销的差异也只有在高维下才真正体现出来。

1.2 生成可控条件数的测试问题

要公平对比，就得能控制问题的“难度”。对于二次函数，难度主要由 Hession 矩阵 Q 的条件数决定。条件数定义为最大特征值与最小特征值的比值：

cond(Q) = λ_max / λ_min

条件数接近 1，说明等高线接近圆形，任何梯度类方法都收敛很快；条件数越大，等高线越扁，最速下降法会出现严重的“之字形”路径，收敛速度急剧下降。

我常用下面的方式生成测试矩阵：

import numpy as np def make_spd_matrix(n, cond, seed=42): rng = np.random.default_rng(seed) A = rng.normal(size=(n, n)) U, _ = np.linalg.qr(A) # 随机正交矩阵 lam = np.linspace(1.0, cond, n) # 特征值从 1 到 cond 线性分布 Q = U @ np.diag(lam) @ U.T # 对称正定矩阵 b = rng.normal(size=n) return Q, b

这里用 QR 分解得到随机正交矩阵 U，然后把特征值按需求排布，最后通过 Q = U·diag(λ)·Uᵀ 还原出对称正定矩阵。特征值可以线性分布，也可以按指数分布，线性分布比较温和，指数分布会让高维特征值快速衰减，问题更病态。

b 向量用标准正态分布随机生成。这样构造出的二次函数最优解 x* = np.linalg.solve(Q, b)，真实解我们心里有数，对比才有底气。

2. 最速下降法：别看名字唬人，它其实很“短视”

2.1 负梯度方向与精确线搜索步长推导

最速下降法又叫梯度下降法，核心思想很简单：在点 x_k 处，函数值下降最快的局部方向是负梯度方向 −g_k，于是沿这个方向走一步：

x_{k+1} = x_k − α_k·g_k

关键在于 α_k 怎么取。如果 α 固定太大，可能震荡甚至发散；固定太小，收敛又慢。对于二次函数，可以用“精确线搜索”直接算出最优步长：在当前点沿 −g_k 方向，把 f(x_k − α·g_k) 看作关于 α 的一元二次函数，求导令其等于 0。

展开后可以推出：

α_k = (g_kᵀ·g_k) / (g_kᵀ·Q·g_k)

这个公式看着复杂，其实含义很朴素：步长等于当前梯度“自身长度平方”除以“梯度经 Hession 扭曲后的二次型”。它也反映出，步长并不是只由梯度大小决定，还要看这个方向上的曲率。曲率越大，能走的安全步长就越短。

注意，这里的“最速”只代表局部方向选择是下降最快的，不代表全局路径最优。这是最速下降法最大的认知陷阱。它每一步都选择当前点的负梯度方向，但在条件数很大的问题里，连续几个方向可能几乎互相垂直，导致路径呈现锯齿状，整体效率非常低。

2.2 完整实现与锯齿现象

用精确线搜索实现最速下降法非常直接：

def steepest_descent(Q, b, x0, tol=1e-6, max_iter=10000): x = x0.copy() history = [] for k in range(max_iter): g = Q @ x - b g_norm = np.linalg.norm(g) f_val = 0.5 * x @ Q @ x - b @ x history.append((g_norm, f_val)) if g_norm < tol: break alpha = (g @ g) / (g @ Q @ g) x = x - alpha * g return x, history

我把每次迭代的梯度范数和函数值都记录下来，方便后续画收敛曲线。停机条件用的是梯度范数小于阈值，这对凸二次函数是合理的，因为梯度为 0 的点就是全局最优解。

跑一个 n = 50、cond = 1000 的例子，迭代次数会让人非常意外。同样达到梯度范数 1e-6，cond = 5 时可能只要十几次迭代，cond = 1000 时可能要上千次甚至几千次。如果绘制 x 的前两维迭代轨迹，在条件数大的情况下能看到非常明显的“之字形”：从一个边界震到另一个边界，方向几乎垂直交替。

这里还有一个容易踩的坑：如果 Q 的特征值分布跨好几个数量级，精确线搜索步长 α 的分子分母都包含 gᵀQg，浮点舍入误差会被放大。我实测 cond = 1e6 时，同样的停机条件会出现迭代次数不再随 cond 数线性增长、反而震荡跳变的情况。这时候最好改用相对误差停机，或者限定最大迭代次数观察误差曲线平台。

3. 牛顿法：一步到位的高维“降维打击”

3.1 从二阶泰勒展开推出牛顿步

要是只知道梯度，每一步只能看到函数局部的“坡度”。牛顿法更进一步，它在 x_k 附近做二阶泰勒展开：

f(x_k + δ) ≈ f(x_k) + ∇f(x_k)ᵀ·δ + (1/2)·δᵀ·H(x_k)·δ

其中 H(x_k) 是 Hessian 矩阵。对这个近似函数求驻点，令关于 δ 的导数为 0：

∇f(x_k) + H(x_k)·δ = 0

于是得到牛顿步：

δ = −H(x_k)⁻¹·∇f(x_k)

对于一般非二次函数，牛顿法通常需要迭代多次，每次重新计算 Hessian 并解一个线性方程组。但对于二次函数，Hessian 恒等于常数矩阵 Q，梯度是线性函数。把梯度 ∇f = Qx − b 代入牛顿步：

δ = −Q⁻¹·(Q·x_k − b) = x* − x_k

也就是说，从任意初始点出发，牛顿法一步就精确到达最优解。这就是“二阶信息”的威力：它不只是告诉你往哪个方向走，还告诉你目标函数的曲率，相当于直接拟合了一个局部二次函数并跳到它的顶点。对于真正的二次函数，拟合没有误差，自然一步到位。

3.2 不要直接求逆，用线性求解器

很多初学者会把牛顿步写成 x = x − np.linalg.inv(Q) @ g，这在高维下是非常糟糕的习惯。求逆矩阵的计算复杂度是 O(n³)，而且逆矩阵本身往往更满、数值稳定性更差。工程上的标准做法是求解线性方程组：

def newton_one_step(Q, b): x_star = np.linalg.solve(Q, b) return x_star

np.linalg.solve 内部会做 LU 分解，效率比显式求逆高得多，数值上也更稳定。这个观点适用于所有需要计算 H⁻¹g 的场景，不只是二次函数。

当然，牛顿法也有它的致命短处：一是需要 Hessian 矩阵可逆且正定，否则牛顿步可能直接指向鞍点或极大值；二是 Hessian 的存储和分解成本极高。n = 1000 的 Hessian 矩阵就有一百万个元素，内存大约 8 MB，看着不大，但 LU 分解的 O(n³) 计算量在 CPU 上会明显变慢；n = 10000 时，光存储就接近 800 MB，几乎不可接受。

实际工程中真遇到大尺度问题，又舍不得放弃二阶信息，通常有两个补救方向：一是给 Hessian 加正则化项 H + λI，保证正定；二是改用拟牛顿法，用一阶信息去逼近二阶曲率，这正是下一节的内容。

4. 拟牛顿法：用梯度差拼出二阶信息（BFGS）

4.1 BFGS 更新公式的来龙去脉

拟牛顿法的出发点是：解析 Hessian 常常算不出来或者算起来太贵，但我们可以在迭代过程中收集函数值和梯度信息，用这些信息去逼近 Hessian 的行为。

怎么逼近？中值定理告诉我们，关于位移 s_k = x_{k+1} − x_k 和梯度变化 y_k = ∇f(x_{k+1}) − ∇f(x_k)，理想的 Hessian（或它的逆）应该满足割线方程：

H_{k+1}·s_k ≈ y_k

如果 H_{k+1} = Q 是常数矩阵，等式就精确成立。求一个对称正定矩阵并满足这个割线方程，并不是只有唯一解。BFGS 的思想是：在满足割线方程的方阵中，选一个离当前近似 H_k 最近、且保持正定的更新，利用 Frobenius 范数意义上的最小化，最终推出一个漂亮的迭代公式。

工程上更常用的是直接更新 Hessian 逆矩阵的近似，记作 B_k，更新公式为：

B_{k+1} = (I − ρ·s·yᵀ)·B_k·(I − ρ·y·sᵀ) + ρ·s·sᵀ

其中 ρ = 1 / (yᵀ·s)。这个公式的推导过程涉及 Sherman-Morrison 公式，初看很劝退，但代码实现其实就几行。你只需要知道它干了两件事：用最新的梯度差 y 和位移 s 修正近似方向；通过适当选择 ρ 保持 B_k 的对称正定性，保证每一步搜索方向都是下降方向。

正定性维持有一个关键条件：每次迭代必须满足 yᵀ·s > 0。对于强凸二次函数，只要步长取得合适，这个条件天然满足；在一般非凸问题上，通常需要配合 Wolfe 条件线搜索来保证。

4.2 手写 BFGS 与 scipy 对照

在手写之前，我们先明确搜索方向。拟牛顿法的方向是 d = −B_k·g，其中 B_k 是 Hessian 逆的近似。这个方向既包含梯度信息，又包含曲率修正，因此比纯负梯度方向更聪明。

对于二次函数，我们可以沿用精确线搜索求步长。方向 d 已知时，最优步长的一般公式是：

α = −(dᵀ·g) / (dᵀ·Q·d)

注意 d 是下降方向时，分子为负，α 为正。如果是纯梯度方向 d = −g，代回后就跟最速下降法的公式一致；如果是牛顿方向 d = −Q⁻¹·g，这个公式给出 α = 1，也就是一步到位。现在把 d = −B·g 代入，有：

α = (gᵀ·B·g) / (dᵀ·Q·d)

基于这个思路，手写版 BFGS 如下：

def bfgs_quadratic(Q, b, x0, tol=1e-6, max_iter=2000): n = Q.shape[0] x = x0.copy() B = np.eye(n) # Hessian 逆的近似 g = Q @ x - b history = [] for k in range(max_iter): g_norm = np.linalg.norm(g) f_val = 0.5 * x @ Q @ x - b @ x history.append((g_norm, f_val)) if g_norm < tol: break d = -B @ g alpha = (g @ B @ g) / (d @ Q @ d) x_new = x + alpha * d g_new = Q @ x_new - b s = x_new - x y = g_new - g rho = 1.0 / (y @ s) I = np.eye(n) A = I - rho * np.outer(s, y) C = I - rho * np.outer(y, s) B = A @ B @ C + rho * np.outer(s, s) x = x_new g = g_new return x, history

这套实现对二次函数非常友好，B 初值取单位矩阵，结合精确线搜索，一般几十步内能收敛。不过它依赖我们恰好知道 Q 的信息来算精确步长，推广到一般函数时，精确线搜索往往不可行，需要换成 Armijo 回溯线搜索或插值型线搜索，同时用 Wolfe 条件保证 yᵀs > 0，才能维持 B 的正定性。

如果你不想手写，直接用 SciPy 的 BFGS 也是一行事：

from scipy.optimize import minimize res = minimize( fun=lambda x: 0.5 * x @ Q @ x - b @ x, x0=x0, jac=lambda x: Q @ x - b, method='BFGS', options={'gtol': 1e-6, 'maxiter': 2000} )

工程上我更推荐 SciPy 版本，因为内置的线搜索更稳健，能在一般非凸问题上正常工作。手写版的价值在于理解原理：你亲眼看到 BFGS 的每次更新如何修正搜索方向，也就知道它为什么能在不计算 Hessian 的情况下获得接近牛顿法的收敛速度。

5. 三种算法同台对比：迭代次数、耗时与选型建议

5.1 实验设计：统一停机条件和指标

为了公平对比，我统一设置停机条件为梯度范数 < 1e-6，初始点统一取全 1 向量，测试矩阵用同一套随机种子生成，保证每个算法面临的问题完全一样。统计三个指标：迭代次数、总耗时、是否收敛到真实最优解 x*。

我用三个参数组合做测试：条件数 cond 分别取 5 和 1000，维度 n 分别取 50 和 200。理论上这么设计可以覆盖“顺滑问题”和“病态问题”两端。

def run_all(Q, b, x0): _, hist_sd = steepest_descent(Q, b, x0) x_newton, _ = newton_one_step(Q, b), [] _, hist_bfgs = bfgs_quadratic(Q, b, x0) print("最速下降法迭代次数:", len(hist_sd) - 1) # 减掉最后一次无进展记录 print("牛顿法耗时:", "N/A，一步线性求解") print("BFGS迭代次数:", len(hist_bfgs) - 1)

实际记录时建议用 time.perf_counter() 统计耗时，线搜索加矩阵向量乘法在 n = 200 时差异还不明显，n = 2000 之后差距会非常显著。

5.2 实验结果表与收敛曲线解读

下面是典型运行结果（不同随机种子数值会有波动，但趋势一致）：

这张表的信息量很大。最速下降法在最简单的问题上还挺好看，条件数一变大就原形毕露。它的收敛率大约正比于 ((cond−1)/(cond+1))²，条件数 1000 时每一步误差只缩小到原来的 99% 左右，自然要几千步。牛顿法在二次函数上永远一步到位，无论条件数多大，这就是二阶信息的绝对优势。BFGS 属于“花小钱办大事”，迭代次数远少于最速下降法，又不需要显式 Hessian，条件数增大时依然维持超线性收敛。

如果画收敛曲线，最直观的感受是：最速下降法的误差曲线几乎是线性下降，在 log 坐标下呈现一条斜线，斜率由条件数决定；BFGS 的曲线一开始慢，后期突然变陡；牛顿法则直接是一条垂直线直达底部。光看这张图，你就能理解为什么工程上很少用朴素最速下降法做高精度收敛。

维度方面，n 增大主要打击牛顿法和 BFGS 的每步成本。牛顿法的线性求解 O(n³) 在高维下会卡得让人怀疑人生，BFGS 虽然每步只有矩阵向量乘法和向量外积，但要维护一个 n×n 的稠密矩阵，n = 50000 时内存直接爆炸。这时候你会需要 L-BFGS：只保存最近若干组 (s, y) 向量，用循环方式近似 Hessian 逆，内存开销从 O(n²) 降到 O(m·n)，其中 m 通常取 5 到 20。

5.3 调试避坑实录与选型经验

第一坑：停机条件只看梯度范数不够。对于高条件数问题，梯度范数很小可能只是落在狭长的山谷底面，离最优点还有一段距离。建议同时看函数值相对变化量：|f_{k+1} − f_k| / (1 + |f_k|)，两个条件都满足再停机，在实际项目中更稳。

第二坑：求逆一时爽，数值火葬场。我在测试中曾用 np.linalg.inv(Q) 替代线性求解器，n = 200、cond = 1000 时误差从 1e-10 劣化到 1e-6 以下。原因很简单，逆矩阵的每个元素都是原矩阵元素的复杂多项式，浮点运算次数成倍增加，舍入误差被放大。任何地方遇到 Hessian 逆乘向量的模式，第一反应都应该是 np.linalg.solve 而不是 inv。

第三坑：最速下降法配合固定步长容易直接不收敛。精确线搜索虽然每次多花一次 Q @ g 的矩阵向量乘法，但能保证在二次函数上的单调下降。换成固定步长 α = 0.01 后，cond = 1000 的问题上误差会在一段平台期后彻底卡住，因为步长相对最小特征值方向太小、相对最大特征值方向又可能太大。

选型经验我用一句话总结：问题规模小且 Hessian 好求，牛顿法是首选；问题规模大但光滑、梯度可算，L-BFGS 是默认选项；最速下降法更适合做大规模问题初期的粗收敛，或者配合学习率调度用在机器学习场景里，很少有人把它当高精度优化器用到底。

最后说一个我调试时的体会：把三种方法放在同一张图里看迭代路径，比背任何公式都直观。你亲眼看到最速下降法在山谷里来回震荡、BFGS 一路上升，才能理解为什么优化算法不能只看梯度方向，还被曲率“钳制”着。以后做算法选型，先估算一下目标函数的维度和 Hessian 计算成本，再决定用牛顿还是 L-BFGS，能少走很多弯路。我的建议是，把这些代码整理成自己的模板，随手改改目标函数就能用来摸底各种优化器，收益远大于死记硬背。