非线性优化（一）：优化方法

JamesB · 发表于 2021-12-29 18:16

讲完了李群和李代数，就轮到最小二乘优化了。和之前一样，应该会是一系列文章。这一章就先讲讲一些基本的概念和常见的最小二乘算法。
1.问题的定义

一个最小二乘问题的定义如下：
寻找一个关于下列函数的局部最小值  ：

$F(x) = \frac{1}{2} \sum^m_{i=1}(f_i (\mathbf{x}) )^2 \quad f_i: \mathbf{R}^n \to \mathbf{R}, \quad m \ge n$

这里， $f_i ( \mathbf{x} )$ 是误差方程。

而局部最小值的定义为：

给定 $F: \mathbf{R}^n \to \mathbf{R}$ ，找到一个  ，使得

$F(\mathbf{x}^*) \le F(\mathbf{x}) \quad \text{for} \ || \mathbf{x} - \mathbf{x}^* || < \delta$

这里， $\delta$ 是一个给定的小正数。

假设损失函数F是光滑、可微的，就可以对其进行泰勒展开:

$F(\mathbf{x} + \mathbf{h}) = F(\mathbf{x}) + \mathbf{h}^T \mathbf{g} + \frac{1}{2} \mathbf{h}^T \mathbf{H} \mathbf{h} + O(|| \mathbf {h} ||^3)$

其中，g是一阶导数（gradient）：

$g = F' (\mathbf{x}) = \begin{bmatrix} \frac{\partial F}{\partial x_1} (\mathbf{x}) \\ \\ \frac{\partial F}{\partial x_2} (\mathbf{x}) \\ \\ \vdots \\ \frac{\partial F}{\partial x_n} (\mathbf{x}) \\ \end{bmatrix}$

而H为海塞矩阵（Hessian）

$\mathbf{H} = F''(\mathbf{x}) = \begin{bmatrix} \frac{ \partial^2 F }{\partial x_i \partial x_j } (\mathbf{x}) \end{bmatrix}$

显然，如果  是一个局部最小值且 $||\mathbf{h}||$ 足够小，那么我们不可能再找到一个更小的 $F(\mathbf{x}^* + \mathbf{h})$ 。

如此，可以先给出一个局部极小值的必要条件：

如果  是一个局部极小值，那么 $\mathbf{g}^* = F'(\mathbf{x}^*) = 0$ 。

我们知道，满足上述条件的点包括：局部极小值，局部极大值和鞍点，统称为驻点（stationary point）。为了判断一个驻点 $\mathbf{x}_s$ 是否是一个局部极小值，将之反代回F的泰勒展开式：

$F( \mathbf{x}_s + \mathbf{h} ) = F( \mathbf{x}_s ) + \frac{1}{2} \mathbf{h}^T \mathbf{H}_s \mathbf{h} + h.o.t$

由海塞矩阵的定义可知，H是一个对称矩阵。如果它同时是一个正定矩阵的话，那么其特征值会大于某个数值 $\delta > 0$ ，并且

$\frac{1}{2} \mathbf{h}^T \mathbf{H}_s \mathbf{h} > \delta ||\mathbf{h}|| ^2$

显然，此时 $F(\mathbf{x}_s + \mathbf{h}) > F(\mathbf{x}_s)$ 。

由此，可以得到一个局部最小值的充分条件:

假定  是一个驻点，且 $F''(\mathbf{x}^*)$ 是一个正定矩阵，那么  是一个局部极小值。

显然，如果H是一个负定矩阵，那么x是一个局部极大值；如果H既有正的特征值，又有负的特征值，那么s就是个鞍点。

2. 迭代梯度下降方法

对于非线性优化方法，其解法基本都是迭代方法：给定一个起始点 $\mathbf{x}_0$ ，不断迭代产生新的 $\mathbf{x}_1, \mathbf{x}_2 \dots$ ，并最终收敛至  。为此，引入下降条件（decending condition）：

$F(\mathbf{x}_{k+1}) < F(\mathbf{x}_k)$

回到F的泰勒展开式，我们有

$F(\mathbf{x} + \alpha \mathbf{h}) \approx F(\mathbf{x}) + \alpha \mathbf{h}^T \mathbf{F}' (\mathbf{x})$

可以看到，如果 $\mathbf{h}^T \mathbf{F}' (\mathbf{x}) < 0$ ，那么h就是F的一个下降方向。

很显然，基本所有下降方法都可以概括为以下两步：

