最小二乘估计（Least squares estimation）
最小二乘估计面向的问题
现有一色环模糊的电阻，不知道其真实电阻值，但是手头有一个万用表。由于测量方法和万用表精度问题，测量误差不可避免。这就需要我们从具有加性噪声的量测中估计出电阻的真值。电阻真值是个未知的恒定标量。为了使问题描述更具有普遍性，将电阻真值看作一恒定向量，也就是说，待估计量的数值是不随着时间改变的。
下面用数学术语来描述这个估计问题，假设 x x x是一个 n n n维待估计的未知恒定向量，现在有 k k k个量测数据 y 1 , y 2 , . . . , y k y_1,y_2,...,y_k y1,y2,...,yk，如何得到真值 x x x的最优估计 x ^ \hat{x} x^，便是最小二乘估计解决的问题。假设每个量测数据是真值的线性组合并加上噪声 y 1 = H 11 x 1 + H 12 x 2 + . . . + H 1 n x n + ν 1 ⋮ y k = H k 1 x 1 + H k 2 x 2 + . . . + H k n x n + ν k y_{1}=H_{11}x_{1}+H_{12}x_{2}+...+H_{1n}x_{n}+\nu_{1}\\ \vdots\\ y_{k}=H_{k1}x_{1}+H_{k2}x_{2}+...+H_{kn}x_{n}+\nu_{k} y1=H11x1+H12x2+...+H1nxn+ν1⋮yk=Hk1x1+Hk2x2+...+Hknxn+νk上式可以看做由 k k k个方程组解得 n n n个未知数，方程有且只有唯一解则满足 k ≥ n k\ge n k≥n，量测数据是线性独立的，将其写成矩阵的形式为 y = H x + ν y=Hx+\nu y=Hx+ν待估量 x x x和量测 y y y之间的观测矩阵 H H H假定是已知的。
最小二乘估计（Least squares estimation）
最优准则的选取不同决定了不同的估计方法。其中最为直观的一种方法为：最优估计 x ^ \hat{x} x^可以使所有量测的误差总体达到最小，这便是最小二乘估计采取的最优准则。第 i i i个量测残差为 ε i = y i − H i x ^ \varepsilon_{i}=y_{i}-H_{i}\hat{x} εi=yi−Hix^所有量测误差之和为 J = ε 1 2 + ε 2 2 + . . . + ε k 2 = ε T ε J=\varepsilon_{1}^2+\varepsilon_{2}^2+...+\varepsilon_{k}^2=\varepsilon^{T}\varepsilon J=ε12+ε22+...+εk2=εTε其中 ε = y − H x ^ \varepsilon=y-H\hat{x} ε=y−Hx^，将 J J J改写为 J = ( y − H x ^ ) T ( y − H x ^ ) = y T y − y T H x ^ − x ^ T H T y + x ^ T H T H x ^ J=(y-H\hat{x})^{T}(y-H\hat{x})=y^Ty-y^TH\hat{x}-\hat{x}^TH^Ty+\hat{x}^TH^TH\hat{x} J=(y−Hx^)T(y−Hx^)=yTy−yTHx^−x^THTy+x^THTHx^让 J J J对 x ^ \hat{x} x^的一阶偏导等于0，找到使 J J J最小的值 x ^ \hat{x} x^ ∂ J ∂ x ^ = − y T H − y T H + 2 x ^ T H T H = 0 \dfrac{\partial J}{\partial \hat{x}}=-y^TH-y^TH+2\hat{x}^TH^TH=0 ∂x^∂J=−yTH−yTH+2x^THTH=0得到估计值 x ^ = ( H T H ) − 1 H T y \hat{x}=(H^TH)^{-1}H^Ty x^=(HTH)−1HTy并且二阶导数为 ∂ 2 J ∂ x ^ 2 = 2 H T H \dfrac{\partial^2J}{\partial\hat{x}^2}=2H^TH ∂x^2∂2J=2HTH，这是一个正半定阵，那么 x ^ \hat{x} x^为最小值点。
加权最小二乘估计（Weighted least squares estimation）
让我们回到最初电阻测量的问题，不同的是，现在有多个价位的万用表可供使用，不同价位的万用表测量精度不同，表现在噪声的方差大小不一。这使得不同精度万用表测得数据的可信度也是不一样的，精度差的万用表测量的可信度较低，但无论精度如何，量测数据都不应丢弃，因为即使可信度较低的量测仍包含少量信息。
假定每次测量的噪声是零均值且独立的，噪声协方差矩阵为 R = E ( ν ν T ) = ( σ 1 2 ⋯ 0 ⋮ ⋱ ⋮ 0 ⋯ σ k 2 ) R=E(\nu\nu^T)=\begin{pmatrix} \sigma_1^2&\cdots&0\\ \vdots&\ddots&\vdots\\ 0&\cdots& \sigma_k^2 \end{pmatrix} R=E(ννT)=⎝⎜⎛σ12⋮0⋯⋱⋯0⋮σk2⎠⎟⎞其中 σ i 2 = E [ ν i 2 ] \sigma_i^2=E[\nu_i^2] σi2=E[νi2]，方差越大可信度越低。因此，对量测方差进行加权，加权后的总体误差为 J = ε 1 2 σ 1 2 + ε 2 2 σ 2 2 + . . . + ε k 2 σ k 2 = ε T R − 1 ε J=\dfrac{\varepsilon_{1}^2}{\sigma_1^2}+\dfrac{\varepsilon_{2}^2}{\sigma_2^2}+...+\dfrac{\varepsilon_{k}^2}{\sigma_k^2}=\varepsilon^{T}R^{-1}\varepsilon J=σ12ε12+σ22ε22+...+σk2εk2=εTR−1ε将 ε = y − H x ^ \varepsilon=y-H\hat{x} ε=y−Hx^代入展开得到 J = y T R − 1 y − y T R − 1 H x ^ − x ^ T H T R − 1 y + x ^ T H T R − 1 H x ^ J=y^TR^{-1}y-y^TR^{-1}H\hat{x}-\hat{x}^TH^TR^{-1}y+\hat{x}^TH^TR^{-1}H\hat{x} J=yTR−1y−yTR−1Hx^−x^THTR−1y+x^THTR−1Hx^对 x ^ \hat{x} x^求一阶导为 ∂ J ∂ x ^ = − y T R − 1 H − y T R − 1 H + 2 x ^ T H T R − 1 H = 0 \dfrac{\partial J}{\partial \hat{x}}=-y^TR^{-1}H-y^TR^{-1}H+2\hat{x}^TH^TR^{-1}H=0 ∂x^∂J=−yTR−1H−yTR−1H+2x^THTR−1H=0得到估计值 x ^ = ( H T R − 1 H ) − 1 H T R − 1 y \hat{x}=(H^TR^{-1}H)^{-1}H^TR^{-1}y x^=(HTR−1H)−1HTR−1y注意量测噪声存在逆矩阵， R R R是非奇异矩阵，也就是说加权最小二乘要求每个测量数据被噪声干扰。 J J J对 x ^ \hat{x} x^的二阶偏导为 2 H T R − 1 H 2H^TR^{-1}H 2HTR−1H，它是一个正半定阵。
递推最小二乘估计（Recursive least squares estimation）
递推最小二乘估计
上述最小二乘估计是建立在已获得全部量测数据的基础之上，若想在测量数据的同时，实现对真值的最小二乘估计则需要用到递推形式。已知 k − 1 k-1 k−1时刻的估计值 x ^ k − 1 \hat{x}_{k-1} x^k−1和 k k k时刻的量测 y k y_k yk，实现估计 x ^ k \hat{x}_k x^k。 k k k时刻的观测数据获得根据 y k = H k x + ν k y_k=H_kx+\nu_k yk=Hkx+νk k k k时刻估计和 x ^ k − 1 \hat{x}_{k-1} x^k−1以及 y k y_k yk的关系为 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+Kk(yk−Hkx^k−1)括号里被称作修正项， K k K_k Kk是修正增益，首先考虑估计的无偏性（估计的无偏性是指估计值 x ^ \hat{x} x^的期望等于真值 x x x，i.e. E [ x − x ^ ] = 0 E[x-\hat{x}]=0 E[x−x^]=0） E [ x − x ^ k ] = E [ x − x ^ k − 1 − K k ( H k x + ν k − H k x ^ k − 1 ) ] = E [ ε k − 1 − K k H k ε k − 1 − K k ν k ] = ( I − K k H k ) E [ ε k − 1 ] − K k E [ ν k ] \begin{aligned}E[x-\hat{x}_k]&=E[x-\hat{x}_{k-1}-K_{k}(H_kx+\nu_k-H_{k}\hat{x}_{k-1})]\\&=E[\varepsilon_{k-1}-K_kH_k\varepsilon_{k-1}-K_k\nu_k]\\&=(I-K_kH_k)E[\varepsilon_{k-1}]-K_kE[\nu_k]\end{aligned} E[x−x^k]=E[x−x^k−1−Kk(Hkx+νk−Hkx^k−1)]=E[εk−1−KkHkεk−1−Kkνk]=(I−KkHk)E[εk−1]−KkE[νk]如果 E [ ν k ] = 0 E[\nu_k]=0 E[νk]=0和 E [ ε k − 1 ] = 0 E[\varepsilon_{k-1}]=0 E[εk−1]=0，则 E [ ε k ] = 0 E[\varepsilon_{k}]=0 E[εk]=0。若初值设置等于真值，那么所有时刻的估计值都等于真值，并且不受增益 K k K_k Kk的影响。
最优准则
递归最小二乘估计的最优准则为，使 k k k时刻估计误差方差之和最小，数学表述为 J k = E [ ( x 1 − x ^ 1 ) 2 ] + E [ ( x 2 − x ^ 2 ) 2 ] + ⋯ + E [ ( x n − x ^ n ) 2 ] = E [ ε x 1 , k 2 + ε x 2 , k 2 + ⋯ + ε x n , k 2 ] = E [ Tr ( ε k ε k T ) ] = Tr P k \begin{aligned}J_k&=E[(x_1-\hat{x}_1)^2]+E[(x_2-\hat{x}_2)^2]+\cdots+E[(x_n-\hat{x}_n)^2]\\&=E[\varepsilon_{x1,k}^2+\varepsilon_{x2,k}^2+\cdots+\varepsilon_{xn,k}^2]\\&=E[\text{Tr}(\varepsilon_{k}\varepsilon_{k}^T)]\\&=\text{Tr}P_{k}\end{aligned} Jk=E[(x1−x^1)2]+E[(x2−x^2)2]+⋯+E[(xn−x^n)2]=E[εx1,k2+εx2,k2+⋯+εxn,k2]=E[Tr(εkεkT)]=TrPk将 P k P_k Pk展开得到 P k = E [ ε k ε k T ] = ( I − K k H k ) E [ ε k − 1 ε k − 1 T ] ( I − K k H k ) T + K k E [ ν k ν k T ] K k T = ( I − K k H k ) P k − 1 ( I − K k H k ) T + K k R k K k T \begin{aligned} P_k&=E[\varepsilon_{k}\varepsilon_{k}^T]\\&=(I-K_kH_k)E[\varepsilon_{k-1}\varepsilon_{k-1}^T](I-K_kH_k)^T+K_kE[\nu_k\nu_k^T]K_k^T\\&=(I-K_kH_k)P_{k-1}(I-K_kH_k)^T+K_kR_kK_k^T\end{aligned} Pk=E[εkεkT]=(I−KkHk)E[εk−1εk−1T](I−KkHk)T+KkE[νkνkT]KkT=(I−KkHk)Pk−1(I−KkHk)T+KkRkKkT该式由于对称性保证了良好协方差矩阵的正定性，后面将会举例说明在计算精度有限的系统上，这种形式能保证数值计算。有矩阵求导公式 ∂ Tr ( A B A T ) ∂ A = 2 A B \dfrac{\partial \text{Tr}(ABA^T)}{\partial A}=2AB ∂A∂Tr(ABAT)=2AB（要求 B B B是对称阵）。 J k J_k Jk对 K k K_k Kk求偏导为 ∂ J k ∂ K k = 2 ( I − K k H k ) P k − 1 ( − H k T ) + 2 K k R k \dfrac{\partial J_k}{\partial K_k}=2(I-K_kH_k)P_{k-1}(-H_k^T)+2K_kR_k ∂Kk∂Jk=2(I−KkHk)Pk−1(−HkT)+2KkRk一阶偏导数等于零，计算得到增益 K k K_k Kk为 K k = P k − 1 H k T ( H k P k − 1 H k T + R k ) − 1 K_k=P_{k-1}H_k^T(H_kP_{k-1}H_k^T+R_k)^{-1} Kk=Pk−1HkT(HkPk−1HkT+Rk)−1
递推最小二乘估计步骤
计算增益
K
k
=
P
k
−
1
H
k
T
(
H
k
P
k
−
1
H
k
T
+
R
k
)
−
1
K_k=P_{k-1}H_k^T(H_kP_{k-1}H_k^T+R_k)^{-1}
Kk=Pk−1HkT(HkPk−1HkT+Rk)−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+Kk(yk−Hkx^k−1)
协方差更新
P
k
=
(
I
−
K
k
H
k
)
P
k
−
1
(
I
−
K
k
H
k
)
T
+
K
k
R
k
K
k
T
P_k=(I-K_kH_k)P_{k-1}(I-K_kH_k)^T+K_kR_kK_k^T
Pk=(I−KkHk)Pk−1(I−KkHk)T+KkRkKkT
协方差更新公式的三种形式
协方差更新公式有三种形式，在数学上是等价的，但三者在数值计算上有所不同，分别为
P
k
=
(
I
−
K
k
H
k
)
P
k
−
1
(
I
−
K
k
H
k
)
T
+
K
k
R
k
K
k
T
P
k
=
(
I
−
K
k
H
k
)
P
k
−
1
P
k
=
[
P
k
−
1
−
1
+
H
k
T
R
k
−
1
H
k
]
−
1
\begin{aligned}&P_k=(I-K_kH_k)P_{k-1}(I-K_kH_k)^T+K_kR_kK_k^T\\&P_k=(I-K_kH_k)P_{k-1}\\&P_k=[P_{k-1}^{-1}+H_k^TR_k^{-1}H_k]^{-1}\end{aligned}
Pk=(I−KkHk)Pk−1(I−KkHk)T+KkRkKkTPk=(I−KkHk)Pk−1Pk=[Pk−1−1+HkTRk−1Hk]−1第二种形式计算更为简便，但无法确保协方差的正定性；相比之下，第一种形式在结构上对称，能保证协方差矩阵的正定。
例： 假设
H
1
=
1
,
R
1
=
0
H_1=1,R_1=0
H1=1,R1=0，初始协方差设为
P
0
=
6
P_0=6
P0=6，计算精度只能达到
0.001
0.001
0.001，计算增益
K
1
K_1
K1为
K
1
=
P
0
(
P
0
)
−
1
=
6
(
1
6
)
=
(
6
)
(
0.167
)
=
1.002
K_1=P_0(P_0)^{-1}=6(\dfrac{1}{6})=(6)(0.167)=1.002
K1=P0(P0)−1=6(61)=(6)(0.167)=1.002若用第二种协方差更新形式得到
P
1
P_1
P1为
P
1
=
(
1
−
K
1
)
P
0
=
(
−
0.002
)
(
6
)
=
−
0.012
P_1=(1-K_1)P_0=(-0.002)(6)=-0.012
P1=(1−K1)P0=(−0.002)(6)=−0.012这样计算的方差是个负值，显然是错误的。若采用第一种对称形式计算
P
1
=
(
1
−
1.002
)
P
0
(
1
−
1.002
)
=
0
P_1=(1-1.002)P_0(1-1.002)=0
P1=(1−1.002)P0(1−1.002)=0因为超出计算机计算精度，所得结果为0，这在理论上是个合理的值。第一种协方差更新公式在计算精度有限情况下保证了协方差矩阵的正定。
2万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
