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

目录

背景

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

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

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

数学描述

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

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

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

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

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

以上是数学描述。

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

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

正交投影引理

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

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

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

OP[A1θ1+A2θ2∣X]=A1OP[θ1∣X]+A2OP[θ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)xk]Tx(k)=[x(k-1) \; x_k]^Tx(k)=[x(k−1)xk​]T，这里面x(k)，x(k−1)，xkx(k)，x(k-1)，x_kx(k)，x(k−1)，xk​都是随机矢量。那么，对于随机矢量sss，有：

OP[s∣x(k)]=OP[s∣x(k−1)]+OP[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∣xk​]

其中，s~=s−OP[s∣s(k−1)],x~k=xk−OP[xk∣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)],xk​=xk​−OP[xk​∣x(k−1)]

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

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

可以察看参考文献。

LMSE算法

先验条件;

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

XXX:均值XθX_\thetaXθ​,协方差矩阵CXC_XCX​已知

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

解析法

要想使θ^\hat{\theta}θ^与θ\thetaθ尽可能接近，只需求解函数E[(θ−θ^)T(θ−θ^)]E[(\theta- \hat{\theta})^T(\theta-\hat{\theta})]E[(θ−θ^)T(θ−θ^)]的最小值，这里面θ\thetaθ是真实值，是一

个常数。函数可以理解为内积后取平均。

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

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

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

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

得到：

aL=μθ−CθXCX−1μXa_L=\mu_\theta-C_{\theta X}C_{X}^{-1}\mu_XaL​=μθ​−CθX​CX−1​μX​

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

代入原函数中得到;

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

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

θ^=μθ+CθHT(HCθHT+Cn)−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​)

迭代法

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

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

OP[θ(k)∣x(k)]=OP[θ(k−1)∣x(k−1)]+OP[θ~(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)=θ−OP[θ∣X(k−1)],x~k=xk−OP[xk∣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)],xk​=xk​−OP[xk​∣X(k−1)] ,这一步是更新的增量，相当于把第kkk次观测中与前k−1k-1k−1次观测相关的信息去掉，留下新的信息。

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

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

M0=CθM_0=C_\thetaM0​=Cθ​

for1:Kfor \ 1:Kfor1:K

Kk=Mk−1HkT(HkMk−1HkT+Cnk)−1K_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+Kk(Xk−Hkθ^k−1−μnk)\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​​)

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

endendend

LS算法

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

解析法

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

求导=0得:

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

迭代法

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

递推原理如下:

OP[Xr∣X(k)]=OP[Xr∣X(k−1)]+OP[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=Xr−OP[Xr∣X(k−1)]M=X_r-OP[X_r|X(k-1)]M=Xr​−OP[Xr​∣X(k−1)]

N=xk−OP[xk∣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矩阵是模型矩阵,里面不可能有观测量XH=[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_ntheta=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]. 电子工业出版社, .

如果觉得《线性最小均方误差算法（LMSE） 最小二乘法（LS）》对你有帮助，请点赞、收藏，并留下你的观点哦！