从最小二乘法到卡尔曼滤波

2021 年 12 月 22 日 图与推荐


© 作者 | 李承蒙
学校 | 桂林电子科技大学
研究方向 | 计算机视觉


本文旨在梳理总结学习到的一些知识。由于笔者水平有限,文中难免存在一些不严谨和错误之处,诚请各位批评指正。

最近看了一篇文章,文章从最小二乘法的角度推导了卡尔曼滤波的公式(链接在文末)。看完后震惊不已,很受启发,于是写了这篇文章。一是为了倒逼输出从而达到知新的效果,二是为了记录一下自己的理解以便日后自己翻阅。



最小二乘法


最小二乘法是解决曲线拟合问题最常用的方法,本质是进行参数估计。实现方式是先使用待定系数法设出拟合函数,然后使用拟合函数和观测数据构造出目标函数来评估拟合效果,并取拟合效果最好时的拟合函数。拟合函数由参数和事先选定的一组线性无关的函数(这里我们姑且叫它基函数)构成,基函数已事先确定,求解拟合函数的问题便转变为参数估计的问题。

是一组观测数据:


假定其噪声(随机变量) 分别为:

并且 ,噪声之间互不相关。

为拟合函数,其中 参数(也就是我们求解的重点) 为事先选定的一组线性无关的函数:


为目标函数:


其中除以 的目的是加权,噪声方差越大的观测数据其权值越小,保证整体结果不会被个别噪声较大的数据所影响。同时也满足了高斯-马尔可夫定理的要求:

高斯-马尔可夫定理

在线性回归模型中,如果误差满足零均值、同方差且互不相关,则回归系数的最佳线性无偏估计BLUE, Best Linear Unbiased Estimator)就是普通最小二乘法估计。

我们的目的便是使 取到最小值,将此时得到的 作为最佳的拟合函数。

举个例子:


如马同学的这张图为例,图中的 便是观测数据并且每个都对应 轴上的一个坐标。通过观察我们可以取拟合函数为 ,其中 便是我们要估计的参数。

构造目标函数:


分别对 求偏导,当满足下式时目标函数 取到最小值,可求出参数


从而得到拟合函数 。当然如果你对拟合出的函数不满意,可以再取其他的函数作为拟合函数并去估计其参数。



一个简单的状态估计问题



还是这张图,不过这次我们给坐标轴一些物理意义。横坐标代表时间,纵坐标代表位置,我们将其视为一维的匀速运动目标的位置-时间图。你可以想象成一个小车沿着直线匀速前进,你每隔一段时间观测它一次并记录位置和时间数据。

分我们定义状态量 ,其中 代表目标所处的位置, 代表目标当前的速度。

观测数据


自然还有我们刚刚得出来的拟合函数 

自然还有我们刚刚得出来的拟合函数

不过在被赋予了物理意义之后它现在长这样: ,其中 代表函数与纵轴的截距,也就是 时刻时的位置, 代表目标的运动速度。不过我们对目标 0 时刻的位置不感兴趣,我们对它当前的状态估计值感兴趣,于是有:


其中 代表对当前目标位置的估计值, 代表对当前目标速度的估计值。这样的拟合函数 才是我们想要的,它包含了目标在 时刻的位置与速度。换句话说,我们通过这五个观测数据得到了目标在 时刻的状态估计值

于是,如果我们能继续获得更多的观测数据(一直到 ),那么有:


从而有:


于是:



其中 时刻下的观测矩阵,它将目标 时刻下的状态 转化观测值

这是个很怪的拟合函数,非常反直觉。它的观测方式似乎是使用最新的状态估计值,去获得其他时间节点上的观测值;如此一来它既是在对状态进行观测,又是实现了不同时刻间状态的转移。但我们根本不关心拟合函数本身,我们只关心构成它的参数。它的参数包含了目标当前时刻的状态估计值 。如果我们不断地继续获得数据,我们也能不断地对拟合函数的参数进行估计,从而得出目标最新的状态估计值。



铺垫


但是我们发现这个方法 太 慢 了。随着迭代推进,每一次迭代都需要用到历史的所有数据来估计目标当前的状态。也就是说随着运行时间增长,积累的历史数据越多,计算出目标当前状态估计值所需要的时间也就越长。

我们先把最小二乘法拓展为矩阵形式来看一看。

将历史观测量合起来:


观测噪声:


同理,也把对应的观测矩阵合起来:


对于目标函数:




虽然变成了矩阵形式,但最小二乘法的思想没有变,通过类比一维的情况可以快速的理解矩阵形式。 

我们不妨来实实在在的求一下矩阵形式下拟合函数的参数。

在当前我们已经有 个观测数据的情况下,我们对当前目标的状态估计量求导,并令其为零:


于是就有:


上式中的 便是我们目前获得到的观测量数量,可见随着迭代的进行我们的运算量会越来越大。

