递推最小二乘法的推导和理解

  • 最小二乘法
    • 快速回顾最小二乘法的推导
      • 建立误差平方
      • 将其最小化
    • 一种对最小二乘法理解的视角
  • 递推最小二乘法
    • 在线实时预测问题
    • 推导思路与详细过程
      • 将k时刻的表达式写成k-1时刻表达式加某一个量
      • 写出k-1时刻满足的最小二乘表达式
      • 将前两步的公式带入第k时刻的最小二乘表达式中
    • 公式的简单理解
      • 角度一:回归在线预测问题
      • 角度二:梯度下降视角
      • 角度三:状态方程视角下的(XkTXk)−1(X_k^{T}X_k)^{-1}(XkT​Xk​)−1:
  • 数据量太大: 矩阵求逆公式

本文的框架如下:

  • 首先回忆一些最小二乘法的概念,如果很熟悉可以直接跳到递推最小二乘法,评判标准就是可以理解(XkTXk)−1XkTYk(X_k^{T}X_k)^{-1}X_k^{T}Y_k(XkT​Xk​)−1XkT​Yk​这个公式的推导。
  • 之后介绍在线实时预测问题,引出递推最小二乘法并进行一些简单的理解。

最小二乘法

快速回顾最小二乘法的推导

最小二乘法解决的是给定一组输入数据和输出数据(xi,yi)(x_i,y_i)(xi​,yi​),对其进行参数估计的问题。达到的效果是估计出的参数可以使得这个方程很好拟合当前的数据。其实就是拟合。这里面最重要的思想就是误差平方最小。根据这个思想,推导的思路总体可以概括为:

  1. 根据误差平方生成损失函数
  2. 要最小化损失函数

建立误差平方

yiy_iyi​是现实中真正测量的值,而我们希望通过f(xi)f(x^i)f(xi)求出的yyy能够很好的和yiy_iyi​吻合(这里xix^ixi可以表示多个值)。现在已知有一组数据{(x1,y1),(x2,y2),...,(xn,yn)(x^1,y_1),(x^2,y_2),...,(x^n,y_n)(x1,y1​),(x2,y2​),...,(xn,yn​)},假设对于每一组点都希望满足:
yi=f(xi)=(x1i⋯xni)(θ1⋮θn)y_i=f(x^i)= \left( \begin{array}{ccc} x_{1}^i & \cdots & x_{n}^i \end{array} \right)\left( \begin{array}{c} θ_{1} \\ \vdots \\ θ_{n} \end{array} \right) yi​=f(xi)=(x1i​​⋯​xni​​)⎝⎜⎛​θ1​⋮θn​​⎠⎟⎞​
但问题是,不可能所有的f(xi)f(x_i)f(xi​)都和测量值相等。那我们自然想到,如果计算每个测量点和理论计算值的误差,让他们总和积累量最小,就可以近似认为满足需求,但是这个误差有正有负,所以我们对其进行平方处理。这里直接写成矩阵的形式:
J(θ)=12(Xθ−Y)T(Xθ−Y)J(θ) =\frac{1}{2}(Xθ-Y)^T(Xθ-Y) J(θ)=21​(Xθ−Y)T(Xθ−Y)
注意这里X,YX,YX,Y都是测量出来的数据。未知量只有θθθ。

将其最小化

让损失对于θθθ最小,自然就是对其进行求导等于零:
∂∂θJ(θ)=XT(Xθ−Y)=0\frac{\partial}{\partial θ}J(θ)=X^T(Xθ-Y)=0 ∂θ∂​J(θ)=XT(Xθ−Y)=0
整理得出:
θ=(XTX)−1XTYθ=(X^TX)^{-1}X^TY θ=(XTX)−1XTY

一种对最小二乘法理解的视角

对于最小二乘问题,其实是一个求解方程组的问题:
Xθ=YXθ=Y Xθ=Y
最理想的情况是XXX满秩,那么这样θθθ就可以直接解出,为:
θ=X−1Yθ=X^{-1}Y θ=X−1Y
但是一般情况下,XXX是一个长条的矩阵。这在数学上有个处理,等式左右同时乘XTX^TXT(不做解释,MIT线性代数中老爷子讲的很好,大家如果想深入理解,可以去看,其实是最小二乘的思想,在Chapter4),这就变成了最小二乘法:
θ=(XTX)−1XTYθ=(X^TX)^{-1}X^TY θ=(XTX)−1XTY
但是其本身想求解的还是:
θ=X−1Yθ=X^{-1}Y θ=X−1Y
所以之后对递推最小二乘法的结果分析中,我们也可以进行这种简化的理解,即就是要方程组的解
在进行了简单的回顾之后,下面我们引出递推最小二乘法。

