优化理论之非线性方法-高斯牛顿法与莱文贝格-马夸特方法

Nov 4, 2019 · 优化理论 ·

分享到:

高斯-牛顿法(Guass-Newton Algorithm)与莱文贝格-马夸特方法(Levenberg–Marquardt algorithm)求解非线性最小二乘问题

众所周知，最小二乘法通过最小化误差平方和获得最佳函数。有时候你可能产生疑问，为什么不能通过其他方式获得最优函数，比如说最小化误差的绝对值的和？本文中，我将会从概率的角度解释最小二乘法的依据（参考自andrew ng 《机器学习课程》第三讲）。最小二乘问题可以分为线性最小二乘和非线性最小二乘两类，本文的目标是介绍两种经典的最小二乘问题解法：高斯牛顿法与莱文贝格-马夸特方法。实际上，后者是对前者以及梯度下降法的综合。

最小二乘法的概率解释(probabilistic interpretation)

以线性回归为例，假设最佳函数为 \(y=\mathbf{\theta}^T\mathbf{x}\) (\(\theta,x\)为向量), 对于每对观测结果\((x^{(i)},y^{(i)})\)，都有 \[y^{(i)}=\theta^Tx^{(i)}+\epsilon^{(i)}\] 其中 \(\epsilon\)为误差，基于一种合理的假设（中心极限定理），我们可以认为误差的分布服从正态分布(又称高斯分布)，即 \(\epsilon \sim N(0,\sigma^2)\) ，那么，我们可以认为\(y^{(i)} \sim N(\theta^Tx^{(i)},\sigma^2)\),根据正态分布的概率公式 \[P(y^{(i)}|x^{(i)};\theta)=\frac{1}{\sqrt{2\pi}\sigma}exp(-\frac{(y^{(i)}-\theta^Tx^{(i)})^2}{2\sigma^2})\] 在统计学中，将所有的\(P(y|x)\)累乘作为\(\theta\)的似然函数,用以衡量\(\theta\)的或然性(likelihood) \[L(\theta)=P(y|x;\theta)=\prod _{i=0}^m \frac{1}{\sqrt{2\pi}\sigma}exp(-\frac{(y^{(i)}-\theta^Tx^{(i)})^2}{2\sigma^2})\] 最佳的参数\(\theta\)应该是使得所有数据出现的概率最大的那个，这个过程称之为最大似然估计(maximum likelihood estimation)。为了数学计算上的便利，采用单调函数:log函数\(l(\theta)=log(L(\theta))\)，称为对数似然函数代表\(L(\theta)\): \[l(\theta)=log\prod_{i=0}^m \frac{1}{\sqrt{2\pi}\sigma}exp(-\frac{(y^{(i)}-\theta^Tx^{(i)})^2}{2\sigma^2}) \\ =\sum_{i=0}^mlog\frac{1}{\sqrt{2\pi}\sigma}exp(-\frac{(y^{(i)}-\theta^Tx^{(i)})^2}{2\sigma^2})\\=mlog\frac{1}{\sqrt{2\pi}\sigma} – \frac{1}{2\sigma^2}\sum_{i=0}^m(y^{(i)}-\theta^Tx^{(i)})^2\] 需要指出的是，在Andrew ng的讲解认为\(\sigma\)不影响\(\theta\)的决定，因此，最大化\(l(\theta)\)等同于最小化\(\sum_{i=0}^m(y^{(i)}-\theta^Tx^{(i)})^2\)。最小二乘法通过最小化误差平方和得到最佳函数的方法存在概率论方面的基础。

高斯-牛顿法(Guass-Newton Algorithm)

高斯-牛顿法是在牛顿法基础上进行修改得到的，用来(仅用于)解决非线性最小二乘问题。高斯-牛顿法相较牛顿法的最大优点是不需要计算二阶导数矩阵(Hessian矩阵)，当然，这项好处的代价是其仅适用于最小二乘问题。如下是其推导过程：

