递推最小二乘法(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}}Y∈Rm×1，mmm为样本数，ddd为特征数，考虑最小二乘解
θ0=(XTX)−1XTY=Σ0−1XTY(1)\begin{aligned}{\theta_0} = {\left( {{X^{\rm{T}}}X} \right)^{ - 1}}{X^{\rm{T}}}Y = {\Sigma_0}^{-1}{X^{\rm{T}}}Y \tag{1}\end{aligned}θ0=(XTX)−1XTY=Σ0−1XTY(1)
⇒Σ0θ0=XTY(2)\Rightarrow {\Sigma_0}{\theta_0} = {X^{\rm{T}}}Y \tag{2}⇒Σ0θ0=XTY(2)
当新数据(X1,Y1)\left( {{X_1},{Y_1}} \right)(X1,Y1)到来时，更新模型，得到新的回归系数
θ1=([XX1]T[XX1])−1[XX1]T[YY1]=Σ1−1[XX1]T[YY1](3)\begin{aligned} {\theta_1} &= {\left( {{{\left[ {\begin{array}{cc} X\\ {{X_1}} \end{array}} \right]}^{\rm{T}}}\left[ {\begin{array}{cc} X\\ {{X_1}} \end{array}} \right]} \right)^{ - 1}}{\left[ {\begin{array}{cc} X\\ {{X_1}} \end{array}} \right]^{\rm{T}}}\left[ {\begin{array}{cc} Y\\ {{Y_1}} \end{array}} \right] \\ &= {\Sigma _1}^{ - 1}{\left[ {\begin{array}{cc} X\\ {{X_1}} \end{array}} \right]^{\rm{T}}}\left[ {\begin{array}{cc} Y\\ {{Y_1}} \end{array}} \right]\tag{3} \end{aligned}θ1=([XX1]T[XX1])−1[XX1]T[YY1]=Σ1−1[XX1]T[YY1](3)
其中
Σ1=[XX1]T[XX1]=XTX+X1TX1=Σ0+X1TX1(4)\begin{aligned} {\Sigma _1} &= {\left[ {\begin{array}{cc} X\\ {{X_1}} \end{array}} \right]^{\rm{T}}}\left[ {\begin{array}{cc} X\\ {{X_1}} \end{array}} \right] \\ &= {X^{\rm{T}}}X + X_1^{\rm{T}}{X_1} \\ &= {\Sigma _0} + X_1^{\rm{T}}{X_1} \end{aligned} \tag{4}Σ1=[XX1]T[XX1]=XTX+X1TX1=Σ0+X1TX1(4)
⇒Σ0=Σ1−X1TX1(5)\begin{aligned} \Rightarrow {\Sigma _0} = {\Sigma _1} - X_1^{\rm{T}}{X_1} \end{aligned} \tag{5}⇒Σ0=Σ1−X1TX1(5)
根据公式(4)的结果，通过归纳可得
Σk=Σk−1+XkTXk(6)\begin{aligned} {\Sigma _k} = {\Sigma _{k - 1}} + X_k^{\rm{T}}{X_k} \end{aligned} \tag{6}Σk=Σk−1+XkTXk(6)
[XX1]T[YY1]=XTY+X1TY1=Σ0θ0+X1TY1//公式(2)结果替换得到=(Σ1−X1TX1)θ0+X1TY1//公式(5)结果替换得到=Σ1θ0+X1T(Y1−X1θ0)(7)\begin{aligned} {\left[ {\begin{array}{cc} X\\ {{X_1}} \end{array}} \right]^{\rm{T}}}\left[ {\begin{array}{cc} Y\\ {{Y_1}} \end{array}} \right] &= {X^{\rm{T}}}Y + X_1^{\rm{T}}{Y_1}\\ &= {\Sigma _0}{\theta_0} + X_1^{\rm{T}}{Y_1} \quad //公式(2)结果替换得到\\ &= \left( {{\Sigma _1} - X_1^{\rm{T}}{X_1}} \right){\theta_0} + X_1^{\rm{T}}{Y_1} \quad //公式(5)结果替换得到\\ &= {\Sigma _1}{\theta_0} + X_1^{\rm{T}}\left( {{Y_1} - {X_1}{\theta_0}} \right) \end{aligned} \tag{7}[XX1]T[YY1]=XTY+X1TY1=Σ0θ0+X1TY1//公式(2)结果替换得到=(Σ1−X1TX1)θ0+X1TY1//公式(5)结果替换得到=Σ1θ0+X1T(Y1−X1θ0)(7)
将公式(7)回带到公式(3)：
θ1=Σ1−1(Σ1θ0+X1T(Y1−X1θ0))=θ0+Σ1−1X1T(Y1−X1θ0)(8)\begin{aligned} {\theta_1} &= {\Sigma _1}^{ - 1}\left( {{\Sigma _1}{\theta_0} + X_1^{\rm{T}}\left( {{Y_1} - {X_1}{\theta_0}} \right)} \right) \\ &= {\theta_0} + {\Sigma _1}^{ - 1}X_1^{\rm{T}}\left( {{Y_1} - {X_1}{\theta_0}} \right) \end{aligned} \tag{8}θ1=Σ1−1(Σ1θ0+X1T(Y1−X1θ0))=θ0+Σ1−1X1T(Y1−X1θ0)(8)
根据公式(8)的结果，通过归纳可得
θk=θk−1+Σk−1XkT(Yk−Xkθk−1)(9)\begin{aligned} {\theta_k} = {\theta_{k - 1}} + {\Sigma _k}^{ - 1}X_k^{\rm{T}}\left( {{Y_k} - {X_k}{\theta_{k - 1}}} \right) \end{aligned} \tag{9}θk=θk−1+Σk−1XkT(Yk−Xkθk−1)(9)

