基本最小二乘到递推最小二乘
基本最小二乘(LS)
先导知识:
从函数出发
假设一个函数 y = [ θ 1 θ 2 ⋯ θ n ] [ x 1 x 2 ⋮ x n ] = θ X = ∑ i = 1 n θ i x i y= \begin{bmatrix} \theta_1& \theta_2& \cdots& \theta_n \end{bmatrix} \begin{bmatrix} x_1\\ x_2\\ \vdots\\ x_n \end{bmatrix} =\boldsymbol{\theta}X =\sum_{i=1}^n{\theta_i x_i} y=[θ1θ2⋯θn]⎣⎢⎢⎢⎡x1x2⋮xn⎦⎥⎥⎥⎤=θX=i=1∑nθixi
我们约定粗体和大写字母均表示矩阵或者向量。
残差
现在假设我们有一个这样的函数函数
y
=
θ
X
y=\boldsymbol{\theta}X
y=θX但是我们并不知道
θ
\boldsymbol{\theta}
θ的各个量的具体数值。我们只能用一系列(一共
k
k
k组)的
X
j
=
[
x
1
j
x
2
j
⋮
x
n
j
]
∈
R
n
×
1
(
j
=
1
,
2
,
⋯
,
k
)
X_j=\begin{bmatrix}x_1^j\\x_2^j\\ \vdots\\x_n^j\end{bmatrix}\in\mathbb{R}^{n\times1}\;\;(j=1,2,\cdots,k)
Xj=⎣⎢⎢⎢⎡x1jx2j⋮xnj⎦⎥⎥⎥⎤∈Rn×1(j=1,2,⋯,k)作为输入往这个函数的入口放(这里的
j
j
j并不是次方的意思,而是表示这是第
j
j
j组数据中的某个量)。然后我们得到了
k
k
k组输出
Y
=
[
y
1
y
2
⋯
y
k
]
∈
R
1
×
k
Y=\begin{bmatrix}y_1&y_2&\cdots&y_k\end{bmatrix}\in \mathbb{R}^{1\times k}
Y=[y1y2⋯yk]∈R1×k
也即:
Y
=
θ
[
X
1
X
2
⋯
X
k
]
∈
(
R
1
×
n
∗
R
n
×
k
)
=
R
1
×
k
Y=\boldsymbol{\theta}\begin{bmatrix}X_1&X_2&\cdots &X_k\end{bmatrix}\in(\mathbb{R}^{1\times n}*\mathbb{R}^{n\times k} )=\mathbb{R}^{1\times k}
Y=θ[X1X2⋯Xk]∈(R1×n∗Rn×k)=R1×k
我们可以估计出一个
θ
\boldsymbol{\theta}
θ的取值
θ
^
\hat{\boldsymbol{\theta}}
θ^。我们再把
X
j
(
j
=
1
,
2
,
⋯
,
k
)
X_j(j=1,2,\cdots,k)
Xj(j=1,2,⋯,k)代入我们估计出的方程模型
y
=
θ
^
X
y=\hat{\boldsymbol{\theta}}X
y=θ^X。由于估计的模型不可能完全拟合原始数据,所以我们会得到一组拟合值
Y
^
=
θ
^
[
X
1
X
2
⋯
X
j
]
∈
R
1
×
k
\hat{Y}=\hat{\boldsymbol{\theta}} \begin{bmatrix}X_1&X_2&\cdots&X_j\end{bmatrix} \in\mathbb{R}^{1\times k}
Y^=θ^[X1X2⋯Xj]∈R1×k。
我们把
e
=
Y
−
Y
^
∈
R
1
×
k
\boldsymbol{e}=Y-\hat{Y}\in\mathbb{R}^{1\times k}
e=Y−Y^∈R1×k称为残差
梳理
我们在这里将一些常数和符号做一个约定:
- n——函数的输入变量 x i x_i xi的个数,同时也是 θ i \theta_i θi的总数
- k——测试用的数据总组数
- Φ \boldsymbol\Phi Φ—— [ X 1 X 2 ⋯ X k ] ∈ R n × k \begin{bmatrix}X_1&X_2&\cdots&X_k\end{bmatrix} \in\mathbb{R}^{n\times k} [X1X2⋯Xk]∈Rn×k
- i i i——一般作为输入数据及其参数的下标( θ i x i \theta_i x_i θixi)
- j j j——一般作为输入组数数据的下标( X j X_j Xj)
推导
基本思想:
我们再估计参数
θ
^
\hat{\boldsymbol{\theta}}
θ^时,我们自然是希望与客观存在的值
θ
\boldsymbol{\theta}
θ尽可能的接近。因此,我们需要引入一个评价我们估计的好坏的一个“标准”。前面提到的残差可以作为一个很好的标准。
但是我们发现
e
∈
R
1
×
k
e\in\mathbb{R}^{1\times k}
e∈R1×k作为一个向量,不好作为一个直观的标准。同时它的每一项中既有正也有负,自然是不适合直接作为评价标准来使用。所以我们引入一个指标函数
J
J
J
J
=
∑
k
=
n
+
1
n
+
N
e
2
(
k
)
=
d
e
f
e
⋅
e
T
=
(
Y
−
θ
^
Φ
)
(
Y
−
θ
^
Φ
)
T
J=\sum_{k=n+1}^{n+N}{e^2(k)}\overset{def}{=}\boldsymbol{e}\cdot\boldsymbol{e}^T=(\boldsymbol{Y}-\hat{\boldsymbol{\theta}}\Phi)(\boldsymbol{Y}-\hat{\boldsymbol{\theta}}\Phi)^T
J=k=n+1∑n+Ne2(k)=defe⋅eT=(Y−θ^Φ)(Y−θ^Φ)T
开始推导
所以,最小二乘法就是让指标函数
J
J
J最小的参数估计方法。既有:
θ
^
=
min
θ
J
\hat{\boldsymbol{\theta}}=\min_{\boldsymbol{\theta}}J
θ^=θminJ
而
J
J
J取最小值,我们先讨论J为极值时:
∂
J
∂
θ
^
=
0
\frac{\partial J}{\partial \hat{\boldsymbol{\theta}}}=0
∂θ^∂J=0
⇒ ∂ [ ( Y − θ ^ Φ ) ( Y − θ ^ Φ ) T ] ∂ θ ^ = 0 \Rightarrow\frac{\partial[(\boldsymbol{Y}-\hat{\boldsymbol{\theta}}\boldsymbol\Phi)(\boldsymbol{Y}-\hat{\boldsymbol{\theta}}\boldsymbol\Phi)^T]}{\partial \hat{\boldsymbol{\theta}}}=0 ⇒∂θ^∂[(Y−θ^Φ)(Y−θ^Φ)T]=0
⇒ − 2 Φ ( Y − θ ^ Φ ) T = 0 \Rightarrow -2\boldsymbol\Phi(\boldsymbol Y-\boldsymbol{\hat{\theta}\Phi})^T=0 ⇒−2Φ(Y−θ^Φ)T=0
Φ Φ T θ ^ T = Φ Y T \boldsymbol{\Phi\Phi^T\hat{\theta}^T}=\boldsymbol{\Phi Y^T} ΦΦTθ^T=ΦYT
其中,
Φ
Φ
T
∈
R
n
×
n
\boldsymbol{\Phi\Phi^T}\in\mathbb{R}^{n\times n}
ΦΦT∈Rn×n。
若其逆矩阵存在,则:
θ
^
T
=
(
Φ
Φ
T
)
−
1
Φ
Y
T
\boldsymbol{\hat{\theta}^T}=(\boldsymbol{\Phi\Phi^T})^{-1}\boldsymbol{\Phi Y^T}
θ^T=(ΦΦT)−1ΦYT
θ ^ = Y Φ T ( Φ Φ T ) − 1 \boldsymbol{\hat{\theta}}=\boldsymbol{Y \Phi^T}(\boldsymbol{\Phi\Phi^T})^{-1} θ^=YΦT(ΦΦT)−1
上述结果只是
J
J
J为极值时的结论,
J
J
J可能是极大值也可能是极小值。我们进一步讨论,要使
J
J
J为极小值的条件为
∂
2
J
∂
θ
^
2
>
0
\frac{\partial^2J}{\partial \boldsymbol{\hat{\theta}^2}}>0
∂θ^2∂2J>0
∂ J ∂ θ ^ = − 2 Φ ( Y − Φ θ ^ ) T \frac{\partial J}{\partial \hat{\boldsymbol{\theta}}}=-2\boldsymbol\Phi(\boldsymbol Y-\boldsymbol{\Phi\hat{\theta}})^T ∂θ^∂J=−2Φ(Y−Φθ^)T
⇒ ∂ 2 J ∂ θ ^ 2 = 2 Φ Φ T > 0 \Rightarrow\frac{\partial^2J}{\partial \boldsymbol{\hat{\theta}^2}}=2\boldsymbol{\Phi\Phi^T}>0 ⇒∂θ^2∂2J=2ΦΦT>0
也就是 Φ Φ T \boldsymbol{\Phi\Phi^T} ΦΦT为正定矩阵。
递推最小二乘法
- 暂留2020/11/4,于2020/11/13完成
背景
基本最小二乘法(LS)有诸多缺点,例如对于一组动态的数据,每次接收到新数据,就要全部重算一遍。这种重复的计算的成本很大,导致实用性不好。所以对于一组离线数据,基本最小二乘法是适用的。但是如果是实时统计分析一系列数据,那么基本最小二乘法就会遇到计算困难。
我们希望在获取一个新数据时,可以直接使用该数据和已经计算过的结果进行某种运算,达到“修正”旧结果的目的。
前N个输入输出数据
为此,我们假设在某个时间点已经获取了 N N N组数据:
- Φ N = [ X 1 X 2 ⋯ X N ] ∈ R n × N \boldsymbol{\Phi_N}=\begin{bmatrix}X_1&X_2&\cdots &X_N\end{bmatrix}\in\mathbb{R}^{n\times N} ΦN=[X1X2⋯XN]∈Rn×N(N组输入,每一组有n个分量)
- Y N = [ y 1 , y 2 , ⋯ , y N ] = θ Φ N = θ [ X 1 X 2 ⋯ X N ] ∈ ( R 1 × n ∗ R n × N ) = R 1 × N \boldsymbol{Y_N}=\begin{bmatrix}y_1,y_2,\cdots,y_N\end{bmatrix}=\boldsymbol{\theta\Phi_N}=\boldsymbol{\theta}\begin{bmatrix}X_1&X_2&\cdots &X_N\end{bmatrix}\in(\mathbb{R}^{1\times n}*\mathbb{R}^{n\times N})=\mathbb{R}^{1\times N} YN=[y1,y2,⋯,yN]=θΦN=θ[X1X2⋯XN]∈(R1×n∗Rn×N)=R1×N(就是前N个输入对应的N个输出)
- θ ^ N = Y N Φ N T ( Φ N Φ N T ) − 1 \boldsymbol{\hat{\theta}_N}=\boldsymbol{Y_N}\boldsymbol{\Phi_N}^T(\boldsymbol{\Phi_N\Phi_N^T})^{-1} θ^N=YNΦNT(ΦNΦNT)−1
- θ ~ N = θ − θ ^ N \boldsymbol{\tilde{\theta}_N}=\boldsymbol{\theta}-\boldsymbol{\hat{\theta}_N} θ~N=θ−θ^N
- V a r θ ^ = σ 2 ( Φ N Φ N T ) − 1 Var\boldsymbol{\hat{\theta}}=\sigma^2(\boldsymbol{\Phi_N\Phi_N^T})^{-1} Varθ^=σ2(ΦNΦNT)−1
我们记
P
N
=
(
Φ
N
Φ
N
T
)
−
1
∈
R
n
×
n
\boldsymbol{P_N}=(\boldsymbol{\Phi_N\Phi_N^T})^{-1}\in\mathbb{R}^{n\times n}
PN=(ΦNΦNT)−1∈Rn×n,则
θ
^
N
=
Y
N
Φ
N
T
P
N
\boldsymbol{\hat{\theta}_N}=\boldsymbol{Y_N}\boldsymbol{\Phi_N}^T \boldsymbol{P_N}
θ^N=YNΦNTPN
现在我们已经拥有了 N N N组 I / O I/O I/O数据,现在我们需要结合第 N + 1 N+1 N+1个新数据来修正我们的估计 θ ^ N \boldsymbol{\hat{\theta}_N} θ^N得到 θ ^ N + 1 \boldsymbol{\hat{\theta}_{N+1}} θ^N+1。这个过程就是一种递推,我们需要得到这种递推的通法。我们记之为:
θ ^ N + 1 = f ( θ ^ N , X N + 1 , y N + 1 ) \boldsymbol{\hat{\theta}_{N+1}}=f(\boldsymbol{\hat{\theta}_N},\boldsymbol{X}_{N+1},y_{N+1}) θ^N+1=f(θ^N,XN+1,yN+1)
开始递推
注意到:
θ
^
N
+
1
=
Y
N
+
1
Φ
N
+
1
T
P
N
+
1
\boldsymbol{\hat{\theta}}_{N+1}=\boldsymbol{Y}_{N+1}\boldsymbol{\Phi}_{N+1}^T\boldsymbol{P}_{N+1}
θ^N+1=YN+1ΦN+1TPN+1
先分析
P
N
+
1
∈
R
n
×
n
\boldsymbol{P}_{N+1}\in\mathbb{R}^{n\times n}
PN+1∈Rn×n
P
N
+
1
=
(
Φ
N
+
1
Φ
N
+
1
T
)
−
1
\boldsymbol{P}_{N+1}=(\boldsymbol{\Phi_{N+1}\Phi_{N+1}^T})^{-1}
PN+1=(ΦN+1ΦN+1T)−1
其中
Φ
N
+
1
Φ
N
+
1
T
=
[
X
1
X
2
⋯
X
N
+
1
]
[
X
1
T
X
2
T
⋮
X
N
+
1
T
]
\boldsymbol{\Phi_{N+1}\Phi_{N+1}^T}= \begin{bmatrix}X_1&X_2&\cdots &X_{N+1}\end{bmatrix} \begin{bmatrix}X_1^T\\X_2^T\\ \vdots \\X_{N+1}^T\end{bmatrix}
ΦN+1ΦN+1T=[X1X2⋯XN+1]⎣⎢⎢⎢⎡X1TX2T⋮XN+1T⎦⎥⎥⎥⎤
=
∑
i
=
1
N
+
1
X
i
X
i
T
=
∑
i
=
1
N
X
i
X
i
T
+
X
N
+
1
X
N
+
1
T
=
Φ
N
Φ
N
T
+
X
N
+
1
X
N
+
1
T
=\sum_{i=1}^{N+1}{X_iX_i^T}=\sum_{i=1}^{N}{X_iX_i^T}+X_{N+1}X_{N+1}^T= \boldsymbol{\Phi_{N}\Phi_{N}^T}+X_{N+1}X_{N+1}^T
=i=1∑N+1XiXiT=i=1∑NXiXiT+XN+1XN+1T=ΦNΦNT+XN+1XN+1T
故而
P
N
+
1
−
1
=
P
N
−
1
+
X
N
+
1
X
N
+
1
T
\boldsymbol{P}_{N+1}^{-1}= \boldsymbol{P}_{N}^{-1}+X_{N+1}X_{N+1}^T
PN+1−1=PN−1+XN+1XN+1T
⇒ { P N + 1 = ( P N − 1 + X N + 1 X N + 1 T ) − 1 P N = ( P N + 1 − 1 − X N + 1 X N + 1 T ) − 1 \Rightarrow\left\{\begin{array}{ll} \boldsymbol{P}_{N+1}=(\boldsymbol{P}_{N}^{-1}+X_{N+1}X_{N+1}^T)^{-1}\\ \boldsymbol{P}_{N}= (\boldsymbol{P}_{N+1}^{-1}-X_{N+1}X_{N+1}^T)^{-1} \end{array}\right. ⇒{PN+1=(PN−1+XN+1XN+1T)−1PN=(PN+1−1−XN+1XN+1T)−1
再分析
Y
N
Φ
N
T
∈
R
1
×
n
\boldsymbol{Y}_N\boldsymbol{\Phi}_N^T\in\mathbb{R}^{1\times n}
YNΦNT∈R1×n
Y
N
+
1
Φ
N
+
1
T
=
[
y
1
y
2
⋯
y
N
+
1
]
[
X
1
T
X
2
T
⋮
X
N
+
1
T
]
\boldsymbol{Y}_{N+1}\boldsymbol{\Phi}_{N+1}^T= \begin{bmatrix}y_1& y_2& \cdots & y_{N+1}\end{bmatrix} \begin{bmatrix}X_1^T\\ X_2^T\\ \vdots \\ X_{N+1}^T\end{bmatrix}
YN+1ΦN+1T=[y1y2⋯yN+1]⎣⎢⎢⎢⎡X1TX2T⋮XN+1T⎦⎥⎥⎥⎤
= ∑ i = 1 N + 1 y i X i T = ∑ i = 1 N y i X i T + y N + 1 ⋅ X N + 1 T = Y N Φ N T + y N + 1 ⋅ X N + 1 T =\sum_{i=1}^{N+1}{y_iX_i^T}=\sum_{i=1}^{N}{y_iX_i^T}+y_{N+1}\cdot X_{N+1}^T=\boldsymbol{Y}_{N}\boldsymbol{\Phi}_{N}^T+y_{N+1}\cdot X_{N+1}^T =i=1∑N+1yiXiT=i=1∑NyiXiT+yN+1⋅XN+1T=YNΦNT+yN+1⋅XN+1T
我们得到了关于
θ
^
N
+
1
\boldsymbol{\hat{\theta}}_{N+1}
θ^N+1的两个部分的递推式,我们将其代入到
θ
^
N
+
1
\boldsymbol{\hat{\theta}}_{N+1}
θ^N+1中:
θ
^
N
+
1
=
Y
N
+
1
Φ
N
+
1
T
P
N
+
1
=
(
Y
N
Φ
N
T
+
y
N
+
1
⋅
X
N
+
1
T
)
⋅
(
P
N
−
1
+
X
N
+
1
X
N
+
1
T
)
−
1
\boldsymbol{\hat{\theta}}_{N+1}=\boldsymbol{Y}_{N+1}\boldsymbol{\Phi}_{N+1}^T\boldsymbol{P}_{N+1}=(\boldsymbol{Y}_{N}\boldsymbol{\Phi}_{N}^T+y_{N+1}\cdot X_{N+1}^T)\cdot(\boldsymbol{P}_{N}^{-1}+X_{N+1}X_{N+1}^T)^{-1}
θ^N+1=YN+1ΦN+1TPN+1=(YNΦNT+yN+1⋅XN+1T)⋅(PN−1+XN+1XN+1T)−1
= Y N Φ N T ⋅ P N + 1 + y N + 1 ⋅ X N + 1 T ⋅ P N + 1 =\boldsymbol{Y}_{N}\boldsymbol{\Phi}_{N}^T\cdot\boldsymbol{P}_{N+1}+y_{N+1}\cdot X_{N+1}^T\cdot\boldsymbol{P}_{N+1} =YNΦNT⋅PN+1+yN+1⋅XN+1T⋅PN+1
又因为:
θ
^
N
=
Y
N
Φ
N
T
P
N
⇒
θ
^
N
P
N
−
1
=
Y
N
Φ
N
T
\boldsymbol{\hat{\theta}}_{N}=\boldsymbol{Y}_N\boldsymbol{\Phi}_N^T\boldsymbol{P}_{N} \Rightarrow\boldsymbol{\hat{\theta}}_{N}\boldsymbol{P}_N^{-1}=\boldsymbol{Y}_N\boldsymbol{\Phi}_N^T
θ^N=YNΦNTPN⇒θ^NPN−1=YNΦNT
所以
θ
^
N
+
1
=
θ
^
N
P
N
−
1
⋅
P
N
+
1
+
y
N
+
1
⋅
X
N
+
1
T
⋅
P
N
+
1
\boldsymbol{\hat{\theta}}_{N+1}= \boldsymbol{\hat{\theta}}_{N}\boldsymbol{P}_N^{-1}\cdot \boldsymbol{P}_{N+1} +y_{N+1}\cdot X_{N+1}^T\cdot \boldsymbol{P}_{N+1}
θ^N+1=θ^NPN−1⋅PN+1+yN+1⋅XN+1T⋅PN+1
= θ ^ N ⋅ ( P N + 1 − 1 − X N + 1 X N + 1 T ) P N + 1 + y N + 1 ⋅ X N + 1 T ⋅ P N + 1 =\boldsymbol{\hat{\theta}}_{N}\cdot (\boldsymbol{P}_{N+1}^{-1}-X_{N+1}X_{N+1}^T)\boldsymbol{P}_{N+1} +y_{N+1}\cdot X_{N+1}^T\cdot\boldsymbol{P}_{N+1} =θ^N⋅(PN+1−1−XN+1XN+1T)PN+1+yN+1⋅XN+1T⋅PN+1
= θ ^ N − θ ^ N X N + 1 X N + 1 T P N + 1 + y N + 1 ⋅ X N + 1 T ⋅ P N + 1 =\boldsymbol{\hat{\theta}}_{N}- \boldsymbol{\hat{\theta}}_{N}X_{N+1}X_{N+1}^T\boldsymbol{P}_{N+1} +y_{N+1}\cdot X_{N+1}^T\cdot\boldsymbol{P}_{N+1} =θ^N−θ^NXN+1XN+1TPN+1+yN+1⋅XN+1T⋅PN+1
= θ ^ N + ( y N + 1 − θ ^ N X N + 1 ) ⋅ X N + 1 T P N + 1 =\boldsymbol{\hat{\theta}}_{N}+ (y_{N+1}-\boldsymbol{\hat{\theta}}_{N}X_{N+1}) \cdot X_{N+1}^T\boldsymbol{P}_{N+1} =θ^N+(yN+1−θ^NXN+1)⋅XN+1TPN+1
我们令 K N + 1 = X N + 1 T P N + 1 \boldsymbol{K}_{N+1}=X_{N+1}^T\boldsymbol{P}_{N+1} KN+1=XN+1TPN+1, ε N + 1 = y N + 1 − θ ^ N X N + 1 \varepsilon_{N+1}=y_{N+1}-\boldsymbol{\hat{\theta}}_{N}X_{N+1} εN+1=yN+1−θ^NXN+1。
综上,我们得到了 θ ^ N + 1 = θ ^ N + \boldsymbol{\hat{\theta}}_{N+1}=\boldsymbol{\hat{\theta}}_{N}+ θ^N+1=θ^N+修正量 的形式:
- P N + 1 = ( P N − 1 + X N + 1 X N + 1 T ) − 1 \boldsymbol{P}_{N+1}=(\boldsymbol{P}_{N}^{-1}+X_{N+1}X_{N+1}^T)^{-1} PN+1=(PN−1+XN+1XN+1T)−1
- K N + 1 = X N + 1 T P N + 1 \boldsymbol{K}_{N+1}=X_{N+1}^T\boldsymbol{P}_{N+1} KN+1=XN+1TPN+1
- ε N + 1 = y N + 1 − θ ^ N X N + 1 \varepsilon_{N+1}=y_{N+1}-\boldsymbol{\hat{\theta}}_{N}X_{N+1} εN+1=yN+1−θ^NXN+1
- θ ^ N + 1 = θ ^ N + ε N + 1 K N + 1 \boldsymbol{\hat{\theta}}_{N+1}= \boldsymbol{\hat{\theta}}_{N}+\varepsilon_{N+1} \boldsymbol{K}_{N+1} θ^N+1=θ^N+εN+1KN+1
递推优化
我们发现,
P
N
+
1
=
(
P
N
−
1
+
X
N
+
1
X
N
+
1
T
)
−
1
\boldsymbol{P}_{N+1}=(\boldsymbol{P}_{N}^{-1}+X_{N+1}X_{N+1}^T)^{-1}
PN+1=(PN−1+XN+1XN+1T)−1需要求逆,很麻烦。我们引入公式:
[
A
+
B
C
D
]
−
1
=
A
−
1
−
A
−
1
B
[
C
−
1
+
D
A
−
1
B
]
−
1
D
A
−
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
P N + 1 = ( P N − 1 + X N + 1 X N + 1 T ) − 1 \boldsymbol{P}_{N+1}=(\boldsymbol{P}_{N}^{-1}+X_{N+1}X_{N+1}^T)^{-1} PN+1=(PN−1+XN+1XN+1T)−1
N o t e : A = P N − 1 , B = X N + 1 , C = 1 , D = X N + 1 T Note:A=\boldsymbol{P}_{N}^{-1}, B=\boldsymbol{X}_{N+1}, C=1,D=\boldsymbol{X}_{N+1}^T Note:A=PN−1,B=XN+1,C=1,D=XN+1T
得
P
N
+
1
=
P
N
−
P
N
X
N
+
1
X
N
+
1
T
P
N
1
+
X
N
+
1
T
P
N
X
N
+
1
\boldsymbol{P}_{N+1}=\boldsymbol{P}_{N}-\frac{ \boldsymbol{P}_{N}\boldsymbol{X}_{N+1} \boldsymbol{X}_{N+1}^{T}\boldsymbol{P}_{N}} {1+\boldsymbol{X}_{N+1}^T\boldsymbol{P}_{N}\boldsymbol{X}_{N+1}}
PN+1=PN−1+XN+1TPNXN+1PNXN+1XN+1TPN
结论
对于函数:
y
=
[
θ
1
θ
2
⋯
θ
n
]
[
x
1
x
2
⋮
x
n
]
=
θ
X
=
∑
i
=
1
n
θ
i
x
i
y= \begin{bmatrix} \theta_1& \theta_2& \cdots& \theta_n \end{bmatrix} \begin{bmatrix} x_1\\ x_2\\ \vdots\\ x_n \end{bmatrix} =\boldsymbol{\theta}X =\sum_{i=1}^n{\theta_i x_i}
y=[θ1θ2⋯θn]⎣⎢⎢⎢⎡x1x2⋮xn⎦⎥⎥⎥⎤=θX=i=1∑nθixi
的递推最小二乘估计:
- 1 P N + 1 = P N − P N X N + 1 X N + 1 T P N 1 + X N + 1 T P N X N + 1 \boldsymbol{P}_{N+1}=\boldsymbol{P}_{N}-\frac{ \boldsymbol{P}_{N}\boldsymbol{X}_{N+1} \boldsymbol{X}_{N+1}^{T}\boldsymbol{P}_{N}} {1+\boldsymbol{X}_{N+1}^T\boldsymbol{P}_{N}\boldsymbol{X}_{N+1}} PN+1=PN−1+XN+1TPNXN+1PNXN+1XN+1TPN
- K N + 1 = X N + 1 T P N + 1 \boldsymbol{K}_{N+1}=X_{N+1}^T\boldsymbol{P}_{N+1} KN+1=XN+1TPN+1
- ε N + 1 = y N + 1 − θ ^ N X N + 1 \varepsilon_{N+1}=y_{N+1}-\boldsymbol{\hat{\theta}}_{N}X_{N+1} εN+1=yN+1−θ^NXN+1
- θ ^ N + 1 = θ ^ N + ε N + 1 K N + 1 \boldsymbol{\hat{\theta}}_{N+1}= \boldsymbol{\hat{\theta}}_{N}+\varepsilon_{N+1} \boldsymbol{K}_{N+1} θ^N+1=θ^N+εN+1KN+1
我们需要的是一个估计初值
θ
^
0
\boldsymbol{\hat\theta_0}
θ^0和
P
0
\boldsymbol{P_0}
P0。下面给出常用的初值取值的方法
Matlab 示例
代码部分
为了进一步加深理解,这里附上一段Matlab实现的递推最小二乘代码。
%递推最小二乘示例代码
%完成时间:2020/11/14
%作者:lamphungry
%原创代码,供大家参考
clc;clear;
times=2; %重复4个周期
x1=0;x0=0;u1=0;u0=0;N=50+2; %定义状态参数的4个初值(取0),并定义状态参数的总个数
a1=1.5;a2=-0.7;b1=1;b2=0.5; %初始化系统参数(这是我们需要拟合求的参数))
f=@(x1,x0,u1,u0) a1*x1+a2*x0+b1*u1+b2*u0; %递推函数(这是我们需要拟合的线性多元方程)
x=zeros(1,times*N);x(1:2)=[x0 x1]; %预分配空间并初始化前两个值
%定义输入,使用随机0-1序列,总数为N-1
u=(idinput(N)'+1)/2;u(1:2)=[u0 u1];
u=repmat(u,1,times);
T=1;n=1:times*N;t=n*T;
%定义受噪声的输出和理论输出,引入高斯白噪声,\delta_v^2=0.01^2
delta_v=0.01;
z=zeros(1,times*N);z_ori=zeros(1,times*N);
z(1:2)=x(1:2)+delta_v*randn(1,2); %掺杂噪声
z_ori(1:2)=x(1:2); %无噪声
for i=3:times*N
x(i)=f(x(i-1),x(i-2),u(i-1),u(i-2));
z(i-2)=x(i)+delta_v*randn(1);
z_ori(i-2)=x(i);
end
figure('Name','噪声输出和理论输出');
plot(t,z,t,z_ori);
legend('受噪声的输出','理论输出');
n=4;%待估计的数值为a1,a2,b1,b2.共4个,输入为也为4个分量,共times*N组
in1=x(2:end-1);
in2=x(1:end-2);
in3=u(2:end-1);
in4=u(1:end-2);
P_ori=1e6*eye(4,4); %充分大的数字乘以单位矩阵
theta_ori=[0 0 0 0]; %初始参数估计值设为0
X_ori=[in1(1);in2(1);in3(1);in4(1)];
Phi_ori=[X_ori];
Err=zeros(1,times*N);
Theta=zeros(4,times*N);
for i=2:times*N-2
%求Phi_N
%Phi=zeros(4,i);
X=[in1(i);in2(i);in3(i);in4(i)];
Phi=[Phi_ori X];
%求P_{N+1}
P=P_ori-(P_ori*X*X'*P_ori)/(1+X'*P_ori*X);
%P=(P_ori^-1+X*X')^-1;
%求K_N
K=X'*P;
%求\varepsilon_N
varepsilon=z(i)-theta_ori*X;
Err(i)=varepsilon^2;
%求新值
theta=theta_ori+varepsilon*K;
Theta(:,i)=theta;
%进行下一轮
Phi_ori=Phi;
P_ori=P;
theta_ori=theta;
X_ori=X;
end
% theta
%绘图部分
ff=@(theta,x1,x0,u1,u0) theta*[x1;x0;u1;u0];
xx=zeros(1,times*N);xx(1:2)=[x0 x1];
zz_ori=zeros(1,times*N);
Ori_theta=[1.5 -0.7 1 0.5];
%利用拟合估计的参数再进行递推求值,得到我们的拟合系统
for i=3:times*N
xx(i)=ff(theta,xx(i-1),xx(i-2),u(i-1),u(i-2));
zz_ori(i-2)=xx(i);
end
hold on;plot(t,zz_ori);legend('受噪声的输出','理论输出','仿真结果');
figure('Name','输入');
plot(t,u);
figure('Name','误差值');
plot(t,Err);
figure('Name','参数');hold on;
plot(repmat([0;times*N],1,4),repmat(Ori_theta,2,1),'LineWidth',2);
plot(repmat(t',1,4),repmat(Theta',1,1));
Matlab结果:
可以发现,效果拟合非常好。