1. 找到一个下降方向 $\mathbf{h}_d$ ，然后
2. 找到一个迭代步长使得F的值减小；

其中，这个步长 $\alpha$ 可以通过线性搜索（line search）的方式得到。

2.1 最速下降法

由前面的等式可以得到

$\lim_{\alpha \to 0} \frac{F(\mathbf{x}) - F(\mathbf{x} + \alpha \mathbf{h} ) }{\alpha ||\mathbf{h}||} = \frac{1}{||\mathbf{h}||} \mathbf{h}^T \mathbf{F}'(\mathbf{x}) = -||\mathbf{F}'(\mathbf{x})|| \cos \theta$

该方程表示的是方程F的相对下降（relative gain）。可以发现，它在 $\theta = \pi$ 的时候可以获得最大的下降。

由此，可以得到最速下降法的下降方向

$\mathbf{h}_{sd} = -\mathbf{F}'(\mathbf{x})$

最速下降法的特点是在刚开始的阶段有较好的表现，但是到了最后是线性收敛。

对此，也诞生了混合方法（hybrid methods）：在开始阶段使用一种方法，在最后阶段使用另一种方法。

2.2 牛顿法

我们知道，对于局部极小值有 $\mathbf{F}'(\mathbf{x}^*)= 0$ 。这是一个非线性方程组，为此，对其进行泰勒展开有：

$\mathbf{F}'(\mathbf{x} + \mathbf{h}) \approx \mathbf{F}'(\mathbf{x}) + \mathbf{F}''(\mathbf{x}) \mathbf{h}$

令上式等于0就可以得到牛顿法：

$\mathbf{H} \mathbf{h}_n = -\mathbf{F}'(\mathbf{x})$

下次迭代就有

$\mathbf{x}: = \mathbf{x} + \mathbf{h}_n$

如果H是正定的（说明上式具有唯一解），即对于所有 $\mathbf{u}>0$ 有 $\mathbf{u}^T \mathbf{H} \mathbf{u}>0$ ，则有

$0 < \mathbf{h}^T_n \mathbf{H} \mathbf{h}_n = -\mathbf{h}^T_n \mathbf{F}'(\mathbf{x})$

说明 $\mathbf{h}_n$ 是损失函数F的一个下降方向。

牛顿法在最后的收敛阶段有很好的表现，在一定条件下（局部最小值附近的H是正定的且x已接近局部最小值），其可以达到二次收敛。

问题在于，如果H是负定的，且附近有一个局部最大值，牛顿法可能会收敛至该局部最大值。我们可以要求损失函数每一步都是减小的，或者使用如下的混合方法：

2.3 置信区域法和阻尼法

使用一个模型L来拟合F在其临域内的行为，则有

$F(\mathbf{x} + \mathbf{h}) \approx L(\mathbf{h}) = F(\mathbf{x}) + \mathbf{h}^T \mathbf{c} + \frac{1}{2} \mathbf{h}^T \mathbf{B}\mathbf{h}$

可以看到，L是F的二次泰勒展开式的一个近似。其中，c来近似梯度，而B为近似的海塞矩阵。

对于置信区域法，我们认为这个模型L在一个限定范围  内是足够精确的。如此，就可以得到对应的下降方向

$\mathbf{h} = \mathbf{h}_{tr} = \arg\min_{||h|| < \Delta} \{ L(\mathbf{h}) \}$

而对于阻尼法，有

$\mathbf{h} = \mathbf{h}_{dm} = \arg\min_h \{ L(\mathbf{h} ) + \frac{1}{2} \mu \mathbf{h}^T \mathbf{h}\}$

阻尼因子 $\mu \ge 0$ ，该阻尼是用来惩罚大步长的。
算法2.4就可以变为

其实就是如果h满足下降条件，则 $\alpha = 1$ ；否则 $\alpha = 0$ 。并且通过调整  或  ，算法可以避免被困在某一点。

考虑到L是对F的临域的模拟，算法更新失败可能就是h（步长）太大了，导致近似不太好。为了衡量L对F的拟合度，定义一个gain ratio

$\rho = \frac{F(\mathbf{x}) - F(\mathbf{x} + \mathbf{h})}{L(0) - L(\mathbf{h})}$

表示的是损失函数真实下降值和预计下降值之间的比例。由于分母是正数，当分子很小或者为负时，说明此时损失函数没有下降多少甚至还上升来，我们就知道当前的h太大了，导致L不是一个好的近似。

对于置信区域法和阻尼法，常见的更新策略为

