从零发明数值最优化 6：高斯-牛顿法与列文伯格-马夸尔特法

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

我们之前说过很多普适的优化方法，大概就是给定一个任意的连续函数 ϕ(x)\phi(x)，然后我们盯着它优化就可以了。不过对于一些形式比较好的函数，比如说 ϕ(x)=x⊤Ax−b⊤x\phi(x) = x^\top A x - b^\top x，我们就不必用 GD 或是什么 L-BFGS，而是可以用共轭梯度法针对性地进行超级加速。

今天我们要来研发的另外一类形式也不错的函数，即最小二乘问题：

min⁡x∈Rnϕ(x):=12∥r(x)∥22=12∑i=1mri(x)2(m≥n)\min_{x \in \mathbb{R}^n} \phi(x) := \frac{1}{2} \Vert r(x) \Vert _2^2 = \frac{1}{2}\sum_{i=1}^m r_i(x)^2 \quad (m \ge n)

高斯-牛顿法（GN 法）

在高斯的年代，拟牛顿法还没有出生呢！但当时的人们已经发现，梯度下降法总是会遇到不收敛之类的麻烦问题，所以牛顿法就成为了当时的首选。那么我们不妨先用牛顿法试一下。牛顿法的核心更新公式是：

∇2ϕ(x)⋅p=−∇ϕ(x)\nabla^2 \phi(x) \cdot p = -\nabla \phi(x)

我们先求梯度：

dϕ(x)=r⊤(x)d(r(x))=r⊤(x)∂r∂xdx=tr⁡(((∂r∂x)⊤r(x))⊤dx)d \phi (x) = r^\top(x) d(r(x)) = r^\top(x) \frac{\partial r}{\partial x} dx = \operatorname{tr} \left(\left(\left(\frac{\partial r}{\partial x}\right)^\top r(x)\right)^\top dx\right)

所以 ∇ϕ(x)=(∂r∂x)⊤r(x)\nabla \phi(x) = \left(\frac{\partial r}{\partial x}\right)^\top r(x)。其中，∂r∂x\frac{\partial r}{\partial x} 就是我们熟知的 Jacobian 矩阵 J(x)=[∇r1(x)⊤⋮∇rm(x)⊤]J(x) = \begin{bmatrix} \nabla r_1(x)^\top \\ \vdots \\ \nabla r_m(x)^\top \end{bmatrix}。为了方便书写，我们之后就写成 J(x)J(x) 的形式了。

总而言之，梯度就是 ∇ϕ(x)=J(x)⊤r(x)\nabla \phi(x) = J(x)^\top r(x)。

梯度看上去挺简单，那 Hessian 矩阵呢？

d2ϕ(x)=d(dϕ(x))=d(r⊤Jdx)=(dr)⊤Jdx+r⊤(dJ)dxd^2\phi(x) = d(d\phi(x)) = d\big(r^\top J dx\big) = (dr)^\top J dx + r^\top (dJ) dx

(dr)⊤Jdx=(Jdx)⊤Jdx=dx⊤J⊤Jdx(dr)^\top J dx = (J dx)^\top J dx = dx^\top J^\top J dx

dJ=[d(∇r1(x))⊤⋮d(∇rm(x))⊤]dJ = \begin{bmatrix} d(\nabla r_1(x))^\top \\ \vdots \\ d(\nabla r_m(x))^\top \end{bmatrix}

因为 d(∇ri(x))=∂∇ri(x)∂xdx=∇2ri(x)dxd(\nabla r_i(x)) = \frac{\partial \nabla r_i(x)}{\partial x} dx = \nabla^2 r_i(x) dx，且 ∇2ri(x)\nabla^2 r_i(x) 是对称矩阵，所以

dJ=[(∇2r1(x)dx)⊤⋮(∇2rm(x)dx)⊤]=[dx⊤∇2r1(x)⋮dx⊤∇2rm(x)]dJ = \begin{bmatrix} (\nabla^2 r_1(x) dx)^\top \\ \vdots \\ (\nabla^2 r_m(x) dx)^\top \end{bmatrix} = \begin{bmatrix} dx^\top \nabla^2 r_1(x) \\ \vdots \\ dx^\top \nabla^2 r_m(x) \end{bmatrix}

r⊤(dJ)dx=∑i=1mri(x)(dx⊤∇2ri(x)dx)=dx⊤(∑i=1mri(x)∇2ri(x))dxr^\top (dJ) dx = \sum_{i=1}^m r_i(x) \big(dx^\top \nabla^2 r_i(x) dx\big) = dx^\top \left( \sum_{i=1}^m r_i(x) \nabla^2 r_i(x) \right) dx

