线性最小均方误差算法（LMSE），最小二乘法（LS）

目录

• 背景
• 正交投影引理
• LMSE算法
• LS算法
• 直线拟合

背景

  对于一个系统，在给予一定的输入，那么通常都会产生相对应的输出。在实际的系统中，这样的输出必然伴随着噪声，这样被噪声污染的输出通常是传感器的输出信号，也叫观测信号

  同时，如果系统的模型是清晰的，我们可以通过严格的理论计算来得到真实值，通过作差的方式变把噪声去除了。

  然而在实际的系统中，我们对整个系统的物理模型通常是未知的或者有一些参数未知，也可能是模型不准确，某些参数值有一定的偏差。因此，我们为了得到真实的信号，就需要利用观测信号估计系统的模型，尽可能地去除噪声的影响（其实是信号处理的思想）

数学描述

   X = ( x 1 . . x n ) T X=(x_1 .. x_n)^T X=(x1​..xn​)T是观测信号， f ( u ， θ ) f(u，\theta) f(u，θ)是系统的模型， u , θ u,\theta u,θ 分别是系统的输入与未知参数， n n n是噪声。

  整个观测的过程可以描述为：具有的未知参数的系统在输入 u u u的激励下，产生了伴随噪声 n n n的输出（观测信号） X X X。

  即： X = f ( u , θ ) + n X=f(u,\theta)+n X=f(u,θ)+n

  而LMSE，LS算法要做的工作是利用 X X X,去尽可能地得到一个 θ ^ \hat{\theta} θ^，使 θ ^ \hat{\theta} θ^接近真实的 θ \theta θ，而算法的条件是整个系统模型 f ( u , θ ) f(u,\theta) f(u,θ)得是线性的。

  即： X = H θ + n X=H\theta+n X=Hθ+n,也称为线性观测模型。

以上是数学描述。

  既然是要通过已知的观测信号 X X X来得到真实参数的一个进可能近似 θ ^ \hat{\theta} θ^，那么在LMSE与LS算法中，我们是把 θ ^ \hat{\theta} θ^看作是观测 X X X的线性组合。

即： θ ^ = a + B X \hat{\theta}=a+BX θ^=a+BX,在空间里面表现为真实值在 X X X上的投影。

正交投影引理

  定义：若 θ , X \theta,X θ,X是空间中随机矢量，那么 θ \theta θ在 X X X上的投影定义为 θ , X \theta,X θ,X的内积，记作 O P < θ ∣ X > OP<\theta|X> OP<θ∣X>

引理Ⅰ：若 θ , X \theta,X θ,X是空间中随机矢量， θ \theta θ在 X X X上的投影唯一。

引理Ⅱ：正交投影满足线性性

O P [ A 1 θ 1 + A 2 θ 2 ∣ X ] = A 1 O P [ θ 1 ∣ X ] + A 2 O P [ θ 2 ∣ X ] OP[A_1\theta_1+A_2\theta_2|X]=A_1OP[\theta_1|X]+A_2OP[\theta_2|X] OP[A1​θ1​+A2​θ2​∣X]=A1​OP[θ1​∣X]+A2​OP[θ2​∣X]

以上两个引理比较简单，第三个引理在LMSE递推算法中至关重要。

引理Ⅲ：

记 x ( k ) = [ x ( k − 1 )    x k ] T x(k)=[x(k-1) \; x_k]^T x(k)=[x(k−1)xk​]T，这里面 x ( k ) ， x ( k − 1 ) ， x k x(k)，x(k-1)，x_k x(k)，x(k−1)，xk​都是随机矢量。那么，对于随机矢量 s s s，有：

O P [ s ∣ x ( k ) ] = O P [ s ∣ x ( k − 1 ) ] + O P [ s ~ ∣ x ~ k ] OP[s|x(k)]=OP[s|x(k-1)]+OP[\widetilde{s}|\widetilde{x}_k] OP[s∣x(k)]=OP[s∣x(k−1)]+OP[s ∣x k​]

其中， s ~ = s − O P [ s ∣ s ( k − 1 ) ] , x ~ k = x k − O P [ x k ∣ x ( k − 1 ) ] \widetilde{s}=s-OP[s|s(k-1)],\widetilde{x}_k=x_k-OP[x_k|x(k-1)] s =s−OP[s∣s(k−1)],x k​=xk​−OP[xk​∣x(k−1)]

因为要与随机变量的统计特征联系起来，所以可以推倒出：

O P [ s ~ ∣ x ~ k ] = E ( s ~ x ~ k T ) [ E ( x ~ k x ~ k T ) ] − 1 x ~ k OP[\widetilde{s}|\widetilde{x}_k]=E(\widetilde{s}{\widetilde{x}_k}^T)[E(\widetilde{x}_k{\widetilde{x}_k}^T)]^{-1}\widetilde{x}_k OP[s ∣x k​]=E(s x k​T)[E(x k​x k​T)]−1x k​；具体的推导过程