算法对于其中的常数的设置并不敏感，重要的是选择合适的值使得  或  不震荡。

3. 非线性优化算法

前面的损失函数还可以展开为如下形式：

$F(\mathbf{x}) = \frac{1}{2} \sum_{i=1}^m (f_i (\mathbf{x}) )^2 = \frac{1}{2} || \mathbf{f} ( \mathbf{x} ) ||^2 = \frac{1}{2} \mathbf{f} ( \mathbf{x} )^T \mathbf{f} ( \mathbf{x} )$

误差函数f(x)的泰勒展开式为：

$\mathbf{f} ( \mathbf{x} + \mathbf{h} ) = \mathbf{f} ( \mathbf{x} ) + \mathbf{J} ( \mathbf{x} ) \mathbf{h} + O ( || \mathbf{h} ||^2 )$

其中，雅可以矩阵J可以写为

$( \mathbf{J} ( \mathbf{x} ) )_{ij} = \frac{ \partial f_i }{ \partial x_j } ( \mathbf{x} )$

如此一来，损失函数F的微分就可以写成：

$\frac{ \partial F }{ \partial x_j } ( \mathbf{x} ) = \sum_{i=1}^m f_i ( \mathbf{x} ) \frac{ \partial f_i }{ \partial x_j } ( \mathbf{x} )$

将它写成矩阵的形式就有：

$\mathbf{F}' ( \mathbf{x} ) = \mathbf{J} ( \mathbf{x} )^T \mathbf{f} ( \mathbf{x} )$

最后，F的二阶导（海塞矩阵）为：

$\begin{aligned} \mathbf{F}'' ( \mathbf{x} ) & = \frac{ \partial ^2 F }{ \partial x_j \partial x_k } ( \mathbf{x} ) = \sum_{i=1}^m ( \frac{ \partial f_i }{ \partial x_j } ( \mathbf{x} ) \frac{ \partial f_i }{ \partial x_k } ( \mathbf{x} ) + f_i ( \mathbf{x} ) \frac{ \partial ^2 f_i }{ \partial x_j \partial x_k } ( \mathbf{x} ) ) \\ & = \mathbf{J} ( \mathbf{x} )^T \mathbf{J} ( \mathbf{x} ) + \sum_{i=1}^m f_i ( x ) \mathbf{f}''_i ( \mathbf{x} ) \end{aligned}$

3.1 高斯-牛顿法

GN算法基于对误差函数f的线性近似：

$\mathbf{f} ( \mathbf{x} + \mathbf{h} ) \approx \mathbf{l} ( \mathbf{h} ) = \mathbf{f} ( \mathbf{x} ) + \mathbf{J} ( \mathbf{x} ) \mathbf{h}$

因此，损失函数F就变为

$\begin{aligned} F( \mathbf{x} + \mathbf{h} ) \approx L ( \mathbf{h} ) & = \frac{1}{2} \mathbf{l} ( \mathbf{h} )^T \mathbf{l} ( \mathbf{h} ) \\ & = \frac{1}{2} \mathbf{f}^T \mathbf{f} + \mathbf{h}^T \mathbf{J}^T \mathbf{f} + \frac{1}{2} \mathbf{h}^T \mathbf{J}^T \mathbf{J} \mathbf{h} \\ & = F( \mathbf{x} ) + \mathbf{h}^T \mathbf{J}^T \mathbf{f} + \frac{1}{2} \mathbf{h}^T \mathbf{J}^T \mathbf{J} \mathbf{h} \end{aligned}$

L 的梯度和海塞矩阵就是：

$\mathbf{L}' ( \mathbf{h} ) = \mathbf{J}^T \mathbf{f} + \mathbf{J}^T \mathbf{J} \mathbf{h} \quad \mathbf{L}'' ( \mathbf{h} ) = \mathbf{J}^T \mathbf{J}$

而且，如果梯度J是满秩的，那么L的海塞矩阵H就是正定的，说明L(h)拥有唯一的极小值。

通过令L&#39;(h) = 0就可以得到GN算法的解：

$(\mathbf{J}^T \mathbf{J} ) \mathbf{h}_{gn} = - \mathbf{J}^T \mathbf{f}$

可以证明， $\mathbf{h}_{gn}$ 确为损失函数F的下降方向：

$\mathbf{h}_{gn}^T \mathbf{F}' ( \mathbf{x} ) = \mathbf{h}_{gn}^T ( \mathbf{J}^T \mathbf{f} ) = -\mathbf{h}_{gn}^T ( \mathbf{J}^T \mathbf{J} ) \mathbf{h}_{gn} < 0$

