递推最小二乘法推导
递推最小二乘法（Recursive Least Squares, RLS）是一种时间序列分析方法，它用于在线更新线性回归模型的参数，而不需要重新拟合整个数据集。这种方法特别适用于数据流或实时系统。
1.最小二乘法原理
最小二乘法是一种数学优化技术，用于拟合数据点的最佳函数模型。在线性回归的情况下，我们试图找到一组参数
X
X
X ,使得模型预测 $ \hat{Y} = AX$ 与实际数据
Y
Y
Y 之间的平方误差之和最小。其中，
A
A
A 是特征矩阵，
X
X
X 是参数向量。
Y
=
(
y
1
y
2
y
3
⋮
y
m
)
,
A
=
(
a
11
⋯
a
1
n
⋮
⋮
a
m
1
⋯
a
m
n
)
,
X
=
(
x
1
x
2
y
3
⋮
x
n
)
Y=\begin{pmatrix} y_1\\y_2\\y_3\\\vdots\\y_m \end{pmatrix}, A=\begin{pmatrix} a_{11} & \cdots & a_{1n}\\ \vdots & &\vdots \\a_{m1} & \cdots & a_{mn} \end{pmatrix}, X=\begin{pmatrix} x_1\\x_2\\y_3\\\vdots\\x_n \end{pmatrix}
Y=y1y2y3⋮ym,A=a11⋮am1⋯⋯a1n⋮amn,X=x1x2y3⋮xn
式中，
A
A
A 为输入的数据样本，
X
X
X 为拟合的参数，
Y
Y
Y 为输入的数据样本，
Y
^
\hat{Y}
Y^ 为输出的估计
为了估计出
n
n
n 个参数，要求
m
≥
n
m\ge n
m≥n ,当
m
=
n
m=n
m=n 时，存在唯一解。
m
>
n
m>n
m>n 时，由于数据中存在模型噪声和测量噪声，一般不可能确定一组参数，定义误差矢量
E
=
Y
−
A
X
^
,
E
=
[
e
1
,
⋯
,
e
m
]
T
E=Y-A\hat{X},E=[e_1,\cdots,e_m]^T
E=Y−AX^,E=[e1,⋯,em]T
最小二乘法思想：求测量值与估计值的平方和最小，即：
J
=
E
T
E
=
∑
i
=
1
m
e
i
2
=
(
Y
−
A
X
^
)
T
(
Y
−
A
X
^
)
=
Y
T
Y
−
Y
T
A
X
^
−
X
^
T
A
T
Y
+
X
^
T
A
T
A
X
^
J=E^TE=\sum_{i=1}^{m}e_i^2=(Y-A\hat{X})^T(Y-A\hat{X})=Y^TY-Y^TA\hat{X}-\hat{X}^TA^TY+\hat{X}^TA^TA\hat{X}
J=ETE=∑i=1mei2=(Y−AX^)T(Y−AX^)=YTY−YTAX^−X^TATY+X^TATAX^
对
J
J
J 求偏导：
X
^
T
A
T
A
X
^
\hat{X}^TA^TA\hat{X}
X^TATAX^是一个二次项，其偏导数需要使用矩阵的Frechet导数或利用以下规则：如果
f
(
X
)
=
X
T
B
X
f(X)=X^TBX
f(X)=XTBX，其中
B
B
B是一个常数矩阵，则
d
f
d
X
=
2
B
X
\frac{df}{dX}=2BX
dXdf=2BX在这里，
B
=
A
T
A
B=A^TA
B=ATA。
J
J
J 是标量，X是向量，求导是对X的每一个元素求偏导，
d
J
/
d
X
dJ/dX
dJ/dX 的维度是
n
×
1
n\times1
n×1
所以，
d
J
d
X
^
=
−
2
A
T
Y
+
2
A
T
A
X
^
\frac{dJ}{d\hat{X}}=-2A^TY+2A^TA\hat{X}
dX^dJ=−2ATY+2ATAX^
令偏导等于0，即求出最优的拟合参数
X
^
\hat{X}
X^
−
2
A
T
Y
+
2
A
T
A
X
^
=
0
⇒
X
^
=
(
A
T
A
)
−
1
A
T
Y
-2A^TY+2A^TA\hat{X}=0 \Rightarrow\hat{X}=(A^TA)^{-1}A^TY
−2ATY+2ATAX^=0⇒X^=(ATA)−1ATY
这样，当我们有了足够的样本，可以一次性求解出最优的参数向量X
2.递归最小二乘法
设第
k
−
1
k-1
k−1 时刻的样本矩阵
A
k
−
1
,
Y
k
−
1
A_{k-1},Y_{k-1}
Ak−1,Yk−1 分别为：
A
k
−
1
=
(
a
11
⋯
a
1
n
⋮
⋮
a
(
k
−
1
)
1
⋯
a
(
k
−
1
)
n
)
,
Y
k
−
1
=
(
y
1
y
2
y
3
⋮
y
k
−
1
)
A_{k-1}=\begin{pmatrix} a_{11} & \cdots & a_{1n}\\ \vdots & &\vdots \\a_{(k-1)1} & \cdots & a_{(k-1)n} \end{pmatrix},Y_{k-1}=\begin{pmatrix} y_1\\y_2\\y_3\\\vdots\\y_{k-1} \end{pmatrix}
Ak−1=a11⋮a(k−1)1⋯⋯a1n⋮a(k−1)n,Yk−1=y1y2y3⋮yk−1
那么，第k时刻的样本矩阵可以表示为：
A
k
=
(
A
k
−
1
a
k
)
,
Y
k
=
(
Y
k
−
1
y
k
)
A_k=\begin{pmatrix} A_{k-1} \\ a_k \end{pmatrix},Y_k=\begin{pmatrix} Y_{k-1} \\ y_k \end{pmatrix}
Ak=(Ak−1ak),Yk=(Yk−1yk)
计算中间式：
A
k
T
A
k
=
[
A
k
−
1
T
a
k
T
]
[
A
k
−
1
a
k
]
=
A
k
−
1
T
A
k
−
1
+
a
k
T
a
k
A_k^TA_k=\begin{bmatrix}A_{k-1}^T & a_k^T\end{bmatrix}\begin{bmatrix}A_{k-1} \\ a_k\end{bmatrix}=A_{k-1}^TA_{k-1}+a_k^Ta_k
AkTAk=[Ak−1TakT][Ak−1ak]=Ak−1TAk−1+akTak
A
k
T
Y
k
=
[
A
k
−
1
T
a
k
T
]
[
Y
k
−
1
y
k
]
=
A
k
−
1
T
Y
k
−
1
+
a
k
T
y
k
A_k^TY_k=\begin{bmatrix}A_{k-1}^T & a_k^T\end{bmatrix}\begin{bmatrix}Y_{k-1} \\ y_k\end{bmatrix}=A_{k-1}^TY_{k-1}+a_k^Ty_k
AkTYk=[Ak−1TakT][Yk−1yk]=Ak−1TYk−1+akTyk
由第一节推导出的最优拟合结果公式：
X
^
k
=
(
A
k
T
A
k
)
−
1
A
k
T
Y
k
\hat{X}_{k}=(A_k^TA_k)^{-1}A_k^TY_k
X^k=(AkTAk)−1AkTYk
向前递推得：
X
^
k
−
1
=
(
A
k
−
1
T
A
k
−
1
)
−
1
A
k
−
1
T
Y
k
−
1
⇒
A
k
−
1
T
Y
k
−
1
=
(
A
k
−
1
T
A
k
−
1
)
X
^
k
−
1
\hat{X}_{k-1}=(A_{k-1}^TA_{k-1})^{-1}A_{k-1}^TY_{k-1}\Rightarrow A_{k-1}^TY_{k-1}=(A_{k-1}^TA_{k-1})\hat{X}_{k-1}
X^k−1=(Ak−1TAk−1)−1Ak−1TYk−1⇒Ak−1TYk−1=(Ak−1TAk−1)X^k−1
那么：
X
^
k
=
(
A
k
T
A
k
)
−
1
A
k
T
Y
k
=
(
A
k
T
A
k
)
−
1
(
A
k
−
1
T
Y
k
−
1
+
a
k
T
y
k
)
=
(
A
k
T
A
k
)
−
1
(
A
k
−
1
T
A
k
−
1
X
^
k
−
1
+
a
k
T
y
k
)
\hat{X}_{k}=(A_k^TA_k)^{-1}A_k^TY_k=(A_k^TA_k)^{-1}(A_{k-1}^TY_{k-1}+a_k^Ty_k)=(A_k^TA_k)^{-1}(A_{k-1}^TA_{k-1}\hat{X}_{k-1}+a_k^Ty_k)
X^k=(AkTAk)−1AkTYk=(AkTAk)−1(Ak−1TYk−1+akTyk)=(AkTAk)−1(Ak−1TAk−1X^k−1+akTyk)
将
A
k
−
1
T
A
k
−
1
A_{k-1}^TA_{k-1}
Ak−1TAk−1 用中间式1进行替换：
X
^
k
=
(
A
k
T
A
k
)
−
1
(
(
A
k
T
A
k
−
a
k
T
a
k
)
X
^
k
−
1
+
a
k
T
y
k
)
=
[
I
−
(
A
k
T
A
k
)
−
1
a
k
T
a
k
]
X
^
k
−
1
+
(
A
k
T
A
k
)
−
1
a
k
T
y
k
=
X
^
k
−
1
+
(
A
k
T
A
k
)
−
1
a
k
T
(
y
k
−
a
k
X
^
k
−
1
)
\hat{X}_{k}=(A_k^TA_k)^{-1}((A_k^TA_k-a_k^Ta_k)\hat{X}_{k-1}+a_k^Ty_k)=[I-(A_k^TA_k)^{-1}a_k^Ta_k]\hat{X}_{k-1}+(A_k^TA_k)^{-1}a_k^Ty_k=\hat{X}_{k-1}+(A_k^TA_k)^{-1}a_k^T(y_k-a_k\hat{X}_{k-1})
X^k=(AkTAk)−1((AkTAk−akTak)X^k−1+akTyk)=[I−(AkTAk)−1akTak]X^k−1+(AkTAk)−1akTyk=X^k−1+(AkTAk)−1akT(yk−akX^k−1)
即迭代公式初步化简为：
X
^
k
=
X
^
k
−
1
+
(
A
k
T
A
k
)
−
1
a
k
T
(
y
k
−
a
k
X
^
k
−
1
)
\hat{X}_{k}=\hat{X}_{k-1}+(A_k^TA_k)^{-1}a_k^T(y_k-a_k\hat{X}_{k-1})
X^k=X^k−1+(AkTAk)−1akT(yk−akX^k−1)
考虑到
(
A
k
T
A
k
)
−
1
(A_k^TA_k)^{-1}
(AkTAk)−1 计算复杂，使用谢尔曼公式对矩阵逆运算进行化简。
谢尔曼公式为：
(
A
+
u
v
)
−
1
=
A
−
1
−
(
A
−
1
u
)
(
I
+
v
A
−
1
u
)
−
1
(
v
A
−
1
)
(A+uv)^{-1}=A^{-1}-(A^{-1}u)(I+vA^{-1}u)^{-1}(vA^{-1})
(A+uv)−1=A−1−(A−1u)(I+vA−1u)−1(vA−1)
结合中间式1求解
(
A
k
T
A
k
)
−
1
(A_k^TA_k)^{-1}
(AkTAk)−1
令
G
k
−
1
=
(
A
k
T
A
k
)
−
1
G_k^{-1}=(A_k^TA_k)^{-1}
Gk−1=(AkTAk)−1 ,带入中间式1则
G
k
−
1
=
(
G
k
−
1
+
a
k
T
a
k
)
−
1
G_k^{-1}=(G_{k-1}+a_k^Ta_k)^{-1}
Gk−1=(Gk−1+akTak)−1
套用谢尔曼公式：
G
k
−
1
=
G
k
−
1
−
1
−
（
G
k
−
1
−
1
a
k
T
）
(
I
+
a
k
G
k
−
1
−
1
a
k
T
)
−
1
(
a
k
G
k
−
1
−
1
)
G_k^{-1}=G_{k-1}^{-1}-（G_{k-1}^{-1}a_k^T）(I+a_kG_{k-1}^{-1}a_k^T)^{-1}(a_kG_{k-1}^{-1})
Gk−1=Gk−1−1−（Gk−1−1akT）(I+akGk−1−1akT)−1(akGk−1−1)
不难看出：
(
I
+
a
k
G
k
−
1
−
1
a
k
T
)
−
1
=
I
1
+
a
k
G
k
−
1
−
1
a
k
T
(I+a_kG_{k-1}^{-1}a_k^T)^{-1}=\frac{I}{1+a_kG_{k-1}^{-1}a_k^T}
(I+akGk−1−1akT)−1=1+akGk−1−1akTI
并且令
P
k
=
G
k
−
1
P_k=G_k^{-1}
Pk=Gk−1
可以得出另一个递推公式：
P
k
=
P
k
−
1
−
P
k
−
1
a
k
T
a
k
P
k
−
1
1
+
a
k
P
k
−
1
a
k
T
P_k=P_{k-1}-\frac{P_{k-1}a_k^Ta_kP_{k-1}}{1+a_kP_{k-1}a_k^T}
Pk=Pk−1−1+akPk−1akTPk−1akTakPk−1
所以最小二乘法递推形式表示为：
{
P
k
=
P
k
−
1
−
P
k
−
1
a
k
T
a
k
P
k
−
1
1
+
a
k
P
k
−
1
a
k
T
X
^
k
=
X
^
k
−
1
+
P
k
a
k
T
(
y
k
−
a
k
X
^
k
−
1
)
\begin{cases} P_k=P_{k-1}-\frac{P_{k-1}a_k^Ta_kP_{k-1}}{1+a_kP_{k-1}a_k^T} \\ \hat{X}_{k}=\hat{X}_{k-1}+P_ka_k^T(y_k-a_k\hat{X}_{k-1}) \end{cases}
{Pk=Pk−1−1+akPk−1akTPk−1akTakPk−1X^k=X^k−1+PkakT(yk−akX^k−1)
本文转载自: https://blog.csdn.net/weixin_53084505/article/details/140746005
版权归原作者 认真听讲的学渣 所有， 如有侵权，请联系我们删除。
版权归原作者 认真听讲的学渣 所有， 如有侵权，请联系我们删除。