可以察看参考文献。

LMSE算法

先验条件;

θ \theta θ：均值 μ θ \mu_\theta μθ​,协方差矩阵 C θ C_\theta Cθ​已知

X X X:均值 X θ X_\theta Xθ​,协方差矩阵 C X C_X CX​已知

互协方差矩阵 C θ X C_{\theta X} CθX​已知。

解析法

要想使 θ ^ \hat{\theta} θ^与 θ \theta θ尽可能接近，只需求解函数 E [ ( θ − θ ^ ) T ( θ − θ ^ ) ] E[(\theta- \hat{\theta})^T(\theta-\hat{\theta})] E[(θ−θ^)T(θ−θ^)]的最小值，这里面 θ \theta θ是真实值，是一
个常数。函数可以理解为内积后取平均。

解析法思想很简单，求导数取极值即可。

令 θ ^ = a + B X \hat{\theta}=a+BX θ^=a+BX

E [ ( θ − θ ^ ) T ( θ − θ ^ ) ] E[(\theta- \hat{\theta})^T(\theta-\hat{\theta})] E[(θ−θ^)T(θ−θ^)]= E [ ( θ − a − B X ) T ( θ − a − B X ) ] E[(\theta- a-BX)^T(\theta-a-BX)] E[(θ−a−BX)T(θ−a−BX)]

分别对 a , B a,B a,B求偏导=0

得到：
a L = μ θ − C θ X C X − 1 μ X a_L=\mu_\theta-C_{\theta X}C_{X}^{-1}\mu_X aL​=μθ​−CθX​CX−1​μX​

B L = C θ X C X − 1 B_L=C_{\theta X}C_{ X}^{-1} BL​=CθX​CX−1​

代入原函数中得到;

θ ^ = μ θ + C θ X C X − 1 ( X − μ X ) \hat{\theta}=\mu_\theta+C_{\theta X}C_{X}^{-1}(X-\mu_X) θ^=μθ​+CθX​CX−1​(X−μX​)

若观测与噪声独立，且噪声的统计特征已知，解析公式可以进一步简化为：

θ ^ = μ θ + C θ H T ( H C θ H T + C n ) − 1 ( X − H μ θ − μ n ) \hat{\theta}=\mu_\theta+C_{\theta}H^T(HC_{\theta}H^T+C_n)^{-1}(X-H{\mu_{\theta}}-\mu_n) θ^=μθ​+Cθ​HT(HCθ​HT+Cn​)−1(X−Hμθ​−μn​)

迭代法

基本思想：需要迭代 k k k次，且 θ ^ ( k ) = O P [ θ ( k ) ∣ x ( k ) ] \hat{\theta}(k)=OP[\theta(k)|x(k)] θ^(k)=OP[θ(k)∣x(k)]，且 X ( k ) = [ x ( k − 1 )     x k ] T X(k)=[x(k-1) \ \ \ x_k]^T X(k)=[x(k−1)   xk​]T，在这里， x k x_k xk​是新的一个观测值。

那么，根据正交投影引理Ⅲ，有：

O P [ θ ( k ) ∣ x ( k ) ] = O P [ θ ( k − 1 ) ∣ x ( k − 1 ) ] + O P [ θ ~ ( k ) ∣ x ~ ( k ) ] OP[\theta(k)|x(k)]=OP[\theta(k-1)|x(k-1)]+OP[\widetilde{\theta}(k)|\widetilde{x}(k)] OP[θ(k)∣x(k)]=OP[θ(k−1)∣x(k−1)]+OP[θ (k)∣x (k)]

同理，根据正交投影引理Ⅲ：

θ ~ ( k ) = θ − O P [ θ ∣ X ( k − 1 ) ] , x ~ k = x k − O P [ x k ∣ X ( k − 1 ) ] \widetilde{\theta}(k)=\theta-OP[\theta|X(k-1)],\widetilde{x}_k=x_k-OP[x_k|X(k-1)] θ (k)=θ−OP[θ∣X(k−1)],x k​=xk​−OP[xk​∣X(k−1)] ,这一步是更新的增量，相当于把第 k k k次观测中与前 k − 1 k-1 k−1次观测相关的信息去掉，留下新的信息。

进一步推倒可得到递推算法;

θ ^ 0 = μ θ \hat{\theta}_0=\mu_\theta θ^0​=μθ​

M 0 = C θ M_0=C_\theta M0​=Cθ​

f o r   1 : K for \ 1:K for 1:K

K k = M k − 1 H k T ( H k M k − 1 H k T + C n k ) − 1 K_k=M_{k-1}{H_k}^T(H_kM_{k-1}{H_k}^T+C_{n_k})^{-1} Kk​=Mk−1​Hk​T(Hk​Mk−1​Hk​T+Cnk​​)−1

