最速下降法、牛顿法与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 − bHessian 矩阵就是 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, seed42): 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(sizen) return Q, b这里用 QR 分解得到随机正交矩阵 U然后把特征值按需求排布最后通过 Q U·diag(λ)·Uᵀ 还原出对称正定矩阵。特征值可以线性分布也可以按指数分布线性分布比较温和指数分布会让高维特征值快速衰减问题更病态。b 向量用标准正态分布随机生成。这样构造出的二次函数最优解 x* np.linalg.solve(Q, b)真实解我们心里有数对比才有底气。2. 最速下降法别看名字唬人它其实很“短视”2.1 负梯度方向与精确线搜索步长推导最速下降法又叫梯度下降法核心思想很简单在点 x_k 处函数值下降最快的局部方向是负梯度方向 −g_k于是沿这个方向走一步x_{k1} 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, tol1e-6, max_iter10000): 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-6cond 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_starnp.linalg.solve 内部会做 LU 分解效率比显式求逆高得多数值上也更稳定。这个观点适用于所有需要计算 H⁻¹g 的场景不只是二次函数。当然牛顿法也有它的致命短处一是需要 Hessian 矩阵可逆且正定否则牛顿步可能直接指向鞍点或极大值二是 Hessian 的存储和分解成本极高。n 1000 的 Hessian 矩阵就有一百万个元素内存大约 8 MB看着不大但 LU 分解的 O(n³) 计算量在 CPU 上会明显变慢n 10000 时光存储就接近 800 MB几乎不可接受。实际工程中真遇到大尺度问题又舍不得放弃二阶信息通常有两个补救方向一是给 Hessian 加正则化项 H λI保证正定二是改用拟牛顿法用一阶信息去逼近二阶曲率这正是下一节的内容。4. 拟牛顿法用梯度差拼出二阶信息BFGS4.1 BFGS 更新公式的来龙去脉拟牛顿法的出发点是解析 Hessian 常常算不出来或者算起来太贵但我们可以在迭代过程中收集函数值和梯度信息用这些信息去逼近 Hessian 的行为。怎么逼近中值定理告诉我们关于位移 s_k x_{k1} − x_k 和梯度变化 y_k ∇f(x_{k1}) − ∇f(x_k)理想的 Hessian或它的逆应该满足割线方程H_{k1}·s_k ≈ y_k如果 H_{k1} Q 是常数矩阵等式就精确成立。求一个对称正定矩阵并满足这个割线方程并不是只有唯一解。BFGS 的思想是在满足割线方程的方阵中选一个离当前近似 H_k 最近、且保持正定的更新利用 Frobenius 范数意义上的最小化最终推出一个漂亮的迭代公式。工程上更常用的是直接更新 Hessian 逆矩阵的近似记作 B_k更新公式为B_{k1} (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, tol1e-6, max_iter2000): 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( funlambda x: 0.5 * x Q x - b x, x0x0, jaclambda x: Q x - b, methodBFGS, 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 实验结果表与收敛曲线解读下面是典型运行结果不同随机种子数值会有波动但趋势一致测试条件最速下降法牛顿法BFGSn50, cond519 次迭代1 次求解6 次迭代n50, cond1000约 3200 次迭代1 次求解约 28 次迭代n200, cond1000收敛极慢超 2 万次1 次求解约 40 次迭代这张表的信息量很大。最速下降法在最简单的问题上还挺好看条件数一变大就原形毕露。它的收敛率大约正比于 ((cond−1)/(cond1))²条件数 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_{k1} − 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能少走很多弯路。我的建议是把这些代码整理成自己的模板随手改改目标函数就能用来摸底各种优化器收益远大于死记硬背。