我们肯定无法接受这一点!

所以我们还需要更进一步。



递推最小二乘法

我们的目的就是避免运算量随着时间而增长,所以必须想方设法将 改为递推形式。换句话说,只使用 时刻的各种数据来推算 时刻的状态,而不是像现在这样将所有的历史信息全部用上。但是忘记历史意味着背叛,所以我们选一个折中的办法——“只送大脑”,我们可以递推地求解目标每个时刻状态的协方差矩阵。这样可以将历史信息蕴含在协方差矩阵中,达到我们的目的。
对于:


有以下性质:





于是求解 的协方差矩阵  :



其中:


所以:


并且我们发现

然后我们把 拆开可以看到:


于是有状态协方差矩阵递推公式:


同理把 拆开可得:


并且我们知道 ,于是推理可得:


于是:


至此我们得到了以下公式:



的形式非常有趣,它以 为基准,并通过 的方式来估计 时刻观测量的误差( 表示对 时刻观测量的估计),最后乘以 并补偿到 上。我们一般称 这样的系数为观测增益。这种综合了观测数据与状态数据的计算方式已经有些接近卡尔曼滤波了,但是只有更新过程,几乎没有预测过程。比如在 式中的 式中的 完全可以替换成由他们自己经过预测后得出的 时刻的估计量,然后再进行更新,这也是卡尔曼滤波的思想之一。

问题的根源在于我们怪异而又反直觉的观测矩阵。它虽然叫观测矩阵,但它耦合了观测与状态转移两个功能。观测指的是将状态量转化为观测量的过程,状态转移指的是将某时刻的状态转化为另一个时刻的状态。与卡尔曼滤波相比,递推最小二乘法在整个递推过程中缺少了程噪声。所以我们需要对其进行解耦合,重新定义观测矩阵与状态转移矩阵。



从最小二乘法到卡尔曼滤波


让我们先用正常的方式来描述目标状态的递推过程:


状态转移矩阵,它可以将 时刻的状态转变为 时刻的。其中 为目标在 时刻的先验状态量,也就是说它是我们使用 时刻的状态预测得来的。 时刻的后验状态量,也就是经过上轮递推后得到的 时刻的最优状态估计值。

同时我们用 来表示过程噪声,其中 代表位置噪声, 代表速度噪声,且有:


注意,我们将 视为一个随机变量,且有:


同时设


并称 为过程噪声协方差矩阵。


现在我们再去求一次状态协方差矩阵:


其中:


但要注意的是此时( 时刻)计算的状态协方差矩阵为先验状态协方差矩阵 ,也就是我们根据 时刻的状态协方差矩阵推算得出的。

现在我们获得了预测过程的公式:



从而对应的更新过程公式应该改写为:



发现式 中只出现过 ,也就是观测矩阵的下标只有 。观测矩阵在下标为 (也就是最新的时刻)时只具有观测意义,不存在状态转移。因为在定义时观测矩阵 的意义是将最新时刻 目标的状态转变为 时刻的观测量(可以看二、中举的例子),当 时便只有观测意义。于是我们不需要重新定义观测矩阵,当前的式子符合我们的要求。



总结


最小二乘法求解最小目标函数的过程与卡尔曼滤波中求解协方差矩阵迹的最小值的过程非常相似,只不过最小二乘法是使用所有数据来进行优化,而卡尔曼滤波通过对迹求导来获得最优卡尔曼增益。我们可以通过构造特殊拟合函数的方式来让最小二乘法来递推运作,从而只传递状态噪声。然后定义状态转移过程,并加入过程噪声来描述状态转移误差。

相比使用贝叶斯公式与正态分布假设来推导卡尔曼滤波,最小二乘法推导的条件更为宽松。使用贝叶斯公式来推导时需要假设所有的噪声服从正态分布,而最小二乘法推导仅需要满足高斯-马尔可夫定理(噪声零均值、同方差、不相关)即可。当然推导卡尔曼滤波的方法并不只这两种,但了解的越多理解越深刻,对学习也更有帮助。

从递推最小二乘法推导出卡尔曼滤波的过程并不严谨,直接使用先验估计来替换掉递推最小二乘法中的一些项。如果以后有更加严谨的方式我会补上证明。


参考文献

[1] https://zhuanlan.zhihu.com/p/67250500

[2] https://zhuanlan.zhihu.com/p/339118204

[3] https://www.zhihu.com/question/37031188/answer/411760828



登录查看更多
1

相关内容

