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

目录

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

背景

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

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

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

数学描述

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

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

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

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

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

以上是数学描述。

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

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

正交投影引理

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

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

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

$OP[A_1{\theta }_1+A_2{\theta }_2|X]=A_{1}OP[{\theta }_1|X]+A_{2}OP[{\theta }_2|X]$

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

引理Ⅲ：

记$x(k)=[x(k-1)x_k{]}^T$，这里面$x(k)，x(k-1)，x_k$都是随机矢量。那么，对于随机矢量$s$，有：

$OP[s|x(k)]=OP[s|x(k-1)]+OP[\tilde{s}|{\tilde{x}}_k]$

其中，$\tilde{s}=s-OP[s|s(k-1)],{\tilde{x}}_k=x_k-OP[x_k|x(k-1)]$

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

$OP[\tilde{s}|{\tilde{x}}_k]=E(\tilde{s}{{\tilde{x}}_k}^T)[E({\tilde{x}}_k{{\tilde{x}}_k}^T){]}^{-1}{\tilde{x}}_k$；具体的推导过程

可以察看参考文献。

LMSE算法

先验条件;

$\theta$：均值${\mu }_{\theta }$,协方差矩阵$C_{\theta }$已知

$X$:均值$X_{\theta }$,协方差矩阵$C_X$已知

互协方差矩阵$C_{\theta X}$已知。

解析法

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

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

令$\hat{\theta }=a+BX$

$E[(\theta -\hat{\theta }{)}^T(\theta -\hat{\theta })]$=$E[(\theta -a-BX{)}^T(\theta -a-BX)]$

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

得到：
$a_L={\mu }_{\theta }-C_{\theta X}{C}_X^{-1}{\mu }_X$

$B_L=C_{\theta X}{C}_X^{-1}$

代入原函数中得到;

$\hat{\theta }={\mu }_{\theta }+C_{\theta X}{C}_X^{-1}(X-{\mu }_X)$

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

$\hat{\theta }={\mu }_{\theta }+C_{\theta }{H}^T(H{C}_{\theta }{H}^T+C_n{)}^{-1}(X-H{\mu }_{\theta }-{\mu }_n)$

迭代法

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

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

$OP[\theta (k)|x(k)]=OP[\theta (k-1)|x(k-1)]+OP[\tilde{\theta }(k)|\tilde{x}(k)]$

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

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

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

${\hat{\theta }}_0={\mu }_{\theta }$

$M_0=C_{\theta }$

$for1:K$

$K_k=M_{k-1}{H_k}^T(H_k{M}_{k-1}{H_k}^T+C_{n_k}{)}^{-1}$

${\hat{\theta }}_k={\hat{\theta }}_{k-1}+K_k(X_k-H_k{\hat{\theta }}_{k-1}-{\mu }_{n_k})$

$M_k=(I-K_k{H}_k)M_{k-1}$

$end$

LS算法

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

解析法

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

求导=0得:

$\hat{\theta }=(H^{T}H{)}^{-1}{H}^{T}X$

迭代法

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

递推原理如下:

$OP[X_r|X(k)]=OP[X_r|X(k-1)]+OP[M|N]$

$M=X_r-OP[X_r|X(k-1)]$

$N=x_k-OP[x_k|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.