3.1.1 收敛速度

上一小节提到过，牛顿法可以达到二次收敛。但GN法却不行。为此，我们可以对比这两种方法：

$\mathbf{F}''(x) \mathbf{h}_n = -\mathbf{F}'(x) \quad \text{and} \quad \mathbf{L}'' (\mathbf{h}) \mathbf{h}_{gn} = -\mathbf{L}'(0)$

已知这两个等式的右边是一样的，但是等式左边却是不一样的：

$\mathbf{F}'' ( \mathbf{x} ) = \mathbf{L}'' ( \mathbf{h} ) + \sum_{i=1}^m f_i ( x ) \mathbf{f}''_i ( \mathbf{x} )$

显然，如果 $\mathbf{f} ( \mathbf{x}^* ) = 0$ ，那么当x接近x*时， $\mathbf{L}'' (\mathbf{h}) \backsimeq \mathbf{F}''(x)$ 。此时，GN算法可以有二次收敛；但更多时候，GN算法只能达到线性收敛。

3.2 Levenberg-Marquardt 法

LM算法是阻尼法的一种，其步长可以这么算出：

$( \mathbf{J}^T \mathbf{J} + \mu \mathbf{I} ) \mathbf{h}_{lm} = -\mathbf{g} \quad \text{with} \ \mathbf{g} = - \mathbf{J}^T \mathbf{f} \ \text{and} \ \mu \ge 0$

其中，u是拉格朗日算子，它有以下作用：

1. 对于u > 0，可以保证系数矩阵是正定的。这就保证了 $\mathbf{h}_{lm}$ 必然是一个下降方向；
2. 如果u值很大，那么 $\mathbf{h}_{lm} \backsimeq - \frac{1}{\mu} \mathbf{g} = -\frac{1}{\mu} \mathbf{F}' ( \mathbf{x} )$ . 这是最速下降法下的一小步。这应用在当前迭代距离结果很远的时候；
3. 如果u值很小，那么 $\mathbf{h}_{lm} \simeq \mathbf{h}_{gn}$ ，说明当前的线性近似比较准确。

因此，阻尼因子同时影响了当前迭代的方向和步长。其更新方式在前面已经介绍过了，而初始值可以根据 $\mathbf{A}_0 = \mathbf{J}^T(\mathbf{x}_0) \mathbf{J}(\mathbf{x}_0)$ 中元素的大小给出：

$\mu_0 = \tau \cdot \max_i \{ a^{(0)}_{ii} \}$

对于  的选择，如果起始点x0是对结果的一个比较好的近似，  应该选取一个比较小的值，如 $10^{-6}$ ；否则，可以选择 $\tau = 10^{-3}$ 或者 $\tau = 1$ .

其更新策略同样是计算gain ratio：

$\rho = \frac{F(\mathbf{x}) - F(\mathbf{x} + \mathbf{h}_{lm})}{L(0) - L(\mathbf{h}_{lm})}$

其中，

$\begin{aligned} L(0) - L(\mathbf{h}_{lm})& = -\mathbf{h}_{lm}^T \mathbf{J}^T \mathbf{f} - \frac{1}{2}\mathbf{h}_{lm}^T \mathbf{J}^T \mathbf{J} \mathbf{h}_{lm} \\ & = -\frac{1}{2} \mathbf{h}_{lm}^T ( 2 \mathbf{g} + (\mathbf{J}^T \mathbf{J} + \mu \mathbf{I} - \mu \mathbf{I}) \mathbf{h}_{lm}) \\ & = \frac{1}{2} \mathbf{h}_{lm}^T ( \mu \mathbf{h}_{lm} - \mathbf{g} ) \end{aligned}$

注意到等式最后两项必然是正的，因此，就像前面说的，  的值取决于F的增量。

最后，LM算法的流程如下：

3.2.1 停止条件

首先，当算法到达一个最小值时，有 $\mathbf{F}' ( \mathbf{x}^* ) = \mathbf{g} ( \mathbf{x}^* ) = 0$ 。因此，第一种停止条件可以为

$|| \mathbf{g}||_\infty \le \epsilon_1$

这里， $\epsilon_1$ 是一个很小的正数。

此外，第二种停止条件为在x变化很小时停止：

$|| \mathbf{x}_{new} - \mathbf{x} || \le \epsilon_2 ( || \mathbf{x} || + \epsilon_2 )$