到这里，已经能够实现对最小二乘的递推，其过程可概括如下，我们称为算法1:

根据公式(5)更新Σk=Σk−1+XkTXk{\Sigma _k} = {\Sigma _{k - 1}} + X_k^{\rm{T}}{X_k}Σk=Σk−1+XkTXk；
根据公式(9)更新θk=θk−1+Σk−1XkT(Yk−Xkθk−1){\theta_k} = {\theta_{k - 1}} + {\Sigma _k}^{ - 1}X_k^{\rm{T}}\left( {{Y_k} - {X_k}{\theta_{k - 1}}} \right)θk=θk−1+Σk−1XkT(Yk−Xkθk−1)。

但以上过程存在两个问题：

对矩阵Σk\Sigma_kΣk的求逆计算复杂度比较高，我们能否在递推过程中避免对Σk\Sigma_kΣk的求逆计算，而直接更新它的逆矩阵；
矩阵Σk\Sigma_kΣk中的元素会随着数据量的增加不断增大，可能会发生数值溢出的问题。

针对以上问题，我们要对公式进一步改造，根据Sherman-Morrison-Woodbury公式：
(A+UVT)−1=A−1−A−1U(I+VTA−1U)−1VTA−1{\left( {A + U{V^{\rm{T}}}} \right)^{ - 1}} = {A^{ - 1}} - {A^{ - 1}}U{\left( {I + {V^{\rm{T}}}{A^{ - 1}}U} \right)^{ - 1}}{V^{\rm{T}}}{A^{ - 1}}(A+UVT)−1=A−1−A−1U(I+VTA−1U)−1VTA−1
公式(6)的逆可写成如下形式
Σk−1=(Σk−1+XkTXk)−1=Σk−1−1−Σk−1−1XkT(I+XkΣk−1−1XkT)−1XkΣk−1−1(10)\begin{aligned} {\Sigma _k}^{ - 1} &= {\left( {{\Sigma _{k - 1}} + X_k^{\rm{T}}{X_k}} \right)^{ - 1}} \\ &= \Sigma _{k - 1}^{ - 1} - \Sigma _{k - 1}^{ - 1}X_k^{\rm{T}}{\left( {I + {X_k}\Sigma _{k - 1}^{ - 1}X_k^{\rm{T}}} \right)^{ - 1}}{X_k}\Sigma _{k - 1}^{ - 1} \end{aligned} \tag{10}Σk−1=(Σk−1+XkTXk)−1=Σk−1−1−Σk−1−1XkT(I+XkΣk−1−1XkT)−1XkΣk−1−1(10)
令Pk=∑k−1{P_k} = {\sum _k}^{ - 1}Pk=∑k−1，公式(10)变为：
Pk=Pk−1−Pk−1XkT(I+XkPk−1XkT)−1XkPk−1(11)\begin{aligned} {P_k} = {P_{k - 1}} - {P_{k - 1}}X_k^{\rm{T}}{\left( {I + {X_k}{P_{k - 1}}X_k^{\rm{T}}} \right)^{ - 1}}{X_k}{P_{k - 1}} \end{aligned} \tag{11}Pk=Pk−1−Pk−1XkT(I+XkPk−1XkT)−1XkPk−1(11)
公式(9)变为：
θk=θk−1+PkXkT(Yk−Xkθk−1)(12)\begin{aligned} {\theta_k} = {\theta_{k - 1}} + {P_k}X_k^{\rm{T}}\left( {{Y_k} - {X_k}{\theta_{k - 1}}} \right) \end{aligned} \tag{12}θk=θk−1+PkXkT(Yk−Xkθk−1)(12)
注意到，公式(11)依然存在对I+XkPk−1XkT{I + {X_k}{P_{k - 1}}X_k^{\rm{T}}}I+XkPk−1XkT的求逆运算，这似乎依然没有解决上述问题1，我们避免了对Σk\Sigma_kΣk的求逆，但却又引入了一个新的逆。事实上，如果数据是逐个到达的，则XkX_kXk为一个行向量(在本文中，一个样本我们用行向量表示，这主要是因为本文规定数据矩阵中每一行代表一个样本)，因此I+XkPk−1XkT{I + {X_k}{P_{k - 1}}X_k^{\rm{T}}}I+XkPk−1XkT最终得到结果为一个数值，我们无需矩阵求逆计算，只需要取它的倒数就好了，即
Pk=Pk−1−Pk−1XkTXkPk−11+XkPk−1XkT(13)\begin{aligned} {P_k} = {P_{k - 1}} - \frac{{{P_{k - 1}}X_k^{\rm{T}}{X_k}{P_{k - 1}}}}{{1 + {X_k}{P_{k - 1}}X_k^{\rm{T}}}} \end{aligned} \tag{13}Pk=Pk−1−1+XkPk−1XkTPk−1XkTXkPk−1(13)
于是我们得到了新的递推算法如下,我们称为算法2：