递推最小二乘法

在线实时预测问题

现在我们改变一下需求,假如上述的数据不是一次性给出的,而是隔一段时间给一个数据,且需要根据之前的数据和现在多加进来的一个数据重新进行最小二乘的预测。或者换句话说,如果数据是实时在线给出的,我们需要怎么进行求解?

一种简单的想法就是每一时刻都进行一次最小二乘法的计算,这个当然可以,但是这是相当消耗内存和时间的。还有一种想法是采用迭代:既然我们在上一时刻已经计算出过一组参数,那么下一时刻能否只用上一时刻计算出来的参数,加上这一时刻得到的新的数据,计算出这一时刻最小二乘的结果。这样不仅不用重复计算,而且也不需要记忆数据,大大降低了数据的存储量。主要的思想其实就是找到一种迭代格式,让其满足:
θk+1=θk+ε(θk,xk+1)θ_{k+1} = θ_{k}+ ε(θ_{k},x_{k+1}) θk+1​=θk​+ε(θk​,xk+1​)

推导思路与详细过程

首先需要知道,一组数据是一个列向量,在第kkk时刻,一共有kkk组数据,这个时候这些列向量共同组成了XkX_kXk​。这里的下标可以理解成有kkk组数据。

将k时刻的表达式写成k-1时刻表达式加某一个量

对于最小二乘法的解中的XkTXkX_k^TX_kXkT​Xk​:
XkTXk=[x1...xk][x1T⋮xkT]=[x1...xk−−1][x1T⋮xk−1T]+xkxkT=Xk−1TXk−1+xkxkTX_k^TX_k=\begin{bmatrix} x_1 & ... & x_k \end{bmatrix}\begin{bmatrix} x_1^T\\\\ \vdots \\\\ x_k^T \end{bmatrix}=\begin{bmatrix} x_1 & ... & x_{k--1} \end{bmatrix}\begin{bmatrix} x_1^T\\\\ \vdots \\\\ x_{k-1}^T \end{bmatrix}+x_kx_k^T=X_{k-1}^TX_{k-1}+x_kx_k^T XkT​Xk​=[x1​​...​xk​​]⎣⎢⎢⎢⎢⎢⎡​x1T​⋮xkT​​⎦⎥⎥⎥⎥⎥⎤​=[x1​​...​xk−−1​​]⎣⎢⎢⎢⎢⎢⎡​x1T​⋮xk−1T​​⎦⎥⎥⎥⎥⎥⎤​+xk​xkT​=Xk−1T​Xk−1​+xk​xkT​
同理,对于XkTYkX_k^TY_kXkT​Yk​做相同的处理,得出:
XkTYk=Xk−1TYk−1+xkykX_k^TY_k=X_{k-1}^TY_{k-1}+x_ky_k XkT​Yk​=Xk−1T​Yk−1​+xk​yk​

写出k-1时刻满足的最小二乘表达式

在k−1k-1k−1时刻,也应该满足最小二乘的表达式:
θk−1=(Xk−1TXk−1)−1Xk−1TYk−1θ_{k-1}=(X_{k-1}^{T}X_{k-1})^{-1}X_{k-1}^{T}Y_{k-1} θk−1​=(Xk−1T​Xk−1​)−1Xk−1T​Yk−1​
做一步变换:
(Xk−1TXk−1)θk−1=Xk−1TYk−1(X_{k-1}^{T}X_{k-1})θ_{k-1}=X_{k-1}^{T}Y_{k-1} (Xk−1T​Xk−1​)θk−1​=Xk−1T​Yk−1​

将前两步的公式带入第k时刻的最小二乘表达式中