这说明了在x数值较大时，使用前者（相对量）来说为停止条件；在x数值较小时，使用后者（绝对量）来作为停止条件。

最后，还可以限制迭代次数：

$k \ge k_{max}$
3.3 狗腿法（Powell&#39;s dog-leg method)

我们已经知道，GN算法通过求解如下等式（normal equation）得到迭代方向和步长：

$\mathbf{J}^T \mathbf{J} \mathbf{h}_{gn} = - \mathbf{J}^T \mathbf{f}$

而最速下降法则用逆梯度方向作为迭代方向：

$\mathbf{h}_{sd} = - \mathbf{g} = - \mathbf{J}^T \mathbf{f}$

其步长可以通过以下方式得到：

$\mathbf{f} ( \mathbf{x} + \alpha \mathbf{h}_{sd} ) \simeq \mathbf{f} ( \mathbf{x} ) + \alpha \mathbf{J} \mathbf{h}_{sd}$

对损失函数则有

$F( \mathbf{x} + \alpha \mathbf{h}_{sd} ) \simeq \frac{1}{2} || \mathbf{f} ( \mathbf{x} ) + \alpha \mathbf{J} \mathbf{h}_{sd} || ^2 = F( \mathbf{x} ) + \alpha \mathbf{h}_{sd}^T \mathbf{J}^T \mathbf{J} \mathbf{f} + \frac{1}{2} \alpha^2 || \mathbf{J} \mathbf{h}_{sd} ||^2$

通过对上式进行求导，并令导数为0可得

$\alpha = -\frac{ \mathbf{h}_{sd}^T \mathbf{J}^T \mathbf{f}}{ || \mathbf{J} \mathbf{h}_{sd} || ^2 } = \frac{|| \mathbf{g}||^2 }{ || \mathbf{J} \mathbf{g} ||^2 }$

这样一来，在更新x的时候，我们就有两个可选的解了： $\mathbf{a} = \alpha \mathbf{h}_{sd}$ 和 $\mathbf{b} = \mathbf{h}_{gn}$ 。

而狗腿法是一种置信区域方法。其选取迭代方向和步长的方法如下图所示：

通过上图可以发现，还有一个参数 $\beta$ 没有确定。令 $c = \mathbf{a}^T ( \mathbf{b} - \mathbf{a} )$ ，和

$\psi(\beta) = || \mathbf{a} + \beta ( \mathbf{b} - \mathbf{a}) ||^2 - \Delta^2 = || \mathbf{b} - \mathbf{a} ||^2 \beta^2 + 2 c \beta + || \mathbf{a} || - \Delta^2$

则有

$\begin{cases} \beta = (-c + \sqrt{c^2 + ||\mathbf{b} - \mathbf{a}||^2 (\Delta^2 - ||\mathbf{a}||^2)}) / || \mathbf{b} - \mathbf{a} ||^2 & \quad c \le 0 \\ \beta =(\Delta^2 - ||\mathbf{a}||^2) / (c + \sqrt{c^2 + ||\mathbf{b} - \mathbf{a}||^2 (\Delta^2 - ||\mathbf{a}||^2)}) & \quad c > 0 \end{cases}$

最终的算法如下图所示

<hr/>参考文献：METHODS FOR NON-LINEAR LEAST SQUARES PROBLEMS

如果觉得有用，还请点个赞: )

KaaPexei · 发表于 2021-12-29 18:21

以李群为代表的微分流行方法，为众多高维问题提供了精致完备简洁的描述；非线性优化（高斯-牛顿，LM等）和状态观测器（贝叶斯，卡尔曼滤波等）提供了递归迭代求解最优值的途径。难点就在于，如何将欧式空间的优化方法用在流行（manifold）空间，对实际问题建立合理的manifold模型。欢迎讨论！正研究相关问题。

kyuskoj · 发表于 2021-12-29 18:27

解释的真清楚

mastertravels77 · 发表于 2021-12-29 18:33

最速下降法第一个式子的中间部分少了一个负号

acecase · 发表于 2021-12-29 18:36

非常棒[赞同]

ainatipen · 发表于 2021-12-29 18:36

谢谢博主，解释得很清楚！

zt3ff3n · 发表于 2021-12-29 18:39

为什么lm部分，L(0)-L(hlm)的最后两项一定是正的呀？

NoiseFloor · 发表于 2021-12-29 18:47

这是我们期望的误差减小量，公式代入就可以得到

		自动登录	找回密码
密码			立即注册

非线性优化（一）：优化方法

本帖子中包含更多资源