递推最小二乘遗忘因子法(Recursive Forgetting Factor, RFF)

本文探讨了数据饱和现象在递推最小二乘算法中的影响,介绍了带遗忘因子的方法,通过遗忘因子调整数据权重,避免模型过饱和。通过Sherman-Morrison-Woodbury公式,提出了两种改进算法:算法2利用逆矩阵递推和算法3简化为增益与新息的更新。
摘要由CSDN通过智能技术生成

在普通的递推最小二乘算法中,随着数据的不断到来,显然矩阵 X T X X^TX XTX中的元素会变得越来越大,而矩阵 P k P_k Pk作为 X T X X^TX XTX的逆矩阵,则会逐渐趋于零,这时,模型将无法继续更新或者更新极其缓慢,这就是数据饱和现象。我们也可以更加直观的理解为,由于过去累加了很多的数据,当前到来的数据与之前累加的数据相比就如一滴水掉落到大海,不会惊起任何波澜。针对数据饱和问题,其中一个解决方案就是带遗忘因子的递推最小二乘法。

假设有数据 ( X , Y ) (X,Y) (X,Y),其中 X ∈ R m × d X \in {\mathbb{R}^{m \times d}} XRm×d Y ∈ R m × 1 Y \in {\mathbb{R}^{m \times 1}} YRm×1 m m m为样本数, d d d为特征数,考虑最小二乘解