θ ^ k = θ ^ k − 1 + K k ( X k − H k θ ^ k − 1 − μ n k ) \hat{\theta}_k=\hat{\theta}_{k-1}+K_k(X_k-H_k\hat{\theta}_{k-1}-\mu_{n_k}) θ^k​=θ^k−1​+Kk​(Xk​−Hk​θ^k−1​−μnk​​)

M k = ( I − K k H k ) M k − 1 M_k=(I-K_kH_k)M_{k-1} Mk​=(I−Kk​Hk​)Mk−1​

e n d end end

LS算法

LS算法的特点是不需要各类参数以及噪声的先验信息，相当于直接把真实值看成参数，并且直接在 X X X上面进行投影。然后误差均方最小。

解析法

解析法与LSME算法一样，只是目标函数变成了 ( X − H θ ^ ) T ( X − H θ ^ ) (X-H\hat{\theta})^T(X-H\hat{\theta}) (X−Hθ^)T(X−Hθ^)

求导=0得:

θ ^ = ( H T H ) − 1 H T X \hat{\theta}=(H^TH)^{-1}H^TX θ^=(HTH)−1HTX

迭代法

令系统输出的真实值为 X r , X_r, Xr​,，观测值为 X X X

递推原理如下:

O P [ X r ∣ X ( k ) ] = O P [ X r ∣ X ( k − 1 ) ] + O P [ M ∣ N ] OP[X_r|X(k)]=OP[X_r|X(k-1)]+OP[M|N] OP[Xr​∣X(k)]=OP[Xr​∣X(k−1)]+OP[M∣N]

M = X r − O P [ X r ∣ X ( k − 1 ) ] M=X_r-OP[X_r|X(k-1)] M=Xr​−OP[Xr​∣X(k−1)]

N = x k − O P [ x k ∣ X ( k − 1 ) ] N=x_k-OP[x_k|X(k-1)] N=xk​−OP[xk​∣X(k−1)]

直线拟合

%%%本程序基于LMSE算法，利用20个点的数据拟合直线。

%%产生数据
clear;clc;
N=50;
x=1:1:N; 
k=1;
b=0.5;

n=normrnd(0.6,1.5,[N,1]);  %N(0.6,0.3) 20*1的噪声
y=k*x'+b+n;


%%%把数据转换为观测方程
X=y;
AT=1:N;
%H=[X,ones(N,1)];
%问题一 H矩阵是模型矩阵,里面不可能有观测量X
H=[AT',ones(N,1)];

%%给定初始值
u_theta=[0.5;0.8];     %theta均值矩阵
C_theta=[0.1,0;0,0.4]';%theta协方差矩阵

u_n=ones(N,1);           %n的均值矩阵
C_n=diag(0.3*ones(N,1)); %n的协方差矩阵

%%LMSE计算解析解：
%假设观测与噪声独立的

%theta=u_theta+C_theta*H'/(H*C_theta*H'+C_n)*(X-H*u_theta-u_n);
%问题二 多减了u_n
theta=u_theta+C_theta*H'/(H*C_theta*H'+C_n)*(X-H*u_theta);


%%%下面展示利用递推算法;
%%给定初始条件

M=C_theta;
theta_1=u_theta;
%以上相当于迭代第一次

for i=2:N
  
    %K=M*H(1:i,:)'/(H(1:i,:)*M*H(1:i,:)'+C_n(1:i,1:i));%计算K矩阵
    %问题三 递推过程中，H矩阵一直都是1*2的，元素值是随采样时刻变化的
    %C_n是2*1的
    H_1=[i;1]';%1行2列的递推H矩阵
    K=M*H_1'/(H_1*M*H_1'+C_n(i,i));%2行1列的矩阵
    theta_1=theta_1+K*(X(i,1)-H_1*theta_1);%更新theta矩阵
    
    %更新M矩阵
    [k,l]=size(M);
    M=(eye(k)-K*H_1)*M;
       
end
   
   
%%下面利用最小二乘法
theta_2=(H'*H)\H'*X;

%%作图对比
figure(1);
plot(x',y,'*');  %观测点

hold on;
x=1:0.1:N;
H_=[x',ones(length(x),1)];
y_=H_*theta;

plot(x',y_,'r');   %LMSE解析法

hold on;

y_1=H_*theta_1;
plot(x',y_1,'g');  %LMSE迭代法 

hold on;
y_2=H_*theta_2;
plot(x',y_2,'b');  %LS算法

legend('测量曲线','LMSE解析法','LMSE迭代法','LS');
hold off

效果：

三种算法拟合效果几乎一样。

参考文献：

赵树杰, 赵建勋. 信号检测与估计理论[M]. 电子工业出版社, 2013.