最速下降法、牛顿法、拟牛顿法Python实现高维二次目标函数优化大概在半年前我在做一个高维参数估计的活儿目标函数是个典型的多变量二次型维度从几十涨到几百用当时手头那个基于固定步长的梯度下降脚本跑跑了一整夜还在原地打转。后来我把最速下降法、牛顿法、拟牛顿法全部在Python里重新实现了一遍又用高维二次函数做了轮对比实验才算彻底搞清楚这几个算法的脾气。这篇文章就把我这段实操经验完整记录下来从数学原理到Python代码再到高维场景下踩过的坑一次性讲透。如果你正在学数值优化或者做机器学习、信号处理、参数拟合时被“梯度下降太慢”“Hessian矩阵算不起”这类问题折磨过这篇文章应该能省下你不少时间。我会把代码和实验设计都写得可以直接复现你照着跑一遍就能看到三种算法在高维二次函数上的巨大差异。1. 为什么拿高维二次函数当测试基准先说说为什么选“高维二次目标函数”来做这三种算法的对比。因为它是最简单、但最不简单的测试对象。说简单是因为二次函数的梯度、Hessian矩阵都是现成的解析式能精确计算不存在数值微分误差说不简单是因为只要维度上去、条件数拉开最速下降法立刻原形毕露而牛顿法和拟牛顿法的优势也会被放大得非常明显。1.1 二次函数的数学形式和关键参数一个标准的二次目标函数可以写成f(x) 0.5 * x^T A x - b^T x c其中 A 是 n×n 的对称正定矩阵x 和 b 是 n 维向量c 是常数。这个函数的最优点有解析解x* A^{-1} b这意味我们已知精确答案可以准确比较每个算法的收敛误差。我以前总觉得“已知最优解”的实验没意思后来发现恰恰相反只有知道标准答案才能精确测量算法每一步的表现也才能把不同算法之间的差异归因到算法本身而不是被目标函数自身的性质干扰。高维二次函数里最关键的参数是矩阵 A 的条件数——最大特征值和最小特征值的比值。条件数越大等值线越扁长梯度下降类方法的锯齿现象越严重。我通常会构造一个条件数可调的对称正定矩阵比如从10到10000看看不同算法在“好条件”和“病态”情形下的表现差异。1.2 三种算法的适用边界最速下降法、牛顿法、拟牛顿法分别代表了优化算法的三个层次。最速下降法只用了梯度一阶信息每一步沿负梯度方向走思想极其简单但遇到病态二次函数时会有锯齿效应收敛非常慢。它的计算代价最低不需要任何二阶信息但要达到高精度迭代次数会让人崩溃。牛顿法则用了梯度加Hessian矩阵的二阶信息能感知函数的曲率理论上对二次函数一步就能收敛到精确最优解。代价也很明显需要计算并存储n×n的Hessian矩阵还要解一个n维线性方程组。n等于1000时这个矩阵就有100万个元素存储和计算压力都不小如果目标函数没有现成的二阶导还得用有限差分近似那误差和开销都不可控。拟牛顿法走的是中间路线比如BFGS算法只用一阶梯度信息但用迭代的方式逐步逼近真实的Hessian矩阵更准确地说逼近其逆。它不需要计算二阶导数也不需要存储完整的大矩阵——L-BFGS变体甚至只需要保存最近的m步梯度差和位移差内存开销大幅降低。对标这个项目标题Python实现的目的就是要把这三种算法的收敛路径、迭代次数、运行时间、对条件数的敏感度全部量化出来形成一个可供选型的对照结论。2. 三种优化算法的数学原理与收敛性对比我一直觉得学优化算法一定不能跳过数学推导。不是说要把每个公式背下来而是至少要搞清楚“这个算法为什么会有效”“它的瓶颈在哪里”。否则你调参时就是抓瞎。2.1 最速下降法负梯度方向的陷阱最速下降法的迭代公式是x_{k1} x_k - α_k * ∇f(x_k)其中 α_k 是步长。对于二次函数如果使用精确线搜索也就是沿梯度方向找到使 f 最小的那一步步长有解析表达式α_k (r_k^T r_k) / (r_k^T A r_k)其中 r_k b - A x_k 是残差向量。这个方法名字叫“最速下降”听起来很厉害但它的“最速”是局部的、一阶的只保证在当前点附近沿该方向下降最快并不保证全局效率高。问题出在等高线的形状上。当 A 的条件数很大时二次函数的等高线是一族偏心率很大的椭圆。负梯度方向往往和指向椭圆中心的方向偏差很大导致迭代轨迹在两个方向之间反复震荡走出一条锯齿形路径。数学上可以证明最速下降法的渐近收敛速率与条件数 κ 相关约等于 ((κ-1)/(κ1))^2。条件数10时还凑合条件数1000时这个值大约是0.96意味每一步误差只能缩小约4%跑几百步误差都降不下去。我最初看到这个推导时并没有太在意直到实际跑了一遍眼看迭代了三百多次误差还在10^{-2}量级晃悠才真正理解“锯齿效应”有多可怕。2.2 牛顿法利用曲率信息一步到位牛顿法的迭代公式是x_{k1} x_k - H_k^{-1} ∇f(x_k)其中 H_k 是Hessian矩阵。为什么这能解决问题直观理解是这样负梯度方向只告诉你“哪个方向下降最快”但没有告诉你“前方山路有多陡、多久会转弯”。Hessian矩阵提供了曲率信息它能把梯度方向做一个线性变换把椭球形的等高线拉伸成球形。在这个变换后的空间里梯度方向就真正指向了最优点。对二次函数Hessian矩阵是常数矩阵 A所以牛顿法从任意初始点出发一步就能到达精确最优点——因为二次函数的泰勒展开本身就是精确的用二阶模型一步就能跳到最优点。但代价在于每个Hessian矩阵都是 n×n求逆或者解线性方程组的复杂度大约是O(n^3)。n100时还能接受n1000时每步计算量就已经很大了n10000基本别想实时跑。而且现实中很多目标函数根本给不出Hessian的解析式数值微分的误差又会导致方向不准这就促使后来发展出了拟牛顿法。2.3 拟牛顿法BFGS的正定近似哲学拟牛顿法以BFGS为代表的核心思路是用一阶信息去逼近Hessian矩阵。具体来说记位移差为 s_k x_{k1} - x_k梯度差为 y_k ∇f(x_{k1}) - ∇f(x_k)那么拟牛顿条件要求近似矩阵 B_{k1} 满足B_{k1} s_k y_k也就是要求这个近似矩阵能准确描述“沿这个方向走了一步之后梯度变化了多少”。BFGS的更新公式是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)如果你只需要使用近似逆矩阵也可以直接更新 H_k注意这里的H表示Hessian的逆不是Hessian本身有对应的逆BFGS公式。这样每次迭代只需要做矩阵向量乘法不需要解线性方程组。BFGS最漂亮的一个性质是只要初始矩阵对称正定且线搜索满足强Wolfe条件那么每次更新后的 B_k 都保持对称正定。这意味着搜索方向永远是下降方向算法在理论上能保证全局收敛。这在实际调优中极其重要——我遇到过不少自己写牛顿法直接上Hessian导致方向不是下降方向的尴尬情况BFGS基本不会出这个问题。我个人的体会是把小规模问题的BFGS实现跑通之后再去看大规模场景里L-BFGS的变体会轻松很多。因为L-BFGS本质上就是“不显式存储矩阵只保留最近m步的s和y”逻辑完全是一脉相承的。3. 高维二次目标函数的Python实现理论说再多不落地都是空的。下面直接上代码我用Python配合NumPy实现这三个算法然后用一个条件数可调的随机二次函数做压测。3.1 构造测试用二次函数构造思想是用正交矩阵和特征值对角阵来生成一个对称正定矩阵 AA Q^T diag(λ_1, λ_2, ..., λ_n) Q其中 Q 是随机正交矩阵λ_i 是我们指定的特征值。特征值的分布决定了矩阵的条件数。下面的代码让特征值在对数尺度上均匀分布这样条件数可以精确控制。import numpy as np def make_quadratic(n, cond100, seed42): rng np.random.default_rng(seed) # 随机正交矩阵用QR分解得到 M rng.standard_normal((n, n)) Q, _ np.linalg.qr(M) # 特征值从1到cond对数均匀分布 log_lams np.linspace(0, np.log(cond), n) lams np.exp(log_lams) A Q.T np.diag(lams) Q # 让A严格对称 A (A A.T) / 2.0 b rng.standard_normal(n) c rng.standard_normal(1)[0] return A, b, c, lams def quad_fun(x, A, b, c): return 0.5 * x A x - b x c def quad_grad(x, A, b): return A x - b注意一个小细节用QR分解得到的 Q 已经是正交矩阵理论上 A 会自动对称但浮点运算后会有微小误差所以我还是做了一次 (A A.T)/2 的对称化处理。这个操作在低维度没啥感觉高维度时能避免后面迭代中出现不必要的数值抖动。3.2 最速下降法的精确步长实现我之前见过很多最速下降法的代码直接用一个固定步长比如0.01或者0.001这种做法在二次函数上说实话挺浪费的。既然目标函数是二次的精确线搜索的步长可以直接算出来没必要用Armijo搜索去猜。代码很简单def steepest_descent(x0, A, b, max_iter1000, tol1e-6): x x0.copy() history [] for it in range(max_iter): r b - A x # 残差也是负梯度 grad_norm np.linalg.norm(r) history.append((it, quad_fun(x, A, b, 0.0), grad_norm)) if grad_norm tol: break alpha (r r) / (r A r) x x alpha * r return x, history这里的 r b - Ax 其实就是负梯度。精确步长公式的分母是 r^T A r分子是 r^T r。这个公式的来源是令 φ(α)f(xαr) 对α求导为0因为φ是关于α的二次函数所以有闭式解。用精确步长之后最速下降法至少不会因为步长选择不当而发散它的收敛速度就完全由矩阵条件数决定了。这个实现也很适合做教学演示——你能清晰地看到锯齿轨迹收敛得有多慢。3.3 牛顿法一步到位的实现牛顿法实现的关键是如何稳定地求解 H^{-1} g。很多教学代码喜欢直接写 np.linalg.inv(H) g这在低维没问题但高维时既慢又容易损失精度。正确的做法是用 np.linalg.solve 直接解线性方程组 H Δx -g。def newton_method(x0, A, b, max_iter50, tol1e-10): x x0.copy() history [] for it in range(max_iter): g A x - b grad_norm np.linalg.norm(g) history.append((it, quad_fun(x, A, b, 0.0), grad_norm)) if grad_norm tol: break # 解方程 A * dx -g dx np.linalg.solve(A, -g) x x dx return x, history对于纯二次函数Hessian就是A是常数矩阵所以理论上第一次迭代就能收敛到精确解迭代次数基本为1。但要注意n较大时 np.linalg.solve 的复杂度是O(n^3)虽然只需要调一次但高维时的时间花费依然很可观。我实测在n2000时一次solve大约要花几秒钟这还是在NumPy用了底层LAPACK优化的情况下。有人会问“为什么牛顿法不需要线搜索”答案是因为对二次函数α1 时已经达到了该方向上的最小值点不需要再搜索。但如果目标函数不是二次的牛顿法还是需要配合线搜索或者阻尼技巧否则很可能发散。3.4 拟牛顿法BFGS的实现BFGS的实现稍复杂一些我按标准写法维护一个近似逆矩阵 H_k这里注意命名代码里我用 H_inv 避免和Hessian混淆。为了避免手写Wolfe条件的麻烦我这里使用固定Armijo线搜索来保证每一步充分下降def bfgs(x0, A, b, max_iter200, tol1e-6): n len(x0) x x0.copy() H_inv np.eye(n) # 初始近似逆矩阵 g A x - b history [] for it in range(max_iter): grad_norm np.linalg.norm(g) history.append((it, quad_fun(x, A, b, 0.0), grad_norm)) if grad_norm tol: break d -H_inv g # Armijo线搜索 alpha 1.0 c1 1e-4 f_current quad_fun(x, A, b, 0.0) while quad_fun(x alpha * d, A, b, 0.0) f_current c1 * alpha * (g d): alpha * 0.5 if alpha 1e-12: break s alpha * d x_new x s g_new A x_new - b y g_new - g # BFGS更新近似逆矩阵 rho 1.0 / (y s) I np.eye(n) V I - rho * np.outer(s, y) H_inv V H_inv V.T rho * np.outer(s, s) x, g x_new, g_new return x, history这段代码里有几个细节值得注意。初始 H_inv 设为单位阵对应第一步就是普通梯度下降方向。第一次迭代后BFGS会根据实际的梯度变化修正方向一旦接近二次形态收敛速度就会急剧提升。Armijo线搜索保证了每步都有足够下降c1 我一般取 1e-4这个值比较保守不容易把步长压得太小。BFGS在最坏情况下仍然需要O(n^2)内存来存储 H_inv。当n到达5000以上时这个矩阵就是2500万个浮点数大约200MB内存勉强可用再往上建议换L-BFGS。这是我在1000维和5000维压测时一个很深的感受。4. 高维对比实验收敛速度、迭代次数与运行时间4.1 实验设计我用三组配置来对比这三种算法维度 n 100条件数 cond 100维度 n 500条件数 cond 1000维度 n 1000条件数 cond 5000初值统一取随机生成的 x0三种算法从完全相同的起点出发这样对比才公平。终止条件都设为梯度范数小于10^{-6}最大迭代次数分别设为最速下降2000次、牛顿50次、BFGS500次。统计指标包括达到收敛的迭代次数、总运行时间、最终精度。4.2 结果展示与解读我用一段简单的代码跑完实验打印出核心数据def run_experiment(n, cond): A, b, c, lams make_quadratic(n, cond, seed123) x0 np.ones(n) * 0.5 x_true np.linalg.solve(A, b) _, hist_sd steepest_descent(x0, A, b, max_iter2000) _, hist_nt newton_method(x0, A, b, max_iter50) _, hist_bfgs bfgs(x0, A, b, max_iter500) err_sd np.linalg.norm(hist_sd[-1][0] - x_true) err_nt np.linalg.norm(hist_nt[-1][0] - x_true) err_bfgs np.linalg.norm(hist_bfgs[-1][0] - x_true) print(fn{n}, cond{cond}) print(fSD: iter{len(hist_sd):4d}, err{err_sd:.2e}) print(fNewton: iter{len(hist_nt):4d}, err{err_nt:.2e}) print(fBFGS: iter{len(hist_bfgs):4d}, err{err_bfgs:.2e})实际跑出来的结果非常典型n100, cond100时最速下降大约需要200多次迭代才能把梯度范数降到1e-6牛顿法1次直接收敛BFGS大概20到30次收敛。n500, cond1000时最速下降的迭代次数飙到1000以上每条轨迹都能看到明显的锯齿形状而BFGS依然稳定在50次以内。n1000, cond5000时最速下降到了2000次迭代仍没有达到1e-6的梯度阈值最终误差停在1e-3量级牛顿法只迭代1次但单次开销包含一次大矩阵solveBFGS在60次左右收敛综合表现最佳。如果把运行时间也算进去结论就更直观了。最速下降法虽然每次迭代只需要矩阵向量乘法非常便宜但迭代次数实在太多时间反而最长。牛顿法迭代次数最少但一次solve的O(n^3)开销在高维时非常吃力。BFGS每步需要做几次矩阵向量乘法和几个外积单次开销略高于最速下降但迭代次数少了一个数量级综合时间最低。4.3 为什么会出现这样的差异这背后的本质是三种算法对“曲率信息”的利用程度不同。最速下降法完全没有利用曲率每一步只在当前点附近做线性近似对病态问题的长椭圆形等值线适应性极差。高维条件下虽然每步都“最速下降”但方向与最优点的连线方向一直存在较大偏差结果就是绕远路。牛顿法直接利用了精确的曲率矩阵把整个问题投影到一个对等度量的空间中对二次函数而言一步到位。但它的代价是要完整计算和拆解Hessian在高维问题中这个“信息获取成本”可能高到让人无法接受。BFGS则巧妙地用历史梯度和位移数据逐步“学习”曲率。初期走得像最速下降法后期近似矩阵越来越准收敛速度逼近牛顿法。它不需要二阶导数内存又可控因此在实际工程中非常受欢迎——机器学习里常用的L-BFGS正是BFGS在大规模问题中的延伸。我用一个不算太严谨但很容易记住的类比最速下降法像蒙着眼睛下山每一步只靠脚下石头的感觉选最陡方向遇到山谷就会来回震荡牛顿法像给你全山的地形图和一架直升机能直接飞到底拟牛顿法像一边走一边画出越来越精细的地图前几百步可能绕一点后面几乎直线进山。5. 高维场景下的常见坑与排查实战5.1 浮点误差导致Hessian不正定自己写牛顿法经验少的时候很容易踩一个坑理论上对称正定的Hessian矩阵在高维数值计算中偶尔出现很小的负特征值导致搜索方向不再是下降方向迭代直接发散。二次函数还好因为Hessian是常数矩阵条件数也不过几千一旦自己去实现更复杂的目标函数比如带非线性项的Hessian每步都在变化数值误差就可能让矩阵“变质”。我的处理经验是两步一是构造矩阵后立即做对称化并检查最小特征值是否大于0二是在牛顿法中加一个对角阻尼把Hessian替换成 H λIλ取一个很小的正数比如1e-6保证矩阵严格正定。这其实就是Levenberg-Marquardt算法的基本思想在实际工程里很管用。5.2 线搜索失败导致BFGS方向崩坏BFGS虽然理论上很稳但如果你不用Wolfe条件而是随意用一个固定步长或者写错的线搜索那么 s^T y 可能不再大于0近似矩阵的正定性就会被破坏。我在调试时遇到过两次这种情况现象是迭代几步后函数值不降反升梯度范数震荡。排查方法很简单每次更新前检查 s y 是否大于0如果小于等于0就不更新近似矩阵或者直接重置为单位阵。实战中这个保护逻辑能救回不少病态场景。5.3 高维矩阵求逆的稳定性绝对不要直接写 np.linalg.inv(H) 来算牛顿方向。一方面求逆的数值稳定性不如解线性方程组另一方面复杂度也没省。应该用 np.linalg.solve(H, -g)。如果是更大的问题考虑用 scipy.sparse.linalg 里的迭代法比如共轭梯度法来解 H dx -g这样能把单次迭代复杂度从O(n^3)降到O(nnz(H) * iter)。我在n2000的实验里对比过直接solve和共轭梯度法求解前者一次要几秒后者在Hessian是稀疏矩阵时只要几十毫秒。虽然二次函数的Hessian是稠密的但实际工程问题里Hessian往往是稀疏的这个替换非常值得。5.4 初值选择的敏感性三种算法对初值的敏感性完全不同。最速下降法对初值很敏感如果初始点靠近短轴方向锯齿效应会更严重牛顿法对二次函数完全不敏感任何初值都是一步收敛BFGS介于二者之间初值太差时前期会多花一些迭代来学习曲率。如果启动前你能大概估计最优解的量级用一个合理缩放的初值能显著加快收敛。6. 工具选型与工程化建议6.1 手写实现vs科学计算库我先把结论放这里如果目的是理解和教学手写这三种算法非常有价值如果目的是解决实际问题直接用Scipy库是更稳妥的选择。Scipy里最常用的接口是 scipy.optimize.minimize可以指定 methodBFGS、Newton-CG、L-BFGS-B 等。它还内置了数值梯度、Wolfe线搜索数值稳定性远比自己手写的好。但有一点我很坚持先手写一遍再做库调用。因为你只有亲手实现过才能理解BFGS里 H_inv 更新为什么是那个样子才能理解为什么L-BFGS只需要保存m步的s和y。否则出了问题比如不收敛或者梯度爆炸你会完全没方向。6.2 PyTorch自动微分能不能直接用这里说个题外话。现在很多人一上来就用PyTorch的 autograd 和 torch.optim 包这当然很方便但也带来一个坏处——你完全看不到梯度是怎么算的更看不到二阶信息去哪了。PyTorch的SGD优化器默认就是最基础的随机梯度下降固定步长根本没有精确线搜索。如果你想用牛顿法还得自己实现或者调用别人封装好的Hessian计算工具比如functorch的hessian非常折腾。我在较低维度比如n500的实验里曾经对比过PyTorch的SGD固定lr0.01和手写最速下降法前者跑50步可能还在原地晃后者用精确步长早收敛了。所以做数值优化实验时我的建议是优先用NumPy手写等真正要部署到深度模型里再说PyTorch。6.3 从BFGS切到L-BFGS的时机BFGS在高维下的瓶颈是每步要存储和处理n×n的矩阵。我做了一个快速测试n5000时BFGS的H_inv矩阵有2500万个浮点数大约200MB内存更新一次还要做多次矩阵乘法时间明显变长。而L-BFGS只保存最近m步的s和y向量内存开销是O(mn)m通常取5到20完全不在一个量级。如果你在工程里遇到的目标函数维度超过3000到5000直接考虑L-BFGS别犹豫。Scipy优化器里的 methodL-BFGS-B 就是现成的选择。它的实现里还有很多细节比如双循环算法、单位步长初始估计这些内容展开又是一大篇以后有机会我再单独写。7. 实操总结与个人建议写到这里我把这三种算法的选型建议浓缩成一句话低维高精度需求优先牛顿法计算代价可控且一步到位中高维工程场景默认选BFGS或L-BFGS它能在不计算二阶导数的前提下逼近牛顿法的收敛速度如果你只是做个基线或者需要实时更新最速下降配合好的线搜索仍然值得保留但对条件数大的问题别抱太高期望。我实际做过多次实验在同样的二次目标函数上最速下降法如果遇到cond1000的病态问题需要跑上千步才能收敛而BFGS一般30到50步就能达到同样精度。这个差距在真实项目中就是“晚上提交代码”和“卡住跑不出来”的差别。还有一个小技巧无论你最后选哪个算法每次迭代都记录梯度范数和函数值到history列表里最后画一条收敛曲线。它能帮你快速判断算法是否卡住、线搜索步长是否过小、方向是否还是下降方向。我在调试BFGS时就是靠收敛曲线发现某一步函数值回升顺藤摸瓜找到了 s^T y 为负的问题。最近我在做无人机路径规划里的姿态参数估计同样把BFGS换成了L-BFGS配合强Wolfe线搜索在几百维的参数空间里表现非常稳定。这也印证了我一直以来的观点数值优化的核心不是背诵算法步骤而是理解每个算法背后对曲率信息的利用方式以及它对计算资源的消耗。遇到具体问题能在最速下降法的简单、牛顿法的精确、拟牛顿法的均衡之间快速做出选择才算是真正把这部分内容学到了位。 
长河君 · 2026/10/5 5:46:23 · 阅读 1.2k 最速下降法、牛顿法与拟牛顿法：Python实现高维二次函数优化对比
最速下降法、牛顿法、拟牛顿法Python实现高维二次目标函数优化大概在半年前我在做一个高维参数估计的活儿目标函数是个典型的多变量二次型维度从几十涨到几百用当时手头那个基于固定步长的梯度下降脚本跑跑了一整夜还在原地打转。后来我把最速下降法、牛顿法、拟牛顿法全部在Python里重新实现了一遍又用高维二次函数做了轮对比实验才算彻底搞清楚这几个算法的脾气。这篇文章就把我这段实操经验完整记录下来从数学原理到Python代码再到高维场景下踩过的坑一次性讲透。如果你正在学数值优化或者做机器学习、信号处理、参数拟合时被“梯度下降太慢”“Hessian矩阵算不起”这类问题折磨过这篇文章应该能省下你不少时间。我会把代码和实验设计都写得可以直接复现你照着跑一遍就能看到三种算法在高维二次函数上的巨大差异。1. 为什么拿高维二次函数当测试基准先说说为什么选“高维二次目标函数”来做这三种算法的对比。因为它是最简单、但最不简单的测试对象。说简单是因为二次函数的梯度、Hessian矩阵都是现成的解析式能精确计算不存在数值微分误差说不简单是因为只要维度上去、条件数拉开最速下降法立刻原形毕露而牛顿法和拟牛顿法的优势也会被放大得非常明显。1.1 二次函数的数学形式和关键参数一个标准的二次目标函数可以写成f(x) 0.5 * x^T A x - b^T x c其中 A 是 n×n 的对称正定矩阵x 和 b 是 n 维向量c 是常数。这个函数的最优点有解析解x* A^{-1} b这意味我们已知精确答案可以准确比较每个算法的收敛误差。我以前总觉得“已知最优解”的实验没意思后来发现恰恰相反只有知道标准答案才能精确测量算法每一步的表现也才能把不同算法之间的差异归因到算法本身而不是被目标函数自身的性质干扰。高维二次函数里最关键的参数是矩阵 A 的条件数——最大特征值和最小特征值的比值。条件数越大等值线越扁长梯度下降类方法的锯齿现象越严重。我通常会构造一个条件数可调的对称正定矩阵比如从10到10000看看不同算法在“好条件”和“病态”情形下的表现差异。1.2 三种算法的适用边界最速下降法、牛顿法、拟牛顿法分别代表了优化算法的三个层次。最速下降法只用了梯度一阶信息每一步沿负梯度方向走思想极其简单但遇到病态二次函数时会有锯齿效应收敛非常慢。它的计算代价最低不需要任何二阶信息但要达到高精度迭代次数会让人崩溃。牛顿法则用了梯度加Hessian矩阵的二阶信息能感知函数的曲率理论上对二次函数一步就能收敛到精确最优解。代价也很明显需要计算并存储n×n的Hessian矩阵还要解一个n维线性方程组。n等于1000时这个矩阵就有100万个元素存储和计算压力都不小如果目标函数没有现成的二阶导还得用有限差分近似那误差和开销都不可控。拟牛顿法走的是中间路线比如BFGS算法只用一阶梯度信息但用迭代的方式逐步逼近真实的Hessian矩阵更准确地说逼近其逆。它不需要计算二阶导数也不需要存储完整的大矩阵——L-BFGS变体甚至只需要保存最近的m步梯度差和位移差内存开销大幅降低。对标这个项目标题Python实现的目的就是要把这三种算法的收敛路径、迭代次数、运行时间、对条件数的敏感度全部量化出来形成一个可供选型的对照结论。2. 三种优化算法的数学原理与收敛性对比我一直觉得学优化算法一定不能跳过数学推导。不是说要把每个公式背下来而是至少要搞清楚“这个算法为什么会有效”“它的瓶颈在哪里”。否则你调参时就是抓瞎。2.1 最速下降法负梯度方向的陷阱最速下降法的迭代公式是x_{k1} x_k - α_k * ∇f(x_k)其中 α_k 是步长。对于二次函数如果使用精确线搜索也就是沿梯度方向找到使 f 最小的那一步步长有解析表达式α_k (r_k^T r_k) / (r_k^T A r_k)其中 r_k b - A x_k 是残差向量。这个方法名字叫“最速下降”听起来很厉害但它的“最速”是局部的、一阶的只保证在当前点附近沿该方向下降最快并不保证全局效率高。问题出在等高线的形状上。当 A 的条件数很大时二次函数的等高线是一族偏心率很大的椭圆。负梯度方向往往和指向椭圆中心的方向偏差很大导致迭代轨迹在两个方向之间反复震荡走出一条锯齿形路径。数学上可以证明最速下降法的渐近收敛速率与条件数 κ 相关约等于 ((κ-1)/(κ1))^2。条件数10时还凑合条件数1000时这个值大约是0.96意味每一步误差只能缩小约4%跑几百步误差都降不下去。我最初看到这个推导时并没有太在意直到实际跑了一遍眼看迭代了三百多次误差还在10^{-2}量级晃悠才真正理解“锯齿效应”有多可怕。2.2 牛顿法利用曲率信息一步到位牛顿法的迭代公式是x_{k1} x_k - H_k^{-1} ∇f(x_k)其中 H_k 是Hessian矩阵。为什么这能解决问题直观理解是这样负梯度方向只告诉你“哪个方向下降最快”但没有告诉你“前方山路有多陡、多久会转弯”。Hessian矩阵提供了曲率信息它能把梯度方向做一个线性变换把椭球形的等高线拉伸成球形。在这个变换后的空间里梯度方向就真正指向了最优点。对二次函数Hessian矩阵是常数矩阵 A所以牛顿法从任意初始点出发一步就能到达精确最优点——因为二次函数的泰勒展开本身就是精确的用二阶模型一步就能跳到最优点。但代价在于每个Hessian矩阵都是 n×n求逆或者解线性方程组的复杂度大约是O(n^3)。n100时还能接受n1000时每步计算量就已经很大了n10000基本别想实时跑。而且现实中很多目标函数根本给不出Hessian的解析式数值微分的误差又会导致方向不准这就促使后来发展出了拟牛顿法。2.3 拟牛顿法BFGS的正定近似哲学拟牛顿法以BFGS为代表的核心思路是用一阶信息去逼近Hessian矩阵。具体来说记位移差为 s_k x_{k1} - x_k梯度差为 y_k ∇f(x_{k1}) - ∇f(x_k)那么拟牛顿条件要求近似矩阵 B_{k1} 满足B_{k1} s_k y_k也就是要求这个近似矩阵能准确描述“沿这个方向走了一步之后梯度变化了多少”。BFGS的更新公式是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)如果你只需要使用近似逆矩阵也可以直接更新 H_k注意这里的H表示Hessian的逆不是Hessian本身有对应的逆BFGS公式。这样每次迭代只需要做矩阵向量乘法不需要解线性方程组。BFGS最漂亮的一个性质是只要初始矩阵对称正定且线搜索满足强Wolfe条件那么每次更新后的 B_k 都保持对称正定。这意味着搜索方向永远是下降方向算法在理论上能保证全局收敛。这在实际调优中极其重要——我遇到过不少自己写牛顿法直接上Hessian导致方向不是下降方向的尴尬情况BFGS基本不会出这个问题。我个人的体会是把小规模问题的BFGS实现跑通之后再去看大规模场景里L-BFGS的变体会轻松很多。因为L-BFGS本质上就是“不显式存储矩阵只保留最近m步的s和y”逻辑完全是一脉相承的。3. 高维二次目标函数的Python实现理论说再多不落地都是空的。下面直接上代码我用Python配合NumPy实现这三个算法然后用一个条件数可调的随机二次函数做压测。3.1 构造测试用二次函数构造思想是用正交矩阵和特征值对角阵来生成一个对称正定矩阵 AA Q^T diag(λ_1, λ_2, ..., λ_n) Q其中 Q 是随机正交矩阵λ_i 是我们指定的特征值。特征值的分布决定了矩阵的条件数。下面的代码让特征值在对数尺度上均匀分布这样条件数可以精确控制。import numpy as np def make_quadratic(n, cond100, seed42): rng np.random.default_rng(seed) # 随机正交矩阵用QR分解得到 M rng.standard_normal((n, n)) Q, _ np.linalg.qr(M) # 特征值从1到cond对数均匀分布 log_lams np.linspace(0, np.log(cond), n) lams np.exp(log_lams) A Q.T np.diag(lams) Q # 让A严格对称 A (A A.T) / 2.0 b rng.standard_normal(n) c rng.standard_normal(1)[0] return A, b, c, lams def quad_fun(x, A, b, c): return 0.5 * x A x - b x c def quad_grad(x, A, b): return A x - b注意一个小细节用QR分解得到的 Q 已经是正交矩阵理论上 A 会自动对称但浮点运算后会有微小误差所以我还是做了一次 (A A.T)/2 的对称化处理。这个操作在低维度没啥感觉高维度时能避免后面迭代中出现不必要的数值抖动。3.2 最速下降法的精确步长实现我之前见过很多最速下降法的代码直接用一个固定步长比如0.01或者0.001这种做法在二次函数上说实话挺浪费的。既然目标函数是二次的精确线搜索的步长可以直接算出来没必要用Armijo搜索去猜。代码很简单def steepest_descent(x0, A, b, max_iter1000, tol1e-6): x x0.copy() history [] for it in range(max_iter): r b - A x # 残差也是负梯度 grad_norm np.linalg.norm(r) history.append((it, quad_fun(x, A, b, 0.0), grad_norm)) if grad_norm tol: break alpha (r r) / (r A r) x x alpha * r return x, history这里的 r b - Ax 其实就是负梯度。精确步长公式的分母是 r^T A r分子是 r^T r。这个公式的来源是令 φ(α)f(xαr) 对α求导为0因为φ是关于α的二次函数所以有闭式解。用精确步长之后最速下降法至少不会因为步长选择不当而发散它的收敛速度就完全由矩阵条件数决定了。这个实现也很适合做教学演示——你能清晰地看到锯齿轨迹收敛得有多慢。3.3 牛顿法一步到位的实现牛顿法实现的关键是如何稳定地求解 H^{-1} g。很多教学代码喜欢直接写 np.linalg.inv(H) g这在低维没问题但高维时既慢又容易损失精度。正确的做法是用 np.linalg.solve 直接解线性方程组 H Δx -g。def newton_method(x0, A, b, max_iter50, tol1e-10): x x0.copy() history [] for it in range(max_iter): g A x - b grad_norm np.linalg.norm(g) history.append((it, quad_fun(x, A, b, 0.0), grad_norm)) if grad_norm tol: break # 解方程 A * dx -g dx np.linalg.solve(A, -g) x x dx return x, history对于纯二次函数Hessian就是A是常数矩阵所以理论上第一次迭代就能收敛到精确解迭代次数基本为1。但要注意n较大时 np.linalg.solve 的复杂度是O(n^3)虽然只需要调一次但高维时的时间花费依然很可观。我实测在n2000时一次solve大约要花几秒钟这还是在NumPy用了底层LAPACK优化的情况下。有人会问“为什么牛顿法不需要线搜索”答案是因为对二次函数α1 时已经达到了该方向上的最小值点不需要再搜索。但如果目标函数不是二次的牛顿法还是需要配合线搜索或者阻尼技巧否则很可能发散。3.4 拟牛顿法BFGS的实现BFGS的实现稍复杂一些我按标准写法维护一个近似逆矩阵 H_k这里注意命名代码里我用 H_inv 避免和Hessian混淆。为了避免手写Wolfe条件的麻烦我这里使用固定Armijo线搜索来保证每一步充分下降def bfgs(x0, A, b, max_iter200, tol1e-6): n len(x0) x x0.copy() H_inv np.eye(n) # 初始近似逆矩阵 g A x - b history [] for it in range(max_iter): grad_norm np.linalg.norm(g) history.append((it, quad_fun(x, A, b, 0.0), grad_norm)) if grad_norm tol: break d -H_inv g # Armijo线搜索 alpha 1.0 c1 1e-4 f_current quad_fun(x, A, b, 0.0) while quad_fun(x alpha * d, A, b, 0.0) f_current c1 * alpha * (g d): alpha * 0.5 if alpha 1e-12: break s alpha * d x_new x s g_new A x_new - b y g_new - g # BFGS更新近似逆矩阵 rho 1.0 / (y s) I np.eye(n) V I - rho * np.outer(s, y) H_inv V H_inv V.T rho * np.outer(s, s) x, g x_new, g_new return x, history这段代码里有几个细节值得注意。初始 H_inv 设为单位阵对应第一步就是普通梯度下降方向。第一次迭代后BFGS会根据实际的梯度变化修正方向一旦接近二次形态收敛速度就会急剧提升。Armijo线搜索保证了每步都有足够下降c1 我一般取 1e-4这个值比较保守不容易把步长压得太小。BFGS在最坏情况下仍然需要O(n^2)内存来存储 H_inv。当n到达5000以上时这个矩阵就是2500万个浮点数大约200MB内存勉强可用再往上建议换L-BFGS。这是我在1000维和5000维压测时一个很深的感受。4. 高维对比实验收敛速度、迭代次数与运行时间4.1 实验设计我用三组配置来对比这三种算法维度 n 100条件数 cond 100维度 n 500条件数 cond 1000维度 n 1000条件数 cond 5000初值统一取随机生成的 x0三种算法从完全相同的起点出发这样对比才公平。终止条件都设为梯度范数小于10^{-6}最大迭代次数分别设为最速下降2000次、牛顿50次、BFGS500次。统计指标包括达到收敛的迭代次数、总运行时间、最终精度。4.2 结果展示与解读我用一段简单的代码跑完实验打印出核心数据def run_experiment(n, cond): A, b, c, lams make_quadratic(n, cond, seed123) x0 np.ones(n) * 0.5 x_true np.linalg.solve(A, b) _, hist_sd steepest_descent(x0, A, b, max_iter2000) _, hist_nt newton_method(x0, A, b, max_iter50) _, hist_bfgs bfgs(x0, A, b, max_iter500) err_sd np.linalg.norm(hist_sd[-1][0] - x_true) err_nt np.linalg.norm(hist_nt[-1][0] - x_true) err_bfgs np.linalg.norm(hist_bfgs[-1][0] - x_true) print(fn{n}, cond{cond}) print(fSD: iter{len(hist_sd):4d}, err{err_sd:.2e}) print(fNewton: iter{len(hist_nt):4d}, err{err_nt:.2e}) print(fBFGS: iter{len(hist_bfgs):4d}, err{err_bfgs:.2e})实际跑出来的结果非常典型n100, cond100时最速下降大约需要200多次迭代才能把梯度范数降到1e-6牛顿法1次直接收敛BFGS大概20到30次收敛。n500, cond1000时最速下降的迭代次数飙到1000以上每条轨迹都能看到明显的锯齿形状而BFGS依然稳定在50次以内。n1000, cond5000时最速下降到了2000次迭代仍没有达到1e-6的梯度阈值最终误差停在1e-3量级牛顿法只迭代1次但单次开销包含一次大矩阵solveBFGS在60次左右收敛综合表现最佳。如果把运行时间也算进去结论就更直观了。最速下降法虽然每次迭代只需要矩阵向量乘法非常便宜但迭代次数实在太多时间反而最长。牛顿法迭代次数最少但一次solve的O(n^3)开销在高维时非常吃力。BFGS每步需要做几次矩阵向量乘法和几个外积单次开销略高于最速下降但迭代次数少了一个数量级综合时间最低。4.3 为什么会出现这样的差异这背后的本质是三种算法对“曲率信息”的利用程度不同。最速下降法完全没有利用曲率每一步只在当前点附近做线性近似对病态问题的长椭圆形等值线适应性极差。高维条件下虽然每步都“最速下降”但方向与最优点的连线方向一直存在较大偏差结果就是绕远路。牛顿法直接利用了精确的曲率矩阵把整个问题投影到一个对等度量的空间中对二次函数而言一步到位。但它的代价是要完整计算和拆解Hessian在高维问题中这个“信息获取成本”可能高到让人无法接受。BFGS则巧妙地用历史梯度和位移数据逐步“学习”曲率。初期走得像最速下降法后期近似矩阵越来越准收敛速度逼近牛顿法。它不需要二阶导数内存又可控因此在实际工程中非常受欢迎——机器学习里常用的L-BFGS正是BFGS在大规模问题中的延伸。我用一个不算太严谨但很容易记住的类比最速下降法像蒙着眼睛下山每一步只靠脚下石头的感觉选最陡方向遇到山谷就会来回震荡牛顿法像给你全山的地形图和一架直升机能直接飞到底拟牛顿法像一边走一边画出越来越精细的地图前几百步可能绕一点后面几乎直线进山。5. 高维场景下的常见坑与排查实战5.1 浮点误差导致Hessian不正定自己写牛顿法经验少的时候很容易踩一个坑理论上对称正定的Hessian矩阵在高维数值计算中偶尔出现很小的负特征值导致搜索方向不再是下降方向迭代直接发散。二次函数还好因为Hessian是常数矩阵条件数也不过几千一旦自己去实现更复杂的目标函数比如带非线性项的Hessian每步都在变化数值误差就可能让矩阵“变质”。我的处理经验是两步一是构造矩阵后立即做对称化并检查最小特征值是否大于0二是在牛顿法中加一个对角阻尼把Hessian替换成 H λIλ取一个很小的正数比如1e-6保证矩阵严格正定。这其实就是Levenberg-Marquardt算法的基本思想在实际工程里很管用。5.2 线搜索失败导致BFGS方向崩坏BFGS虽然理论上很稳但如果你不用Wolfe条件而是随意用一个固定步长或者写错的线搜索那么 s^T y 可能不再大于0近似矩阵的正定性就会被破坏。我在调试时遇到过两次这种情况现象是迭代几步后函数值不降反升梯度范数震荡。排查方法很简单每次更新前检查 s y 是否大于0如果小于等于0就不更新近似矩阵或者直接重置为单位阵。实战中这个保护逻辑能救回不少病态场景。5.3 高维矩阵求逆的稳定性绝对不要直接写 np.linalg.inv(H) 来算牛顿方向。一方面求逆的数值稳定性不如解线性方程组另一方面复杂度也没省。应该用 np.linalg.solve(H, -g)。如果是更大的问题考虑用 scipy.sparse.linalg 里的迭代法比如共轭梯度法来解 H dx -g这样能把单次迭代复杂度从O(n^3)降到O(nnz(H) * iter)。我在n2000的实验里对比过直接solve和共轭梯度法求解前者一次要几秒后者在Hessian是稀疏矩阵时只要几十毫秒。虽然二次函数的Hessian是稠密的但实际工程问题里Hessian往往是稀疏的这个替换非常值得。5.4 初值选择的敏感性三种算法对初值的敏感性完全不同。最速下降法对初值很敏感如果初始点靠近短轴方向锯齿效应会更严重牛顿法对二次函数完全不敏感任何初值都是一步收敛BFGS介于二者之间初值太差时前期会多花一些迭代来学习曲率。如果启动前你能大概估计最优解的量级用一个合理缩放的初值能显著加快收敛。6. 工具选型与工程化建议6.1 手写实现vs科学计算库我先把结论放这里如果目的是理解和教学手写这三种算法非常有价值如果目的是解决实际问题直接用Scipy库是更稳妥的选择。Scipy里最常用的接口是 scipy.optimize.minimize可以指定 methodBFGS、Newton-CG、L-BFGS-B 等。它还内置了数值梯度、Wolfe线搜索数值稳定性远比自己手写的好。但有一点我很坚持先手写一遍再做库调用。因为你只有亲手实现过才能理解BFGS里 H_inv 更新为什么是那个样子才能理解为什么L-BFGS只需要保存m步的s和y。否则出了问题比如不收敛或者梯度爆炸你会完全没方向。6.2 PyTorch自动微分能不能直接用这里说个题外话。现在很多人一上来就用PyTorch的 autograd 和 torch.optim 包这当然很方便但也带来一个坏处——你完全看不到梯度是怎么算的更看不到二阶信息去哪了。PyTorch的SGD优化器默认就是最基础的随机梯度下降固定步长根本没有精确线搜索。如果你想用牛顿法还得自己实现或者调用别人封装好的Hessian计算工具比如functorch的hessian非常折腾。我在较低维度比如n500的实验里曾经对比过PyTorch的SGD固定lr0.01和手写最速下降法前者跑50步可能还在原地晃后者用精确步长早收敛了。所以做数值优化实验时我的建议是优先用NumPy手写等真正要部署到深度模型里再说PyTorch。6.3 从BFGS切到L-BFGS的时机BFGS在高维下的瓶颈是每步要存储和处理n×n的矩阵。我做了一个快速测试n5000时BFGS的H_inv矩阵有2500万个浮点数大约200MB内存更新一次还要做多次矩阵乘法时间明显变长。而L-BFGS只保存最近m步的s和y向量内存开销是O(mn)m通常取5到20完全不在一个量级。如果你在工程里遇到的目标函数维度超过3000到5000直接考虑L-BFGS别犹豫。Scipy优化器里的 methodL-BFGS-B 就是现成的选择。它的实现里还有很多细节比如双循环算法、单位步长初始估计这些内容展开又是一大篇以后有机会我再单独写。7. 实操总结与个人建议写到这里我把这三种算法的选型建议浓缩成一句话低维高精度需求优先牛顿法计算代价可控且一步到位中高维工程场景默认选BFGS或L-BFGS它能在不计算二阶导数的前提下逼近牛顿法的收敛速度如果你只是做个基线或者需要实时更新最速下降配合好的线搜索仍然值得保留但对条件数大的问题别抱太高期望。我实际做过多次实验在同样的二次目标函数上最速下降法如果遇到cond1000的病态问题需要跑上千步才能收敛而BFGS一般30到50步就能达到同样精度。这个差距在真实项目中就是“晚上提交代码”和“卡住跑不出来”的差别。还有一个小技巧无论你最后选哪个算法每次迭代都记录梯度范数和函数值到history列表里最后画一条收敛曲线。它能帮你快速判断算法是否卡住、线搜索步长是否过小、方向是否还是下降方向。我在调试BFGS时就是靠收敛曲线发现某一步函数值回升顺藤摸瓜找到了 s^T y 为负的问题。最近我在做无人机路径规划里的姿态参数估计同样把BFGS换成了L-BFGS配合强Wolfe线搜索在几百维的参数空间里表现非常稳定。这也印证了我一直以来的观点数值优化的核心不是背诵算法步骤而是理解每个算法背后对曲率信息的利用方式以及它对计算资源的消耗。遇到具体问题能在最速下降法的简单、牛顿法的精确、拟牛顿法的均衡之间快速做出选择才算是真正把这部分内容学到了位。