卡尔曼滤波是一种高效率的递归滤波器(自回归滤波器),它能够从一系列的不完全及包含噪声的测量中,估计动态系统的状态。
【经典书】时间序列分析与应用R语言,第二版,501页pdf
【经典书】时间序列分析:预测与控制(原书第5版),709页pdf
NeurIPS 2021 Spotlight | 针对有缺失坐标的聚类问题的核心集
专知会员服务
14+阅读 · 2021年11月27日
【NeurIPS 2021】设置多智能体策略梯度的方差
专知会员服务
20+阅读 · 2021年10月24日
【经典书】信息论原理,774页pdf
专知会员服务
254+阅读 · 2021年3月22日
【TAMU】最新《时间序列分析》课程笔记,527页pdf
专知会员服务
179+阅读 · 2020年9月12日
手写实现李航《统计学习方法》书中全部算法
专知会员服务
48+阅读 · 2020年8月2日
交替方向乘子法(ADMM)算法原理详解
PaperWeekly
3+阅读 · 2022年1月21日
【夯实基础】卡尔曼滤波
极市平台
1+阅读 · 2021年11月3日
经典重温:卡尔曼滤波器介绍与理论分析
极市平台
0+阅读 · 2021年10月25日
从泰勒展开来看梯度下降算法
深度学习每日摘要
13+阅读 · 2019年4月9日
面试时让你手推公式不在害怕 | 梯度下降
计算机视觉life
14+阅读 · 2019年3月27日
R语言时间序列分析
R语言中文社区
12+阅读 · 2018年11月19日
零基础概率论入门:最大似然估计
论智
12+阅读 · 2018年1月18日
[有意思的数学] 参数估计
机器学习和数学
15+阅读 · 2017年6月4日
回归预测&时间序列预测
GBASE数据工程部数据团队
43+阅读 · 2017年5月17日
国家自然科学基金
1+阅读 · 2014年12月31日
国家自然科学基金
1+阅读 · 2014年12月31日
国家自然科学基金
0+阅读 · 2014年12月31日
国家自然科学基金
0+阅读 · 2013年12月31日
国家自然科学基金
0+阅读 · 2013年12月31日
国家自然科学基金
2+阅读 · 2013年12月31日
国家自然科学基金
0+阅读 · 2012年12月31日
国家自然科学基金
1+阅读 · 2012年12月31日
国家自然科学基金
0+阅读 · 2011年12月31日
国家自然科学基金
0+阅读 · 2011年12月31日
Arxiv
0+阅读 · 2022年4月18日
Arxiv
15+阅读 · 2021年2月19日
Arxiv
31+阅读 · 2020年9月21日
Arxiv
17+阅读 · 2018年4月2日
VIP会员
相关VIP内容
【经典书】时间序列分析与应用R语言,第二版,501页pdf
【经典书】时间序列分析:预测与控制(原书第5版),709页pdf
NeurIPS 2021 Spotlight | 针对有缺失坐标的聚类问题的核心集
专知会员服务
14+阅读 · 2021年11月27日
【NeurIPS 2021】设置多智能体策略梯度的方差
专知会员服务
20+阅读 · 2021年10月24日
【经典书】信息论原理,774页pdf
专知会员服务
254+阅读 · 2021年3月22日
【TAMU】最新《时间序列分析》课程笔记,527页pdf
专知会员服务
179+阅读 · 2020年9月12日
手写实现李航《统计学习方法》书中全部算法
专知会员服务
48+阅读 · 2020年8月2日
相关资讯
交替方向乘子法(ADMM)算法原理详解
PaperWeekly
3+阅读 · 2022年1月21日
【夯实基础】卡尔曼滤波
极市平台
1+阅读 · 2021年11月3日
经典重温:卡尔曼滤波器介绍与理论分析
极市平台
0+阅读 · 2021年10月25日
从泰勒展开来看梯度下降算法
深度学习每日摘要
13+阅读 · 2019年4月9日
面试时让你手推公式不在害怕 | 梯度下降
计算机视觉life
14+阅读 · 2019年3月27日
R语言时间序列分析
R语言中文社区
12+阅读 · 2018年11月19日
零基础概率论入门:最大似然估计
论智
12+阅读 · 2018年1月18日
[有意思的数学] 参数估计
机器学习和数学
15+阅读 · 2017年6月4日
回归预测&时间序列预测
GBASE数据工程部数据团队
43+阅读 · 2017年5月17日
相关基金
国家自然科学基金
1+阅读 · 2014年12月31日
国家自然科学基金
1+阅读 · 2014年12月31日
国家自然科学基金
0+阅读 · 2014年12月31日
国家自然科学基金
0+阅读 · 2013年12月31日
国家自然科学基金
0+阅读 · 2013年12月31日
国家自然科学基金
2+阅读 · 2013年12月31日
国家自然科学基金
0+阅读 · 2012年12月31日
国家自然科学基金
1+阅读 · 2012年12月31日
国家自然科学基金
0+阅读 · 2011年12月31日
国家自然科学基金
0+阅读 · 2011年12月31日
Top
微信扫码咨询专知VIP会员