Skip to content

Latest commit

 

History

History
184 lines (137 loc) · 6.86 KB

LeastSquare.md

File metadata and controls

184 lines (137 loc) · 6.86 KB

最小二乘法


说明:网上关于最小二乘法的介绍也挺多了,这里只是做一些简单的总结,包括一些常用的求解方式方法,并不赘述其推导过程。

最小二乘法与最大似然估计

对于一批数据$(\mathbb x_i ,\mathbb y_i)$,通常会试图寻找一个模型拟合这批数据,也即 $$ \mathbb y_i = \mathbb f(\mathbb x_i) $$ 但是实际过程中数据是存在噪声的,所以模型也表示为: $$ \mathbb y_i = \mathbb f(\mathbb x_i) + \eta_i $$ 通常假设$\eta_i$是符合$N(0, \Sigma)$的高斯分布,也即 $$ \mathbb y_i - \mathbb f(\mathbb x_i) \sim N(0, \Sigma) $$ 为了求解使得概率最大的$\mathbb y_i - \mathbb f(\mathbb x_i)$,实际也就是使得 $$ (\mathbb y_i -\mathbb f(\mathbb x_i))^T \Sigma^{-1} (\mathbb y_i - \mathbb f(\mathbb x_i)) $$ 最小。

$J(\mathbb x_i, \mathbb y_i)$是代价函数,有: $$ J(\mathbb x_i, \mathbb y_i) =(\mathbb y_i - \mathbb f(\mathbb x_i))^T\Sigma^{-1}(\mathbb y_i - \mathbb f(\mathbb x_i)) $$ 如果$\Sigma$是个对角矩阵,也就表示,这里只关心不同变量/状态/特征的权重。如果是$\Sigma$是个一般的协方差矩阵,也就表示,不仅考虑变量/状态/特征的权重,还考虑两者之间的相关性。

特别地,当$\Sigma$是单位矩阵时, $$ J(\mathbb x_i, \mathbb y_i) = (\mathbb y_i - \mathbb f(\mathbb x_i))^T (\mathbb y_i - \mathbb f(\mathbb x_i)) $$ 如果$\mathbb y_i$是实值,则有 $$ J(\mathbb x_i, y_i) = (y_i - f(\mathbb x_i))^2 $$ 对$ J(\mathbb x_i, \mathbb y_i)$求导有: $$ \nabla J(\mathbb x_i, \mathbb y_i) = 2H^T\Sigma^{-1}(\mathbb y_i - \mathbb f(\mathbb x_i)) $$ 其中,$H$是$\mathbb y_i - \mathbb f(\mathbb x_i)$的一阶导,即雅克比矩阵。

最小二乘法

一个简单的最小二乘法定义如下: $$ \min_{\mathbb x} \frac{1}{2} \left |\mathbb f(\mathbb x) \right |_{2}^2 $$ 如果$\mathbb f$是数学形式简单的函数,可直接使用解析形式求解,即求一阶导数,并令一阶导数为0,求得极值,选取其中使得代价函数最小的一组解,作为最优解。

但是通常$\mathbb f$比较复杂,需要使用迭代的方式求解。

对于一个更一般的实值代价函数$F(\mathbb x)$,为了求解最优解,通常做泰勒展开处理: $$ F(\mathbb x + \Delta \mathbb x) \approx F(\mathbb x) + J^T(\mathbb x) \Delta \mathbb x +\frac{1}{2} (\Delta \mathbb x)^TH(\mathbb x) \Delta \mathbb x $$ 其中,$J(x)$是$F(\mathbb x)$的雅克比矩阵,$H(\mathbb x)$是海塞矩阵。

对于一般的向量值函数$\mathbb f(\mathbb x)$通常做一阶泰勒展开(二阶导矩阵是个高阶张量,一般不讨论) $$ \mathbb f(\mathbb x + \Delta \mathbb x) \approx \mathbb f(\mathbb x) + J(\mathbb x) \Delta \mathbb x $$ 其中,$J(\mathbb x)$是$\mathbb f(\mathbb x)$的雅克比矩阵。

考虑公式(10),直接求解$\mathbb x$比较复杂,我们考虑求解$\Delta \mathbb x$,然后让$\mathbb x$沿着梯度下降的方向走。

对$\Delta \mathbb x$求导有 $$ \nabla _{\Delta x}F(\mathbb x + \Delta \mathbb x) \approx J(\mathbb x) + \frac {1} {2} (H(\mathbb x) + H^T(\mathbb x)) \Delta x \= J(\mathbb x) + H(\mathbb x) \Delta \mathbb x $$ 令导数为0,则有: $$ H(\mathbb x) \Delta x = - J(\mathbb x) $$ 求解步骤如下:

  1. 给定初值$\mathbb x_0$

  2. while True do

    ​ 根据式(13)计算$\Delta x$

    ​ 若$\Delta x$足够小,即$(\Delta \mathbb x)^T\Delta \mathbb x = 0$,则停止

    ​ 否则,$\mathbb x = \mathbb x + \Delta \mathbb x$

