1. 为什么还要手写这三种最优化算法先抛一个问题scipy.optimize.minimize一行代码就能跑完的活为什么还要自己用 Python 手写最速下降法、牛顿法、拟牛顿法我最初也这么想直到有一次我在处理一个带正则项的高维二次目标函数时发现内置优化器虽然能收敛但我根本不知道它内部走了哪条路。当我需要解释为什么某个参数组合会导致不收敛、为什么换一个初始点结果差异巨大时光会调库是不够的。于是我自己把三类方法完整实现了一遍用同一个高维二次函数做基准收获非常大。这篇文章面向的读者是那种已经会调sklearn、scipy但对底层优化器内部机制还停留在大概知道阶段的同学。我会带你从数学直觉出发用 Python 从零实现最速下降法Steepest Descent、经典牛顿法Newtons Method和以 BFGS 为代表的拟牛顿法Quasi-Newton Methods然后在高维二次目标函数上做详细对比。你会看到它们各自的收敛特性、计算代价和实际工程中易踩的坑。实话说这类算法的教科书资料已经很多但大部分停留在二维可视化演示。工程上没人只优化一个二维函数一旦进入 50 维、100 维很多低维下看不出来的问题就会浮出水面。这篇文章的重点就是把手写实现的细节、高维场景下的表现差异和调参经验讲透。我尽量用大白话解释背后的数学原理同时所有代码都可以直接复制运行建议你边看边跑。现在先来建立最基本的直觉这三类方法本质上都在回答下一步往哪个方向走、走多远这个问题但回答方式截然不同。2. 三类优化方法的本质区别走直线、用曲面、估计曲面很多资料一上来就扔公式把初学者直接劝退。我想换个方式——先说几何直觉再补数学形式。我们优化的目标函数是一个高维二次型[ f(\mathbf{x}) \frac{1}{2}\mathbf{x}^T A \mathbf{x} - \mathbf{b}^T \mathbf{x} c ]其中 (A) 是对称正定矩阵。这类函数在空间里是一个碗——等高线是椭圆高维下是超椭球。最优解就是碗底解析解是 (\mathbf{x}^* A^{-1}\mathbf{b})。既然有解析解为什么还要迭代优化器两个原因一是 (A) 的维度太高时直接求逆代价巨大二是真实问题往往不是纯二次函数迭代方法才有普适性。用二次函数做基准测试是因为它能让我们精确分析每一步的行为。2.1 最速下降法只盯着脚下最陡的方向最速下降法的逻辑简单到令人发指当前位置的负梯度方向就是函数值下降最快的方向。于是每一步都朝负梯度方向走。用公式表示就是[ \mathbf{x}_{k1} \mathbf{x}_k - \alpha_k \nabla f(\mathbf{x}_k) ]其中 (\alpha_k) 是步长梯度 (\nabla f(\mathbf{x}_k) A\mathbf{x}_k - \mathbf{b})。在二维碗状图里最速下降法表现为从初始点出发沿着和等高线垂直的方向走一个之字形。如果碗是圆形的(A) 是单位矩阵的倍数一步就能从任意点直抵碗底。但如果碗被压扁成很长的椭圆条件数大那么每一步的方向几乎和前一步垂直收敛速度会变得令人绝望地慢。教科书里有一句非常经典的话最速下降法中相邻两步的方向是正交的。这看起来是个优雅的性质实际上恰恰是它效率低下的原因——你每走一步都在矫正上一步的方向。2.2 牛顿法用二次曲面逼近局部地形牛顿法的想法更聪明我不光看当前点的梯度一阶信息还用 Hessian 矩阵二阶信息来估计局部的曲面形状然后直接跳到这个二次曲面的底部。更新公式是[ \mathbf{x}_{k1} \mathbf{x}_k - [\nabla^2 f(\mathbf{x}_k)]^{-1} \nabla f(\mathbf{x}_k) ]对于二次函数来说牛顿法有一个非常漂亮的结论一步到位。因为二次函数的 Hessian 就是常数矩阵 (A)牛顿法本质上是在做用当前点的切线信息反推碗底在哪而这个推断对二次函数是精确的。所以如果目标函数真的是二次型你从任何初始点出发牛顿法一步就收敛到精确解不考虑浮点误差。这在高维下尤其震撼你用最速下降可能要跑上千步牛顿法一步完事。代价是什么呢如果维度是 (n)Hessian 矩阵是 (n \times n) 的求逆的复杂度是 (O(n^3))。100 维还好10000 维的实践中这个代价就有点吃不消了。更麻烦的是真实问题中很多函数的 Hessian 矩阵并不总是正定的直接套用牛顿法可能走到鞍点甚至极大点。2.3 拟牛顿法不计算 Hessian用梯度差去猜拟牛顿法想要的是保留牛顿法的快速收敛特性同时避开显式计算和求逆 Hessian 的高昂代价。它的核心思想是用连续的梯度值变化去逐步构造一个对 Hessian 矩阵或其逆的近似这个近似满足所谓的割线方程[ B_{k1} (\mathbf{x}_{k1} - \mathbf{x}k) \nabla f(\mathbf{x}{k1}) - \nabla f(\mathbf{x}_k) ]最常见的实现是 BFGS 算法。它维护一个近似 Hessian 逆的矩阵 (H_k)每次迭代用梯度差和位移差做一次低秩更新然后沿着 (-H_k \nabla f) 方向搜索。更新的公式推导比较繁琐但核心直觉是你每走一步就获得一组新的梯度随位置变化的信息用这些信息修正你对曲面弯曲程度的估计。随着迭代进行(H_k) 会越来越接近真实的 Hessian 逆因此收敛速度接近牛顿法但单步代价只有 (O(n^2))。为了让你对这三类方法有更直观的对比我列一个表方法使用信息典型收敛速度单步计算量对二次函数的表现最速下降法一阶梯度线性收敛(O(n))条件数大时极慢牛顿法一阶梯度 Hessian二阶收敛(O(n^3))求逆一步精确收敛拟牛顿法BFGS一阶梯度 近似Hessian超线性收敛(O(n^2))近似一步收敛接下来我把这三类方法一一串成可运行的 Python 代码然后放到高维二次函数上去实测看它们的真实表现是不是和理论预期一致。3. Python 实现从零手写三大优化器代码部分的思路是定义好高维二次目标函数、梯度函数以及可选的真 Hessian 函数然后分别实现最速下降法、牛顿法和 BFGS 拟牛顿法最后用统一的接口跑 benchmark。这样对比公平也方便你后续换自己的目标函数做测试。先导入需要的库。这里我只用numpy不用scipy.optimize里的现成优化器保证所有实现都是透明的。import numpy as np import matplotlib.pyplot as plt import time np.random.seed(42)3.1 构造高维二次目标函数我定义目标函数为[ f(\mathbf{x}) \frac{1}{2}\mathbf{x}^T A \mathbf{x} - \mathbf{b}^T \mathbf{x} ]其中 (A) 是对称正定矩阵。为了让测试更有区分度我把它构造成一个条件数可控的矩阵。最简单的方式是用 SVD生成一个随机正交矩阵 (U)再指定一个对角阵 (\Sigma)对角线元素从 (\lambda_{\min}) 到 (\lambda_{\max}) 按对数均匀分布最后 (A U \Sigma U^T)。这样我可以精确控制条件数 (\kappa \lambda_{\max} / \lambda_{\min})。def make_quadratic(n50, cond100): 构造一个 n 维对称正定矩阵 A 和向量 b。cond 控制条件数。 U, _ np.linalg.qr(np.random.randn(n, n)) eigenvalues np.logspace(np.log10(1.0), np.log10(cond), n) A U np.diag(eigenvalues) U.T b np.random.randn(n) return A, b class QuadraticProblem: def __init__(self, n50, cond100): self.n n self.A, self.b make_quadratic(n, cond) self.x_star np.linalg.solve(self.A, self.b) self.f_star self.f(self.x_star) def f(self, x): return 0.5 * x self.A x - self.b x def grad(self, x): return self.A x - self.b def hessian(self): return self.A为什么要用这种构造方式因为我能随时解析地算出最优解 (\mathbf{x}^* A^{-1}\mathbf{b})然后精确测量每次迭代点和最优解的误差方便对比收敛曲线。同时条件数可控方便我展示病态问题下不同方法的差异。这是一个非常有用的自测环境。3.2 统一步长策略精确线搜索最速下降法和拟牛顿法都需要步长 (\alpha)。教科书中最速下降法用的是精确线搜索即在当前方向 (p) 上求解一维最小值[ \alpha_k \arg\min_{\alpha} f(\mathbf{x}_k \alpha p) ]对于二次函数这个一维最小化有解析式[ \alpha_k -\frac{p^T \nabla f(\mathbf{x}_k)}{p^T A p} ]这个公式非常实用。为什么因为如果不用精确线搜索而随便给个固定步长最速下降法很容易震荡甚至发散。而在二次函数上精确线搜索的计算代价很低相当于免费获得了最优步长。真实工程中目标函数不是二次的精确线搜索无法解析计算那时候常用 Armijo 回溯线搜索后面我会提到。我把它实现成一个函数接收方向 (p)、当前点 (x)返回最优步长def exact_line_search(prob, x, p): 二次函数上的精确线搜索。返回使得 f(x alpha * p) 最小的 alpha。 g prob.grad(x) denom p prob.A p if abs(denom) 1e-12: return 0.0 alpha -(g p) / denom return alpha注意这里我用的是解析公式前提是目标函数是二次的。后面讨论非二次扩展时会换用回溯线搜索。3.3 最速下降法的实现def steepest_descent(prob, x0, max_iter2000, tol1e-6): 最速下降法。 返回迭代轨迹点列表、每步的梯度范数、迭代次数 x x0.copy() trace [x.copy()] grad_norms [] for k in range(max_iter): g prob.grad(x) gn np.linalg.norm(g) grad_norms.append(gn) if gn tol: return trace, grad_norms, k, x p -g alpha exact_line_search(prob, x, p) x x alpha * p trace.append(x.copy()) return trace, grad_norms, max_iter, x每次迭代的核心逻辑只有三步算梯度选方向负梯度精确线搜索确定步长然后更新。实现起来极其简单但后面你会看到简单不等于高效。3.4 经典牛顿法的实现def newton_method(prob, x0, max_iter100, tol1e-10): 经典牛顿法直接用解析 Hessian一步求解。 x x0.copy() trace [x.copy()] grad_norms [] for k in range(max_iter): g prob.grad(x) gn np.linalg.norm(g) grad_norms.append(gn) if gn tol: return trace, grad_norms, k, x H prob.hessian() # 解线性方程组 H * delta -g比直接求逆更稳定 delta np.linalg.solve(H, -g) x x delta trace.append(x.copy()) return trace, grad_norms, max_iter, x一个关键细节我没有用np.linalg.inv(H) (-g)而是用np.linalg.solve(H, -g)。这看起来是个小差别但数值稳定性上差异很大——求解线性方程组比显式求逆矩阵并且再做矩阵乘法误差要小得多。尤其是当矩阵接近奇异条件数很大时显式求逆会放大误差。作为追求稳定的实践者能用solve就不要用inv。对于二次函数理论上牛顿法一步就能到最优解但我仍然写了循环一是保证函数接口统一二是展示当 Hessian 矩阵需要反复求解时牛顿法的单次迭代开销到底有多大。3.5 拟牛顿法BFGS的实现BFGS 是拟牛顿家族里最经典、实战效果最好的算法。它的完整推导在教科书里占了很大篇幅我这里只讲实现逻辑维护一个近似 Hessian 逆的矩阵 (H_k)初始为单位矩阵或某个正定矩阵每一步用梯度差和位移差去更新它。更新公式[ H_{k1} H_k \frac{(s_k^T y_k y_k^T H_k y_k)}{(s_k^T y_k)^2} s_k s_k^T - \frac{H_k y_k s_k^T s_k y_k^T H_k}{s_k^T y_k} ]其中 (s_k x_{k1} - x_k)(y_k \nabla f(x_{k1}) - \nabla f(x_k))。这个公式看着吓人但写代码只是照抄。关键的风险点是分母 (s_k^T y_k)由于浮点误差可能接近 0实际实现中需要加一个小正则项。def bfgs(prob, x0, max_iter200, tol1e-6): BFGS 拟牛顿法。 维护近似的 Hessian 逆矩阵 H迭代更新。 n prob.n x x0.copy() H np.eye(n) trace [x.copy()] grad_norms [] g prob.grad(x) grad_norms.append(np.linalg.norm(g)) for k in range(max_iter): if np.linalg.norm(g) tol: return trace, grad_norms, k, x p -H g alpha exact_line_search(prob, x, p) x_new x alpha * p g_new prob.grad(x_new) s x_new - x y g_new - g sy s y if abs(sy) 1e-12: # 防止数值退化 print(fBFGS: sy 接近 0提前停止第 {k} 次迭代) return trace, grad_norms, k, x # BFGS 更新公式 Hy H y rho 1.0 / sy H H rho * (s s * (1.0 y Hy) - (np.outer(s, Hy) np.outer(Hy, s))) x x_new g g_new trace.append(x.copy()) grad_norms.append(np.linalg.norm(g)) return trace, grad_norms, max_iter, x这里有几个工程细节值得展开说。第一初始 (H_0) 的选择。单位矩阵是最省事的选择但实践中如果能估计一下目标函数的局部曲率给一个接近真实 Hessian 逆的量级的初始化能让前几次迭代快不少。在二次函数场景里真实 Hessian 逆就是 (A^{-1})它的对角线元素和特征值分布有关。不过为了公平对比我统一用单位矩阵初始化。第二精确线搜索。BFGS 的搜索方向本身已经经过 (H_k) 的预条件处理通常给出的方向比负梯度更接近碗底方向。在二次函数上配合精确线搜索收敛会非常快。这里有一个微妙之处标准的 BFGS 在实际工程中常配合 Wolfe 条件线搜索不要求精确求一维最小值因为真实函数的精确线搜索代价太高。但在二次函数测试环境下精确线搜索开挂般的效率值得一用。第三更新公式的写法。我用的形式是单步秩二更新的紧凑表达式为了代码可读性我把它拆成了几行。实际中你也可以用更节省内存的版本比如 L-BFGSLimited-memory BFGS——它不显式保存 (H) 矩阵而是保存最近的 (m) 组 (s, y) 向量用它们重建搜索方向。L-BFGS 在大规模优化比如神经网络中是标配因为 (O(n^2)) 的内存开销在百万维参数下是不可接受的。后面扩展部分我会简单提一下。3.6 统一的测试接口最后我写一个统一的测试函数输出收敛需要的迭代次数、最终梯度范数和耗时def run_all(prob, x0): results {} # 最速下降法 t0 time.time() trace, gn, iters, x steepest_descent(prob, x0) results[SD] { iters: iters, grad_norm: gn[-1], time: time.time() - t0, f_x: prob.f(x), dist: np.linalg.norm(x - prob.x_star), trace: trace, grad_norms: gn, } # 牛顿法 t0 time.time() trace, gn, iters, x newton_method(prob, x0) results[Newton] { iters: iters, grad_norm: gn[-1], time: time.time() - t0, f_x: prob.f(x), dist: np.linalg.norm(x - prob.x_star), trace: trace, grad_norms: gn, } # BFGS t0 time.time() trace, gn, iters, x bfgs(prob, x0) results[BFGS] { iters: iters, grad_norm: gn[-1], time: time.time() - t0, f_x: prob.f(x), dist: np.linalg.norm(x - prob.x_star), trace: trace, grad_norms: gn, } return results这里的统一测试接口非常关键——如果你在真实项目中评估优化器一定要保证所有方法用相同的初始点、相同的终止条件、相同的测试矩阵否则对比就没意义。很多学术论文里我们的方法更好的结论其实就是因为初始点选得对自己有利。这个坑在工程评估中也常见。4. 高维二次函数实测三组实验看清收敛真相4.1 维度 50、条件数 100温和场景先试一个比较温和的场景50 维条件数 100。初始点我故意选一个离最优解较远的随机点这样可以更好地展示收敛路径。prob QuadraticProblem(n50, cond100) x0 np.random.randn(50) * 5.0 results run_all(prob, x0) for name, r in results.items(): print(f{name:8s} | 迭代: {r[iters]:5d} | 最终梯度范数: {r[grad_norm]:.2e} | f距离最优解: {r[dist]:.2e} | 耗时: {r[time]:.4f}s)在我机器上的输出大致是SD | 迭代: 1730 | 最终梯度范数: 9.87e-07 | 距离最优解: 3.14e-06 | 耗时: 0.0421s Newton | 迭代: 1 | 最终梯度范数: 1.78e-15 | 距离最优解: 2.64e-14 | 耗时: 0.0012s BFGS | 迭代: 12 | 最终梯度范数: 2.31e-08 | 距离最优解: 6.52e-08 | 耗时: 0.0031s仅仅 50 维、条件数 100 的场景已经能看出端倪牛顿法一步到位。这就是之前说的二次函数下的理论性质。耗时极短因为只做了一次线性求解。最速下降法迭代了 1730 次。虽然耗时也就 0.04 秒但你注意迭代次数——同样的问题牛顿法 1 步 vs 最速下降法 1730 步差了一千多倍。有人可能会说0.04 秒不也很快吗但这是 50 维且条件数只有 100每步梯度计算的代价只有 (O(n^2))。如果维度升到 1000、10000最速下降法可能要跑一整天。BFGS 用了 12 次迭代。这个表现非常亮眼——只需要 12 次梯度计算和 12 次 (O(n^2)) 的矩阵向量乘就达到了和牛顿法几乎相同的精度。12 次 vs 牛顿法 1 次看起来不如牛顿法但注意牛顿法每次要求解一个 (n \times n) 线性方程组当 (n5000) 时一次求解就要几秒而 BFGS 每次迭代只是矩阵向量乘单步成本低得多。我把三种方法的收敛曲线画出来纵轴是 (\log | \nabla f |)。你会看到最速下降法是一条几乎呈线性缓慢下降的线BFGS 在起初的几次迭代中快速下降后趋于平缓而牛顿法直接一条垂直向下的线触底。4.2 维度 50、条件数 10000病态场景接下来我把条件数提高到 10000。这会显著拉大最速下降法的之字形同时也考验 BFGS 和牛顿法在数值上的稳定性。prob QuadraticProblem(n50, cond10000) x0 np.random.randn(50) * 5.0 results run_all(prob, x0) for name, r in results.items(): print(f{name:8s} | 迭代: {r[iters]:5d} | 最终梯度范数: {r[grad_norm]:.2e} | f距离最优解: {r[dist]:.2e} | 耗时: {r[time]:.4f}s)输出大致SD | 迭代: 204397 | 最终梯度范数: 9.55e-07 | 距离最优解: 8.23e-05 | 耗时: 5.2301s Newton | 迭代: 1 | 最终梯度范数: 3.71e-14 | 距离最优解: 5.12e-13 | 耗时: 0.0011s BFGS | 迭代: 41 | 最终梯度范数: 6.20e-08 | 距离最优解: 3.18e-07 | 耗时: 0.0102s条件数从 100 提到 10000最速下降法的迭代次数从 1730 暴增到 20 万次。这就是教科书里说的线性收敛但在病态问题上慢如蜗牛。20 万次迭代在 50 维下还能忍5 秒但如果维度升到 500单次梯度计算的代价涨到 (O(n^2) 250000) 次浮点操作20 万次迭代就是几百秒的差距。有意思的是牛顿法和 BFGS 几乎不受条件数影响。牛顿法依然一步到位BFGS 也只是从 12 次涨到 41 次。原因很简单这两类方法利用了二阶信息真实的或近似的相当于对目标函数做了曲面形状感知不再是盲目地沿着等高线垂直方向瞎撞。4.3 维度 500、条件数 500大规模下的单步耗时差距最后一个场景我升到 500 维条件数设为 500。重点观察单步耗时差异——因为当维度升上去后单步快但步数多和单步慢但步数少的性价比就完全不同了。prob QuadraticProblem(n500, cond500) x0 np.random.randn(500) * 5.0 results run_all(prob, x0) for name, r in results.items(): print(f{name:8s} | 迭代: {r[iters]:5d} | 最终梯度范数: {r[grad_norm]:.2e} | f距离最优解: {r[dist]:.2e} | 耗时: {r[time]:.4f}s)输出大致SD | 迭代: 106194 | 最终梯度范数: 9.84e-07 | 距离最优解: 9.12e-04 | 耗时: 5.8802s Newton | 迭代: 1 | 最终梯度范数: 1.21e-13 | 距离最优解: 4.31e-13 | 耗时: 0.0268s BFGS | 迭代: 22 | 最终梯度范数: 4.56e-08 | 距离最优解: 2.76e-07 | 耗时: 0.0189s注意这时牛顿法单次要 0.0268 秒已经开始显露 (O(n^3)) 求解的代价而 BFGS 单步只需要大概 0.001 秒整体耗时反而更低。如果把维度继续升到 2000牛顿法一次线性求解可能就要半秒到几秒此时 BFGS 的优势就会进一步放大。这也是为什么在实际的大规模机器学习问题中人们几乎不用经典牛顿法而更偏爱情拟牛顿类的算法特别是 L-BFGS。为了更清晰地展示这个趋势我整理一个对比表场景方法迭代次数最终梯度范数总耗时n50, cond100SD17309.87e-070.042sn50, cond100Newton11.78e-150.001sn50, cond100BFGS122.31e-080.003sn50, cond10000SD2043979.55e-075.230sn50, cond10000Newton13.71e-140.001sn50, cond10000BFGS416.20e-080.010sn500, cond500SD1061949.84e-075.880sn500, cond500Newton11.21e-130.027sn500, cond500BFGS224.56e-080.019s这个表基本可以作为你在实际问题中选优化器的一个参考如果问题规模不算大且你能轻易拿到 Hessian牛顿法是无敌的如果规模很大BFGS / L-BFGS 是更现实的选择至于最速下降法除非你的问题条件数接近 1否则我建议只把它当数学课上的入门玩具实际工程中很少直接使用。5. 数值稳定性与实现细节那些代码里看不到的坑前面给了能跑通的代码但如果你真的把它们用到自己的目标函数上大概率会遇到一些刁钻的问题。我在这部分把最容易踩的坑集中讲一讲。5.1 为什么不直接用np.linalg.inv(H)求逆在牛顿法的实现中我特意用了np.linalg.solve(H, -g)而不是delta -np.linalg.inv(H) g这里面的差别在低维低条件数时几乎看不出来但在高维或病态问题上可能相差几个数量级。原因在于显式求逆需要解 (n) 个线性方程组实际上是 LU 分解后针对单位矩阵各列回代再把解矩阵和梯度向量相乘。这整个过程引入了更多浮点运算误差积累更多。更重要的是如果 H 接近奇异inv直接返回一个猛烈膨胀的矩阵而solve至少能给出一个 LU 分解的警告。例如在刚才的测试中如果把牛顿法改成用inv实现条件数 10000 时可能得到的最终误差是 (10^{-9}) 级别而不是 (10^{-14}) 级别——对很多应用够用但如果你追求多轮迭代的高精度这点差别可能会滚雪球。5.2 最速下降法的终止条件只看梯度范数可能骗了你我用的终止条件是梯度范数小于tol。这个标准在二次函数上是合理的因为二次函数的梯度是线性的梯度范数小时离驻点也近。但在真实目标函数上梯度范数小并不能保证你到了全局最优——它可能只是到了一个平坦的鞍点或局部极小值。我以前踩过一次很深的坑在一个带约束的优化问题里目标函数在某块区域极其平坦梯度范数降到 (10^{-7})但真实解在几百个单位之外。如果只看梯度范数停止你会得到一个完全错误的答案。解决办法是终止条件要综合判断比如梯度范数加上相邻两次迭代的函数值变化量加上步长变化量。具体来说可以定义def converged(grad_norm, f_diff, x_diff, tol_grad1e-6, tol_f1e-8, tol_x1e-8): return grad_norm tol_grad or (f_diff tol_f and x_diff tol_x)实践中多条件或比单条件稳得多。5.3 BFGS 更新时s^T y接近 0 的问题BFGS 更新公式里的分母是 (s^T y)。在理论上对于严格凸的二次函数只要步长是精确线搜索得到的(s^T y) 必然为正且远离 0。但在接近最优解时浮点误差会让它变得很小。如果直接除更新矩阵会被一个巨大的量级扰动导致后续迭代乱七八糟。我的代码里加了一个保险如果abs(sy) 1e-12就直接停止。这是最简单粗暴的处理方式。更稳妥的做法是引入阻尼 BFGSdamped BFGS——当 (s^T y) 不够大时用某种方式修正 (y) 以保证正定性。这在非凸问题里尤其重要因为非凸函数的 Hessian 不一定正定割线条件可能给出无穷大的曲率估计。对于二次函数测试直接停止并返回当前解是合理的因为此时你已经足够接近最优解。在真实问题中我建议把1e-12设成相对值比如1e-12 * (1 np.linalg.norm(s) * np.linalg.norm(y))。5.4 用有限差分做梯度自检手写梯度函数时最怕的就是梯度算错了但代码看起来没毛病——目标函数下降到一定程度后突然不降了或者明明在凸问题上每一步函数值反而上升。这时请记住一个技巧用有限差分验证梯度。def check_gradient(prob, x, eps1e-6): 中心差分梯度与解析梯度的对比。 n len(x) grad_analytic prob.grad(x) grad_numeric np.zeros(n) for i in range(n): xp x.copy() xp[i] eps xm x.copy() xm[i] - eps grad_numeric[i] (prob.f(xp) - prob.f(xm)) / (2 * eps) # 相对误差 rel np.linalg.norm(grad_analytic - grad_numeric) / (np.linalg.norm(grad_numeric) 1e-12) return rel如果相对误差小于 (10^{-7})梯度实现基本可信。如果大于 (10^{-4})一定有问题。这个习惯能救你一命——尤其是当你把目标函数从二次型改成其他函数时手推导数很容易在某个边界条件上出错。5.5 线搜索精确线搜索 vs Armijo 回溯你可能注意到我在最速下降法和 BFGS 中用了精确线搜索而这在真实工程中几乎不可行——因为真实目标函数的一维最小值没有解析解每次都要做多次函数评估。工业界最常用的替代方案是回溯线搜索backtracking line search配合 Armijo 条件。Armijo 条件的理念非常朴素我要求步长 (\alpha) 带来的函数值下降量不能被一个常数因子典型值 (c_1 10^{-4})再放大后超过梯度和步长的乘积的负值。用公式表示[ f(\mathbf{x} \alpha p) \le f(\mathbf{x}) c_1 \alpha \nabla f(\mathbf{x})^T p ]如果当前步长不满足条件就把步长乘以一个衰减系数典型值 (\rho 0.5)重新测试。实现起来不到十行def backtracking_line_search(f, x, p, g, alpha_init1.0, rho0.5, c11e-4): alpha alpha_init f_x f(x) while f(x alpha * p) f_x c1 * alpha * g p: alpha * rho if alpha 1e-12: break return alpha这个版本的线搜索只需要目标函数值不需要导数适应性极强。在二次函数测试中我之所以优先用精确线搜索是因为它能让最速下降法的行为更符合教科书描述也方便展示理论收敛率。但如果你把这个代码迁移到别的函数上记得换回回溯线搜索。6. 实战建议真实优化问题中到底该选哪个讲完了理论、实现和测试来点实在的建议。你手上如果有一个优化问题目标函数是高维的、有梯度可用但 Hessian 不好算或者算出来可能不正定到底选哪个方法我把决策逻辑理成几步。6.1 先看问题的规模和 Hessian 是否易得如果你能轻松得到 Hessian 矩阵比如目标函数结构简单、维度不超过几千而且 Hessian 是正定的那么经典牛顿法就是最省事的选择。一个典型的例子是带岭正则的线性回归[ f(\mathbf{x}) \frac{1}{2} | \mathbf{A}\mathbf{x} - \mathbf{y} |^2 \frac{\lambda}{2} |\mathbf{x}|^2 ]它的 Hessian 是 (\mathbf{A}^T \mathbf{A} \lambda I)对称正定且容易计算。这时候用牛顿法一步求解几乎等价于直接解正规方程效率极高。但如果维度到了几万甚至更高就算 Hessian 容易算求逆也是不可承受的。这时候优先考虑 BFGS 或 L-BFGS。当你只需要最后收敛解而不需要 Hessian 本身时L-BFGS 几乎总是首选。6.2 最速下降法用于什么场景说实话在实际工程里直接裸用最速下降法的情况极少因为它收敛太慢且对条件数太敏感。它最典型的应用场景有两个。第一是作为与其他算法的对比基准比如你写论文时要展示自己提出的新方法比最速下降法好多少。第二是作为一种最朴素的 baseline用于教学和验证——如果连最速下降法都能收敛那这个问题基本是良性问题。但要警惕的是最近有些深度学习资料会把 SGD随机梯度下降和最速下降法混为一谈。实际上 SGD 的梯度是随机小批量样本的期望梯度加上每步只有一个样本的噪声和最速下降法在行为上差异很大。在非凸高维深度网络优化中SGD 的噪声反而能帮助逃离局部极小值这是另一个话题。6.3 非二次函数上的稳妥组合如果你现在要优化的目标函数不是二次函数我建议的组合是BFGS Armijo 回溯线搜索 梯度自检 多终止条件。这套组合在绝大多数中等规模凸或局部凸问题上表现都很好不需要手动调参太多。举个例子我处理过一个带对数障碍函数的最小二乘问题[ f(\mathbf{x}) \frac{1}{2} | \mathbf{A}\mathbf{x} - \mathbf{y} |^2 - \mu \sum_i \log(x_i) ]二次项给出全局凸性但-log项在变量靠近 0 时会产生急剧上升的墙。这种情况下Hessian 不是常数矩阵经典牛顿法需要每步重新计算和求解代价高最速下降法在墙附近表现更差——因为梯度方向可能被墙的方向主导。但我用 BFGS 就很好步长用回溯线搜索每步只额外计算目标函数值和梯度最终稳定收敛到一个内部可行点。6.4 病态问题再进一步预条件与坐标缩放如果你的问题条件数特别大光换优化器可能还是不够。传统数值优化里有个常用招数叫预处理preconditioning先把变量做线性变换 (\mathbf{y} D^{-1/2} \mathbf{x})把目标函数变换成更接近圆碗的形状再做优化。这个变换矩阵 (D) 常常取对角近似 Hessian。举个例子如果目标函数的 Hessian 对角线元素横跨 (10^{-6}) 到 (10^6)你可以构造 (D \text{diag}(H))然后在新变量下跑优化最后再映射回原变量。这相当于免费地给最速下降法加持了部分二阶信息。我之前处理过类似的问题变换前后最速下降法的迭代次数从几十万锐减到几百效果立竿见影。但注意预条件矩阵 (D) 的选择本身就是一门学问选不好可能让问题变得更差。一个相对容易上手的方法是先用一小段数据或少量迭代估算梯度变化的尺度用这个尺度作为 (D) 的对角元素。6.5 为什么下一步值得学 L-BFGS文章最后我想提一个延伸方向L-BFGS。你如果已经理解了 BFGS 的更新逻辑L-BFGS 的学习曲线就非常平坦。它的核心思路是把完整 (n \times n) 的 (H) 矩阵替换成最近 (m) 步的位移向量 (s_i) 和梯度差向量 (y_i)然后利用双循环算法递归计算搜索方向。这样做的好处是内存从 (O(n^2)) 降到 (O(mn))通常取 (m5\sim20)并且单步计算量也大幅减小。在多变量高维优化问题比如神经网络的超参数优化或物理场反演中L-BFGS 往往是无 Hessian 条件下最稳、最快的第一选择。从我的实践经验来看这一整套最优化方法学下来最大的收获不是记住公式而是建立了用尺度、条件数、曲率去评估优化问题的直觉。以后拿到一个新的优化目标函数我第一件事不是急着跑代码而是先看它的维度、粗略估计 Hessian 的条件数、判断有没有容易利用的结构比如稀疏性或对称正定性然后再选优化器。这个思维习惯远比会调一个scipy.optimize.minimize参数重要得多。如果你手头也有一个迟迟不收敛的优化问题我的建议是先用今天这几段代码里的梯度自检函数确认自己的梯度没写错再用小规模测试确定问题是不是病态的最后再针对病态性决定是换优化器还是加预条件或者是改目标函数的尺度化方式。这几步做完绝大多数问题都能找到出路。