根据公式(13)更新Pk=Pk−1−Pk−1XkTXkPk−11+XkPk−1XkT；{P_k} = {P_{k - 1}} - \frac{{{P_{k - 1}}X_k^{\rm{T}}{X_k}{P_{k - 1}}}}{{1 + {X_k}{P_{k - 1}}X_k^{\rm{T}}}}；Pk=Pk−1−1+XkPk−1XkTPk−1XkTXkPk−1；
根据公式(12)更新θk=θk−1+PkXkT(Yk−Xkθk−1){\theta_k} = {\theta_{k - 1}} + {P_k}X_k^{\rm{T}}\left( {{Y_k} - {X_k}{\theta_{k - 1}}} \right)θk=θk−1+PkXkT(Yk−Xkθk−1)。

一些书上的递推算法可能并非这样的形式，我们可以进一步对上述过程进行一些整理。在一些书中，Kk=PkXkT{K_k} = {P_k}X_k^{\rm{T}}Kk=PkXkT也被称为增益，Yk−Xkθk−1{Y_k} - {X_k}{\theta_{k - 1}}Yk−Xkθk−1被称为新息，顾名思义，就是引入的新信息。
Kk=PkXkT=(Pk−1−Pk−1XkT(I+XkPk−1XkT)−1XkPk−1)XkT//公式(11)结果替换得到=Pk−1XkT(I−(I+XkPk−1XkT)−1XkPk−1XkT)=Pk−1XkT(I+XkPk−1XkT)−1((I+XkPk−1XkT)−XkPk−1XkT)=Pk−1XkT(I+XkPk−1XkT)−1(14)\begin{aligned} {K_k} &= {P_k}X_k^{\rm{T}}\\ &= \left( {{P_{k - 1}} - {P_{k - 1}}X_k^{\rm{T}}{{\left( {I + {X_k}{P_{k - 1}}X_k^{\rm{T}}} \right)}^{ - 1}}{X_k}{P_{k - 1}}} \right)X_k^{\rm{T}} \quad //公式(11)结果替换得到\\ &= {P_{k - 1}}X_k^{\rm{T}}\left( {I - {{\left( {I + {X_k}{P_{k - 1}}X_k^{\rm{T}}} \right)}^{ - 1}}{X_k}{P_{k - 1}}X_k^{\rm{T}}} \right)\\ &= {P_{k - 1}}X_k^{\rm{T}}{\left( {I + {X_k}{P_{k - 1}}X_k^{\rm{T}}} \right)^{ - 1}}\left( {\left( {I + {X_k}{P_{k - 1}}X_k^{\rm{T}}} \right) - {X_k}{P_{k - 1}}X_k^{\rm{T}}} \right)\\ &= {P_{k - 1}}X_k^{\rm{T}}{\left( {I + {X_k}{P_{k - 1}}X_k^{\rm{T}}} \right)^{ - 1}} \end{aligned} \tag{14}Kk=PkXkT=(Pk−1−Pk−1XkT(I+XkPk−1XkT)−1XkPk−1)XkT//公式(11)结果替换得到=Pk−1XkT(I−(I+XkPk−1XkT)−1XkPk−1XkT)=Pk−1XkT(I+XkPk−1XkT)−1((I+XkPk−1XkT)−XkPk−1XkT)=Pk−1XkT(I+XkPk−1XkT)−1(14)
将公式(14)的结果代入到公式(11)可得
Pk=Pk−1−KkXkPk−1=(I−KkXk)Pk−1(15)\begin{aligned} {P_k} = {P_{k - 1}} - {K_k}{X_k}{P_{k - 1}} = \left( {I - {K_k}{X_k}} \right){P_{k - 1}} \end{aligned} \tag{15}Pk=Pk−1−KkXkPk−1=(I−KkXk)Pk−1(15)
于是，算法2可进一步的写为如下形式，我们称为算法3：