最小二乘方法的目标是令残差的平方和最小： \[f(\theta) = \frac{1}{2}\sum_{i=0}^mr(\mathbf{x_i})^2\] 采用牛顿法求解\(f(\theta)\)的最小值，需要计算其梯度向量与Hessian矩阵。 \[\nabla_\theta f =\frac{\partial f}{\partial \theta}=\sum r_i\frac{\partial r_i}{\partial \theta}\\ =\begin{bmatrix} r(x_1)&r(x_2)&\dots&r(r_m)\end{bmatrix} \begin{bmatrix} \nabla_\theta r(x_1)^T \\ \nabla_\theta r(x_2)^T \\ \vdots\\ \nabla_\theta r(x_m)^T \\ \end{bmatrix}\]

其中 \[J_r(\theta)=\begin{bmatrix}\frac{\partial r_j}{\partial\theta_i}\end{bmatrix}_{j=1,\dots,m;i=1,\dots,n}=\begin{bmatrix} \nabla_\theta r(x_1)^T \\ \nabla_\theta r(x_2)^T \\ \vdots\\ \nabla_\theta r(x_m)^T \\\ \end{bmatrix}\] 称为\(r\)的雅各比(Jacobian)矩阵。因此上式可以写作 \[\nabla_\theta f = r^TJ_r=J_r^Tr\] 其中\(r=\begin{bmatrix}r(x_1)&r(x_2)&\dots&r_m\end{bmatrix}^T\)。再看Hessian矩阵的计算： \[H = \begin{bmatrix}\frac{\partial^2f}{\partial \theta^2} \end{bmatrix}=\sum\begin{bmatrix}r_i\frac{\partial^2r_i}{\partial\theta^2}+(\frac{\partial r_i}{\partial \theta})(\frac{\partial r_i}{\partial \theta})^T \end{bmatrix}\] 观察二阶导数项\(r_i\frac{\partial^2r_i}{\partial\theta^2}\)，因为残差 \(r_i\approx 0\),因此我们可以认为此项接近于0而舍去。所以Hessian矩阵可以近似写成： \[H\approx\sum\begin{bmatrix}(\frac{\partial r_i}{\partial\theta})(\frac{\partial r_i}{\partial\theta})^T\end{bmatrix}=J_r^TJ_r\] 这里我们可以看到高斯-牛顿法相对于牛顿法的不同就是在于采用了近似的Hessian矩阵降低了计算的难度，但是同时，舍去项仅适用于最小二乘问题中残差较小的情形。

将梯度向量，Hessian矩阵(近似)带入牛顿法求根公式，得到高斯-牛顿法的迭代式： \[\theta_i = \theta_{i-1}-{(J_r^TJ_r)}^{-1}J_r^Tr\] 只需要计算出\(m\times n\)的Jacobian矩阵便可以进行高斯-牛顿法的迭代，计算已经算是非常简便的了。

Levenberg-Marquart 算法

与牛顿法一样，当初始值距离最小值较远时，高斯-牛顿法的并不能保证收敛。并且当\(J_r^TJ_r\)近似奇异的时候，高斯牛顿法也不能正确收敛。Levenberg-Marquart 算法是对上述缺点的改进。L-M方法是对梯度下降法与高斯-牛顿法进行线性组合以充分利用两种算法的优势。通过在Hessian矩阵中加入阻尼系数\(\lambda\)来控制每一步迭代的步长以及方向： \[(H+\lambda I)\epsilon = -J_r^Tr\]

当\(\lambda\)增大时，\(H+\lambda I\)趋向于\(\lambda I\)，因此\(\epsilon\)趋向于 \(-\lambda J_r^Tr\)，也就是梯度下降法给出的迭代方向；
当\(\lambda\)减小时，\(H+\lambda I\)趋向于\(H\)，\(\epsilon\)趋向于\(-H^{-1}J_r^Tr\)，也就是高斯-牛顿法给出的方向。