θk=(XkTXk)−1XkTYk=(XkTXk)−1(Xk−1TYk−1+xkyk)=(XkTXk)−1((Xk−1TXk−1)θk−1+xkyk)=(XkTXk)−1((XkTXk−xkxkT)θk−1+xkyk)=θk−1−(XkTXk)−1xkxkTθk−1+(XkTXk)−1xkyk=θk−1+(XkTXk)−1xk(yk−xkTθk−1)\begin{aligned} θ_{k}&=(X_k^{T}X_k)^{-1}X_k^{T}Y_k\\[2mm] &=(X_k^{T}X_k)^{-1}(X_{k-1}^TY_{k-1}+x_ky_k)\\[2mm] &=(X_k^{T}X_k)^{-1}((X_{k-1}^{T}X_{k-1})θ_{k-1}+x_ky_k)\\[2mm] &=(X_k^{T}X_k)^{-1}((X_{k}^{T}X_{k}-x_kx_k^T)θ_{k-1}+x_ky_k)\\[2mm] &=θ_{k-1}-(X_k^{T}X_k)^{-1}x_kx_k^Tθ_{k-1}+(X_k^{T}X_k)^{-1}x_ky_k\\[2mm] &=θ_{k-1}+(X_k^{T}X_k)^{-1}x_k(y_k-x_k^Tθ_{k-1}) \end{aligned} θk​​=(XkT​Xk​)−1XkT​Yk​=(XkT​Xk​)−1(Xk−1T​Yk−1​+xk​yk​)=(XkT​Xk​)−1((Xk−1T​Xk−1​)θk−1​+xk​yk​)=(XkT​Xk​)−1((XkT​Xk​−xk​xkT​)θk−1​+xk​yk​)=θk−1​−(XkT​Xk​)−1xk​xkT​θk−1​+(XkT​Xk​)−1xk​yk​=θk−1​+(XkT​Xk​)−1xk​(yk​−xkT​θk−1​)​
至此其实递推最小二乘法的算法已经推到结束了。不过这里其实还有值也可以通过迭代来,就是(XkTXk)−1(X_k^{T}X_k)^{-1}(XkT​Xk​)−1这个量。这个我们之前已经推到过了,所以个人感觉最好是写成:
θk=θk−1+(Xk−1TXk−1+xkxkT)−1xk(yk−xkTθk−1)θ_{k}=θ_{k-1}+(X_{k-1}^TX_{k-1}+x_kx_k^T)^{-1}x_k(y_k-x_k^Tθ_{k-1}) θk​=θk−1​+(Xk−1T​Xk−1​+xk​xkT​)−1xk​(yk​−xkT​θk−1​)
因为这样,就可以只用上一时刻已经求出来的参数和当前得到的数据计算当前时刻的参数了。

公式的简单理解

我们主要看一下多出来的这一项(XkTXk)−1xk(yk−xkTθk−1)(X_k^{T}X_k)^{-1}x_k(y_k-x_k^Tθ_{k-1})(XkT​Xk​)−1xk​(yk​−xkT​θk−1​)是什么:
对于yk−xkTθk−1y_k-x_k^Tθ_{k-1}yk​−xkT​θk−1​这一项,其实可以理解成,使用上一时刻的参数和当前获得的数据计算出来的预测值和当前的测量值之间的误差。即上一时刻参数对当前时刻的误差值。我们写成εεε。此时公式变成:
θk=θk−1+(XkTXk)−1xkεθ_{k}=θ_{k-1}+(X_k^{T}X_k)^{-1}x_kε θk​=θk−1​+(XkT​Xk​)−1xk​ε

角度一:回归在线预测问题

我们将(XkTXk)−1xk(X_k^{T}X_k)^{-1}x_k(XkT​Xk​)−1xk​看成加上当前时刻的信息之后产生的价值,这样可以将其简写成KkK_kKk​。这样就变成了:
θk=θk−1+Kkεθ_{k}=θ_{k-1}+K_kε θk​=θk−1​+Kk​ε
这个公式其实就回到了当时我们提出在线预测问题那里,我们希望使用上一时刻已经计算出来的参数加上某些值就能算出当前的参数,这里加的值就是使用上一时刻参数造成的误差乘上比例,这个比例可以理解为加上当前数据产生的价值

角度二:梯度下降视角

