从零发明数值最优化 3：牛顿法

发表于 更新于 分类于 ， ， ， ，

在梯度下降法中，我们只使用了一阶的泰勒展开信息，那如果使用二阶信息呢？

f(xk+pk)≈f(xk)+∇f(xk)⊤pk+12pk⊤∇2f(xk)pkf(x_k + p_k) \approx f(x_k) + \nabla f(x_k)^\top p_k + \frac 1 2 p_k^\top \nabla^2 f(x_k) p_k

我们将其视作一个关于 pkp_k 的函数，那么我们得到 ψ(pk)=12pk⊤∇2f(xk)pk+∇f(xk)⊤pk+f(xk)\psi(p_k) = \frac 1 2 p_k^\top \nabla^2 f(x_k) p_k + \nabla f(x_k)^\top p_k + f(x_k)。我们希望 ψ(pk)\psi(p_k) 最小，那么就肯定有 ∇ψ(pk)=0\nabla \psi(p_k) = 0，即

∇2f(xk)pk+∇f(xk)=0\nabla^2 f(x_k) p_k + \nabla f(x_k) = 0

另外，如果 ψ(pk)\psi(p_k) 存在唯一最小值，那么它的二次项系数 ∇2f(xk)\nabla^2 f(x_k) 必须是正定的。而正定矩阵总满秩，所以当 ∇2f(xk)\nabla^2 f(x_k) 正定时，搜索方向唯一确定：

pkNewton=−(∇2f(xk))−1∇f(xk)p_k^{\text{Newton}} = - \left( \nabla^2 f(x_k) \right)^{-1} \nabla f(x_k)

收敛速度

利用一阶泰勒展开和积分中值定理，梯度 ∇f(x⋆)\nabla f(x^\star) 在 xkx_k 处的展开式为

0=∇f(x⋆)=∇f(xk)+∫01∇2f(xk+t(x⋆−xk))(x⋆−xk)dt0 = \nabla f(x^\star) = \nabla f(x_k) + \int_0^1 \nabla^2 f(x_k + t(x^\star - x_k)) (x^\star - x_k) dt

将牛顿步 ∇f(xk)=∇2f(xk)(xk−xk+1)\nabla f(x_k) = \nabla^2 f(x_k)(x_k - x_{k+1}) 代入：

∇2f(xk)(xk+1−x⋆)=∫01[∇2f(xk+t(x⋆−xk))−∇2f(xk)](x⋆−xk) dt\nabla^2 f(x_k)(x_{k+1} - x^\star) = \int_0^1 \left[ \nabla^2 f(x_k + t(x^\star - x_k)) - \nabla^2 f(x_k) \right] (x^\star - x_k) \, dt

假设 Hessian 矩阵满足 L2L_2-Lipschitz 连续（即 ∥∇2f(x)−∇2f(y)∥≤L2∥x−y∥\Vert \nabla^2 f(x) - \nabla^2 f(y) \Vert \le L_2 \Vert x - y\Vert），且最优解附近存在逆上界 ∥∇2f(xk)−1∥≤M0\Vert \nabla^2 f(x_k)^{-1} \Vert \le M_0：

∥xk+1−x⋆∥≤M0∫01L2⋅t∥xk−x⋆∥2dt=M0L22∥xk−x⋆∥2\Vert x_{k+1} - x^\star \Vert \le M_0 \int_0^1 L_2 \cdot t \Vert x_k - x^\star \Vert^2 dt = \frac{M_0 L_2}{2} \Vert x_k - x^\star \Vert^2

令 C=M0L22C = \frac{M_0 L_2}{2}，则有

∥xk+1−x∗∥≤C∥xk−x∗∥2\Vert x_{k+1} - x^*\Vert \le C \Vert x_k - x^*\Vert^2

所以，牛顿法的收敛速度是二阶收敛的。

算法瓶颈分析

时间复杂度

Hessian 矩阵有 n(n+1)2\frac{n(n+1)}{2} 个独立的元素，所以计算海森矩阵至少要 O(n2)O(n^2) 时间。计算单个二阶偏导数 ∂2f∂xi∂xj\frac{\partial^2 f}{\partial x_i \partial x_j} 的时间取决于目标函数表达式的规模 TevalT_{\text{eval}}。解析求导或自动微分计算完整 Hessian 矩阵的单步耗时通常为 O(n2⋅Teval)O(n^2 \cdot T_{\text{eval}})；若采用有限差分近似，计算 n2n^2 个元素需要评估 O(n)O(n) 次梯度或 O(n2)O(n^2) 次函数值，计算成本远高于简单的 O(n2)O(n^2)。

另外，求解线性方程组需经过 Cholesky 分解或 LU 分解，单步耗时为 O(n3)O(n^3)。

无论如何，时间复杂度至少是 O(n3)O(n^3) 的。

空间复杂度

Hessian 矩阵有 n(n+1)2\frac{n(n+1)}{2} 个独立的元素，所以存储每次的 Hessian 矩阵要 O(n2)O(n^2) 空间。

非正定性缺陷

当 ∇2f(xk)\nabla^2 f(x_k) 出现负特征值，甚至为不可逆矩阵时，pkp_k 可能指向鞍点或极大值。当 λmin≤0\lambda_\text{min} \leq 0 时，工程上一般将精确 Hessian 矩阵换为 ∇2f(xk)+μI,μ>−λmin+ε\nabla^2 f(x_k) + \mu I, \quad \mu > -\lambda_\text{min} + \varepsilon，通过牺牲拟合的精细度来保证解的唯一性。