根据公式(14)更新模型增益Kk=Pk−1XkT(I+XkPk−1XkT)−1{K_k} = {P_{k - 1}}X_k^{\rm{T}}{\left( {I + {X_k}{P_{k - 1}}X_k^{\rm{T}}} \right)^{ - 1}}Kk=Pk−1XkT(I+XkPk−1XkT)−1；
根据公式(15)更新Pk=(I−KkXk)Pk−1{P_k} = \left( {I - {K_k}{X_k}} \right){P_{k - 1}}Pk=(I−KkXk)Pk−1；
更新回归系数θk=θk−1+Kk(Yk−Xkθk−1){\theta_k} = {\theta_{k - 1}} + {K_k}\left( {{Y_k} - {X_k}{\theta_{k - 1}}} \right)θk=θk−1+Kk(Yk−Xkθk−1)

递推最小二乘法(Recursive least square, RLS)详细推导相关推荐

[线性代数学习笔记] 线性递推数列及 Berlekamp-Massey 算法的详细推导过程
线性递推数列线性递推对于无限数列 {a0,a1,...}\{a_0,a_1,...\}{a0,a1,...} 和有限非空数列 {r0,r1,...,rm−1}\{r_{0},r_1,...,r ...
递推最小二乘法RLS公式详细推导
递推最小二乘法RLS公式详细推导整理递推最小二乘法推导过程自我整理. 递推最小二乘估计(RLS)作为一种估计方式是在最小二乘法(LS)的基础上发展来的. 最小二乘法可以解决的问题是不需要知道先验的概 ...
递推最小二乘法RLS的轮胎侧偏刚度估计(原书缺失代码已补全)
目录 1 参数辨识 1.1 最小二乘法 1.2 递推最小二乘法 RLS 1.3 具有遗忘因子 λ 的递推最小二乘法 2 轮胎线性侧偏刚度估计 2.1 RLS 算法分析 2.2 联合仿真平台的设计 ca ...
递推最小二乘法RLS：推导
递推最小二乘法RLS:推导 I use RLS as a method for parameter identification. Reference: https://blog.csdn.net/w ...
递推最小二乘法的推导和理解
递推最小二乘法的推导和理解最小二乘法快速回顾最小二乘法的推导建立误差平方将其最小化一种对最小二乘法理解的视角递推最小二乘法在线实时预测问题推导思路与详细过程将k时刻的表达式写成k-1 ...
基于遗忘因子递推最小二乘法辨识一阶RC等效电路模型
%% 基于一阶RC等效电路模型实现不同倍率下电模型参数辨识 clear clc%% 载入实验数据 % 导入hppc实验数据 load('hppc_pulse_25deg') temp = hppc_p ...
函数的递推matlab,关于递推最小二乘法辨识参数的matlab编程（含注释）
最近在做关于过热气温的动态建模问题.有现场运行的历史数据,要找出导前区和惰性区的传递函数. 对这类算法了解不多,程序读起来比较吃力,所以就转来一篇完整的辨识程序,在原有基础上进行了简化,并稍加注解一下 ...
一步完成最小二乘法、递推最小二乘法、增广最小二乘法、广义最小二乘法、辅助变量法、二步法辨识（传递函数）
用一步完成最小二乘法.递推最小二乘法.增广最小二乘法.广义最小二乘法.辅助变量法.二步法辨识如下模型的参数: 噪声的成形滤波器采样时间0.01 要求: 1.用matlab 写出程序代码: 2.画出实 ...
【c++递推递归算法】放苹果（详细代码+图解+解题思路）
递推递归算法解决放苹果问题:把M个同样的苹果放在N个同样的盘子里,允许有的盘子空着不放,问共有多少种不同的分法?(用K表示)5,1,1和1,5,1 是同一种分法. 1.[题目描述] 2.[解题思路] ...

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

递推最小二乘法(Recursive least square, RLS)详细推导相关推荐

最新文章

热门文章