之前在推导最小二乘法的时候我们计算过梯度为:
∂∂θJ(θ)=XT(Xθ−Y)=0\frac{\partial}{\partial θ}J(θ)=X^T(Xθ-Y)=0 ∂θ∂​J(θ)=XT(Xθ−Y)=0
如果我们希望更新θθθ,按照传统的梯度下降会怎么做:
θk=θk−1+αXT(Xθ−Y)θ_{k}=θ_{k-1}+αX^T(Xθ-Y) θk​=θk−1​+αXT(Xθ−Y)
但是这里需要进行一个修改,即将XT(Xθ−Y)X^T(Xθ-Y)XT(Xθ−Y)按照乘法对应关系,变成xk(yk−xkTθk−1)x_k(y_k-x_k^Tθ_{k-1})xk​(yk​−xkT​θk−1​)。其实就是将XT(Xθ−Y)X^T(Xθ-Y)XT(Xθ−Y)提取出了一部分,只不过这里的θθθ不是每一个kkk时刻都一样的(当然我们期望是一样的)。然后步长设置为(XkTXk)−1(X_k^{T}X_k)^{-1}(XkT​Xk​)−1,这个时候可能会有疑问,原来梯度下降,梯度是一个列向量,步长是一个标量,为什么这里步长变成了一个矩阵。写到这里我也思考了半天,于是有了第三个视角

角度三:状态方程视角下的(XkTXk)−1(X_k^{T}X_k)^{-1}(XkT​Xk​)−1:

状态方程在控制领域里面很常见,比如一个简单的基本运动学方程为:
x(k+1)=x(k)+Bku(k)x(k+1)=x(k)+B_ku(k) x(k+1)=x(k)+Bk​u(k)
这是一个动态的变化过程,根据运动学方程,我们就可以根据当前的状态和给出的控制,知道下一时刻的状态。我们将状态变量类比乘递推最小二乘中的参数,控制类比成梯度,那么相对应的矩阵BkB_kBk​对应的就是(XkTXk)−1(X_k^{T}X_k)^{-1}(XkT​Xk​)−1。也就是说,从状态方程的角度来讲,(XkTXk)−1(X_k^{T}X_k)^{-1}(XkT​Xk​)−1起到的作用是如何将所谓的梯度,或者控制,作用到参数上面的。这里因为XkX_kXk​一直在变动,所以通俗的解释为:控制的变动是由数据量的变动引起的

数据量太大: 矩阵求逆公式

这里主要是对(Xk−1TXk−1+xkxkT)−1(X_{k-1}^TX_{k-1}+x_kx_k^T)^{-1}(Xk−1T​Xk−1​+xk​xkT​)−1进行一点补充,其实如果求这个东西的逆,如果数据量不大,正常直接求就好。但是如果处理的数据量特别大的话,求逆会变成一件很耗时的事情,所以这里有一个数学公式可以进行处理:
[A+BCD]−1=A−1−A−1B[C−1+DA−1B]−1DA−1[A+BCD]^{-1}=A^{-1}-A^{-1}B[C^{-1}+DA^{-1}B]^{-1}DA^{-1} [A+BCD]−1=A−1−A−1B[C−1+DA−1B]−1DA−1

对应这个公式,可以把(Xk−1TXk−1+xkxkT)−1(X_{k-1}^TX_{k-1}+x_kx_k^T)^{-1}(Xk−1T​Xk−1​+xk​xkT​)−1写成:
(Xk−1TXk−1+xkxkT)−1=Xk−1TXk−1−Xk−1TXk−1xkxkTXk−1TXk−11+xkTXk−1TXk−1xk(X_{k-1}^TX_{k-1}+x_kx_k^T)^{-1}=X_{k-1}^TX_{k-1}-\frac{X_{k-1}^TX_{k-1}x_kx_k^TX_{k-1}^TX_{k-1}}{1+x_k^TX_{k-1}^TX_{k-1}x_k} (Xk−1T​Xk−1​+xk​xkT​)−1=Xk−1T​Xk−1​−1+xkT​Xk−1T​Xk−1​xk​Xk−1T​Xk−1​xk​xkT​Xk−1T​Xk−1​​

