目录
背景正交投影引理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]=A1OP[θ1∣X]+A2OP[θ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(sxkT)[E(xkxkT)]−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θXCX−1μX
BL=CθXCX−1B_L=C_{\theta X}C_{ X}^{-1}BL=CθXCX−1
代入原函数中得到;
θ^=μθ+CθXCX−1(X−μX)\hat{\theta}=\mu_\theta+C_{\theta X}C_{X}^{-1}(X-\mu_X)θ^=μθ+CθXCX−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−1HkT(HkMk−1HkT+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−KkHk)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)》对你有帮助,请点赞、收藏,并留下你的观点哦!