d2ϕ(x)=dx⊤(J(x)⊤J(x)+∑i=1mri(x)∇2ri(x))dxd^2\phi(x) = dx^\top \left( J(x)^\top J(x) + \sum_{i=1}^m r_i(x) \nabla^2 r_i(x) \right) dx

对照定义 d2ϕ(x)=dx⊤(∇2ϕ(x))dxd^2\phi(x) = dx^\top \left(\nabla^2\phi(x)\right) dx，可以求得

∇2ϕ(x)=J(x)⊤J(x)+∑i=1mri(x)∇2ri(x)\nabla^2 \phi(x) = J(x)^\top J(x) + \sum_{i=1}^m r_i(x) \nabla^2 r_i(x)

看上去 Hessian 矩阵的形式非常不妙。问题在两方面：第一，实际问题中 mm 经常成千上万，巨大无比。迭代一次就要算千把个 n×nn \times n 的矩阵加法，这谁受得了？第二，后半项的求和可能让 Hessian 矩阵变成非正定甚至负定，直接向最高点处迭代了。

修剪 Hessian 矩阵

让我们从后人的视角来重新回顾这个问题。让我们来想一想拟牛顿法（虽然它的出现时间在高斯-牛顿法之后，不过我们确实先学的它嘛！）。拟牛顿法其实没有很在意这个 Hessian 矩阵，而是用了一种差不多的方式对它做拟合，用的是自己捣鼓出来的 BB。

让我们学习一下后人的思路（笑），也搞一个对 Hessian 矩阵的拟合出来。我们嫌弃现在看到的这个 Hessian 矩阵太丑了，所以让我们直接把后面那串求和咔吧一下裁掉，只留下一个 J(x)⊤J(x)J(x)^\top J(x)。

J(x)⊤J(x)J(x)^\top J(x) 的性质超级好。首先，它天然就是半正定的，只要列满秩就严格正定，从而保证计算出的方向是一个严格的局部下降方向；其次，我们算它不额外消耗时间，因为梯度本来就需要 J(x)J(x) 参与计算，做个矩阵乘法也是顺手的事。

所以，我们得到了高斯-牛顿法的更新公式：

(J(xk)⊤J(xk))pk=−J(xk)⊤r(xk)(J(x_k)^\top J(x_k)) p_k = -J(x_k)^\top r(x_k)

残差线性化

我们还有第二种认识 GN 法的视角。

先回到我们需要最小化的那个公式 min⁡x∈Rnϕ(x):=12∥r(x)∥22\min_{x \in \mathbb{R}^n} \phi(x) := \frac{1}{2} \Vert r(x) \Vert _2^2。如果它里面的结构是 min⁡p12∥Ap−b∥22\min_p \frac{1}{2} \Vert A p - b \Vert _2^2——这简直是世界上最简单的问题，解一个正规方程 A⊤Ap=A⊤bA^\top A p = A^\top b 直接得到答案了。

那么，其实难点就在于 r(x)r(x) 是非线性的。既然它是非线性的，我们能不能把它转换成线性的呢？

于是想到一阶泰勒展开：

r(xk+pk)≈r(xk)+J(xk)pkr(x_k + p_k) \approx r(x_k) + J(x_k) p_k

原来的非线性问题就变成了局部的线性最小二乘逼近子问题：

min⁡pkϕ~(pk):=12∥r(xk)+J(xk)pk∥22\min_{p_k} \tilde{\phi}(p_k) := \frac{1}{2} \Vert r(x_k) + J(x_k) p_k \Vert _2^2

对步长求导，我们得到的还是上面已经出现过的那个式子：

(J(xk)⊤J(xk))pk=−J(xk)⊤r(xk)(J(x_k)^\top J(x_k)) p_k = -J(x_k)^\top r(x_k)

高斯-牛顿法的问题

然而，看似精巧的 GN 法也存在不少问题。其一，如果雅可比矩阵 J(xk)J(x_k) 不是列满秩的，那么 J(xk)⊤J(xk)J(x_k)^\top J(x_k) 就会变成不可逆的奇异矩阵，使得方程解不出唯一的 pkp_k。其二，线性化假设只在 pkp_k 很小的局部邻域才成立。如果当前点距离最优解很远， GN 法算出来的 pkp_k 可能会巨大无比，下一个迭代点都不知道飞到哪里去了。

列文伯格-马夸尔特法（LM 法）

GN 法的第二个问题出在哪里？出在它笃定地相信这个线性模型对全空间都成立。然而，事实上泰勒展开只在一个很小的邻域内靠谱。一旦超出这个区域，模型就在胡说八道。