θ 0 = ( X T X ) − 1 X T Y = Σ 0 − 1 X T Y (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=Σ01XTY(1)

⇒ Σ 0 θ 0 = X T Y (2) \Rightarrow {\Sigma_0}{\theta_0} = {X^{\rm{T}}}Y \tag{2} Σ0θ0=XTY(2)

当新数据 ( X 1 , Y 1 ) \left( {{X_1},{Y_1}} \right) (X1,Y1)到来时,更新模型。我们希望当前的数据对于回归结果更加重要,而过去数据的重要性随着时间的回溯依次降低,从而使得模型能够更好的适应当前数据的变化,克服数据饱和问题。为了达到上述目的,我们可以给之前数据的损失乘以一个权重 α , 0 < α < 1 \alpha, 0<\alpha<1 α,0<α<1,即
α ∥ Y − X θ ∥ 2 2 + ∥ Y 1 − X 1 θ ∥ 2 2 = ∥ [ α Y Y 1 ] − [ α X X 1 ] θ ∥ 2 2 \alpha \left\| {Y - X\theta } \right\|_2^2 + \left\| {{Y_1} - {X_1}\theta } \right\|_2^2 = \left\| {\left[ {\begin{array}{cc} {\sqrt \alpha Y} \\ {{Y_1}} \end{array}} \right] - \left[ {\begin{array}{cc} {\sqrt \alpha X} \\ {{X_1}} \end{array}} \right]\theta } \right\|_2^2 αYXθ22+Y1X1θ22=[α YY1][α XX1]θ22
于是得到新的回归系数

θ 1 = ( [ α X X 1 ] T [ α X X 1 ] ) − 1 [ α X X 1 ] T [ α Y Y 1 ] = Σ 1 − 1 [ α X X 1 ] T [ α Y Y 1 ] (3) \begin{aligned} {\theta_1} &= {\left( {{{\left[ {\begin{array}{cc} {\sqrt \alpha }X\\ {{X_1}} \end{array}} \right]}^{\rm{T}}}\left[ {\begin{array}{cc} {\sqrt \alpha }X\\ {{X_1}} \end{array}} \right]} \right)^{ - 1}}{\left[ {\begin{array}{cc} {\sqrt \alpha }X\\ {{X_1}} \end{array}} \right]^{\rm{T}}}\left[ {\begin{array}{cc} {\sqrt \alpha }Y\\ {{Y_1}} \end{array}} \right] \\ &= {\Sigma _1}^{ - 1}{\left[ {\begin{array}{cc} {\sqrt \alpha }X\\ {{X_1}} \end{array}} \right]^{\rm{T}}}\left[ {\begin{array}{cc} {\sqrt \alpha }Y\\ {{Y_1}} \end{array}} \right]\tag{3} \end{aligned} θ1=([α XX1]T[α XX1])1[α XX1]T[α YY1]=Σ11[α XX1]T[α YY1](3)

其中

Σ 1 = [ α X X 1 ] T [ α X X 1 ] = α X T X + X 1 T X 1 = α Σ 0 + X 1 T X 1 (4) \begin{aligned} {\Sigma _1} &= {\left[ {\begin{array}{cc} {\sqrt \alpha }X\\ {{X_1}} \end{array}} \right]^{\rm{T}}}\left[ {\begin{array}{cc} {\sqrt \alpha }X\\ {{X_1}} \end{array}} \right] \\ &= \alpha {X^{\rm{T}}}X + X_1^{\rm{T}}{X_1} \\ &= \alpha {\Sigma _0} + X_1^{\rm{T}}{X_1} \end{aligned} \tag{4} Σ1=[α XX1]T[α XX1]=αXTX+X1TX1=αΣ0+X1TX1(4)

⇒ α Σ 0 = Σ 1 − X 1 T X 1 (5) \begin{aligned} \Rightarrow \alpha {\Sigma _0} = {\Sigma _1} - X_1^{\rm{T}}{X_1} \end{aligned} \tag{5} αΣ0=Σ1X1TX1(5)

根据公式(4)的结果,通过归纳可得
Σ k = α Σ k − 1 + X k T X k (6) \begin{aligned} {\Sigma _k} = \alpha {\Sigma _{k - 1}} + X_k^{\rm{T}}{X_k} \end{aligned} \tag{6} Σk=αΣk1+XkTXk(6)
从公式(6)可以很容易看出,这是一个典型的套娃行为,当前的权重为1,上一次的权重为 α \alpha α,通过一次套娃 α Σ k − 1 = α 2 Σ k − 2 + α X k − 1 T X k − 1 \alpha{\Sigma _{k-1}} = \alpha^2 {\Sigma _{k - 2}} + \alpha X_{k-1}^{\rm{T}}{X_{k-1}} αΣk1=α2Σk2+αXk1TXk1,我们发现上上次的权重变成了 α 2 \alpha^2 α2,依次类推,最终我们发现随着时间的回溯,越久的数据权重会越来越低并逐渐趋于0。
[ α X X 1 ] T [ α Y Y 1 ] = α X T Y + X 1 T Y 1 = α Σ 0 θ 0 + X 1 T Y 1 / / 公 式 ( 2 ) 结 果 替 换 得 到 = ( Σ 1 − X 1 T X 1 ) θ 0 + X 1 T Y 1 / / 公 式 ( 5 ) 结 果 替 换 得 到 = Σ 1 θ 0 + X 1 T ( Y 1 − X 1 θ 0 ) (7) \begin{aligned} {\left[ {\begin{array}{cc} {\sqrt \alpha }X\\ {{X_1}} \end{array}} \right]^{\rm{T}}}\left[ {\begin{array}{cc} {\sqrt \alpha }Y\\ {{Y_1}} \end{array}} \right] &= \alpha {X^{\rm{T}}}Y + X_1^{\rm{T}}{Y_1}\\ &= \alpha {\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)=(Σ1X1TX1)θ0+X1TY1//(5)=Σ1θ0+X1T(Y1X1θ0)(7)

将公式(7)回带到公式(3):

θ 1 = Σ 1 − 1 ( Σ 1 θ 0 + X 1 T ( Y 1 − X 1 θ 0 ) ) = θ 0 + Σ 1 − 1 X 1 T ( Y 1 − X 1 θ 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=Σ11(Σ1θ0+X1T(Y1X1θ0))=θ0+Σ11X1T(Y1X1θ0)(8)

根据公式(8)的结果,通过归纳可得
θ k = θ k − 1 + Σ k − 1 X k T ( Y k − X k θ 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=θk1+Σk1XkT(YkXkθk1)(9)

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

  1. 根据公式(5)更新 Σ k = α Σ k − 1 + X k T X k {\Sigma _k} = \alpha {\Sigma _{k - 1}} + X_k^{\rm{T}}{X_k} Σk=αΣk1+XkTXk
  2. 根据公式(9)更新 θ k = θ k − 1 + Σ k − 1 X k T ( Y k − X k θ k − 1 ) {\theta_k} = {\theta_{k - 1}} + {\Sigma _k}^{ - 1}X_k^{\rm{T}}\left( {{Y_k} - {X_k}{\theta_{k - 1}}} \right) θk=θk1+Σk1XkT(YkXkθk1)

但以上过程存在一个问题:

  1. 对矩阵 Σ k \Sigma_k Σk的求逆计算复杂度比较高,我们能否在递推过程中避免对 Σ k \Sigma_k Σk的求逆计算,而直接更新它的逆矩阵;

针对以上问题,我们要对公式进一步改造

根据Sherman-Morrison-Woodbury公式:
( A + U V T ) − 1 = A − 1 − A − 1 U ( I + V T A − 1 U ) − 1 V T A − 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=A1A1U(I+VTA1U)1VTA1
公式(6)的逆可写成如下形式

Σ k − 1 = ( α Σ k − 1 + X k T X k ) − 1 = 1 α Σ k − 1 − 1 − 1 α Σ k − 1 − 1 X k T ( α I + X k Σ k − 1 − 1 X k T ) − 1 X k Σ k − 1 − 1 (10) \begin{aligned} {\Sigma _k}^{ - 1} &= {\left( \alpha {{\Sigma _{k - 1}} + X_k^{\rm{T}}{X_k}} \right)^{ - 1}} \\ &= \frac{1}{\alpha } \Sigma _{k - 1}^{ - 1} - \frac{1}{\alpha }\Sigma _{k - 1}^{ - 1}X_k^{\rm{T}}{\left( {\alpha I + {X_k}\Sigma _{k - 1}^{ - 1}X_k^{\rm{T}}} \right)^{ - 1}}{X_k}\Sigma _{k - 1}^{ - 1} \end{aligned} \tag{10} Σk1=(αΣk1+XkTXk)1=α1Σk11α1Σk11XkT(αI+XkΣk11XkT)1XkΣk11(10)
P k = ∑ k − 1 {P_k} = {\sum _k}^{ - 1} Pk=k1,公式(10)变为:
P k = 1 α P k − 1 − 1 α P k − 1 X k T ( α I + X k P k − 1 X k T ) − 1 X k P k − 1 (11) \begin{aligned} {P_k} = \frac{1}{\alpha }{P_{k - 1}} - \frac{1}{\alpha }{P_{k - 1}}X_k^{\rm{T}}{\left( {\alpha I + {X_k}{P_{k - 1}}X_k^{\rm{T}}} \right)^{ - 1}}{X_k}{P_{k - 1}} \end{aligned} \tag{11} Pk=α1Pk1α1Pk1XkT(αI+XkPk1XkT)1XkPk1(11)
公式(9)变为:
θ k = θ k − 1 + P k X k T ( Y k − X k θ 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=θk1+PkXkT(YkXkθk1)(12)
注意到,公式(11)依然存在对 α I + X k P k − 1 X k T {\alpha I + {X_k}{P_{k - 1}}X_k^{\rm{T}}} αI+XkPk1XkT 的求逆运算,这似乎依然没有解决上述问题1,我们避免了对 Σ k \Sigma_k Σk 的求逆,但却又引入了一个新的逆。事实上,如果数据是逐个到达的,则 X k X_k Xk 为一个行向量(在本文中,一个样本我们用行向量表示,这主要是因为本文规定数据矩阵中每一行代表一个样本),因此 α I + X k P k − 1 X k T {\alpha I + {X_k}{P_{k - 1}}X_k^{\rm{T}}} αI+XkPk1XkT 最终得到结果为一个数值,我们无需矩阵求逆计算,只需要取它的倒数就好了,即
P k = 1 α P k − 1 − 1 α P k − 1 X k T X k P k − 1 α + X k P k − 1 X k T (13) \begin{aligned} {P_k} = \frac{1}{\alpha } {P_{k - 1}} - \frac{1}{\alpha } \frac{{{P_{k - 1}}X_k^{\rm{T}}{X_k}{P_{k - 1}}}}{{\alpha + {X_k}{P_{k - 1}}X_k^{\rm{T}}}} \end{aligned} \tag{13} Pk=α1Pk1α1α+XkPk1XkTPk1XkTXkPk1(13)
于是我们得到了新的递推算法如下,我们称为算法2:

  1. 根据公式(13)更新 P k = 1 α P k − 1 − 1 α P k − 1 X k T X k P k − 1 α + X k P k − 1 X k T ; {P_k} = \frac{1}{\alpha } {P_{k - 1}} - \frac{1}{\alpha } \frac{{{P_{k - 1}}X_k^{\rm{T}}{X_k}{P_{k - 1}}}}{{\alpha + {X_k}{P_{k - 1}}X_k^{\rm{T}}}}; Pk=α1Pk1α1α+XkPk1XkTPk1XkTXkPk1
  2. 根据公式(12)更新 θ k = θ k − 1 + P k X k T ( Y k − X k θ k − 1 ) {\theta_k} = {\theta_{k - 1}} + {P_k}X_k^{\rm{T}}\left( {{Y_k} - {X_k}{\theta_{k - 1}}} \right) θk=θk1+PkXkT(YkXkθk1)

一些书上的递推算法可能并非这样的形式,我们可以进一步对上述过程进行一些整理。在一些书中, K k = P k X k T {K_k} = {P_k}X_k^{\rm{T}} Kk=PkXkT也被称为增益, Y k − X k θ k − 1 {Y_k} - {X_k}{\theta_{k - 1}} YkXkθk1被称为新息,顾名思义,就是引入的新信息。
K k = P k X k T = ( 1 α P k − 1 − 1 α P k − 1 X k T ( α I + X k P k − 1 X k T ) − 1 X k P k − 1 ) X k T / / 公 式 ( 11 ) 结 果 替 换 得 到 = P k − 1 X k T ( 1 α I − 1 α ( α I + X k P k − 1 X k T ) − 1 X k P k − 1 X k T ) = P k − 1 X k T ( α I + X k P k − 1 X k T ) − 1 ( 1 α ( α I + X k P k − 1 X k T ) − 1 α X k P k − 1 X k T ) = P k − 1 X k T ( α I + X k P k − 1 X k T ) − 1 (14) \begin{aligned} {K_k} &= {P_k}X_k^{\rm{T}}\\ &= \left( \frac{1}{\alpha } {{P_{k - 1}} - \frac{1}{\alpha }{P_{k - 1}}X_k^{\rm{T}}{{\left( {\alpha 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( {\frac{1}{\alpha } I - \frac{1}{\alpha } {{\left( {\alpha 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( {\alpha I + {X_k}{P_{k - 1}}X_k^{\rm{T}}} \right)^{ - 1}}\left( \frac{1}{\alpha } {\left( {\alpha I + {X_k}{P_{k - 1}}X_k^{\rm{T}}} \right) - \frac{1}{\alpha } {X_k}{P_{k - 1}}X_k^{\rm{T}}} \right)\\ &= {P_{k - 1}}X_k^{\rm{T}}{\left( {\alpha I + {X_k}{P_{k - 1}}X_k^{\rm{T}}} \right)^{ - 1}} \end{aligned} \tag{14} Kk=PkXkT=(α1Pk1α1Pk1XkT(αI+XkPk1XkT)1XkPk1)XkT//(11)=Pk1XkT(α1Iα1(αI+XkPk1XkT)1XkPk1XkT)=Pk1XkT(αI+XkPk1XkT)1(α1(αI+XkPk1XkT)α1XkPk1XkT)=Pk1XkT(αI+XkPk1XkT)1(14)
将公式(14)的结果代入到公式(11)可得
P k = 1 α P k − 1 − 1 α K k X k P k − 1 = 1 α ( I − K k X k ) P k − 1 (15) \begin{aligned} {P_k} = \frac{1}{\alpha } {P_{k - 1}} - \frac{1}{\alpha } {K_k}{X_k}{P_{k - 1}} = \frac{1}{\alpha } \left( {I - {K_k}{X_k}} \right){P_{k - 1}} \end{aligned} \tag{15} Pk=α1Pk1α1KkXkPk1=α1(IKkXk)Pk1(15)
于是,算法2可进一步的写为如下形式,我们称为算法3:

  1. 根据公式(14)更新模型增益 K k = P k − 1 X k T ( α I + X k P k − 1 X k T ) − 1 {K_k} = {P_{k - 1}}X_k^{\rm{T}}{\left( {\alpha I + {X_k}{P_{k - 1}}X_k^{\rm{T}}} \right)^{ - 1}} Kk=Pk1XkT(αI+XkPk1XkT)1
  2. 根据公式(15)更新 P k = 1 α ( I − K k X k ) P k − 1 {P_k} = \frac{1}{\alpha } \left( {I - {K_k}{X_k}} \right){P_{k - 1}} Pk=α1(IKkXk)Pk1
  3. 更新回归系数 θ k = θ k − 1 + K k ( Y k − X k θ k − 1 ) {\theta_k} = {\theta_{k - 1}} + {K_k}\left( {{Y_k} - {X_k}{\theta_{k - 1}}} \right) θk=θk1+Kk(YkXkθk1)
下面是一个简单的MATLAB脚本,实现了多新息递推最小二乘辨识算法(MI-RLS): ```matlab clear; clf; format short g % 设置参数 p = 1; % 创新长度 p=1,5,10 length1 = 1000; % 数据长度 lambda = 0.98; % 遗忘因子 % 多项式A(z),B(z),D(z)的阶数 na = 2; nb = 2; nc = 2; nd = 2; nl = na + nb; n2 = nc + nd; n = na + nb + nc + nd; % A(z),B(z),D(z)的系数向量 a = [1, -0.56, 0.42]; b = [0, 0.9, 0.6]; c = [1, -0.3, 0.2]; d = [1, 0.3, 0, -0.20]; c1 = conv(a, c); theta = [a(2:na+1); b(2:nb+1); c(2:nc+1); d(2:nd+1)]; % 参数向量 theta n = length(theta); fprintf('\nMI-RLS算法\n'); fprintf('创新长度p=%d\n', p); % 准备输入数据和噪声数据 rand('state', 15); % 设置随机变量的状态 u = (rand(length1, 1) - 0.5) * sqrt(12); % 输入随机序列 randn('state', 15); v = randn(length1, 1); % 噪声随机序列 % 生成输出数据 y = zeros(length1, 1); w = zeros(n, 1); for t = n:length1 w(t) = theta' * [-w(t-1:-1:t-nc); v(t-1:-1:t-nd)] + v(t); y(t) = theta(1:na)' * [-y(t-1:-1:t-na); u(t-1:-1:t-nb)] + w(t); end % MI-RLS算法 Phi = eye(n*p) / lambda; % 初始协方差矩阵 theta_hat = zeros(n*p, 1); % 初始参数估计 for t = p:length1 phi = [-y(t-1:-1:t-na); u(t-1:-1:t-nb); -w(t-1:-1:t-nc); v(t-1:-1:t-nd)]; K = Phi * phi / (lambda + phi' * Phi * phi); % RLS增益 e = y(t) - phi' * theta_hat; % 预测误差 theta_hat = theta_hat + K * e; % 参数更新 Phi = (Phi - K * phi' * Phi) / lambda; % 协方差更新 end % 输出结果 fprintf('真实参数向量:\n'); disp(theta); fprintf('估计参数向量:\n'); disp(theta_hat); ``` 请注意,这只是一个简单的示例,可能需要根据您的具体需求进行修改和优化。希望对您有所帮助!
评论 4
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值