\(\lambda\)的大小通过如下规则调节，也就是L-M算法的流程：

L-Mprocessure.png

初始化 \(\theta_0\)，\(\lambda_0\)。
计算当前点\(\theta_i\)处的残差向量\(r_i\)与雅各比矩阵\(J_r\)。
通过求解\((H_i+\lambda I)\epsilon = -J_r^Tr_i\)求解迭代方向\(\epsilon\)。
计算\(\theta_i^\prime=\theta_i+\epsilon\)点处的残差向量\(r_i^\prime\)。
如果\(\|r_i\prime\|^2>\|r_i\|^2\)?,即残差没有下降，则更新\(\lambda = \beta\lambda\)，增大\(\lambda\)重新回到第三步重新求解新的\(\epsilon\)。如果残差下降，则更新\(\theta_{i+1} = \theta_i+\epsilon\) ，到第二步，并且降低\(\lambda=\alpha\lambda\)，增大迭代步长。

在曲线拟合实践中，\(\alpha\)通常选取 0.1，\(\beta\)选取10。

相比于高斯-牛顿法，L-M算法的优势在于非常的鲁棒，很多情况下即使初始值距离(局部)最优解非常远，仍然可以保证求解成功。作为一种阻尼最小二乘解法，LMA(Levenberg-Marquart Algorithm)的收敛速度要稍微低于GNA(Guass-Newton Algorithm)。L-M算法作为求解非线性最小二乘问题最流行的算法广泛被各类软件包实现，例如google用于求解优化问题的库Ceres Solver。后续，我会通过最小二乘圆拟合的案例给出L-M算法的实现细节。

最小二乘法的改进

传统的最小二乘法是针对线性函数与线性场景，典型的场景如下图：

Linear_regression.svg

其需要优化的方程是： \[ \min\|X\theta-y\|^2 \] 其中，\(\theta\)表示参数向量，参数的个数与特征的个数一样（算上常数项的），\(X\)表示\(m\)个观察到的数据，每个数据有\(n\)个特征（算上常数项的），因此\(X\)是一个\(m\times n\)维矩阵。\(y\)是实际观察到的结果向量。详细表示起来则是： \[ X=\begin{bmatrix} ——(x^{(1)})^T——\\ ——(x^{(2)})^T——\\ \vdots\\ ——(x^{(m)})^T——\\ \end{bmatrix}, x^{(i)}=\begin{bmatrix} x^{(i)}_0(=1)\\ x^{(i)}_1\\ \vdots\\ x^{(i)}_{n-1}\\ \end{bmatrix}, \theta=\begin{bmatrix} \theta_0\\ \theta_1\\ \vdots\\ \theta_{n-1}\\ \end{bmatrix}, y=\begin{bmatrix} y_1\\ y_2\\ \vdots\\ y_m\\ \end{bmatrix} \] 根据矩阵计算结果： \[ \theta=(X^TX)^{-1}X^Ty \] 显然，\(X\theta\)只能表示直线结果，对于曲线结果无能为力。针对曲线场景，通常有两种改进思路，第一个就是将特征\(x_j\)的高阶项加入其中，举个简单例子，比如除了常数项，只有一个特征的场景： \[ y=\theta_0+\theta_0 x_1+\theta_2 x^2_1+\theta_3 x^3_1 \] 这样通过高阶多项式就可以拟合曲线，且根据拉格朗日插值法，只有阶数够多，那么就能够很好的拟合任意曲线。然而，高阶项数都是超参数，而且如果让每个特征都有高阶参数，那么参数数量就是几何数量级增长，难以控制参数规模。比如，有两个特征时，引入二次项就会有\(\theta_0+\theta_1x_1+\theta_2x_2+\theta_3x_1^2+\theta_4x_2^2+\theta_5x_1x_2\)共6项，如果到三次方项，那会有更多参数。我们知道，当参数少时会有欠拟合现象，偏差大；参数过多就是过拟合，方差大，因此引入高次方项后，情况变得难以控制。同时，高次函数拟合会出现龙格现象。

