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

  • 最小二乘法
    • 快速回顾最小二乘法的推导
      • 建立误差平方
      • 将其最小化
    • 一种对最小二乘法理解的视角
  • 递推最小二乘法
    • 在线实时预测问题
    • 推导思路与详细过程
      • 将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. Windows脚本初探之PowerShell流程控制if
  2. TomCat运行struts1的编码问题
  3. 网站长尾词布置需要注意什么事项?
  4. [Head First设计模式]饺子馆(冬至)中的设计模式——工厂模式
  5. 理解一条语句:SELECT difference(sum(value)) FROM mq_enqueue WHERE channel =~ /ActiveMQ_TEST/ AND $tim...
  6. 【IDEA】自动导入无歧义的包
  7. 所有特征在不同分类之间、 train和test之间的列分布差异(图形绘制)
  8. python对比两张图片_用python实现对比两张图片的不同
  9. JS定时器使用,定时定点,固定时刻,循环执行
  10. bash: shasum: command not found
  11. Mac下VirtualBox虚拟机Win7与主机共享文件夹
  12. React基础篇(六)React中绑定事件的注意点
  13. 让px单位自动转换为rem的方法
  14. 交通规划软件功能分析
  15. IT笔试题收集,免费下载
  16. 【Ubuntu 安装】Ubuntu20.04和Win10双系统安装指南
  17. NIST SP 800-108密钥导出函数KDF研究
  18. 新手学Python之学习官网教程(一: Whetting Your Appetite)
  19. 数据结构-malloc申请动态空间-链表的创建
  20. 简单编程---哥德巴赫猜想

热门文章

  1. 2023最新SSM计算机毕业设计选题大全(附源码+LW)之java高校会议室管理系统w169g
  2. 一分钱不用花,一分钟搭建一个微官网
  3. 几种常见的传统汽车总线传输通信技术
  4. MapGIS67编辑线右键结束时,想继续接着编辑线时怎么操作?
  5. java mp3 信息_java读取MP3的信息
  6. access 2000 转换为97
  7. 智能ai文章伪原创工具-智能ai文章原创处理系统
  8. 职称计算机技巧集锦,PowerPoint2003使用技巧集锦(7)
  9. Taskade——Mac最全能的任务管理器
  10. Ecside_项目搭建