贴一篇我在知乎的文章,这里面也阐述了最小二乘法与卡尔曼滤波的关系。两个算法一起说了,而且学习的话,比较轻量级。
https://zhuanlan.zhihu.com/p/67250500
一般最小二乘法
一般最小二乘法是给定若干观测值,计算一个最有可能的估计。
[
y
(
1
)
y
(
2
)
⋮
y
(
k
)
]
=
[
h
(
1
)
1
h
(
1
)
2
⋯
h
(
1
)
n
h
(
2
)
1
h
(
2
)
2
⋯
h
(
2
)
n
⋮
⋮
⋱
⋮
h
(
k
)
1
h
(
k
)
2
⋯
h
(
k
)
n
]
[
x
1
x
2
⋮
x
n
]
\left [\begin{matrix}y(1)\\ y(2)\\ \vdots\\ y(k) \end{matrix}\right]=\left [\begin{matrix}h(1)_1&h(1)_2& \cdots&h(1)_n\\ h(2)_1&h(2)_2& \cdots&h(2)_n\\ \vdots&\vdots&\ddots&\vdots\\ h(k)_1&h(k)_2& \cdots&h(k)_n \end{matrix}\right]\left [\begin{matrix}x_1\\ x_2\\ \vdots\\ x_n \end{matrix}\right]
⎣⎢⎢⎢⎡y(1)y(2)⋮y(k)⎦⎥⎥⎥⎤=⎣⎢⎢⎢⎡h(1)1h(2)1⋮h(k)1h(1)2h(2)2⋮h(k)2⋯⋯⋱⋯h(1)nh(2)n⋮h(k)n⎦⎥⎥⎥⎤⎣⎢⎢⎢⎡x1x2⋮xn⎦⎥⎥⎥⎤
计算上述最有可能的
x
x
x取值,是使得代价函数:
J
(
k
)
=
(
Y
(
k
)
−
H
(
k
)
x
)
T
(
Y
(
k
)
−
H
(
k
)
x
)
=
Y
T
(
k
)
Y
(
k
)
−
Y
T
(
k
)
H
(
k
)
x
−
x
T
H
T
(
x
)
Y
(
k
)
+
x
T
H
(
k
)
T
H
(
k
)
x
J(k)=(Y(k)-H(k)x)^T(Y(k)-H(k)x)\\ =Y^T(k)Y(k)-Y^T(k)H(k)x-x^TH^T(x)Y(k)+x^TH(k)^TH(k)x
J(k)=(Y(k)−H(k)x)T(Y(k)−H(k)x)=YT(k)Y(k)−YT(k)H(k)x−xTHT(x)Y(k)+xTH(k)TH(k)x
最小的
x
x
x的取值。令其为
x
^
(
k
)
\hat x(k)
x^(k),意思是
k
k
k时刻的估计值
上述中,
H
(
k
)
H(k)
H(k)代表
k
k
k时刻以及
k
k
k时刻以前的观测,即为:
H
(
k
)
=
[
h
(
1
)
T
,
h
(
2
)
T
,
⋯
h
(
k
)
T
]
T
H(k)=[h(1)^T,\ h(2)^T,\cdots h(k)^T]^T
H(k)=[h(1)T, h(2)T,⋯h(k)T]T
Y
(
k
)
Y(k)
Y(k)类似:
Y
(
k
)
=
[
y
(
1
)
,
y
(
2
)
,
⋯
y
(
k
)
]
T
Y(k)=[y(1),y(2),\cdots y(k)]^T
Y(k)=[y(1),y(2),⋯y(k)]T
当
J
(
k
)
J(k)
J(k)最小时,有:
∂
J
(
k
)
∂
x
=
−
2
Y
T
(
k
)
H
(
k
)
+
2
x
T
H
T
(
k
)
H
(
k
)
=
0
\frac{\partial J(k)}{\partial x}=-2Y^T(k)H(k)+2x^TH^T(k)H(k)=0
∂x∂J(k)=−2YT(k)H(k)+2xTHT(k)H(k)=0
此时,
x
=
x
^
(
k
)
x = \hat x(k)
x=x^(k)
所以:
x
^
(
k
)
=
(
H
T
(
k
)
H
(
k
)
)
−
H
T
(
k
)
Y
(
k
)
\hat x(k)=(H^T(k)H(k))^-H^T(k)Y(k)
x^(k)=(HT(k)H(k))−HT(k)Y(k)
只有
H
T
(
k
)
H
(
k
)
H^T(k)H(k)
HT(k)H(k) 满秩时,才存在逆矩阵。
所以只有当
k
≥
n
k\geq n
k≥n时,
H
T
(
k
)
H
(
k
)
H^T(k)H(k)
HT(k)H(k)才可能存在逆矩阵。
带权最小二乘法
有时候,测量有好有坏,带权的最小二乘法就比较有必要了。可以人为赋予每个数据的一个置信度
r
(
k
)
r(k)
r(k),有时候也会令
r
2
(
k
)
=
1
D
(
v
(
k
)
)
=
1
σ
k
2
r^2(k)=\frac{1}{D(v(k))}=\frac{1}{\sigma_k^2}
r2(k)=D(v(k))1=σk21
其中,
v
(
k
)
v(k)
v(k) 是第
k
k
k测量的误差。区别于估计误差。上面的式子也比较好理解。测量的不确定性会削弱这次测量的置信度。如果利用测量误差的话,则不能存在
0
0
0误差的情况。改写
J
(
k
)
J(k)
J(k)的形式:
J
(
k
)
=
∑
i
=
1
k
r
2
(
k
)
(
y
(
k
)
−
y
^
(
k
)
)
2
J(k)=\sum_{i=1}^kr^2(k)(y(k)-\hat y(k))^2
J(k)=i=1∑kr2(k)(y(k)−y^(k))2
令:
R
(
k
)
=
[
r
2
(
1
)
0
⋯
0
0
r
2
(
2
)
⋯
0
⋮
⋮
⋱
⋮
0
0
⋯
r
2
(
n
)
]
R(k)=\left [\begin{matrix}r^2(1)&0& \cdots&0\\ 0&r^2(2)& \cdots&0\\ \vdots&\vdots&\ddots&\vdots\\ 0&0& \cdots&r^2(n) \end{matrix}\right]
R(k)=⎣⎢⎢⎢⎡r2(1)0⋮00r2(2)⋮0⋯⋯⋱⋯00⋮r2(n)⎦⎥⎥⎥⎤
那么:
J
(
k
)
=
(
Y
(
k
)
−
H
(
k
)
x
^
)
T
R
(
Y
(
k
)
−
H
(
k
)
x
^
)
=
Y
(
k
)
T
R
(
k
)
Y
(
k
)
−
Y
T
(
k
)
R
(
k
)
H
(
k
)
x
^
−
x
^
T
H
(
k
)
T
R
(
k
)
Y
(
k
)
+
x
^
T
H
T
(
k
)
R
(
K
)
H
T
(
k
)
x
^
J(k)=\Big(Y(k)-H(k)\hat x\Big)^TR\Big(Y(k)-H(k)\hat x\Big)\\=Y(k)^TR(k)Y(k)-Y^T(k)R(k)H(k)\hat x-\hat x^TH(k)^TR(k)Y(k)+\hat x^TH^T(k)R(K)H^T(k)\hat x
J(k)=(Y(k)−H(k)x^)TR(Y(k)−H(k)x^)=Y(k)TR(k)Y(k)−YT(k)R(k)H(k)x^−x^TH(k)TR(k)Y(k)+x^THT(k)R(K)HT(k)x^
令:
∂
J
(
k
)
∂
x
^
=
−
2
Y
T
(
k
)
R
(
k
)
+
2
x
T
H
T
(
k
)
R
(
k
)
H
(
k
)
=
0
\frac{\partial J(k)}{\partial\hat x}=-2Y^T(k)R(k)+2x^TH^T(k)R(k)H(k)=0
∂x^∂J(k)=−2YT(k)R(k)+2xTHT(k)R(k)H(k)=0
则:
x
^
(
k
)
=
(
H
T
(
k
)
R
(
k
)
H
(
k
)
)
−
H
T
(
k
)
R
(
k
)
Y
(
k
)
\hat x(k)= (H^T(k)R(k)H(k))^-H^T(k)R(k)Y(k)
x^(k)=(HT(k)R(k)H(k))−HT(k)R(k)Y(k)
这样设置权重,其实是对噪声的归一化。
R L S RLS RLS 递推算法
有时候,我们不可能一次获得所有的观测结果,有时候随着观测结果的增加。从新修改参数使用最小二乘法计算,会很耗费时间。 R L S RLS RLS 就是等价于一般最小二乘法的递推计算方法。
这里介绍一个精巧的矩阵反演:
在这只要求
A
,
D
A,D
A,D可逆的情况下给定分块矩阵:
[
A
B
C
D
]
\left [\begin{matrix} A& B\\ C&D \end{matrix}\right]
[ACBD]
令:
E
=
D
−
C
A
−
B
,
F
=
A
−
B
D
−
C
E =D-CA^-B,\ \ F = A-BD^-C
E=D−CA−B, F=A−BD−C
计算可得:
[
A
B
C
D
]
[
A
−
+
A
−
B
E
−
C
A
−
−
A
−
B
E
−
−
E
−
C
A
−
E
−
]
=
[
I
O
C
A
−
+
C
A
−
B
E
−
C
A
−
−
D
E
−
C
A
−
−
C
A
−
B
E
−
+
D
E
−
]
=
[
I
O
C
A
−
+
(
C
A
−
B
−
D
)
E
−
C
A
−
(
−
C
A
−
B
+
D
)
E
−
]
=
[
I
O
O
I
]
\left [\begin{matrix} A&B\\ C&D \end{matrix}\right]\left [\begin{matrix} A^-+A^-BE^-CA^-&-A^-BE^-\\ -E^-CA^-&E^- \end{matrix}\right]\\=\left [\begin{matrix} I&O\\ CA^-+CA^-BE^-CA^--DE^-CA^-&-CA^-BE^-+DE^- \end{matrix}\right]\\=\left [\begin{matrix} I&O\\ CA^-+(CA^-B-D)E^-CA^-&(-CA^-B+D)E^- \end{matrix}\right]=\left [\begin{matrix} I&O\\ O&I \end{matrix}\right]
[ACBD][A−+A−BE−CA−−E−CA−−A−BE−E−]=[ICA−+CA−BE−CA−−DE−CA−O−CA−BE−+DE−]=[ICA−+(CA−B−D)E−CA−O(−CA−B+D)E−]=[IOOI]
另一个方向:
[
A
B
C
D
]
[
F
−
−
A
−
B
E
−
−
E
−
C
A
−
E
−
]
=
[
I
O
O
I
]
\left [\begin{matrix} A&B\\ C&D \end{matrix}\right]\left [\begin{matrix}F^-&-A^-BE^-\\-E^-CA^-&E^- \end{matrix}\right]=\left [\begin{matrix} I&O\\ O&I \end{matrix}\right]
[ACBD][F−−E−CA−−A−BE−E−]=[IOOI]
这个方向可以自行计算。
这也就是说,当矩阵
A
,
D
,
E
,
F
A,D,E,F
A,D,E,F可逆时,:
(
A
+
B
D
−
C
)
−
=
A
−
−
A
−
B
(
D
+
C
A
−
B
)
−
C
A
−
(A+BD^-C)^-=A^--A^-B(D+CA^-B)^-CA^-
(A+BD−C)−=A−−A−B(D+CA−B)−CA−
令
P
(
k
)
=
(
H
T
(
k
)
H
(
k
)
)
−
P(k)=(H^T(k)H(k))^-
P(k)=(HT(k)H(k))−
则有:
P
(
k
)
=
(
H
T
(
k
)
H
(
k
)
)
−
=
(
H
T
(
k
−
1
)
H
(
k
−
1
)
+
h
T
(
k
)
h
(
k
)
)
−
=
(
P
−
(
k
−
1
)
+
h
T
(
k
)
h
(
k
)
)
−
=
P
(
k
−
1
)
−
P
(
k
−
1
)
h
T
(
k
)
h
(
k
)
P
(
k
−
1
)
1
+
h
(
k
)
P
(
k
−
1
)
h
T
(
k
)
P(k)=(H^T(k)H(k))^-\\=(H^T(k-1)H(k-1)+h^T(k)h(k))^-\\=(P^-(k-1)+h^T(k)h(k))^-\\=P(k-1)-\frac{P(k-1)h^T(k)h(k)P(k-1)}{1+h(k)P(k-1)h^T(k)}
P(k)=(HT(k)H(k))−=(HT(k−1)H(k−1)+hT(k)h(k))−=(P−(k−1)+hT(k)h(k))−=P(k−1)−1+h(k)P(k−1)hT(k)P(k−1)hT(k)h(k)P(k−1)
又有:
x
^
(
k
)
=
P
(
k
)
H
T
(
k
)
Y
(
k
)
\hat x(k)=P(k)H^T(k)Y(k)
x^(k)=P(k)HT(k)Y(k)
可得:
P
−
(
k
)
x
^
(
k
)
=
H
T
(
k
)
Y
(
k
)
P^-(k)\hat x(k) = H^T(k)Y(k)
P−(k)x^(k)=HT(k)Y(k)
又因为:
P
(
k
)
=
(
P
−
(
k
−
1
)
+
h
T
(
k
)
h
(
k
)
)
−
P(k)=(P^-(k-1)+h^T(k)h(k))^-
P(k)=(P−(k−1)+hT(k)h(k))−
得:
P
−
(
k
−
1
)
=
P
−
(
k
)
−
h
T
(
k
)
h
(
k
)
P^-(k-1)=P^-(k)-h^T(k)h(k)
P−(k−1)=P−(k)−hT(k)h(k)
所以:
x
^
(
k
)
=
P
(
k
)
H
T
(
k
)
Y
(
k
)
=
P
(
k
)
(
P
−
(
k
−
1
)
x
^
(
k
−
1
)
+
h
T
(
k
)
y
(
k
)
)
=
P
(
k
)
(
(
P
−
(
k
)
−
h
T
(
k
)
h
(
k
)
)
x
^
(
k
−
1
)
+
h
T
(
k
)
y
(
k
)
)
=
x
^
(
k
−
1
)
+
P
(
k
)
h
T
(
k
)
(
y
(
k
)
−
h
(
k
)
x
^
(
k
−
1
)
)
\hat x(k) = P(k)H^T(k)Y(k)\\=P(k)\Big(P^-(k-1)\hat x(k-1)+h^T(k)y(k)\Big)\\=P(k)\Big((P^-(k)-h^T(k)h(k))\hat x(k-1)+h^T(k)y(k)\Big)\\=\hat x(k-1)+P(k)h^T(k)\Big(y(k)-h(k)\hat x(k-1)\Big)
x^(k)=P(k)HT(k)Y(k)=P(k)(P−(k−1)x^(k−1)+hT(k)y(k))=P(k)((P−(k)−hT(k)h(k))x^(k−1)+hT(k)y(k))=x^(k−1)+P(k)hT(k)(y(k)−h(k)x^(k−1))
令:
K
(
k
)
=
P
(
k
−
1
)
h
T
(
k
)
1
+
h
(
k
)
P
(
k
−
1
)
h
T
(
k
)
=
P
(
k
)
h
T
(
k
)
K(k)=\frac{P(k-1)h^T(k)}{1+h(k)P(k-1)h^T(k)}=P(k)h^T(k)
K(k)=1+h(k)P(k−1)hT(k)P(k−1)hT(k)=P(k)hT(k)
则:
P
(
k
)
=
(
I
−
K
(
k
)
h
(
k
)
)
P
(
k
−
1
)
x
^
(
k
)
=
x
^
(
k
−
1
)
+
K
(
k
)
(
y
(
k
)
−
h
(
k
)
x
^
(
k
−
1
)
)
P(k)=(I-K(k)h(k))P(k-1)\\\hat x(k)=\hat x(k-1)+K(k)(y(k)-h(k)\hat x(k-1))
P(k)=(I−K(k)h(k))P(k−1)x^(k)=x^(k−1)+K(k)(y(k)−h(k)x^(k−1))
如何确定
P
(
0
)
P(0)
P(0) 呢?
显然:
P
(
k
)
=
(
∑
i
=
1
k
h
T
(
i
)
h
(
i
)
)
−
≈
(
∑
i
=
1
k
h
T
(
i
)
h
(
i
)
+
1
∞
I
)
−
P(k)=\Big(\sum_{i=1}^kh^T(i)h(i)\Big)^- ≈ \Big(\sum_{i=1}^kh^T(i)h(i)+\frac{1}{∞ }I\Big)^-
P(k)=(i=1∑khT(i)h(i))−≈(i=1∑khT(i)h(i)+∞1I)−
令
P
(
0
)
=
∞
I
P(0) = ∞I
P(0)=∞I即可。那么此时
y
(
0
)
=
h
T
(
0
)
=
x
^
(
0
)
=
O
n
×
1
y(0)=h^T(0)=\hat x(0) = O_{n\times 1}
y(0)=hT(0)=x^(0)=On×1
这样定义对代价函数影响忽略不计。
带权RLS 递推计算
当考虑带权递推时
x
^
(
k
)
=
(
H
T
(
k
)
R
(
k
)
H
(
k
)
)
−
H
T
(
k
)
R
(
k
)
Y
(
k
)
\hat x(k)= \Big(H^T(k)R(k)H(k)\Big)^-H^T(k)R(k)Y(k)
x^(k)=(HT(k)R(k)H(k))−HT(k)R(k)Y(k)
一般
R
L
S
RLS
RLS递推的
P
(
k
)
P(k)
P(k)重新定义为:
P
(
k
)
=
(
H
T
(
k
)
R
(
k
)
H
(
k
)
)
−
=
(
H
T
(
k
−
1
)
R
(
k
−
1
)
H
(
k
−
1
)
+
h
T
(
k
)
r
(
k
)
h
(
k
)
)
−
P(k) = (H^T(k)R(k)H(k))^-=\Big(H^T(k-1)R(k-1)H(k-1)+h^T(k)r(k)h(k)\Big)^-
P(k)=(HT(k)R(k)H(k))−=(HT(k−1)R(k−1)H(k−1)+hT(k)r(k)h(k))−
令:
γ
(
k
)
=
1
r
(
k
)
\gamma(k) = \frac{1}{r(k)}
γ(k)=r(k)1
矩阵反演有:
P
(
k
)
=
P
(
k
−
1
)
−
P
(
k
−
1
)
h
T
(
k
)
h
(
k
)
P
(
k
−
1
)
γ
(
k
)
−
h
(
k
)
P
(
k
−
1
)
h
T
(
k
)
P(k)=P(k-1)-\frac{P(k-1)h^T(k)h(k)P(k-1)}{\gamma(k)-h(k)P(k-1)h^T(k)}
P(k)=P(k−1)−γ(k)−h(k)P(k−1)hT(k)P(k−1)hT(k)h(k)P(k−1)
由:
x
^
(
k
)
=
P
(
k
)
H
T
(
k
)
R
(
k
)
Y
(
k
)
\hat x(k)=P(k)H^T(k)R(k)Y(k)
x^(k)=P(k)HT(k)R(k)Y(k)
有:
P
−
(
k
)
x
^
(
k
)
=
H
T
(
k
)
R
(
k
)
Y
(
k
)
P^-(k)\hat x(k)=H^T(k)R(k)Y(k)
P−(k)x^(k)=HT(k)R(k)Y(k)
由:
P
(
k
)
=
(
P
−
(
k
−
1
)
+
h
T
(
k
)
r
(
k
)
h
(
k
)
)
P(k)=\Big(P^-(k-1)+h^T(k)r(k)h(k)\Big)
P(k)=(P−(k−1)+hT(k)r(k)h(k))
有:
P
−
(
k
−
1
)
=
P
−
(
k
)
−
h
T
(
k
)
r
(
k
)
h
(
k
)
P^-(k-1)=P^-(k)-h^T(k)r(k)h(k)
P−(k−1)=P−(k)−hT(k)r(k)h(k)
则:
x
^
(
k
)
=
P
(
k
)
H
T
(
k
)
R
(
k
)
Y
(
k
)
=
P
(
k
)
(
P
−
(
k
−
1
)
x
^
(
k
−
1
)
+
h
T
(
k
)
r
(
k
)
y
(
k
)
)
=
P
(
k
)
(
(
P
−
(
k
)
−
h
T
(
k
)
r
(
k
)
h
(
k
)
)
x
^
(
k
−
1
)
+
h
T
(
k
)
r
(
k
)
y
(
k
)
)
=
x
^
(
k
−
1
)
+
P
(
k
)
h
T
(
k
)
r
(
k
)
(
y
(
k
)
−
h
(
k
)
x
(
k
)
)
\hat x(k)=P(k)H^T(k)R(k)Y(k)\\=P(k)\Big(P^-(k-1)\hat x(k-1)+h^T(k)r(k)y(k)\Big)\\=P(k)\Big((P^-(k)-h^T(k)r(k)h(k))\hat x(k-1)+h^T(k)r(k)y(k)\Big)\\ =\hat x(k-1)+P(k)h^T(k)r(k)\Big(y(k)-h(k)x(k)\Big)
x^(k)=P(k)HT(k)R(k)Y(k)=P(k)(P−(k−1)x^(k−1)+hT(k)r(k)y(k))=P(k)((P−(k)−hT(k)r(k)h(k))x^(k−1)+hT(k)r(k)y(k))=x^(k−1)+P(k)hT(k)r(k)(y(k)−h(k)x(k))
令
K
(
k
)
=
P
(
k
)
h
T
(
k
)
r
(
k
)
=
P
(
k
−
1
)
h
T
(
k
)
γ
(
k
)
−
h
(
k
)
P
(
k
−
1
)
h
T
(
k
)
K(k)=P(k)h^T(k)r(k)=\frac{P(k-1)h^T(k)}{\gamma(k)-h(k)P(k-1)h^T(k)}
K(k)=P(k)hT(k)r(k)=γ(k)−h(k)P(k−1)hT(k)P(k−1)hT(k)
则:
P
(
k
)
=
(
I
−
K
(
k
)
h
(
k
)
)
P
(
k
−
1
)
P(k)=\Big(I-{K(k)h(k)}\Big)P(k-1)
P(k)=(I−K(k)h(k))P(k−1)
x
^
(
k
)
=
x
^
(
k
−
1
)
+
K
(
k
)
(
y
(
k
)
−
h
(
k
)
x
^
(
k
−
1
)
)
\hat x(k) = \hat x(k-1)+K(k)(y(k)-h(k)\hat x(k-1))
x^(k)=x^(k−1)+K(k)(y(k)−h(k)x^(k−1))
同上,当没有很好的初始设置数值时,可以讲,
P
(
0
)
=
∞
I
P(0)=∞I
P(0)=∞I,
x
^
(
0
)
=
h
T
(
0
)
=
O
n
×
1
\hat x(0) = h^T(0)=O_{n\times 1}
x^(0)=hT(0)=On×1
带有遗忘的 R L S RLS RLS
对于带有遗忘版本的
R
L
S
RLS
RLS算法,将对近期数据更为敏感。
设置遗忘因子
λ
∈
(
0
,
1
]
\lambda \in (0,1]
λ∈(0,1]
当
λ
=
1
\lambda = 1
λ=1时,不会遗忘。
每次更新,历史数据的权重会被整体缩小
λ
\lambda
λ
那么此时:
J
(
k
)
=
λ
J
(
k
−
1
)
+
(
y
(
k
)
−
h
(
k
)
x
^
(
k
)
)
2
J(k)=\lambda J(k-1) +(y(k)-h(k)\hat x(k))^2
J(k)=λJ(k−1)+(y(k)−h(k)x^(k))2
则:
J
(
k
)
=
λ
∣
∣
H
(
k
−
1
)
x
^
(
k
)
−
Y
(
k
−
1
)
∣
∣
2
+
(
y
(
k
)
−
h
(
k
)
x
^
(
k
)
)
2
J(k)=\lambda ||H(k-1) \hat x(k)-Y(k-1)||^2+(y(k)-h(k)\hat x(k))^2
J(k)=λ∣∣H(k−1)x^(k)−Y(k−1)∣∣2+(y(k)−h(k)x^(k))2
从新定义
H
(
k
)
=
[
β
H
(
k
−
1
)
,
h
(
k
)
]
H(k)=[\beta H(k-1),\ h(k)]
H(k)=[βH(k−1), h(k)]
Y
(
k
)
=
[
β
Y
(
k
−
1
)
,
y
(
k
)
]
Y(k)=[\beta Y(k-1),\ y(k)]
Y(k)=[βY(k−1), y(k)]
其中,
β
2
=
λ
\beta^2= \lambda
β2=λ
P
(
k
)
=
(
H
T
(
k
)
H
(
k
)
)
−
=
(
H
T
(
k
−
1
)
H
(
k
−
1
)
λ
+
h
T
(
k
)
h
(
k
)
)
−
P(k)=(H^T(k)H(k))^-\\=(H^T(k-1)H(k-1)\lambda+h^T(k)h(k))^-
P(k)=(HT(k)H(k))−=(HT(k−1)H(k−1)λ+hT(k)h(k))−
=
(
P
−
(
k
−
1
)
λ
+
h
T
(
k
)
h
(
k
)
)
−
=(P^-(k-1)\lambda+h^T(k)h(k))^-
=(P−(k−1)λ+hT(k)h(k))−
=
P
(
k
−
1
)
1
λ
−
1
λ
2
P
(
k
−
1
)
h
T
(
k
)
h
(
k
)
P
(
k
−
1
)
1
+
1
λ
h
(
k
)
P
(
k
−
1
)
h
T
(
k
)
=P(k-1)\frac{1}{\lambda}-\frac{\frac{1}{\lambda^2}P(k-1)h^T(k)h(k)P(k-1)}{1+\frac{1}{\lambda}h(k)P(k-1)h^T(k)}
=P(k−1)λ1−1+λ1h(k)P(k−1)hT(k)λ21P(k−1)hT(k)h(k)P(k−1)
=
(
I
−
P
(
k
−
1
)
h
T
(
k
)
h
(
k
)
λ
+
h
(
k
)
P
(
k
−
1
)
h
T
(
k
)
)
P
(
k
−
1
)
1
λ
=\Big(I-\frac{P(k-1)h^T(k)h(k)}{\lambda+h(k)P(k-1)h^T(k)}\Big)P(k-1)\frac{1}{\lambda}
=(I−λ+h(k)P(k−1)hT(k)P(k−1)hT(k)h(k))P(k−1)λ1
另一边:
x
^
(
k
)
=
(
H
T
(
k
)
H
(
k
)
)
−
H
T
(
k
)
Y
(
k
)
=
P
(
k
)
H
T
(
k
)
Y
(
k
)
\hat x(k)=(H^T(k)H(k))^-H^T(k)Y(k)\\ =P(k)H^T(k)Y(k)
x^(k)=(HT(k)H(k))−HT(k)Y(k)=P(k)HT(k)Y(k)
x
^
(
k
)
=
P
(
k
)
H
T
(
k
)
Y
(
k
)
\hat x(k)=P(k)H^T(k)Y(k)
x^(k)=P(k)HT(k)Y(k)
可得:
P
−
(
k
)
x
^
(
k
)
=
H
T
(
k
)
Y
(
k
)
P^-(k)\hat x(k)=H^T(k)Y(k)
P−(k)x^(k)=HT(k)Y(k)
同时:
P
(
k
)
=
(
P
−
(
k
−
1
)
λ
+
h
T
(
k
)
h
(
k
)
)
−
P(k)=(P^-(k-1)\lambda +h^T(k)h(k))^-
P(k)=(P−(k−1)λ+hT(k)h(k))−
P
−
(
k
−
1
)
=
1
λ
(
P
−
(
k
)
−
h
T
(
k
)
h
(
k
)
)
P^-(k-1)=\frac{1}{\lambda}(P^-(k)-h^T(k)h(k))
P−(k−1)=λ1(P−(k)−hT(k)h(k))
x
^
(
k
)
=
P
(
k
)
(
H
T
(
k
−
1
)
Y
(
k
−
1
)
λ
+
h
T
(
k
)
y
(
k
)
)
=
P
(
k
)
(
P
−
(
k
−
1
)
x
^
(
k
−
1
)
λ
+
h
T
(
k
)
y
(
k
)
)
=
P
(
k
)
(
(
P
−
(
k
)
−
h
T
(
k
)
h
(
k
)
)
x
^
(
k
−
1
)
+
h
T
(
k
)
y
(
k
)
)
=
x
^
(
k
−
1
)
+
P
(
k
)
h
T
(
k
)
(
y
(
k
)
−
h
(
k
)
x
^
(
k
−
1
)
)
\hat x(k)=P(k)(H^T(k-1)Y(k-1)\lambda+h^T(k)y(k))\\ =P(k)(P^-(k-1)\hat x(k-1)\lambda +h^T(k)y(k))\\ =P(k)\Big(\big(P^-(k)-h^T(k)h(k)\big)\hat x(k-1)+h^T(k)y(k)\Big)\\ =\hat x(k-1)+P(k)h^T(k)\Big(y(k)-h(k)\hat x(k-1)\Big)
x^(k)=P(k)(HT(k−1)Y(k−1)λ+hT(k)y(k))=P(k)(P−(k−1)x^(k−1)λ+hT(k)y(k))=P(k)((P−(k)−hT(k)h(k))x^(k−1)+hT(k)y(k))=x^(k−1)+P(k)hT(k)(y(k)−h(k)x^(k−1))
令:
K
(
k
)
=
P
(
k
)
h
T
(
k
)
K(k)=P(k)h^T(k)
K(k)=P(k)hT(k)
那么:
x
^
(
k
)
=
x
^
(
k
−
1
)
+
K
(
k
)
(
y
(
k
)
−
h
(
k
)
x
^
(
k
−
1
)
)
\hat x(k)=\hat x(k-1)+K(k)\Big(y(k)-h(k)\hat x(k-1)\Big)
x^(k)=x^(k−1)+K(k)(y(k)−h(k)x^(k−1))
K
(
k
)
=
P
(
k
)
h
(
k
)
=
(
I
−
P
(
k
−
1
)
h
T
(
k
)
h
(
k
)
λ
+
h
(
k
)
P
(
k
−
1
)
h
T
(
k
)
)
P
(
k
−
1
)
1
λ
h
T
(
k
)
=
P
(
k
−
1
)
h
T
(
k
)
λ
+
h
(
k
)
P
(
k
−
1
)
h
T
(
k
)
K(k)=P(k)h(k)=\Big(I-\frac{P(k-1)h^T(k)h(k)}{\lambda+h(k)P(k-1)h^T(k)}\Big)P(k-1)\frac{1}{\lambda}h^T(k)\\ =\frac{P(k-1)h^T(k)}{\lambda +h(k)P(k-1)h^T(k)}
K(k)=P(k)h(k)=(I−λ+h(k)P(k−1)hT(k)P(k−1)hT(k)h(k))P(k−1)λ1hT(k)=λ+h(k)P(k−1)hT(k)P(k−1)hT(k)
P
(
k
)
=
(
I
−
K
(
k
)
h
(
k
)
)
P
(
k
−
1
)
1
λ
P(k)=(I-K(k)h(k))P(k-1)\frac{1}{\lambda}
P(k)=(I−K(k)h(k))P(k−1)λ1
总结
RLS 算法一个比较精巧的在线滤波算法。 还可以增加遗忘因子,让其对近期近期数据更为敏感,比如笔记平滑种可能需要这种形式的滤波器。