递推最小二乘法的推导和理解相关推荐

  1. 递推最小二乘法RLS公式详细推导

    递推最小二乘法RLS公式详细推导 整理递推最小二乘法推导过程自我整理. 递推最小二乘估计(RLS)作为一种估计方式是在最小二乘法(LS)的基础上发展来的. 最小二乘法可以解决的问题是不需要知道先验的概 ...

  2. 递推最小二乘法RLS:推导

    递推最小二乘法RLS:推导 I use RLS as a method for parameter identification. Reference: https://blog.csdn.net/w ...

  3. 递推最小二乘法RLS的轮胎侧偏刚度估计(原书缺失代码已补全)

    目录 1 参数辨识 1.1 最小二乘法 1.2 递推最小二乘法 RLS 1.3 具有遗忘因子 λ 的递推最小二乘法 2 轮胎线性侧偏刚度估计 2.1 RLS 算法分析 2.2 联合仿真平台的设计 ca ...

  4. VINS-Mono之后端非线性优化 (目标函数中视觉残差和IMU残差,及其对状态量的雅克比矩阵、协方差递推方程的推导)

    文章目录 1. 前言 2. 非线性最小二乘 2.1 Guass-Newton 和 Levenberg-Marquardt 2.2 鲁棒核函数下状态量增量方程的构建 3. 局部Bundle Adjust ...

  5. 基于遗忘因子递推最小二乘法辨识一阶RC等效电路模型

    %% 基于一阶RC等效电路模型实现不同倍率下电模型参数辨识 clear clc%% 载入实验数据 % 导入hppc实验数据 load('hppc_pulse_25deg') temp = hppc_p ...

  6. 【递推法】错排问题的递推式和推导过程

    [递推法]错排问题的递推式和推导过程   前言:这篇博客是帮助没有见过错排的新人更好的理解错排问题的递推式和推导过程,各位大佬可自行跳过 题目链接:洛谷P1595信封 一.错排问题的定义:   很多人 ...

  7. 函数的递推matlab,关于递推最小二乘法辨识参数的matlab编程(含注释)

    最近在做关于过热气温的动态建模问题.有现场运行的历史数据,要找出导前区和惰性区的传递函数. 对这类算法了解不多,程序读起来比较吃力,所以就转来一篇完整的辨识程序,在原有基础上进行了简化,并稍加注解一下 ...

  8. 一步完成最小二乘法、递推最小二乘法、增广最小二乘法、广义最小二乘法、辅助变量法、二步法辨识(传递函数)

    用一步完成最小二乘法.递推最小二乘法.增广最小二乘法.广义最小二乘法.辅助变量法.二步法辨识如下模型的参数: 噪声的成形滤波器 采样时间0.01 要求: 1.用matlab 写出程序代码: 2.画出实 ...

  9. 递推最小二乘法(Recursive least square, RLS)详细推导

    假设有数据(X,Y)(X,Y)(X,Y),其中X∈Rm×dX \in {\mathbb{R}^{m \times d}}X∈Rm×d,Y∈Rm×1Y \in {\mathbb{R}^{m \times ...

最新文章

  1. 关于librtmp接收数据(接收网络电视的数据流)
  2. python如何编写爬虫_如何实现一个Python爬虫框架
  3. 39_上下采样、MaxPool2d、AvgPool2d、ReLU案例、二维最大池化层和平均池化层、填充和步幅、多通道
  4. ELK:kibana使用的lucene查询语法
  5. modelsim和matlab联合仿真,Modelsim与Matlab联合仿真
  6. cocosstdio之字体之文本和FNT字体
  7. 计算机用公式找出第一名,用公式查找Excel工作表中重复数据
  8. vue项目 构建 打包 发布 三部曲
  9. 以逗号分隔的正则表达式_再见,正则表达式
  10. Kafka从上手到实践 - Kafka集群:Kafka Listeners | 凌云时刻
  11. UGNX1953~1980版本怎么测量重量
  12. 基于OpenLayers的地图应用中图标汉化
  13. MATLAB语言中int函数
  14. GameofMir引擎架设传奇服务器【3:在服务器上架设引擎】
  15. 数字电位器X9312
  16. MapReduce最佳成绩统计,男生女生比比看
  17. SSD 1306显示屏 adafruit SSD 1306
  18. 分享个被骗10元钱的经历
  19. Traingview MACD自定义指标颜色修改
  20. warning: #546-D: transfer of control bypasses initialization

热门文章

  1. 解决 clean-webpack-plugin www has been removed 问题
  2. excel高效之拆分单元格数据、导入ppt实现修改同步
  3. 公务员考试信息管理系统设计与实现-计算机毕业设计源码+LW文档
  4. 如何用射频接收机测量噪声系数?
  5. C语言|一个简单的文章让你轻松理解猜字小游戏的原理
  6. 云时通 X 九毛九 | SRM系统助力连锁餐饮行业高效管理供应链
  7. Gesture Recognition Dataset: Jester 数据集解压
  8. 告别繁琐的签到,使用 dailycheckin 就可以解决
  9. 知云文献翻译打不开_有了这个英文文献翻译助手,SCI论文阅读不用再复制粘贴...
  10. 数据结构 时间复杂度 空间复杂度 一看就懂版本