前言
最优估计理论中研究的最小二乘估计（LS）为线性最小二乘估计（LLS），包括古典最小二乘估计（CLS）、加权最小二乘估计（WLS）[1]和递推最小二乘估计（RLS）。本文将详细介绍递推最小二乘估计（RLS）。
线性参数估计问题描述
这里重复文章[1]的相关描述。设 X X X为 n n n维未知参数向量， Z k Z_{k} Zk为 k k k维观测向量，表示经过 k k k组实验观测得到的观测值向量，其中元素 z i z_{i} zi表示第i次观测实验得到的观测值，显然其是1维观测标量， V k V_{k} Vk为 k k k维观测噪声向量，其中元素 v i v_{i} vi表示第i次观测实验的观测噪声，显然其是1维噪声标量。一般情况下 k > n k > n k>n且希望 k k k比 n n n大得多。单次观测值为多维的情况将在后文讨论。观测实验依据的自变量为 θ \theta θ，则将观测量 z i z_{i} zi表示为关于 θ \theta θ的未知函数 f ( θ , X ) f(\theta,X) f(θ,X)：
z i = f ( θ , X ) = ∑ j = 1 n [ x j h i , j ( θ ) ] + v i = x 1 h i , 1 ( θ ) + x 2 h i , 2 ( θ ) + ⋯ + x n h i , n ( θ ) + v i \begin{align*} z_{i} &= f(\theta,X) \\ &= \sum_{j=1}^{n} \left [ x_{j}h_{i,j}(\theta) \right ]+ v_{i} \\ &= x_{1}h_{i,1}(\theta)+ x_{2}h_{i,2}(\theta) + \cdots + x_{n}h_{i,n}(\theta) + v_{i} \tag{1} \\ \end{align*} zi=f(θ,X)=j=1∑n[xjhi,j(θ)]+vi=x1hi,1(θ)+x2hi,2(θ)+⋯+xnhi,n(θ)+vi(1)
其中
X = [ x 1 x 2 ⋮ x n ] Z k = [ z 1 z 2 ⋮ z k ] V k = [ v 1 v 2 ⋮ v k ] \begin{align*} X = \begin{bmatrix} x_{1} \\ x_{2} \\ \vdots \\ x_{n} \end{bmatrix} Z_{k} = \begin{bmatrix} z_{1} \\ z_{2} \\ \vdots \\ z_{k} \end{bmatrix} V_{k} = \begin{bmatrix} v_{1} \\ v_{2} \\ \vdots \\ v_{k} \end{bmatrix} \end{align*} X=
x1x2⋮xn
Zk=
z1z2⋮zk
Vk=
v1v2⋮vk
式(1)中 h i , j ( θ ) h_{i,j}(\theta) hi,j(θ)表示第 i i i次观测第 j j j个基函数，常用为多项式、三角函数或自然指数函数形式：
h i , j ( θ ) = θ j − 1 h i , j ( θ ) = s i n ( j θ ) h i , j ( θ ) = e x p ( λ j θ ) \begin{align*} h_{i,j}(\theta) &= \theta ^{j-1} \\ h_{i,j}(\theta) &= sin(j\theta) \\ h_{i,j}(\theta) &= exp(\lambda_{j} \theta) \\ \end{align*} hi,j(θ)hi,j(θ)hi,j(θ)=θj−1=sin(jθ)=exp(λjθ)
其中， λ j \lambda_{j} λj为自然数指数参数。
当观测实验进行，上述基函数均可根据 θ \theta θ求得。令 h i = [ h i , 1 ( θ ) h i , 2 ( θ ) ⋯ h i , n ( θ ) ] h_{i} = \begin{bmatrix} h_{i,1}(\theta) & h_{i,2}(\theta) & \cdots & h_{i,n}(\theta) \\ \end{bmatrix} hi=[hi,1(θ)hi,2(θ)⋯hi,n(θ)]且为已知，其为 n n n维常向量，将式(1)改写为：
Z k = H k X + V k \begin{align*} Z_{k}= H_{k}X+ V_{k} \tag{2} \\ \end{align*} Zk=HkX+Vk(2)
其中， H H H为参数向量 X X X到观测向量 Z Z Z的 k × n k \times n k×n维转移矩阵：
H = [ h 1 h 2 ⋮ h 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 ( θ ) ] \begin{align*} H = \begin{bmatrix} h_{1} \\ h_{2} \\ \vdots \\ h_{k} \end{bmatrix} = \begin{bmatrix} h_{1,1}(\theta) & h_{1,2}(\theta) & \cdots & h_{1,n}(\theta) \\ h_{2,1}(\theta) & h_{2,2}(\theta) & \cdots & h_{2,n}(\theta) \\ \vdots & \vdots & \ddots & \vdots\\ h_{k,1}(\theta) & h_{k,2}(\theta) & \cdots & h_{k,n}(\theta) \end{bmatrix} \\ \end{align*} H=
h1h2⋮hk
=
h1,1(θ)h2,1(θ)⋮hk,1(θ)h1,2(θ)h2,2(θ)⋮hk,2(θ)⋯⋯⋱⋯h1,n(θ)h2,n(θ)⋮hk,n(θ)
显然，观测向量 Z Z Z与被估参数向量 X X X存在线性关系，依据最优准则求对 X X X的估计值 X ^ \hat{X} X^是一个线性参数估计问题，自然对应线性最小二乘估计（LLS）。
这里讨论下超定方程组的矛盾：当 k = n k = n k=n时，线性方程组有唯一精确解，但当 k > n k > n k>n，线性方程数大于未知被估参数向量的维度，线性方程组变成线性超定方程组，其解不唯一。最小二乘法的思想是需求统计意义上的近似解，使线性超定方程组中各方程能得到近似相等。
递推最小二乘估计（Recursive Least Squares Estimation, RLSE）
递推最小二乘估计（RLS） 没有在估计准则方面有新的提法，而是用新的观测信息修正现有估计值计算出新的估计值，实现最小二乘法的递推更新与校正。相比最小二乘估计和加权最小二乘估计的批量处理方法，借助递推滤波的思想，可以在计算空间上有显著的提高[2]。
以加权最小二乘估计中的最优加权估计为例，阐述递推实现过程。
估计准则为：加权残差平方和最小。根据式(3)代价函数改写如下：
J = E k ^ T W E k ^ = ( Z k − H k X k ^ ) T W k ( Z k − H X k ^ ) = ∑ i = 1 k w i e ^ i 2 = ∑ i = 1 k w i ( z i − h i X k ^ ) 2 = m i n \begin{align*} J = \hat{E_{k}}^{T}W\hat{E_{k}} &= (Z_{k}-H_{k}\hat{X_{k}})^{T}W_{k}(Z_{k}-H\hat{X_{k}}) = \sum_{i=1}^{k} w_{i}\hat{e}_{i}^{2} = \sum_{i=1}^{k}w_{i}(z_{i}-h_{i}\hat{X_{k}})^{2}=min \tag{3} \\ \end{align*} J=Ek^TWEk^=(Zk−HkXk^)TWk(Zk−HXk^)=i=1∑kwie^i2=i=1∑kwi(zi−hiXk^)2=min(3)
其中， e ^ i \hat{e}_{i} e^i为第 i i i次观测的残差（Residual Error）， E ^ k \hat{E}_{k} E^k为 k k k维残差向量有：
e ^ i = z i − h i X ^ k E ^ k = Z k − H k X ^ = [ e ^ 1 e ^ 2 ⋮ e ^ k ] \begin{align*} \hat{e}_{i} &= z_{i}-h_{i}\hat{X}_{k} \\ \hat{E}_{k} &= Z_{k}-H_{k}\hat{X} = \begin{bmatrix} \hat{e}_{1} \\ \hat{e}_{2} \\ \vdots \\ \hat{e}_{k} \end{bmatrix} \\ \end{align*} e^iE^k=zi−hiX^k=Zk−HkX^=
e^1e^2⋮e^k
观测噪声向量 V V V为白噪声，即 E [ V ] = 0 E[V] = 0 E[V]=0且方差矩阵为 R k R_{k} Rk：
R k = [ σ 1 2 0 ⋯ 0 0 σ 2 2 ⋯ 0 ⋮ ⋮ ⋱ ⋮ 0 0 ⋯ σ k 2 ] \begin{align*} R_{k} &= \begin{bmatrix} \sigma_{1}^{2} & 0& \cdots& 0\\ 0& \sigma_{2}^{2} & \cdots& 0 \\ \vdots& \vdots& \ddots & \vdots\\ 0& 0& \cdots& \sigma_{k}^{2} \end{bmatrix} \\ \end{align*} Rk=
σ120⋮00σ22⋮0⋯⋯⋱⋯00⋮σk2
W k W_{k} Wk为 k × k k\times k k×k阶对称正定加权矩阵，取 W k = R − 1 W_{k}=R^{-1} Wk=R−1时，加权最小二乘估计为最优加权最小二乘估计，也称马尔可夫（Markov）估计，是没有初始值条件下的线性最小方差无偏估计[1]。
W k = [ w 1 0 ⋯ 0 0 w 2 ⋯ 0 ⋮ ⋮ ⋱ ⋮ 0 0 ⋯ w k ] = [ 1 σ 1 2 0 ⋯ 0 0 1 σ 2 2 ⋯ 0 ⋮ ⋮ ⋱ ⋮ 0 0 ⋯ 1 σ n 2 ] \begin{align*} W_{k} &= \begin{bmatrix} w_{1} & 0& \cdots& 0\\ 0& w_{2} & \cdots& 0 \\ \vdots& \vdots& \ddots & \vdots\\ 0& 0& \cdots& w_{k} \end{bmatrix} = \begin{bmatrix} \frac{1}{\sigma_{1}^{2}} & 0& \cdots& 0\\ 0& \frac{1}{\sigma_{2}^{2}} & \cdots& 0 \\ \vdots& \vdots& \ddots & \vdots\\ 0& 0& \cdots& \frac{1}{\sigma_{n}^{2}} \end{bmatrix} \end{align*} Wk=
w10⋮00w2⋮0⋯⋯⋱⋯00⋮wk
=
σ1210⋮00σ221⋮0⋯⋯⋱⋯00⋮σn21
其估计向量 X ^ k \hat{X}_{k} X^k为:
X ^ k = ( H k T W k H k ) − 1 H k T W k Z k = ( H k T R k − 1 H k ) − 1 H k T R k − 1 Z k \begin{align*} \hat{X}_{k} &= (H_{k}^{T}W_{k}H_{k})^{-1}H_{k}^{T}W_{k}Z_{k} \tag{4} \\ &= (H_{k}^{T}R_{k}^{-1}H_{k})^{-1}H_{k}^{T}R_{k}^{-1}Z_{k} \tag{5} \\ \end{align*} X^k=(HkTWkHk)−1HkTWkZk=(HkTRk−1Hk)−1HkTRk−1Zk(4)(5)
在处理 k k k个观测值之后，又得到第 k + 1 k+1 k+1次观测值 z k + 1 z_{k+1} zk+1：
z k + 1 = ∑ j = 1 n [ x j h k + 1 , j ( θ ) ] + v k + 1 = x 1 h k + 1 , 1 ( θ ) + x 2 h k + 1 , 2 ( θ ) + ⋯ + x n h k + 1 , n ( θ ) + v k + 1 = h k + 1 X + v k + 1 \begin{align*} z_{k+1} &= \sum_{j=1}^{n} \left [ x_{j}h_{k+1,j}(\theta) \right ]+ v_{k+1} \\ &= x_{1}h_{k+1,1}(\theta)+ x_{2}h_{k+1,2}(\theta) + \cdots + x_{n}h_{k+1,n}(\theta) + v_{k+1} \\ &= h_{k+1}X+ v_{k+1} \tag{6} \\ \end{align*} zk+1=j=1∑n[xjhk+1,j(θ)]+vk+1=x1hk+1,1(θ)+x2hk+1,2(θ)+⋯+xnhk+1,n(θ)+vk+1=hk+1X+vk+1(6)
其中，令 h k + 1 = [ h k + 1 , 1 ( θ ) h k + 1 , 2 ( θ ) ⋯ h k + 1 , n ( θ ) ] h_{k+1} = \begin{bmatrix} h_{k+1,1}(\theta) & h_{k+1,2}(\theta) & \cdots & h_{k+1,n}(\theta) \\ \end{bmatrix} hk+1=[hk+1,1(θ)hk+1,2(θ)⋯hk+1,n(θ)]
将式(2)与(5)合并，得
Z k + 1 = H k + 1 X + V k + 1 \begin{align*} Z_{k+1}= H_{k+1}X+ V_{k+1} \tag{7} \\ \end{align*} Zk+1=Hk+1X+Vk+1(7)
其中
Z k + 1 = [ Z k z k + 1 ] H k + 1 = [ H k h k + 1 ] V k + 1 = [ V k v k + 1 ] \begin{align*} Z_{k+1} &= \begin{bmatrix} Z_{k} \\ z_{k+1} \end{bmatrix} H_{k+1}= \begin{bmatrix} H _{k} \\ h_{k+1} \end{bmatrix} V_{k+1}= \begin{bmatrix} V _{k} \\ v_{k+1} \end{bmatrix} \end{align*} Zk+1=[Zkzk+1]Hk+1=[Hkhk+1]Vk+1=[Vkvk+1]
由于加权矩阵 W k + 1 W_{k+1} Wk+1对称正定矩阵，则
W k + 1 = [ W k 0 0 w k + 1 ] \begin{align*} W_{k+1} &= \begin{bmatrix} W_{k}& 0 \\ 0& w_{k+1} \\ \end{bmatrix} \tag{8} \\ \end{align*} Wk+1=[Wk00wk+1](8)
得到估计向量 X ^ \hat{X} X^为：
X ^ k + 1 = ( H k + 1 T W k + 1 H k + 1 ) − 1 H k + 1 T W k + 1 Z k + 1 \begin{align*} \hat{X}_{k+1} &= (H_{k+1}^{T}W_{k+1}H_{k+1})^{-1}H_{k+1}^{T}W_{k+1}Z_{k+1} \tag{9} \\ \end{align*} X^k+1=(Hk+1TWk+1Hk+1)−1Hk+1TWk+1Zk+1(9)
( H k + 1 T W k + 1 H k + 1 ) − 1 = [ [ H k h k + 1 ] T [ W k 0 0 w k + 1 ] [ H k h k + 1 ] ] − 1 = ( H k T W k H k + h k + 1 T w k + 1 h k + 1 ) − 1 \begin{align*} (H_{k+1}^{T}W_{k+1}H_{k+1})^{-1} &= \left [ \begin{bmatrix} H_{k} &h_{k+1} \end{bmatrix}^{T} \begin{bmatrix} W_{k}& 0 \\ 0& w_{k+1} \\ \end{bmatrix} \begin{bmatrix} H_{k} &h_{k+1} \end{bmatrix} \right ]^{-1} \\ &= \left ( H_{k}^{T}W_{k}H_{k} + h_{k+1}^{T}w_{k+1}h_{k+1} \right )^{-1} \tag{10} \end{align*} (Hk+1TWk+1Hk+1)−1=[[Hkhk+1]T[Wk00wk+1][Hkhk+1]]−1=(HkTWkHk+hk+1Twk+1hk+1)−1(10)
令
P k = ( H k T W k H k ) − 1 P k − 1 = H k T W k H k P k + 1 = ( H k + 1 T W k + 1 H k + 1 ) − 1 = ( H k T W k H k + h k + 1 T w k + 1 h k + 1 ) − 1 = ( P k − 1 + h k + 1 T w k + 1 h k + 1 ) − 1 P k + 1 − 1 = P k − 1 + h k + 1 T w k + 1 h k + 1 \begin{align*} P_{k} &= (H_{k}^{T}W_{k}H_{k})^{-1} \tag{11} \\ P_{k}^{-1} &= H_{k}^{T}W_{k}H_{k} \tag{12} \\ P_{k+1} &= (H_{k+1}^{T}W_{k+1}H_{k+1})^{-1} \\ &= ( H_{k}^{T}W_{k}H_{k} + h_{k+1}^{T}w_{k+1}h_{k+1} )^{-1} \\ &= (P_{k}^{-1} + h_{k+1}^{T}w_{k+1}h_{k+1} )^{-1} \tag{13} \\ P_{k+1}^{-1} &= P_{k}^{-1} + h_{k+1}^{T}w_{k+1}h_{k+1} \tag{14} \\ \end{align*} PkPk−1Pk+1Pk+1−1=(HkTWkHk)−1=HkTWkHk=(Hk+1TWk+1Hk+1)−1=(HkTWkHk+hk+1Twk+1hk+1)−1=(Pk−1+hk+1Twk+1hk+1)−1=Pk−1+hk+1Twk+1hk+1(11)(12)(13)(14)
根据矩阵求逆引理，得
P k + 1 = P k − P k h k + 1 T ( h k + 1 P k h k + 1 T + w k + 1 − 1 ) − 1 h k + 1 P k \begin{align*} P_{k+1} &= P_{k} - P_{k}h_{k+1}^{T} (h_{k+1}P_{k}h_{k+1}^{T} + w_{k+1}^{-1})^{-1} h_{k+1}P_{k} \tag{15} \\ \end{align*} Pk+1=Pk−Pkhk+1T(hk+1Pkhk+1T+wk+1−1)−1hk+1Pk(15)
由式(4)(9)(12)(14)，得出估计递推公式：
X ^ k + 1 = P k + 1 H k + 1 T W k + 1 Z k + 1 = P k + 1 ( H k T W k Z k + h k + 1 T w k + 1 z k + 1 ) = P k + 1 ( P k − 1 X ^ k + h k + 1 T w k + 1 z k + 1 ) = P k + 1 ( ( P k + 1 − 1 − h k + 1 T w k + 1 h k + 1 ) X ^ k + h k + 1 T w k + 1 z k + 1 ) = X ^ k + P k + 1 h k + 1 T w k + 1 ( z k + 1 − h k + 1 X ^ k ) = X ^ k + K k + 1 ( z k + 1 − h k + 1 X ^ k ) \begin{align*} \hat{X}_{k+1} &= P_{k+1}H_{k+1}^{T}W_{k+1}Z_{k+1} \\ &= P_{k+1}(H_{k}^{T}W_{k}Z_{k} + h_{k+1}^{T}w_{k+1}z_{k+1} ) \\ &= P_{k+1}(P_{k}^{-1} \hat{X}_{k} + h_{k+1}^{T}w_{k+1}z_{k+1} ) \\ &= P_{k+1}((P_{k+1}^{-1}-h_{k+1}^{T}w_{k+1}h_{k+1}) \hat{X}_{k} + h_{k+1}^{T}w_{k+1}z_{k+1} ) \\ &= \hat{X}_{k} + P_{k+1}h_{k+1}^{T}w_{k+1} (z_{k+1} - h_{k+1} \hat{X}_{k}) \\ &= \hat{X}_{k} + K_{k+1} (z_{k+1} - h_{k+1} \hat{X}_{k}) \tag{16} \\ \end{align*} X^k+1=Pk+1Hk+1TWk+1Zk+1=Pk+1(HkTWkZk+hk+1Twk+1zk+1)=Pk+1(Pk−1X^k+hk+1Twk+1zk+1)=Pk+1((Pk+1−1−hk+1Twk+1hk+1)X^k+hk+1Twk+1zk+1)=X^k+Pk+1hk+1Twk+1(zk+1−hk+1X^k)=X^k+Kk+1(zk+1−hk+1X^k)(16)
其中， K k + 1 K _{k+1} Kk+1为递推增益
K k + 1 = P k + 1 h k + 1 T w k + 1 \begin{align*} K_{k+1} &= P_{k+1}h_{k+1}^{T}w_{k+1} \tag{17} \\ \end{align*} Kk+1=Pk+1hk+1Twk+1(17)
对于加权最小二乘估计，式(15)(16)(17)即为其递推最小二乘估计（RLS）的递推公式：
P k + 1 = P k − P k h k + 1 T ( h k + 1 P k h k + 1 T + w k + 1 − 1 ) − 1 h k + 1 P k K k + 1 = P k + 1 h k + 1 T w k + 1 X ^ k + 1 = X ^ k + K k + 1 ( z k + 1 − h k + 1 X ^ k ) \begin{align*} P_{k+1} &= P_{k} - P_{k}h_{k+1}^{T} (h_{k+1}P_{k}h_{k+1}^{T} + w_{k+1}^{-1})^{-1} h_{k+1}P_{k} \tag{15} \\ K_{k+1} &= P_{k+1}h_{k+1}^{T}w_{k+1} \tag{17} \\ \hat{X}_{k+1} &= \hat{X}_{k} + K_{k+1} (z_{k+1} - h_{k+1} \hat{X}_{k}) \tag{16} \\ \end{align*} Pk+1Kk+1X^k+1=Pk−Pkhk+1T(hk+1Pkhk+1T+wk+1−1)−1hk+1Pk=Pk+1hk+1Twk+1=X^k+Kk+1(zk+1−hk+1X^k)(15)(17)(16)
对于最优加权最小二乘估计，将 w k + 1 = r k + 1 − 1 w_{k+1}=r_{k+1}^{-1} wk+1=rk+1−1代入上式， r k + 1 − 1 r_{k+1}^{-1} rk+1−1为第 k + 1 k+1 k+1次观测噪声 v k + 1 v_{k+1} vk+1的方差，到其递推最小二乘估计（RLS）的递推公式：
P k + 1 = P k − P k h k + 1 T ( h k + 1 P k h k + 1 T + r k + 1 ) − 1 h k + 1 P k K k + 1 = P k + 1 h k + 1 T r k + 1 − 1 X ^ k + 1 = X ^ k + K k + 1 ( z k + 1 − h k + 1 X ^ k ) \begin{align*} P_{k+1} &= P_{k} - P_{k}h_{k+1}^{T} (h_{k+1}P_{k}h_{k+1}^{T} + r_{k+1})^{-1} h_{k+1}P_{k} \tag{18} \\ K_{k+1} &= P_{k+1}h_{k+1}^{T}r_{k+1}^{-1} \tag{19} \\ \hat{X}_{k+1} &= \hat{X}_{k} + K_{k+1} (z_{k+1} - h_{k+1} \hat{X}_{k}) \tag{16}\\ \end{align*} Pk+1Kk+1X^k+1=Pk−Pkhk+1T(hk+1Pkhk+1T+rk+1)−1hk+1Pk=Pk+1hk+1Trk+1−1=X^k+Kk+1(zk+1−hk+1X^k)(18)(19)(16)
关于初值 X ^ 0 \hat{X}_{0} X^0与 P 0 P_{0} P0有两种方法确定[2]：
- 批量计算前n次观测实验数据得到；
- 假设 X ^ 0 = 0 \hat{X}_{0}=0 X^0=0与 P 0 = C T C I P_{0}=C^{T}CI P0=CTCI， C C C为一个充分大的正数。
后者比前者不需要计算 n × n n \times n n×n阶逆矩阵，但是要经过一定次数迭代后才能得到满意估计。
多维观测
假设第 i i i组实验的观测值 Z i Z_{i} Zi为 m m m维观测向量，对应的观测噪声 V i V_{i} Vi为 m m m维观测向量，观测转移矩阵 H i H_{i} Hi为 m × n m \times n m×n阶矩阵，有：
Z i = [ z i , 1 z i , 2 ⋮ z i , m ] , V i = [ v i , 1 v i , 2 ⋮ v i , m ] , V a r ( V i ) = R i , H i = [ h i , 1 , 1 ( θ ) h i , 1 , 2 ( θ ) ⋯ h i , 1 , n ( θ ) h i , 2 , 1 ( θ ) h i , 2 , 2 ( θ ) ⋯ h i , 2 , n ( θ ) ⋮ ⋮ ⋱ ⋮ h m , 2 , 1 ( θ ) h m , 2 , 2 ( θ ) ⋯ h m , 2 , n ( θ ) ] \begin{align*} Z_{i} &= \begin{bmatrix} z_{i,1} \\ z_{i,2} \\ \vdots\\ z_{i,m} \end{bmatrix} ,V_{i} = \begin{bmatrix} v_{i,1} \\ v_{i,2} \\ \vdots\\ v_{i,m} \end{bmatrix}, Var(V_i) = R_{i}, H_{i} = \begin{bmatrix} h_{i,1,1}(\theta)& h_{i,1,2}(\theta)& \cdots& h_{i,1,n}(\theta) \\ h_{i,2,1}(\theta)& h_{i,2,2}(\theta)& \cdots& h_{i,2,n}(\theta) \\ \vdots & \vdots& \ddots& \vdots\\ h_{m,2,1}(\theta)& h_{m,2,2}(\theta)& \cdots& h_{m,2,n}(\theta) \end{bmatrix} \end{align*} Zi=
zi,1zi,2⋮zi,m
,Vi=
vi,1vi,2⋮vi,m
,Var(Vi)=Ri,Hi=
hi,1,1(θ)hi,2,1(θ)⋮hm,2,1(θ)hi,1,2(θ)hi,2,2(θ)⋮hm,2,2(θ)⋯⋯⋱⋯hi,1,n(θ)hi,2,n(θ)⋮hm,2,n(θ)
对于第 i i i组实验，观测值向量为：
Z i = H i X + V i \begin{align*} Z_{i} = H_{i}X+V_{i} \tag{20} \end{align*} Zi=HiX+Vi(20)
对于全部 k k k组实验，观测值向量为：
Z ˉ k = H ˉ k X + V ˉ k (21) \bar{Z}_{k} = \bar{H}_{k}X + \bar{V}_{k} \tag{21} Zˉk=HˉkX+Vˉk(21)
其中
Z ˉ k = [ Z 1 Z 2 ⋮ Z k ] , H ˉ k = [ H 1 H 2 ⋮ H k ] , V ˉ k = [ V 1 V 2 ⋮ V k ] \begin{align*} \bar{Z}_{k} &= \begin{bmatrix} Z_{1} \\ Z_{2} \\ \vdots\\ Z_{k} \end{bmatrix} , \bar{H}_{k} &= \begin{bmatrix} H_{1} \\ H_{2} \\ \vdots\\ H_{k} \end{bmatrix} , \bar{V}_{k} &= \begin{bmatrix} V_{1} \\ V_{2} \\ \vdots\\ V_{k} \end{bmatrix} \end{align*} Zˉk=
Z1Z2⋮Zk
,Hˉk=
H1H2⋮Hk
,Vˉk=
V1V2⋮Vk
根据加权最小二乘估计（WLS）准则：加权残差平方和最小。其代价函数为：
J = ( Z ˉ k − H ˉ k X ^ ) T W ˉ k ( Z ˉ k − H ˉ k X ^ ) \begin{align*} J &= (\bar{Z}_{k} - \bar{H}_{k} \hat{X})^{T} \bar{W}_{k} (\bar{Z}_{k} - \bar{H}_{k} \hat{X}) \tag{22} \end{align*} J=(Zˉk−HˉkX^)TWˉk(Zˉk−HˉkX^)(22)
式中， W ˉ k \bar{W}_{k} Wˉk为对称正定加权矩阵。
可得估计值向量 X ^ k \hat{X}_{k} X^k为：
X ^ k = ( H ˉ k T W ˉ k H ˉ k ) − 1 H ˉ k T W ˉ k Z ˉ k \begin{align*} \hat{X}_{k} &= (\bar{H}_{k}^{T} \bar{W}_{k} \bar{H}_{k} ) ^{-1} \bar{H}_{k}^{T}\bar{W}_{k} \bar{Z}_{k} \tag{23} \end{align*} X^k=(HˉkTWˉkHˉk)−1HˉkTWˉkZˉk(23)
现得到第 k + 1 k+1 k+1次观测值为:
Z k + 1 = H k + 1 X + V k + 1 \begin{align*} Z_{k+1} = H_{k+1}X + V_{k+1} \tag{24} \end{align*} Zk+1=Hk+1X+Vk+1(24)
设第 k + 1 k+1 k+1次加权矩阵为：
W ˉ k + 1 = [ W ˉ k 0 0 W k + 1 ] \begin{align*} \bar{W}_{k+1} &= \begin{bmatrix} \bar{W}_{k}& 0 \\ 0& W_{k+1} \\ \end{bmatrix} \end{align*} Wˉk+1=[Wˉk00Wk+1]
根据式(11)(12)(13)(14)，得
P ˉ k = ( H ˉ k T W ˉ k H ˉ k ) − 1 P ˉ k − 1 = H ˉ k T W ˉ k H ˉ k P ˉ k + 1 = ( P ˉ k − 1 + H k + 1 T W k + 1 H k + 1 ) − 1 P ˉ k + 1 − 1 = P ˉ k − 1 + H k + 1 T W k + 1 H k + 1 \begin{align*} \bar{P}_{k} &= (\bar{H}_{k}^{T}\bar{W}_{k}\bar{H}_{k})^{-1} \tag{25} \\ \bar{P}_{k}^{-1} &= \bar{H}_{k}^{T}\bar{W}_{k}\bar{H}_{k} \tag{26} \\ \bar{P}_{k+1} &= (\bar{P}_{k}^{-1} + H_{k+1}^{T}W_{k+1}H_{k+1} )^{-1} \tag{27} \\ \bar{P}_{k+1}^{-1} &= \bar{P}_{k}^{-1} + H_{k+1}^{T}W_{k+1}H_{k+1} \tag{28} \\ \end{align*} PˉkPˉk−1Pˉk+1Pˉk+1−1=(HˉkTWˉkHˉk)−1=HˉkTWˉkHˉk=(Pˉk−1+Hk+1TWk+1Hk+1)−1=Pˉk−1+Hk+1TWk+1Hk+1(25)(26)(27)(28)
根据式(15)(16)(17)，其加权最小二乘估计的递推最小二乘估计（RLS）的递推公式：
P ˉ k + 1 = P ˉ k − P ˉ k H k + 1 T ( H k + 1 P ˉ k H k + 1 T + W k + 1 − 1 ) − 1 H k + 1 P ˉ k K k + 1 = P ˉ k + 1 H k + 1 T W k + 1 X ^ k + 1 = X ^ k + K k + 1 ( Z k + 1 − H k + 1 X ^ k ) \begin{align*} \bar{P}_{k+1} &= \bar{P}_{k} - \bar{P}_{k}H_{k+1}^{T} (H_{k+1}\bar{P}_{k}H_{k+1}^{T} + W_{k+1}^{-1})^{-1} H_{k+1}\bar{P}_{k} \tag{29} \\ K_{k+1} &= \bar{P}_{k+1}H_{k+1}^{T}W_{k+1} \tag{30} \\ \hat{X}_{k+1} &= \hat{X}_{k} + K_{k+1} (Z_{k+1} - H_{k+1} \hat{X}_{k}) \tag{31} \\ \end{align*} Pˉk+1Kk+1X^k+1=Pˉk−PˉkHk+1T(Hk+1PˉkHk+1T+Wk+1−1)−1Hk+1Pˉk=Pˉk+1Hk+1TWk+1=X^k+Kk+1(Zk+1−Hk+1X^k)(29)(30)(31)
如果取加权矩阵：
W ˉ k = R ˉ k − 1 = [ R 1 0 ⋯ 0 0 R 2 ⋯ 0 ⋮ ⋮ ⋱ ⋮ 0 0 ⋯ R k ] W k + 1 = R k + 1 − 1 = V a r ( V k + 1 ) W ˉ k + 1 = [ W ˉ k 0 0 W k + 1 ] = [ R ˉ k − 1 0 0 R k + 1 − 1 ] \begin{align*} \bar{W}_{k} &= \bar{R}_{k}^{-1} = \begin{bmatrix} R_{1}& 0& \cdots& 0 \\ 0& R_{2}& \cdots& 0 \\ \vdots& \vdots& \ddots& \vdots \\ 0& 0& \cdots& R_{k} \\ \end{bmatrix} \tag{32} \\ W_{k+1} &= R_{k+1}^{-1} = Var(V_{k+1}) \tag{33} \\ \bar{W}_{k+1} &= \begin{bmatrix} \bar{W}_{k}& 0 \\ 0& W_{k+1} \\ \end{bmatrix} = \begin{bmatrix} \bar{R}_{k}^{-1}& 0 \\ 0& R_{k+1}^{-1} \\ \end{bmatrix} \tag{34} \end{align*} WˉkWk+1Wˉk+1=Rˉk−1=
R10⋮00R2⋮0⋯⋯⋱⋯00⋮Rk
=Rk+1−1=Var(Vk+1)=[Wˉk00Wk+1]=[Rˉk−100Rk+1−1](32)(33)(34)
根据式(29-34)，其最优加权最小二乘估计的递推最小二乘估计（RLS）递推公式：
P ˉ k + 1 = P ˉ k − P ˉ k H k + 1 T ( H k + 1 P ˉ k H k + 1 T + R k + 1 ) − 1 H k + 1 P ˉ k K k + 1 = P ˉ k + 1 H k + 1 T R k + 1 − 1 X ^ k + 1 = X ^ k + K k + 1 ( Z k + 1 − H k + 1 X ^ k ) \begin{align*} \bar{P}_{k+1} &= \bar{P}_{k} - \bar{P}_{k}H_{k+1}^{T} (H_{k+1}\bar{P}_{k}H_{k+1}^{T} + R_{k+1})^{-1} H_{k+1}\bar{P}_{k} \tag{35} \\ K_{k+1} &= \bar{P}_{k+1}H_{k+1}^{T}R_{k+1}^{-1} \tag{36} \\ \hat{X}_{k+1} &= \hat{X}_{k} + K_{k+1} (Z_{k+1} - H_{k+1} \hat{X}_{k}) \tag{37} \\ \end{align*} Pˉk+1Kk+1X^k+1=Pˉk−PˉkHk+1T(Hk+1PˉkHk+1T+Rk+1)−1Hk+1Pˉk=Pˉk+1Hk+1TRk+1−1=X^k+Kk+1(Zk+1−Hk+1X^k)(35)(36)(37)
关于初值 X ^ 0 \hat{X}_{0} X^0与 P ˉ 0 \bar{P}_{0} Pˉ0有两种方法确定：
- 批量计算前n次观测实验数据得到；
- 假设 X ^ 0 = 0 \hat{X}_{0}=0 X^0=0与 P ˉ = C T C I \bar{P}=C^{T}CI Pˉ=CTCI， C C C为一个充分大的正数。
后者比前者不需要计算 m n × m n mn \times mn mn×mn阶逆矩阵，但是要经过一定次数迭代后才能得到满意估计。
衰减加权
前文表述的加权矩阵将过去和现在的数据视为同等重要。但有些情况被估参数向量可能会随着时间而缓慢变化， 此时应该重视当前数据而逐渐衰减过去数据的影响，只需如下选择加权矩阵：
W k = [ e − ( k − 1 ) a 0 ⋯ 0 0 0 e − ( k − 2 ) a ⋯ 0 0 ⋮ ⋮ ⋱ ⋮ ⋮ 0 0 ⋯ e − a 0 0 0 ⋯ 0 1 ] = [ W k − 1 e − a 0 0 1 ] \begin{align*} W_{k} &= \begin{bmatrix} e^{-(k-1)a}& 0& \cdots& 0& 0 \\ 0& e^{-(k-2)a}& \cdots& 0& 0 \\ \vdots& \vdots& \ddots& \vdots& \vdots \\ 0& 0& \cdots& e^{-a}& 0 \\ 0& 0& \cdots& 0& 1 \end{bmatrix} = \begin{bmatrix} W_{k-1}e^{-a}& 0 \\ 0& 1 \end{bmatrix} \tag{38} \\ \end{align*} Wk=
e−(k−1)a0⋮000e−(k−2)a⋮00⋯⋯⋱⋯⋯00⋮e−a000⋮01
=[Wk−1e−a001](38)
其中， a a a为一个正数，称为衰减因子。
参考文献
[1] 最优估计准则与方法（5）加权最小二乘估计(WLS)_学习笔记
https://blog.csdn.net/jimmychao1982/article/details/149656746
[2] 《最优估计理论》，周凤歧，2009，高等教育出版社。
[3] 《最优估计理论》，刘胜，张红梅著，2011，科学出版社。