这种方式,通常由于$\mathbb x$的维度比较大,难以计算$H(\mathbb x)$,所以通常不采用这种方式求解。对于一些只知道$F(\mathbb x)$的函数可以采用这种方式求解,激光点云数据配准的NDT方法实际就是采用这种方法求解的。

高斯牛顿法

对于公式(9),可以对$\mathbb f(\mathbb x)$展开,也即式(11),所以有: $$ F(\mathbb x + \Delta \mathbb x) = \frac{1}{2} \left | \mathbb f(\mathbb x)\right|_2^2 \ = \frac{1}{2} \mathbb f^T(\mathbb x) \mathbb f(\mathbb x) \ \approx \frac{1}{2} (\mathbb f(\mathbb x) + J(\mathbb x) \Delta \mathbb x)^T (\mathbb f(\mathbb x) + J(\mathbb x) \Delta \mathbb x) \ = \frac{1}{2} (\mathbb f^T(\mathbb x) \mathbb f(\mathbb x) + 2 \mathbb f^T(\mathbb x) J(\mathbb x)\Delta \mathbb x + \Delta \mathbb x^T J^T(\mathbb x) J(\mathbb x) \Delta \mathbb x) $$ 对$\Delta \mathbb x$求导有: $$ \nabla _{\Delta x} F(\mathbb x + \Delta \mathbb x) \approx J^T(\mathbb x) \mathbb f(\mathbb x) + J^T(\mathbb x)J(\mathbb x) \Delta \mathbb x $$ 令式(15)为0, $$ J^T(\mathbb x) J(\mathbb x) \Delta \mathbb x = - J^T(\mathbb x) \mathbb f(\mathbb x) $$ 式(16)也被称为增量方程,也记做 $$ H \Delta \mathbb x = g $$ 求解步骤如下:

  1. 给定初值$\mathbb x_0$

  2. while True do

    ​ 分别计算$\mathbb f(\mathbb x), J(\mathbb x)$

    ​ 根据式(16)计算$\Delta \mathbb x$

    ​ 若$\Delta \mathbb x$足够小,则停止,

    ​ 否则,$\mathbb x = \mathbb x + \Delta \mathbb x$

由于通常要求$H$是一个正定可逆的矩阵,但是这一点很难得到满足,并且通常为了保证近似的准确性,需要限制$\Delta \mathbb x$的取值。所以通常会给$\Delta \mathbb x$乘上一个系数$\alpha$,即$\mathbb f(\mathbb x + \alpha \Delta \mathbb x)$

列文伯格-马夸尔特方法

求解$\Delta \mathbb x$时限定范围, $$ \mathbb f(\mathbb x + \Delta \mathbb x) - \mathbb f(\mathbb x) = \rho J(\mathbb x) \Delta \mathbb x $$ 其中$\rho$是比例系数,$\rho$接近于1,则近似效果好,如果$\rho$比较小,则说明需要缩小$\Delta \mathbb x$的范围,如果$\rho$比较大,则说明需要扩大$\Delta \mathbb x$的取值范围。

求解算法如下:

  1. 给定初始$\mathbb x_0$

  2. while True do

    ​ 求解优化公式 $$ \min_{\Delta \mathbb x} \frac {1}{2} \left | \mathbb f(\mathbb x_k) + J(\mathbb x_k) \Delta \mathbb x_k \right |_2^2, s.t. \left | D \Delta \mathbb x_k \right |_2^2 \leq \mu $$ ​ 计算$\rho$

    ​ 若$\rho > \frac {3} {4}$,则$\mu = 2 \mu$

    ​ 若$\rho < \frac{1}{4}$,则$\mu=0.5\mu$

    ​ 如果$\rho$大于阈值,则令$\mathbb x_{k+1} = \mathbb x_k + \Delta \mathbb x$

    ​ 判断是否收敛,收敛则跳出。

公式(19)中,方阵D将$\Delta \mathbb x_k$限制在了一个半径为$\mu$的球中,通过拉格朗日算子将式(19)转化为无约束的最优化问题: $$ \min_{\Delta \mathbb x} \frac {1} {2} \left | \mathbb f(\mathbb x_k) + J(\mathbb x_k) \Delta \mathbb x_k \right |_2^2 + \frac {\lambda} {2} \left | D \Delta \mathbb x \right |_2^2 $$ 求导得: $$ (J^T(\mathbb x)J(\mathbb x) + \lambda D^TD) \Delta \mathbb x = - J^T(\mathbb x) \mathbb f(\mathbb x) $$

其他求解算法

求解库或者求解包

ceres

g2o

Reference

  1. 《视觉SLAM十四讲》,第6章