在数值分析领域中，龙格现象是在一组等间插值点上使用具有高次多项式的多项式插值时出现的区间边缘处的振荡问题。

曲线拟合

左图：一阶函数欠拟合；中图：二阶函数适当拟合；右图：5阶函数过拟合

第二个改进思路就是使用局部加权线性回归算法(Locally Weighted Linear Regression, LWLR)。它的核心思想就是弯曲点附近的局部曲线，接近整体曲线。我们知道传统的线性规划的代价函数为 \[\sum_i (y^{(i)} − θ^T x^{(i)})^2\] 而局部加权线性回归算法需要体现加权特性，它在不同线性函数之间进行加权，代价函数修改为： \[\sum_i w^{(i)}(y^{(i)} − θ^T x^{(i)})^2\] 这里\(w^{(i)}\)是非负的加权值，可以令\(w^{(i)}=\exp{(-\frac{(x^{(i)}-x)^2}{2\tau^2})}\)。虽然\(w^{(i)}\)有高斯分布形式，但是和高斯分布没有直接联系。显然，当\(|x^{(i)}-x|\)比较小时，\(|x^{(i)}-x|\rightarrow 0\)，此时\(w^{(i)}\rightarrow 1\)；当\(|x^{(i)}-x|\)比较大时，\(w^{(i)}\rightarrow 0\)，权重较小。

局部加权线性回归的优化函数写成矩阵形式为： \[ J(\theta)=(X\theta-y)^TW(X\theta-y) \] 多了一个权重矩阵\(W\)，它是一个对角矩阵，对角元素\(w_{ii}={1\over 2}w^{(i)}\)。根据矩阵求导公式可得： \[ \begin{aligned} \nabla J(\theta)&=\nabla (X\theta-y)^TW(X\theta-y)\\ &=2\frac{\partial{(X\theta-y)}}{\partial{\theta}}W(X\theta-y)\\ &=2X^TW(X\theta-y) \end{aligned} \] 当取极值时，\(\nabla J(\theta)=0\)，即\(2X^TW(X\theta-y)=0\)，所以有： \[ 2X^TW(X\theta-y)=0\\ X^TWX\theta=X^TWy\\ \theta=(X^TWX)^{-1}X^TWy \] 关于矩阵求导的细节，请查阅网页资料线性代数与矩阵-矩阵求导表

此外，整个代价函数中只有一个超参数\(\tau\)，称为带宽参数。它控制着加权参数随着对\(x\)的远离，衰减的速度。当\(\tau\)值比较大时，\(w^{(i)}\)衰减得慢，容易出现欠拟合；当\(\tau\)值比较小时，\(w^{(i)}\)衰减得快，容易出现过拟合。从词义上来看，带宽参数意味着能够容纳周围参数的多少，容纳的越多函数越平缓，当\(\tau\rightarrow \infty\)时，\(w^{(i)}\rightarrow 1\)此时，带宽无限大，容纳了所有数据点，局部加权线性回归算法退化为传统线性规划算法。当\(\tau\)很小时，\(w^{(i)}\)衰减很快，只有在当前点才有效，此时局部加权线性回归算法退化为插值算法。这个算法的问题在于，对于每一个要预测的点，都要重新依据整个数据集计算一个线性回归模型出来，使得算法代价极高。局部加权线性回归效果如下

局部加权线性规划.png

上图：带宽参数过大，欠拟合；中图：适当拟合；下图：带宽参数过小，过拟合

在局部加权线性回归算法中，参数数量\(w^{(i)}\)是随着数据量增加而线性增加的，这个参数不固定的监督学习算法成为非参数学习算法；相对的参数固定的监督学习算法，比如传统的线性规划就称作参数学习算法。