约束步长

我们希望步子不要跨太大。最简单的想法就是给步长加上一个约束项，不许它迈得太长。所以我们想到机器学习中用来预防过拟合的岭回归。引入阻尼系数 μk>0\mu_k > 0，并重定义最小二乘目标：

min⁡pkΦk(pk):=12∥r(xk)+J(xk)pk∥22+μk2∥pk∥22\min_{p_k} \Phi_k(p_k) := \frac{1}{2} \Vert r(x_k) + J(x_k) p_k \Vert _2^2 + \frac{\mu_k}{2} \Vert p_k \Vert _2^2

对目标函数求导，得到

∇pΦk(pk)=J(xk)⊤(r(xk)+J(xk)pk)+μkpk\nabla_p \Phi_k(p_k) = J(x_k)^\top (r(x_k) + J(x_k) p_k) + \mu_k p_k

令梯度等于零，则有

(J(xk)⊤J(xk)+μkI)pk=−J(xk)⊤r(xk)\big(J(x_k)^\top J(x_k) + \mu_k I\big) p_k = -J(x_k)^\top r(x_k)

这就是列文伯格-马夸尔特法的核心方程。

列文伯格-马夸尔特法的优势

GN 法的两大问题在于 J(xk)⊤J(xk)J(x_k)^\top J(x_k) 不可逆导致多解和 pkp_k 过大导致迭代不保证下降，而这两个问题在 LM 法中都得到了控制。由于 μk>0\mu_k > 0，J(xk)⊤J(xk)+μkIJ(x_k)^\top J(x_k) + \mu_k I 绝对是正定的，于是消除了 GN 法求逆崩溃的可能性。而我们计算迭代方向和梯度的内积 ∇ϕ(x)⊤pk=−∇ϕ(x)⊤(J(xk)⊤J(xk)+μkI)−1∇ϕ(x)\nabla \phi(x)^\top p_k = -\nabla \phi(x)^\top (J(x_k)^\top J(x_k) + \mu_k I)^{-1} \nabla \phi(x)，不难发现只要梯度不为零二次型就严格为负，所以 LM 方向总是指向目标函数下降的方向。

另外，μk\mu_k 带来的连续平滑插值也提供了很好的性质。在 μk→0\mu_k \to 0 时，pk→−(J(xk)⊤J(xk))−1∇ϕ(xk)=pkGNp_k \to -\big(J(x_k)^\top J(x_k))^{-1} \nabla \phi(x_k)= p_k^{\text{GN}}；而在 μk→∞\mu_k \to \infty 时，pk→−1μk∇ϕ(x)p_k \to -\frac{1}{\mu_k} \nabla \phi(x) 退化为梯度下降法，保证不越界。

μk\mu_k 如何自适应更新？

既然 μk\mu_k 是核心参数，我们应该怎么知道让它变大还是变小呢？我们通过分析增益比来控制它。

定义真实下降量 Δϕact=ϕ(xk)−ϕ(xk+pk)\Delta \phi_{\text{act}} = \phi(x_k) - \phi(x_k + p_k)，线性模型预测下降量 Δϕpred=qk(0)−qk(pk)\Delta \phi_{\text{pred}} = q_k(0) - q_k(p_k)（记线性模型 qk(pk)=12∥rk+Jkpk∥22q_k(p_k) = \frac{1}{2} \Vert r_k + J_k p_k \Vert _2^2），增益比即为两者之比：

ρk=ΔϕactΔϕpred=ϕ(xk)−ϕ(xk+pk)12∥rk∥22−12∥rk+Jkpk∥22\rho_k = \frac{\Delta \phi_{\text{act}}}{\Delta \phi_{\text{pred}}} = \frac{\phi(x_k) - \phi(x_k + p_k)}{\frac{1}{2}\Vert r_k \Vert _2^2 - \frac{1}{2}\Vert r_k + J_k p_k \Vert _2^2}

  • 若 ρk>0\rho_k > 0 且较大（如 ρk>0.75\rho_k > 0.75）：说明实际下降与模型很吻合，线性模型极其可信，所以接受步长并减小 μk\mu_k，让算法更大胆，靠近 GN 法以加速迭代。
  • 若 ρk>0\rho_k > 0 但较小（如 ρk<0.25\rho_k < 0.25）：尽管下降，但实际下降落后于预期，这意味着非线性误差影响有点大，于是接受步长并增大 μk\mu_k。
  • 若 ρk≤0\rho_k \le 0：线性近似完全失效，拒绝该步长并大幅增大 μk\mu_k，重新解方程。