在前面几章里,是通过考虑无偏估计的类型和确定最小方差来寻找最优和接近最优估计。
现在我们抛开这些限制,研究另一类估计,这类估计一般不与最优性质相关,但对于很多估计问题却很有意义,这就是最小二乘估计。
最小二乘估计是相当古老的估计方法,早在1795年高斯就使用这种方法来研究行星的运动。
这种方法的突出特点是对观测数据不做任何概率或统计的描述,而仅仅假设一个数学模型,实现容易。所以,应用非常广泛。
8.1 最小二乘估计方法
前面我们确定一个好的估计的焦点就是寻找的估计是无偏的和有最小方差。在选择方差作为“好”的度量中,无疑地是试图使参量的估计值和真值之间的误差(平均)最小。
在最小二乘(Least Square,LS)方法中,是使观测数据和假设信号之间的平方误差最小。
LS法说明:假设取决于未知参量θ的模型产生信号s[n],这个信号s[n]是完全确定的信号,由于观测噪声或模型不准确,观测到的信号是受干扰的信号,用观测数据x[n]表示,θ的最小二乘估计(LSE)就是选择使s[n]最靠近观测数据x[n]的θ值。
图8-1 最小二乘估计方法
靠近程度用LS误差指标来度量
J(θ)=∑(x[n]−s[n]) n=0,1,",N−1
n=0N−12 (8-1)
式中J是θ的函数,是通过s[n]反映出来的。使J(θ)最小的θ就是LSE。注意,这里没有做出关于观测数据x[n]的概率假设。这种方法无论对高斯还是非高斯噪声同样有效。当然,LSE的性能无疑取决于干扰噪声以及模型的性质。
LSE通常应用于数据的准确统计特性未知,或不能找出最优估计,或最优估计在实际应用时太复杂的场合。
例8-1直流电平信号
假设图8-1中的信号是s[n]=A,观测数据为x[n](n=0,1,",N−1),那么,根据LS法,使(8-1)式或下式最小可估计A。
J(A)=N−1
n=02()x[n]−A ∑
上式对θ求导并令结果为零可得
1A=NN−1n=0∑x[n]=
这正是我们熟悉的样本均值估计,然而,这里没有声明在MVU意义上的最优,而仅仅是使LS误差最小。
原因是什么?
根据前面的讨论,如果x[n]=A+w[n],w[n]是零均值的WGN,那么,LSE也就是MVU估计,否则,LSE不是MVU估计。为了理解这一点,考虑噪声不是零均值情况,假设噪声为
w[n]=E(w[n])+w'[n]
式中w'[n]是零均值噪声,这时观测数据则为
x[n]=A+E(w[n])+w'[n]
显然,样本均值估计实际上是A+E(w[n])的估计,这时的LSE是有偏估计,不是MVU估计。
也就是说,必须假设观测数据是由确定信号和零均值噪声组成,在信号参量的正确选择下,误差e[n]=x[n]−s[n]在平均意义下趋于零,则(8-1)式的最小是合理的。
如果假设的直流量信号模型不正确,例如描述观测数据的模型是x[n]=A+
Bn+w[n],则这个模型误差也将引起LSE是有偏的。
例8-2 正弦信号的频率估计
考虑信号模型 s[n]=cos2πf0n
式中f0是待估计参量。f0的LSE通过使下式最小找到。
J(f0)=N−1
n=02()x[n]−cos2πfn ∑0
与直流量信号对比,在直流量信号中,容易找到使LS误差指标最小的参量估计。然而这里的LS误差不是f0的线性函数,不能找到使LS误差最小值的参量估计的闭合形式。
对于信号模型是待估计参量的线性函数的最小二乘情况我们称为线性最小二乘问题,如例8-1。
否则,这个问题称为非线性最小二乘问题,如例8-2。
求解非线性LS问题是通过栅格寻优或迭代最小化方法(前面已介绍)。
应该注意的是线性最小二乘问题是指信号是待估计参量的线性函数,信号本身不必是线性的。
例8-3 正弦信号幅度估计
考虑信号模型
s[n]=Acos2πf0n
式中f0是已知的,A是待估计参量。则LSE在A上通过使下式最小找到。
J(A)=∑(x[n]−Acos2πf0n)
n=0N−12
因为J(A)是A的二次型函数,通过求导可以容易求出最小值。这里信号是待估计参量的线性函数,而信号本身不是线性的。实际中,线性最小二乘问题是我们希望的。
如果A 和f0均为待估计参量,
J(A,f0)=∑(x[n]−Acos2πf0n)
n=0N−12
则是可分离的最小二乘问题。
8.2 线性最小二乘估计
8.2.1 标量最小二乘
对于标量参数估计,如应用线性最小二乘方法,假设
s[n]=θh[n] (8-2)
式中h[n]是已知序列,那么LS误差准则为
J(θN−1
n=0)=∑(x[n]−θh[n])2 (8-3)
N−1使上式最小可得LSE如下
θˆ=n=0
N−1
n=0∑x[n]h[n] (8-4)
h2∑[n]
将(8-4)式代入(8-3)式得LS误差的最小值为
Jmin =Jθˆ=()∑(x[n]−θˆh[n])(x[n]−θˆh[n])
n=0
N−1N−1
=
n=0∑x[n](x[n]−θˆh[n])−θˆ∑h[n](x[n]−θˆh[n]) n=0N−1
N−1
=
n=0∑x[n]−θˆ∑x[n]h[n] s
2n=0N−1
(8-5)
代入θˆ得
⎛N−1⎞⎜∑x[n]h[n]⎟⎜⎟N−1n=0⎝⎠=∑x2[n]−N−1n=0
n=02Jmin (8-6) ∑h2[n]
N−12显然,最小的LS误差是观测数据的原始能量或∑nx[n]减去信=0
号拟合产生的能量。
ˆ=,θ=A,对于例8-1,则h[n]=1,代入(8-4)式得A根据(8-5)
式有
N−1
Jmin=
n=0∑x2[n]−N2
如果观测数据是无噪声的,即x[n]=A,则Jmin=0,或者说,对于这些观测数据有一个理想的LS拟合。
另一方面,如果x[n]=A+w[n],且E(w2[n])>>A2,那么,∑n=0x2[n]N−1N>>2,则最小的LS误差为
N−1
Jmin≈
n=0∑x2[n]
与原始能量相差不多。LS误差的最小值总是满足下式
N−1
0≤Jmin≤
n=0∑x2[n] (8-7)
8.2.2 矢量最小二乘
上述标量LSE的结果可直接扩展到矢量参量LSE,且有很大的实际应用价值。假设矢量参量θ是p×1维的,信号s=[s[0]s[1]"s[N−1]]T是待估计参量的线性函数,假设
s=Hθ (8-8)
式中H是一个满秩的N×p(N>p)矩阵,H称为观测矩阵。这是一个线性模型,没有通常的噪声PDF假设。在第4章很多信号满足这个模型。θ的LSE可使下式最小容易找到,即
J(θ)=∑(x[n]−s[n]])=(x−Hθ)
n=0N−12T(x−Hθ) (8-9)
进一步得
J(θ)=xTx−xTHθ−θTHTx+θTHTHθ
=xTx−2xTHθ+θTHTHθ
(注意xTHθ是标量),J对θ的导数(或梯度)
∂J(θ)=−2HTx+2HTHθ ∂θ
令这个导数为零,得LSE
ˆ=HTHθ()−1HTx (8-10)
ˆ的方程HTHθ=HTx称为标准方程。H的满秩保证了用来求解θ
H
TH的逆阵存在。
我们发现,这个LSE与线性模型的有效估计以及BLUE有相
ˆ是对观测数据没有做任何假设同的函数形式,式(8-10)中的θ
下产生的估计,而对于BLUE要求E(x)=Hθ和Cx=σ2I(见第6章),对于线性模型的有效估计还另外要求x是高斯分布的(见第4章)。
如果LSE满足这些假设,则容易确定LSE的统计特性,也就是第4,6章已经给出的形式。否则,确定LSE的统计特性相当困难。
将(8-10)代入(8-9)式可得LS误差最小值为
ˆ=x−HθˆJmin=Jθ
T=⎛⎜x−HHH⎝
TT()()T(x−Hθˆ) ()−1⎛THx⎞⎟⎜x−HHH⎠⎝T−1TTT()−1HTx⎞⎟ ⎠T=xI−H(HH)H
TT−1)(I−H(HH)H)x ( =x(I−H(HH)H)x (8-11)−1T
上式推导中,利用了I−H(HTH)HT是幂等矩阵的性质,即A2=A。−1
Jmin也可写成另一种形式
Jmin=xx−xHHHHx (8-12)
ˆ (8-13)
=xTx−HθTT()T−1T()
8.2.3 加权最小二乘
线性LS问题可扩展到加权LS,即在(8-9)式中增加一个N×N维的正定加权矩阵W如下
J(θ)=(x−Hθ)TW(x−Hθ) (8-14)
=wi>0的对角矩阵,那在误差准则中引进加权因子是为了强调那些可靠观测数据的重要性。例如,在例8-1中,如果W是[W]ii
么,LS误差为 J(A)=N−1
n=0
2如果w[n]是零均值和具有方差σn的不相关噪声,那么,加权因子∑wn(x[n]−A)2
的合理选择是wn2=n。
这个选择将产生A的LSE为
N−1
2n=0σn
N−1
n=0ˆ=A∑x[n]1 (8-15) ∑σ2n
这正是我们熟悉的BLUE,这是由于w[n]是不相关的,所以有W=C−1。
对于加权LSE的一般形式容易求得
ˆ=HTWHθ()−1HTWx (8-16)
它的LS误差的最小值为
Jmin=xT(W−WH(HWH)T−1HTWx (8-17)
)
8.3 LSE的几何解释
下面用几何的观点来解释线性LS法,通过几何解释能更清楚地揭示这种方法的本质。
回顾一般的线性信号模型s=Hθ,如果用hi表示H的各列,则有
s=h1[h2⎡θ1⎤⎢θ⎥1⎥⎢"hp=⎢#⎥⎢⎥⎢⎦⎣θp⎥]∑θhii=1pi
因此,信号模型可以看作是“信号”矢量{h1,h2,",hp}的线性组合。
例8-4 傅里叶分析 在例4-2中,假设M
=1,则信号模型为
s[n]=acos2πf0n+bsin2πf0n n=0,1,",N−1
式中f0是已知频率,θ=[a
b]
T
是待估计量。写成矢量形式为
10⎤⎡s[0]⎤⎡
⎥a⎢s[1]⎥⎢cos2πfπsin2f00⎥⎡⎤⎥=⎢⎢
⎥ (8-18) ⎥⎢⎢#⎥⎢##b⎣⎦⎥⎥⎢⎢
⎣s[N−1]⎦⎣cos2πf0(N−1)sin2πf0(N−1)⎦
上式看出,H的两列分别由余弦和正弦序列组成,即
h1=[1cos2πf0"cos2πf0(N−1)]h2=[0sin2πf0"sin2πf0(N−1)]
T
T
所以有 s=
ah
1
+bh
2
由LS误差定义可得
J(θ)=(x−Hθ)(x−Hθ)
T
T
如果定义N×1维矢量ξ
=[ξ1ξ2"ξN]的欧几里几何长度为
ξ=
=
那么,LS误差也可写为
J(θ)=x−Hθ
2
=x−∑θihi (8-19)
i=1
p
p
2
因此,线性LS法就是使数据矢量x与信号矢量∑i=1θihi之间的距离平方最小,其中信号矢量必须是H各列的线性组合。
数据矢量可以位于N维空间RN中的任意位置上,由p
p
上,记为Sp(由于H的满秩确保其各列是线性独立的,因此由各列所确定的子空间是真正的p维)。现以N=3和p=2为例,对所有可能的θ1和θ2(假设−∞
2上,而数据x通常不在这个子空间上。
图8-2 在R3空间的线性最小二乘的几何解释
ˆ应该是x很明显,在子空间S2上且在几何意义上最靠近x的s
ˆ是x在S2上的正交投影。这意味着在S2上的分量,换句话说,s
ˆ与S2上的所有矢量一定是正交的。误差矢量x−s(如果xTy=0,
ˆ,我则定义在RN上的这两个矢量正交。)在图8-2中,为了确定s
们使用正交条件,也就是误差矢量与信号子空间正交,即
ˆ)⊥S2 (x−s
式中⊥表示正交或垂直。如果上式成立,则一定有
ˆ)⊥h2 ˆ)⊥h1 (x−s(x−s
那么,误差矢量也将与h1和h2任意线性组合正交,根据正交的定义,则有
ˆ)h1=0 (x−sˆ)h=0 (x−s
T
T2
ˆ=θ1h1+θ2h2,则 令s
(x−θ1h1−θ2h2)Th1
=0
(x−θ1h1−θ2h2)
其矩阵形式
T
h2=0
(x−Hθ)[h1
Th2]=0T
或 (x−Hθ)
T
H=
0T (8-20)
由此得LSE为
ˆ=HTHθ
()
−1
HTx
注意到如果e=x−Hθ表示误差矢量,那么由(8-20)式得
eTH=0T (8-21)
即误差矢量一定与H的各列正交,这就是著名的正交原理。实际上,这个误差表示了不能被信号模型描述的x的那部分分量。 根据图8-2,LS误差的最小值为
ˆJmin=x−s
2
ˆ=x−Hθ
2
ˆ=x−Hθ
()(x−Hθˆ)
T
代入(8-21)式的结果,则
ˆJmin=x−Hθ
()(x−Hθˆ)
T
ˆ−θˆTHTe =xTx−Hθˆ =xTx−xTHθ
()
=xTI−HHTH
(
()
−1
HTx
)
(8-22)
ˆ拟合或接近数据矢量x的综上,LS法可看作是用另一个矢量sˆ是在RN的p维子空间上矢量问题,其中x是RN空间的矢量,s
{h,h
1
2
ˆ
为x在子空间的这个问题的解就是选择s,",hp}的线性组合。
正交投影。
如果H的各列是标准正交的,那么
(HH)
T
−1
=(I)=I
−1
则
ˆ=HTHθ
()
−1
HTx=HTx
显然,不需要求逆阵。 例8-4 傅里叶分析(继续) 继续考虑例8-4问题,如果f0
意整数,由式(4-13)容易得出
k⎞⎛k⎞⎛
hh2=∑cos⎜2πn⎟sin⎜2πn⎟=0
N⎠⎝N⎠⎝n=0
T
1
N−1
=kN,k是k=1,2,",N2−1中的任
另外,
Th1h1=
N
2
hT2h2=
N
2
因此,h1和h2是正交的,但不是标准正交。
将上面3个等式进行组合,可表示为HTH=(N2)I,所以有
ˆ⎤⎡aˆθ=⎢ˆ⎥=HTH⎣b⎦
()
−1
HTx
=
2THx N
⎡2N−1⎛
x[n]cos⎜2π⎢N∑⎝=0
=⎢n
N−1
⎛⎢2x[n]sin⎜2π∑⎢⎝⎣Nn=0k⎞⎤
n⎟N⎠⎥⎥ k⎞⎥n⎟N⎠⎥⎦
如果用下式代替原来的信号
s[n]=a′
22k⎞k⎞⎛⎛
cos⎜2πn⎟+b′sin⎜2πn⎟ NN⎠NN⎠⎝⎝
那么,H的列矢量是标准正交(单位正交)的。
然而,一般情况下,H的各列不是正交的,所以,信号矢量估计为
ˆ=HHTHˆ=Hθs
()
−1
HTx
这个信号矢量的估计是x在p维子空间的正交投影。如果令N×N维矩阵P=H(HTH)下性质:
1. PT=P,所以P是对称阵。 2. P2=P,所以P是幂等阵。
−1
HT
,则P称为正交投影矩阵或投影阵,P有如
ˆ中恢投影阵是奇异矩阵,如果不是奇异的,那么x就能从s
复出来。这显然是不可能的,因为如果那样,许多x将会有相同的投影。
ˆ=(I−P)x是x在信号矢量子空间的互补子同样,误差矢量e=x−s
空间或正交子空间的投影,矩阵P⊥=I−P也是投影阵,容易验证,它满足上述投影阵的性质。
根据(8-22)式,LS误差的最小值为
Jmin=xT(I−P)x
=xTP⊥x =xTP⊥P⊥x
=P⊥x
2
T
这恰恰是
e
2
。
8.4 按阶递推最小二乘估计
前面介绍的LSE是在信号模型已知的情况下进行的,但有很多情况是信号模型未知,这时,必须假设信号模型。例如,图8-5a给出了实验观测数据,信号模型可以分别假设为
s1(t)=A1 s2(t)=A2+B2t
式中0≤t≤T。如果在时刻t=nΔ时对x(t)采样得x[n],其中Δ=1和
n=0,1,",N−1,则相应的离散信号模型为
s1[n]=A1
s2[n]=A2+B2n
根据信号模型有
⎡1⎤⎢1⎥H1=⎢⎥
⎢#⎥⎢⎥⎣1⎦
⎡10⎤⎢11⎥
⎥ H2=⎢
⎢##⎥⎢⎥
N1−1⎣⎦
利用LSE得截矩和斜率的估计为
ˆ= A1
(8-23)
和
N=1N=1()221N1−ˆ=Anx[n] ∑x[n]−NN+1∑2
NN+1n=0n=0
N=1N=1
112
nx[n] (8-24)
∑x[n]+NN+1∑NN+1n=0n=0
ˆ=−B2
图8-5用最小二乘模拟实验数据
ˆ和sˆ+Bˆt,ˆ1(t)=Aˆ2(t)=A图8-5b和图8-5c分别画出了s其中T=100。122
从图中明显看出,使用两个参量的拟合效果较好。
后面将看到,随着参量的增加,LS误差的最小值将减小。然而在实际应用中,在足以描述观测数据的情况下,应选择最简单的信号模型。
具体实现就是增加多项式的阶次直到增加阶次后使LS误差的最小值仅稍微减小为止。对于图8-5的观测数据,LS误差的最小值Jmin
与参量数目的关系示于图8-6。
图8-6 最小值Jmin与参量数目的关系
模型s1(t)和s2(t)分别对应于k=1和k=2。从图中可以看出,Jmin从一个参量变化到两个参量有大幅度的下降,当再增大阶次时,
Jmin减小很少。
在这个例子中,信号的实际模型是
s(t)=1+0.03t
噪声w[n]是σ2
=0.1的
WGN。如果A和B估计的很准确,那么,LS
误差的最小值为
ˆ,Bˆ=J(A,B)=∑(x[n]−s[n])=∑w2[n]≈Nσ2Jmin=JA
2
n=0
n=0
()
N−1N−1
所以,当达到实际的阶次时,Jmin≈10,从图8-6可验证该值,通过该值也可增加我们对所选择模型的可信度。
在前面的例子中我们看到,当需要确定几个信号模型的LSE时,一个直接的方法就是对每个模型使用(8-10)式
ˆ=HTHθ
()
−1
HTx
另一种方法就是使用按阶递推LS法,后种方法可以减少计算量。
按阶递推LS法是按阶次逐次递推更新LSE,即,它能够根据基于N×k维的H矩阵的LSE计算基于N×(k+1)维的H矩阵的LSE。为了看出其中的原因,我们假设信号模型为
s1[n]=A1 s2[n]=A2+B2n
其中,−M
−M⎤⎡1
⎢1−(M−1)⎥
⎥H2=⎢
⎢#⎥#⎢⎥1M⎣⎦
⎡2M+1
⎢=HTH22⎢0⎢⎣
⎤
⎥ M
2⎥n∑⎥n=−M⎦
容易求得LSE解为
M
1ˆ=Ax[n]∑1
2M+1n=−M
M
1ˆ=Ax[n] ∑ 2
2M+1n=−M
ˆ=−B2
n=−M
M
∑
M
nx[n]n
2
n=−M
∑
可以看出,在这种情况下,当模型中增加一个参数时,A的LSE不会发生变化,这是由H2TH
2的对角性质得出的结论。
h
i
的正交性使得能把x分别投影到h1和h
2
方向,然后把
结果相加。
一般而言,列矢量不是正交的,因此每个投影不是独立的,但能够利用另一组正交的p矢量来代替(Gram-Schmidt正交化)。新的列h
2
被投影到与h1正交的子空间h
'
2
上。那么,LSE
信号估计变为
ˆ=h1θˆ1+h2'θˆ2 s
其中h1θˆ1也是信号仅根据H=h1 的LSE。
因此,这种方法在信号模型中可以用递推的形式将每一项加到信号模型上来表示阶数的更新。
按阶递推LS法可经过纯粹的代数推导获得(附录8A),这里
ˆk表示基于Hk的仅给出最后结果。令Hk表示N×k维观测矩阵,θ
LSE,即
ˆk=HTθkHk
基于Hk的LS误差的最小值为
()
−1
HTkx (8-25)
ˆJmink=x−Hkθk
(ˆ))(x−Hθ (8-26)
T
kk
现将阶次增加为k+1,这将产生一个新的观测矩阵Hk+1,其分块矩阵形式如下
Hk+1=[Hkhk+1]=[N×kN×1] (8-27)
那么,参量的按阶递推LSE为
−1T⊥⎤TT⎡HHHhhPkkkk+1k+1kxˆ⎢θk−⎥T⊥⎡k×1⎤hPh⎢⎥k+1kk+1==⎢T⊥
⎢⎥⎣1×1⎥hk+1Pkx⎦ (8-28) ⎢⎥T⊥
hk+1Pkhk+1⎣⎦
−1
()
ˆk+1θ
T⊥
式中Pk⊥=I−Hk(HT这个kHk)Hk。Pk是x在子空间上的投影矩阵,
子空间与Hk的各列张成的子空间正交。 为了避开计算HTkHk的逆矩阵,令
Dk=HHk
并使用下面的递推公式
(
T
k
)
−1
(8-29)
T
⎡⎤DkHTDkHTkhk+1hk+1HkDkkhk+1
−T⊥⎢Dk+⎥⎡k×kk×1⎤T⊥
hPhhPhk+1kk+1k+1kk+1⎥=⎢Dk+1=⎢T (8-30) hk+1HkDk1⎢⎥⎣1×k1×1⎥⎦−T⊥
⎢⎥T⊥
hPhhPhk+1kk+1k+1kk+1⎦⎣
式中Pk⊥=I−HkDkHTk。
LS误差最小值的阶次递推式为
Jmink+1=Jmink
(h
−
PxT
hk+1Ph
Tk+1
2⊥k⊥kk+1
)
(8-31)
整套算法不需要求逆阵。递推的初值θˆ1、Jmin1和D1
分别由式(8-25)、(8-26)和(8-29)确定,式(8-28)、(8-30)和(8-31)称为按阶递推最小二乘。
下面将这个方法应用于前面的直线拟合举例中。 例8-5 直线拟合
由于s1[n]=A1和s2[n]=A2+B2n(n=0,1,",N−1),所以有
⎡1⎤
⎢1⎥
H1=⎢⎥=1
⎢#⎥⎢⎥⎣1⎦
H2=[H1h2]
⎡0⎤⎢1⎥
⎥。利用式(8-25)和(8-26)计算递推初值 式中h2=⎢
⎢#⎥⎢⎥N−1⎣⎦
−1TTˆˆA1=θ1=(H1H1)H1x=
ˆJ=x−Hθ11 min1
()(x−Hθˆ)=∑(x[n]−)
T
11
n=0
N−1
2
ˆBˆˆ2=A下一步,利用(8-28)计算θ22
[
],即
T
−1TTT⊥⎤⎡HHHhhP111221xˆ⎥⎢θ1−T⊥h2P1h2⎥ˆ2=⎢θ⊥ ⎥⎢hTPx21
⎥⎢T⊥
hPh212⎦⎣
()
其中所需要的项
(H
T
1
H1
)
−1
=
1N
T
P1⊥=I−H1H1H1
()
−1
T
=I−H1
1T
11 N
P1⊥x=x−
111N
T
x=x−
⊥TT
hTPx=hx−21221=
∑nx[n]−∑n
n=0
n=0
N−1
N−1
2
N−1N−1
1T21⎛⎞2hP
h2=hh2−h21=∑n−⎜∑n⎟
NN⎝n=0⎠n=0
T
2
⊥1
T2
()
N−1
⎡1N−1⎡N−1⎤⎤
nnxnn[]−∑⎢⎢∑⎥⎥Nn=0⎣n=0n=0⎦⎥⎢−
2N−1N−1⎥⎢1⎛⎞2N−1
1⎤⎡nn−⎥⎢⎜⎟∑∑ˆN⎝n=0⎠n=0⎥=⎢−N∑nB2⎥ˆ2=⎢θn=0N−1N−1代入得 ⎥ ⎥⎢⎢
ˆnx[n]−nB2⎥⎥⎢⎢∑⎦⎣
n=0n=0⎥⎢
2N−1⎥⎢1⎛N−1⎞2
n−⎜∑n⎟⎥⎢∑N⎝n=0⎠⎥⎢n=0⎦⎣
其中
ˆ=B2
∑∑
N−1n=0
N−1n=0
nx[n]−∑
N−1n=0
n
2
1⎛2n−
N⎜⎝
∑
N−1n=0
⎞
n⎟⎠
化简得
ˆ=B2
∑
N−1n=0
nx[n]−N
N
(N−1)
N
2
2
−12
N−1
6
x[n]+=−∑NN+1n=0N
N−1
12
nx[n] ∑2
N−1n=0
和
N−11N−1ˆˆˆA2=−∑nB2=−B2Nn=02
N−12(2N−1)N−16
x[n]−nx[n] =∑∑NN+1n=0NN+1n=0
这个结果与(
8-24)式相同。
根据(8-31)式LS误差的最小值为
2
Jmin2=Jmin1
h(−
N−1
T
2
Px)
⊥1
⊥
hTP21h2
N−1
2
=Jmin1
⎛⎞
[]nxn−n∑⎜∑⎟
n=0−n=0
2N−1N−1
1⎛⎞2
n−⎜∑n⎟∑N⎝n=0⎠n=0
N−1
N−1
=Jmin1
⎛⎞[]−nxnn∑⎜∑⎟
n=0⎠−⎝n=0
NN2−12
通常,递推过程是连续求解LS问题,直到获得理想阶次为止。
这样不仅获得了理想模型的LSE,同时,也获得了所有低阶次模型的LSE。正如前面所述,这对于模型阶次未知的情况很有用。
8.5 序贯(递推)最小二乘估计
在多数信号处理应用中,都是对连续时间波形进行采样而获得数据,显然,这些数据是按时间顺序获得的,当需要很多的数据时,我们有两种选择,一种是在获得所有有用的数据后一次处理,另一种就是这一节将介绍的,及时处理每一个数据。 特别是,假设已经确定了基于{x[0],x[1],",x[N−1]}的LSE
ˆ,θ现在又观测到了新的数据x[N],能否在不计算
LSE式—(8-10)
ˆ获得基于{x[0],x[1],",x[N]}的LSE? 式,通过修改θ
回答是肯定的。为了区别起见,这种方法称为序贯(递推)最小二乘(sequential least squares)法,前面的最小二乘法称为批量(batch)法。
8.5.1 标量递推最小二乘
仍考虑例8-1直流量信号的估计。已知基于当前观测到的数据的LSE为
ˆ[N−1]=1A
N
ˆ[N]=A
∑x[n]
n=0
N−1
现在又观测到了一个新的数据x[N],则这时的LSE应为
N
1
x[n] ∑N+1n=0
ˆ[N−1]和x[N]求出Aˆ[N]。 现在我们不重新计算这个求和,而通过A
ˆ[NA
1⎛N−1⎞+xnxN[]]=[]∑⎟N+1⎜⎝n=0⎠
N1ˆ[N−1]+=Ax[N] (8-35)N+1N+1ˆ[N−1]+1x[N]−Aˆ[N−1]=A
N+1 (8-36)
()
由此看出,新的估计等于原来的估计加上一个修正项。
ˆ[N−1]是基于很多的 修正项随着N的增大而减小,这反映了A
ˆ[N−1]可以看作数据,因此,应该给以更多的权重。同时,x[N]−A
是预测误差,如果这个误差是零,则不进行修正。否则,修正新的估计。
LS误差的最小值也可得出递推形式。基于截止N−1时刻数据的LS误差的最小值为
ˆ[N−1]2 Jmin[N−1]=∑x[n]−A
n=0N−1
((
)
增加新的数据后的这个误差
Jmin[N]=∑
n=0N
2
ˆx[n]−A[N]
)
整理后得
J
min[N]=Jmin[N−1]+
Nˆ[N−1]2 x[N]−AN+1
()
由此看出,增加新的数据使LS误差的最小值增大,这很容
易解释,这是因为对于每一个新的数据,都使平方误差项增加。
8.5.2 加权递推最小二乘
如果对线性最小二乘问题进行加权,假设加权矩阵W是
[W]ii
=2的对角阵,由(8-15)式得加权LSE为
ˆ[N−1]=A
∑
N−1n=0
N−1
x[n]
σσ
2n2n
∑
1
n=0
不难推出,其递推形式为
ˆ[N]=A
∑σ
n=0
Nn=0
N
x[n]
2n2n
∑σ
1
=
∑σ
n=0
N−1
x[n]
2nN
+1
x[N]
2σN
∑σ
n=0
⎛N−11⎞ˆx[N]
−AN1[]⎜∑2⎟2
σNn=0σn=+N
N
11
2
n
∑σ
n=0
2n
∑σ
n=0
2n
ˆ[N−1]−=A
σN
1ˆ
A[N−1]2
x[N]+
2
σN
∑σ
n=0
N
1
2n
∑σ
n=0
N
1
2n
1
ˆ[N−1]+=A
2
σN
∑σ
n=0
N
1
n
ˆ[N−1])(x[N]−A
(8-37)
正如所料,如果σn2=
σ2,则与前面的结果(8-36)式相同。
式(8-37)中修正项的增益因子取决于对新的采样值的可信
2
度,如果新的采样值是噪声或σN
→∞,则不改变前面的LSE,即
2ˆ[N]=Aˆ[N−1]。相反,如果新的采样值是无噪声的或σN那么,A→0,
ˆ[N]→x[N],即忽略所有前面的采样值。即采样因子表示对新的A
数据采样相对于前面采样值的可信度。
现在我们进一步解释这个结果。如果x[n]=A+w[n],w[n]是具有方差σn2的零均值不相关噪声。那么,根据例6-2,这个LSE就是BLUE,所以有
ˆ[N−1]=varA
()
1
∑
N−1n=0
1
σ
2
n
根据(8-37)式,第N次修正项的增益因子为
1
K[N]=
σN
1
1
2n
σ
n=0
N
=
σN
σ
2N
+
[N−1]varA
=
ˆ[N−1]varA
)
ˆ[N−1]+σ2 (8-38) varAN
(
)
ˆ[N−1]很大,则修正项大,否则,0≤K[N]≤1,由此看出,如果varA
ˆ[N−1]),可推导得出K[N]的递推形式:相反。由于K[N]取决于var(A
(
ˆ[NvarA
(])=
1
∑
N
1
=
1
n=0
σ
2
n
∑
N−1n=0
1
σ
2n
+
2N2N
1
=
1
[N−1]varA
σ
2
N
+
σ
2N
=
ˆ[N−1]σvarA
[N−1]+σvarA
(
)
ˆ[N−1]⎛varA
=⎜1−
[N−1]+
σ⎜varA
⎝
(
)
2N
⎞⎟var⎟⎠
(
ˆ[N−1]A
)
或
ˆ[N]=(1−K[N])varAˆ[N−1]varA (8-39)
()()
利用(8-38)和(8-39)式可递推K[N]。这样以来,在递推增益的过程中,同时也递推了LSE的方差。现将这些结果总结如下: 估计量的递推式:
ˆ[N]=Aˆ[N−1]+K[N]x[N]−Aˆ[N−1] (8-40) A
其中
ˆ[N−1]varA
K[N]=
N−1+σ2varA
()
(
)
(8-41)
N
方差的递推式:
ˆ[N]=(1−K[N])varAˆ[N−1]varA (8-42)
()()
具体递推过程如下: 1. 在开始递推时,令
ˆ[0]=x[0]A
ˆ[0]=σ2 varA0
()
ˆ[1]。ˆ[1]2. 由(8-41)式确定K[1],(8-40)式确定A下次迭代的varA
由(8-42)式确定。 3. 重复第2步。 图8-9示出了A=10和σ
n
2
仿真结果。
=1时的递推LS法的Monte-Carlo计算机
图8-9 WGN中的直流量的递推LS估计
图中看到,随着N的增大,方差和增益趋于零,这是由于
11ˆ[N])= var(A=N
∑
1N+1
n=0
σ
2
n
和
K[N]=
1=1
N+1
+1N
同时也看到,估计近似收敛于真值A=10。最后,给出LS误差的最小值的递推式
2
ˆx[N]−A[N−1]Jmin[N]=Jmin[N−1]+
N−1+σ2 (8-43)
varA
(
)
N
8.5.3 矢量加权递推最小二乘
现将获得的标量参量的递推LSE扩展到矢量参量(附录8C)。考虑W=C−1的加权LS误差准则,其中C表示零均值噪声的协方差矩阵,或
J=(x−Hθ)C−1(x−Hθ)
T
这与BLUE有同样的假设条件。根据(8-16)式这个加权LSE
T−1−1T−1ˆ为 θ=HCHHCx
()
ˆ的协方差,根由于C是噪声的协方差矩阵,所以我们用Cθˆ表示θ
据(6-17)式有
Cθˆ=HCH
(
T−1
)
−1
ˆ的递推 如果C是对角阵或噪声是不相关的,那么可以得到θ
LSE,否则,则不能。现假设满足条件,令
222
C[n]=diagσ0,σ1,",σn
()
⎡H[n−1]⎤⎡n×p⎤H[n]=⎢T=⎢⎥⎥ []nph1×⎣⎦⎣⎦
x[n]=[x[0]x[1]"n]T
ˆ的加权ˆ[n]表示基于x[n]或(n+1)个采样数据的θθ
LSE,那么,批量
LSE估计为
ˆ[n]=HT[n]C−1[n]H[n]θ
ˆ[n]的协方差阵为 θ
()
−1
HT[n]C−1[n]x[n] (8-44)
Cθˆ=Σ[n]=HT[n]C−1[n]H[n]
()
−1
(8-45)
再次强调,C[n]是噪声的协方差,∑[n]是LSE的协方差。
递推LSE为: 估计量的递推式:
ˆ[n]=θˆ[n−1]+K[n]x[n]−hT[n]θˆ[n−1] (8-46) θ
()
式中
K[n]=
σn2
Σ[n−1]h[n]
(8-47)
+hT[n]Σ[n−1]h[n]
方差的递推式:
Σ[n]=I−K[n]hT[n]Σ[n−1] (8-48)
()
这里的增益因子K[n]是p×1矢量,协方差阵Σ[n]是p×p维的。从上述等式中发现,矢量参量的递推LSE不需要求逆阵。
估计量的递推过程总结于图8-10,图中粗箭头表示矢量。递
ˆ[n−1]和推过程:首先使用批量估计(8-44)和(8-45)式确定θ
Σ[n−1],然后,对于n≥p利用递推等式(8-46)-(8-48)式进行
递推。需要指出的是,在计算 ˆ[n−1]θ和Σ[n−1]时,要求
HT[n−1]C−1[n−1]H[n−1]可逆。为保证这一点,H[n−1]的秩必须大于
或等于p,又由于H[n−1]是n×p维矩阵,所以必须n≥p(假设它的所有列是线性独立的)。
图8-10 序贯LSE估计
最后,给出LS误差的最小值的递推式
2
ˆx[n]−h[n]θ[n−1]Jmin[n]=Jmin[n−1]+2 (8-49)
σn+hT[n]Σ[n−1]h[n]
T
()
下面讨论关于LSE的应用问题。LSE通常在下述情况下使用: 1. 噪声的统计特性未知,或噪声的统计特性已知,但不能找到最优的MVU估计。
2. 对于已知噪声统计特性,但渐近最优的MLE实现起来太复杂。
8.6 约束最小二乘估计
有时,我们会遇到未知参数受到限制的LS问题。例如:我们希望估计几个信号的幅度,但预先知道这些信号的幅度中有部分是相等的。于是,利用先验知识,总的被估计参数的数目应当有所减少。对于这个例子,参数是线性相关的,从而引出线性约束的线性最小二乘估计问题。
假设参数θ受到r
Aθ=b (8-50)
其中A是一个r×p矩阵,而且已知b是一个r×1矢量。
举例来说,如果p=2而且已知一个参数是另一个参数的负数,于是约束可以表示为θ1+θ2=0。于是A=[1 1]和b=0。
矩阵A总是假设为满秩的(等于r),这对于约束条件相互独立是必要的。应当认识到在约束LS问题中,实际上只有p−r个独立参数。
为了求出线性约束的LSE,我们使用拉格朗日乘因子。即通过使如下的拉格朗日乘因子最小来确定θˆc(c表示约束LSE),
Jc=(x−Hθ)
T
T
x−Hθ+λ(Aθ−b) ()
其中λ是拉格朗日乘因子的r×1矢量。
展开这个表达式,得
Jc=xTx−2θTHTx+θTHTHx+λTAθ−λTb
取关于θ的梯度并利用(4-3)式得
∂Jc
Jc=−2HTx+2HTHx+ATλ
∂θ
令上式等于零,得
1TT−1ˆθc=(HH)Hx−(HTH)−1ATλ
2
ˆ−(HTH)−1ATλ/2=θ (8-51)
其中θ是无约束的LSE。
ˆ
为了求λ,我们利用(8-50)的约束条件,于是
ˆ=Aθˆ−A(HTH)−1ATλ/2=b Aθc
因此,
ˆ−b) λ/2=[A(HTH)−1AT]-1(Aθ
代入(8-51)式,得解为
ˆ=θˆ−(HTH)−1AT[A(HTH)−1AT]-1(Aθˆ−b)(8-52) θc
ˆ=(HTH)−1HTx。 其中,θ
约束LSE是无约束lSE的修正形式。如θˆ恰好满足约束条件,即Aθ=b,则约束LSE和无约束lSE具有相同的估计量。
例8.8 约束信号 如果信号模型为
⎧θ1⎪
s[n]=⎨θ2
⎪⎩01
n=0n=1n=2
我们观测到的数据为{x[0],x[1],x[2]},那么观测矩阵为
⎡1H=⎢⎢0
⎢⎣0
10⎤⎥⎥⎥⎦
观测到的信号矢量
⎡θ1⎤⎥θs=Hθ=⎢⎢2⎥⎢⎣0⎥⎦
则无约束lSE为
⎡x[0]⎤T−1Tˆθ=(HH)Hx=⎢⎥⎣x[1]⎦
信号估计为
⎡x[0]⎤ˆ=⎢x[1]⎥ˆ=Hθs⎢⎥
⎢⎣0⎥⎦
现假设我们预先知道θ1
=θ2
.根据(8-50)式,有
−1]θ=0
[1
则
A=[1−1] b=0
根据(8-52)式,得到约束lSE为
ˆ=θˆ−(HTH)−1AT[A(HTH)−1Aθc
ˆ−AT(AAT)-1Aθˆ=θ
=[I−A⎡⎢=⎢⎢⎢⎣
T
T
ˆ−b)]-1(Aθ
(AA
T
ˆ)-1A]θ
1⎤
(x[0]+x[1])⎥2
⎥
1
(x[0]+x[1])⎥
⎥⎦2
约束信号估计为
⎡1⎤([0][1])xx+⎢2⎥⎢⎥1ˆ=⎢(x[0]+x[1])⎥ˆc=Hθsc
⎢2⎥⎢⎥0⎢⎥⎢⎥⎣⎦
因为θ1=θ2,正好取两个观测的平均。
8.7 非线性最小二乘估计
LS方法是通过使下式最小来估计模型参数θ
J=(x−s(θ))
T
(x−s(θ))
其中s(θ)是x的信号模型,显然与θ有关。在线性LSE问题中,信号表示为特殊的形式s(θ)=Hθ,从而引出简单的线性LSE。 然而,一般情况下,s(θ)不能表示为这种形式,而是θ的一个N维非线性函数。这种非线性LS问题的形式在统计学上称为非线性回归问题。非线性LS问题的求解必须基于迭代方法。
在讨论非线性LS问题求解方法之前,我们先讲述两个能降低这种问题复杂度的方法,即参数变换和参数分离方法。 8.7.1 参数变换方法
该方法中,首先寻找一个θ的一对一变换,从而得到新空间中的一个线性信号模型。为此,我们令
α=g(θ)
其中g是θ的一个p维函数,它的反函数存在。如果能找到这样一个g满足
s(θ(α))=s(g(α))=Hα
则信号模型与α呈线性关系。
-1
因此,很容易求得α的线性LSE,的非线性LSE可以通过下式
ˆ=g-1(αˆ) θ
ˆ=(HH)Hx。 求得。其中,α
例8.9 正弦信号参数估计 设一个正弦信号模型为
s[n]=Acos(2πf0n+φ)n=0,1,...,N-1
式中f0是已知的,幅度A和相位φ是待估计参量。则LSE在A和φ上通过使下式最小找到。
J(A)=∑(x[n]−Acos(2πf0n+φ))
n=0N−1
2
T−1T
这是一个非线性最小二乘问题。
然而,因为
Acos(2πf0n+φ)=Acosφcos2πf0n−Asinφsin2πf0n
α1=Acosφ
如果令 α=Asinφ
2
信号模型变为
s[n]=α1cos2πf0n−α2sin2πf0n
利用矩阵形式表示,上式为 s=Hα 其中
10⎡⎤
⎢cos2πf⎥sin2fπ00⎥H=⎢⎢⎥##
⎢⎥ cos2f(N1)sin2f(N1)ππ−−00⎣⎦
现在,它与新参数呈线性关系,α的LSE为
ˆ=(HTH)−1HTx α
ˆ,必须求出逆变换g-1(α
ˆ) 为求出θ
φ=arctan(
−α2
α1
)
则非线性的LSE为
⎡⎤
ˆ⎡A⎤⎢⎥ˆ=⎢⎥=θ⎢arctan(−α2⎥ˆφ⎢⎦⎥⎢⎣α1⎥⎣⎦
然而,只有一小部分非线性LSE问题能用本方法解决。
8.7.2 参数分离方法
假设正弦信号模型为
s[n]=Acos2πf0n n=0,1,...,N-1
如果A 和f0均为待估计参量,
J(A,f0)=∑(x[n]−Acos2πf0n)
n=0N−1
2
则是可分离的最小二乘问题。一般而言,可分离的信号模型具有如下的形式 s=H(α)β
⎡α⎤⎡(p-q)×1⎤θ=⎢⎥=⎢⎥×βq1其中, ⎣⎦⎣⎦
这里,H(α)是一个与α有关的N×q矩阵,模型与β呈线性关系但是与α呈非线性关系。因此,LS误差可以求关于β的最小值,于是化简成仅为α的函数。
J(α,β)=(x-H(α)β)T(x-H(α)β)
对于给定的α,使J达到最小的β为
ˆ=(HT(α)H(α))−1HT(α)xβ (8-53)
根据式(8-22),得到LS误差为
ˆ)=xT[I-H(α)(HT(α)H(α))−1HT(α)]xJ(α,β
则问题化为在上求下式的最大值
xTH(α)(HT(α)H(α))−1HT(α)x
例8.10 阻尼指数信号
假设有如下的信号模型
s[n]=A1rn+A2r2n+A3r3n n=0,1,...,N-1
其中未知参数为{A1,A2,A3,r}。已知0
β=[A1,A2,A3]T呈线性关系但是与α=r呈非线性关系。利用式
(8-54),通过在0
xTH(r)(HT(r)H(r))−1HT(r)x
其中
1⎡1
⎢r2
r
H(r)=⎢
⎢##⎢N−12(N−1)
r⎣r
ˆ
一旦求得r,就可以求出幅度的LSE
1⎤
r3⎥⎥#⎥3(N−1)⎥r⎦
ˆ=(HT(rˆ)H(rˆ))−1HT(rˆ)x β
8.8 小结
1. 最小二乘法是选择θ使(8-1)式最小的参量估计,式中的信号取决于θ。
2. 在8.2节提出了线性LS问题和非线性LS问题。
3. 一般的线性LS问题是使(8-9)式最小,从而导出(8-10)式的LSE和(8-11)-(8-13)式的LS误差的最小值。
4. 式(8-14)是加权最小二乘误差准则,由它产生了(8-16)式的加权LSE和(8-17)加权LS误差的最小值。
5. 在8.4节描述了LS的几何解释,从而导出了重要的正交原理。
6. 当矢量参量的维数未知时,或信号模型未知时,阶次递推最小二乘法是非常有用的。它是随着未知参量数目的增加,递推计算最小二乘估计。递推过程总结于(8-28)-(8-31)式。
7. 如果希望随着新的数据的增加及时修正LSE,那么可利用递推法,它是基于前次的估计和新的数据确定LSE。等式(8-46)-(8-48)给出了所需要的计算等式。
在前面几章里,是通过考虑无偏估计的类型和确定最小方差来寻找最优和接近最优估计。
现在我们抛开这些限制,研究另一类估计,这类估计一般不与最优性质相关,但对于很多估计问题却很有意义,这就是最小二乘估计。
最小二乘估计是相当古老的估计方法,早在1795年高斯就使用这种方法来研究行星的运动。
这种方法的突出特点是对观测数据不做任何概率或统计的描述,而仅仅假设一个数学模型,实现容易。所以,应用非常广泛。
8.1 最小二乘估计方法
前面我们确定一个好的估计的焦点就是寻找的估计是无偏的和有最小方差。在选择方差作为“好”的度量中,无疑地是试图使参量的估计值和真值之间的误差(平均)最小。
在最小二乘(Least Square,LS)方法中,是使观测数据和假设信号之间的平方误差最小。
LS法说明:假设取决于未知参量θ的模型产生信号s[n],这个信号s[n]是完全确定的信号,由于观测噪声或模型不准确,观测到的信号是受干扰的信号,用观测数据x[n]表示,θ的最小二乘估计(LSE)就是选择使s[n]最靠近观测数据x[n]的θ值。
图8-1 最小二乘估计方法
靠近程度用LS误差指标来度量
J(θ)=∑(x[n]−s[n]) n=0,1,",N−1
n=0N−12 (8-1)
式中J是θ的函数,是通过s[n]反映出来的。使J(θ)最小的θ就是LSE。注意,这里没有做出关于观测数据x[n]的概率假设。这种方法无论对高斯还是非高斯噪声同样有效。当然,LSE的性能无疑取决于干扰噪声以及模型的性质。
LSE通常应用于数据的准确统计特性未知,或不能找出最优估计,或最优估计在实际应用时太复杂的场合。
例8-1直流电平信号
假设图8-1中的信号是s[n]=A,观测数据为x[n](n=0,1,",N−1),那么,根据LS法,使(8-1)式或下式最小可估计A。
J(A)=N−1
n=02()x[n]−A ∑
上式对θ求导并令结果为零可得
1A=NN−1n=0∑x[n]=
这正是我们熟悉的样本均值估计,然而,这里没有声明在MVU意义上的最优,而仅仅是使LS误差最小。
原因是什么?
根据前面的讨论,如果x[n]=A+w[n],w[n]是零均值的WGN,那么,LSE也就是MVU估计,否则,LSE不是MVU估计。为了理解这一点,考虑噪声不是零均值情况,假设噪声为
w[n]=E(w[n])+w'[n]
式中w'[n]是零均值噪声,这时观测数据则为
x[n]=A+E(w[n])+w'[n]
显然,样本均值估计实际上是A+E(w[n])的估计,这时的LSE是有偏估计,不是MVU估计。
也就是说,必须假设观测数据是由确定信号和零均值噪声组成,在信号参量的正确选择下,误差e[n]=x[n]−s[n]在平均意义下趋于零,则(8-1)式的最小是合理的。
如果假设的直流量信号模型不正确,例如描述观测数据的模型是x[n]=A+
Bn+w[n],则这个模型误差也将引起LSE是有偏的。
例8-2 正弦信号的频率估计
考虑信号模型 s[n]=cos2πf0n
式中f0是待估计参量。f0的LSE通过使下式最小找到。
J(f0)=N−1
n=02()x[n]−cos2πfn ∑0
与直流量信号对比,在直流量信号中,容易找到使LS误差指标最小的参量估计。然而这里的LS误差不是f0的线性函数,不能找到使LS误差最小值的参量估计的闭合形式。
对于信号模型是待估计参量的线性函数的最小二乘情况我们称为线性最小二乘问题,如例8-1。
否则,这个问题称为非线性最小二乘问题,如例8-2。
求解非线性LS问题是通过栅格寻优或迭代最小化方法(前面已介绍)。
应该注意的是线性最小二乘问题是指信号是待估计参量的线性函数,信号本身不必是线性的。
例8-3 正弦信号幅度估计
考虑信号模型
s[n]=Acos2πf0n
式中f0是已知的,A是待估计参量。则LSE在A上通过使下式最小找到。
J(A)=∑(x[n]−Acos2πf0n)
n=0N−12
因为J(A)是A的二次型函数,通过求导可以容易求出最小值。这里信号是待估计参量的线性函数,而信号本身不是线性的。实际中,线性最小二乘问题是我们希望的。
如果A 和f0均为待估计参量,
J(A,f0)=∑(x[n]−Acos2πf0n)
n=0N−12
则是可分离的最小二乘问题。
8.2 线性最小二乘估计
8.2.1 标量最小二乘
对于标量参数估计,如应用线性最小二乘方法,假设
s[n]=θh[n] (8-2)
式中h[n]是已知序列,那么LS误差准则为
J(θN−1
n=0)=∑(x[n]−θh[n])2 (8-3)
N−1使上式最小可得LSE如下
θˆ=n=0
N−1
n=0∑x[n]h[n] (8-4)
h2∑[n]
将(8-4)式代入(8-3)式得LS误差的最小值为
Jmin =Jθˆ=()∑(x[n]−θˆh[n])(x[n]−θˆh[n])
n=0
N−1N−1
=
n=0∑x[n](x[n]−θˆh[n])−θˆ∑h[n](x[n]−θˆh[n]) n=0N−1
N−1
=
n=0∑x[n]−θˆ∑x[n]h[n] s
2n=0N−1
(8-5)
代入θˆ得
⎛N−1⎞⎜∑x[n]h[n]⎟⎜⎟N−1n=0⎝⎠=∑x2[n]−N−1n=0
n=02Jmin (8-6) ∑h2[n]
N−12显然,最小的LS误差是观测数据的原始能量或∑nx[n]减去信=0
号拟合产生的能量。
ˆ=,θ=A,对于例8-1,则h[n]=1,代入(8-4)式得A根据(8-5)
式有
N−1
Jmin=
n=0∑x2[n]−N2
如果观测数据是无噪声的,即x[n]=A,则Jmin=0,或者说,对于这些观测数据有一个理想的LS拟合。
另一方面,如果x[n]=A+w[n],且E(w2[n])>>A2,那么,∑n=0x2[n]N−1N>>2,则最小的LS误差为
N−1
Jmin≈
n=0∑x2[n]
与原始能量相差不多。LS误差的最小值总是满足下式
N−1
0≤Jmin≤
n=0∑x2[n] (8-7)
8.2.2 矢量最小二乘
上述标量LSE的结果可直接扩展到矢量参量LSE,且有很大的实际应用价值。假设矢量参量θ是p×1维的,信号s=[s[0]s[1]"s[N−1]]T是待估计参量的线性函数,假设
s=Hθ (8-8)
式中H是一个满秩的N×p(N>p)矩阵,H称为观测矩阵。这是一个线性模型,没有通常的噪声PDF假设。在第4章很多信号满足这个模型。θ的LSE可使下式最小容易找到,即
J(θ)=∑(x[n]−s[n]])=(x−Hθ)
n=0N−12T(x−Hθ) (8-9)
进一步得
J(θ)=xTx−xTHθ−θTHTx+θTHTHθ
=xTx−2xTHθ+θTHTHθ
(注意xTHθ是标量),J对θ的导数(或梯度)
∂J(θ)=−2HTx+2HTHθ ∂θ
令这个导数为零,得LSE
ˆ=HTHθ()−1HTx (8-10)
ˆ的方程HTHθ=HTx称为标准方程。H的满秩保证了用来求解θ
H
TH的逆阵存在。
我们发现,这个LSE与线性模型的有效估计以及BLUE有相
ˆ是对观测数据没有做任何假设同的函数形式,式(8-10)中的θ
下产生的估计,而对于BLUE要求E(x)=Hθ和Cx=σ2I(见第6章),对于线性模型的有效估计还另外要求x是高斯分布的(见第4章)。
如果LSE满足这些假设,则容易确定LSE的统计特性,也就是第4,6章已经给出的形式。否则,确定LSE的统计特性相当困难。
将(8-10)代入(8-9)式可得LS误差最小值为
ˆ=x−HθˆJmin=Jθ
T=⎛⎜x−HHH⎝
TT()()T(x−Hθˆ) ()−1⎛THx⎞⎟⎜x−HHH⎠⎝T−1TTT()−1HTx⎞⎟ ⎠T=xI−H(HH)H
TT−1)(I−H(HH)H)x ( =x(I−H(HH)H)x (8-11)−1T
上式推导中,利用了I−H(HTH)HT是幂等矩阵的性质,即A2=A。−1
Jmin也可写成另一种形式
Jmin=xx−xHHHHx (8-12)
ˆ (8-13)
=xTx−HθTT()T−1T()
8.2.3 加权最小二乘
线性LS问题可扩展到加权LS,即在(8-9)式中增加一个N×N维的正定加权矩阵W如下
J(θ)=(x−Hθ)TW(x−Hθ) (8-14)
=wi>0的对角矩阵,那在误差准则中引进加权因子是为了强调那些可靠观测数据的重要性。例如,在例8-1中,如果W是[W]ii
么,LS误差为 J(A)=N−1
n=0
2如果w[n]是零均值和具有方差σn的不相关噪声,那么,加权因子∑wn(x[n]−A)2
的合理选择是wn2=n。
这个选择将产生A的LSE为
N−1
2n=0σn
N−1
n=0ˆ=A∑x[n]1 (8-15) ∑σ2n
这正是我们熟悉的BLUE,这是由于w[n]是不相关的,所以有W=C−1。
对于加权LSE的一般形式容易求得
ˆ=HTWHθ()−1HTWx (8-16)
它的LS误差的最小值为
Jmin=xT(W−WH(HWH)T−1HTWx (8-17)
)
8.3 LSE的几何解释
下面用几何的观点来解释线性LS法,通过几何解释能更清楚地揭示这种方法的本质。
回顾一般的线性信号模型s=Hθ,如果用hi表示H的各列,则有
s=h1[h2⎡θ1⎤⎢θ⎥1⎥⎢"hp=⎢#⎥⎢⎥⎢⎦⎣θp⎥]∑θhii=1pi
因此,信号模型可以看作是“信号”矢量{h1,h2,",hp}的线性组合。
例8-4 傅里叶分析 在例4-2中,假设M
=1,则信号模型为
s[n]=acos2πf0n+bsin2πf0n n=0,1,",N−1
式中f0是已知频率,θ=[a
b]
T
是待估计量。写成矢量形式为
10⎤⎡s[0]⎤⎡
⎥a⎢s[1]⎥⎢cos2πfπsin2f00⎥⎡⎤⎥=⎢⎢
⎥ (8-18) ⎥⎢⎢#⎥⎢##b⎣⎦⎥⎥⎢⎢
⎣s[N−1]⎦⎣cos2πf0(N−1)sin2πf0(N−1)⎦
上式看出,H的两列分别由余弦和正弦序列组成,即
h1=[1cos2πf0"cos2πf0(N−1)]h2=[0sin2πf0"sin2πf0(N−1)]
T
T
所以有 s=
ah
1
+bh
2
由LS误差定义可得
J(θ)=(x−Hθ)(x−Hθ)
T
T
如果定义N×1维矢量ξ
=[ξ1ξ2"ξN]的欧几里几何长度为
ξ=
=
那么,LS误差也可写为
J(θ)=x−Hθ
2
=x−∑θihi (8-19)
i=1
p
p
2
因此,线性LS法就是使数据矢量x与信号矢量∑i=1θihi之间的距离平方最小,其中信号矢量必须是H各列的线性组合。
数据矢量可以位于N维空间RN中的任意位置上,由p
p
上,记为Sp(由于H的满秩确保其各列是线性独立的,因此由各列所确定的子空间是真正的p维)。现以N=3和p=2为例,对所有可能的θ1和θ2(假设−∞
2上,而数据x通常不在这个子空间上。
图8-2 在R3空间的线性最小二乘的几何解释
ˆ应该是x很明显,在子空间S2上且在几何意义上最靠近x的s
ˆ是x在S2上的正交投影。这意味着在S2上的分量,换句话说,s
ˆ与S2上的所有矢量一定是正交的。误差矢量x−s(如果xTy=0,
ˆ,我则定义在RN上的这两个矢量正交。)在图8-2中,为了确定s
们使用正交条件,也就是误差矢量与信号子空间正交,即
ˆ)⊥S2 (x−s
式中⊥表示正交或垂直。如果上式成立,则一定有
ˆ)⊥h2 ˆ)⊥h1 (x−s(x−s
那么,误差矢量也将与h1和h2任意线性组合正交,根据正交的定义,则有
ˆ)h1=0 (x−sˆ)h=0 (x−s
T
T2
ˆ=θ1h1+θ2h2,则 令s
(x−θ1h1−θ2h2)Th1
=0
(x−θ1h1−θ2h2)
其矩阵形式
T
h2=0
(x−Hθ)[h1
Th2]=0T
或 (x−Hθ)
T
H=
0T (8-20)
由此得LSE为
ˆ=HTHθ
()
−1
HTx
注意到如果e=x−Hθ表示误差矢量,那么由(8-20)式得
eTH=0T (8-21)
即误差矢量一定与H的各列正交,这就是著名的正交原理。实际上,这个误差表示了不能被信号模型描述的x的那部分分量。 根据图8-2,LS误差的最小值为
ˆJmin=x−s
2
ˆ=x−Hθ
2
ˆ=x−Hθ
()(x−Hθˆ)
T
代入(8-21)式的结果,则
ˆJmin=x−Hθ
()(x−Hθˆ)
T
ˆ−θˆTHTe =xTx−Hθˆ =xTx−xTHθ
()
=xTI−HHTH
(
()
−1
HTx
)
(8-22)
ˆ拟合或接近数据矢量x的综上,LS法可看作是用另一个矢量sˆ是在RN的p维子空间上矢量问题,其中x是RN空间的矢量,s
{h,h
1
2
ˆ
为x在子空间的这个问题的解就是选择s,",hp}的线性组合。
正交投影。
如果H的各列是标准正交的,那么
(HH)
T
−1
=(I)=I
−1
则
ˆ=HTHθ
()
−1
HTx=HTx
显然,不需要求逆阵。 例8-4 傅里叶分析(继续) 继续考虑例8-4问题,如果f0
意整数,由式(4-13)容易得出
k⎞⎛k⎞⎛
hh2=∑cos⎜2πn⎟sin⎜2πn⎟=0
N⎠⎝N⎠⎝n=0
T
1
N−1
=kN,k是k=1,2,",N2−1中的任
另外,
Th1h1=
N
2
hT2h2=
N
2
因此,h1和h2是正交的,但不是标准正交。
将上面3个等式进行组合,可表示为HTH=(N2)I,所以有
ˆ⎤⎡aˆθ=⎢ˆ⎥=HTH⎣b⎦
()
−1
HTx
=
2THx N
⎡2N−1⎛
x[n]cos⎜2π⎢N∑⎝=0
=⎢n
N−1
⎛⎢2x[n]sin⎜2π∑⎢⎝⎣Nn=0k⎞⎤
n⎟N⎠⎥⎥ k⎞⎥n⎟N⎠⎥⎦
如果用下式代替原来的信号
s[n]=a′
22k⎞k⎞⎛⎛
cos⎜2πn⎟+b′sin⎜2πn⎟ NN⎠NN⎠⎝⎝
那么,H的列矢量是标准正交(单位正交)的。
然而,一般情况下,H的各列不是正交的,所以,信号矢量估计为
ˆ=HHTHˆ=Hθs
()
−1
HTx
这个信号矢量的估计是x在p维子空间的正交投影。如果令N×N维矩阵P=H(HTH)下性质:
1. PT=P,所以P是对称阵。 2. P2=P,所以P是幂等阵。
−1
HT
,则P称为正交投影矩阵或投影阵,P有如
ˆ中恢投影阵是奇异矩阵,如果不是奇异的,那么x就能从s
复出来。这显然是不可能的,因为如果那样,许多x将会有相同的投影。
ˆ=(I−P)x是x在信号矢量子空间的互补子同样,误差矢量e=x−s
空间或正交子空间的投影,矩阵P⊥=I−P也是投影阵,容易验证,它满足上述投影阵的性质。
根据(8-22)式,LS误差的最小值为
Jmin=xT(I−P)x
=xTP⊥x =xTP⊥P⊥x
=P⊥x
2
T
这恰恰是
e
2
。
8.4 按阶递推最小二乘估计
前面介绍的LSE是在信号模型已知的情况下进行的,但有很多情况是信号模型未知,这时,必须假设信号模型。例如,图8-5a给出了实验观测数据,信号模型可以分别假设为
s1(t)=A1 s2(t)=A2+B2t
式中0≤t≤T。如果在时刻t=nΔ时对x(t)采样得x[n],其中Δ=1和
n=0,1,",N−1,则相应的离散信号模型为
s1[n]=A1
s2[n]=A2+B2n
根据信号模型有
⎡1⎤⎢1⎥H1=⎢⎥
⎢#⎥⎢⎥⎣1⎦
⎡10⎤⎢11⎥
⎥ H2=⎢
⎢##⎥⎢⎥
N1−1⎣⎦
利用LSE得截矩和斜率的估计为
ˆ= A1
(8-23)
和
N=1N=1()221N1−ˆ=Anx[n] ∑x[n]−NN+1∑2
NN+1n=0n=0
N=1N=1
112
nx[n] (8-24)
∑x[n]+NN+1∑NN+1n=0n=0
ˆ=−B2
图8-5用最小二乘模拟实验数据
ˆ和sˆ+Bˆt,ˆ1(t)=Aˆ2(t)=A图8-5b和图8-5c分别画出了s其中T=100。122
从图中明显看出,使用两个参量的拟合效果较好。
后面将看到,随着参量的增加,LS误差的最小值将减小。然而在实际应用中,在足以描述观测数据的情况下,应选择最简单的信号模型。
具体实现就是增加多项式的阶次直到增加阶次后使LS误差的最小值仅稍微减小为止。对于图8-5的观测数据,LS误差的最小值Jmin
与参量数目的关系示于图8-6。
图8-6 最小值Jmin与参量数目的关系
模型s1(t)和s2(t)分别对应于k=1和k=2。从图中可以看出,Jmin从一个参量变化到两个参量有大幅度的下降,当再增大阶次时,
Jmin减小很少。
在这个例子中,信号的实际模型是
s(t)=1+0.03t
噪声w[n]是σ2
=0.1的
WGN。如果A和B估计的很准确,那么,LS
误差的最小值为
ˆ,Bˆ=J(A,B)=∑(x[n]−s[n])=∑w2[n]≈Nσ2Jmin=JA
2
n=0
n=0
()
N−1N−1
所以,当达到实际的阶次时,Jmin≈10,从图8-6可验证该值,通过该值也可增加我们对所选择模型的可信度。
在前面的例子中我们看到,当需要确定几个信号模型的LSE时,一个直接的方法就是对每个模型使用(8-10)式
ˆ=HTHθ
()
−1
HTx
另一种方法就是使用按阶递推LS法,后种方法可以减少计算量。
按阶递推LS法是按阶次逐次递推更新LSE,即,它能够根据基于N×k维的H矩阵的LSE计算基于N×(k+1)维的H矩阵的LSE。为了看出其中的原因,我们假设信号模型为
s1[n]=A1 s2[n]=A2+B2n
其中,−M
−M⎤⎡1
⎢1−(M−1)⎥
⎥H2=⎢
⎢#⎥#⎢⎥1M⎣⎦
⎡2M+1
⎢=HTH22⎢0⎢⎣
⎤
⎥ M
2⎥n∑⎥n=−M⎦
容易求得LSE解为
M
1ˆ=Ax[n]∑1
2M+1n=−M
M
1ˆ=Ax[n] ∑ 2
2M+1n=−M
ˆ=−B2
n=−M
M
∑
M
nx[n]n
2
n=−M
∑
可以看出,在这种情况下,当模型中增加一个参数时,A的LSE不会发生变化,这是由H2TH
2的对角性质得出的结论。
h
i
的正交性使得能把x分别投影到h1和h
2
方向,然后把
结果相加。
一般而言,列矢量不是正交的,因此每个投影不是独立的,但能够利用另一组正交的p矢量来代替(Gram-Schmidt正交化)。新的列h
2
被投影到与h1正交的子空间h
'
2
上。那么,LSE
信号估计变为
ˆ=h1θˆ1+h2'θˆ2 s
其中h1θˆ1也是信号仅根据H=h1 的LSE。
因此,这种方法在信号模型中可以用递推的形式将每一项加到信号模型上来表示阶数的更新。
按阶递推LS法可经过纯粹的代数推导获得(附录8A),这里
ˆk表示基于Hk的仅给出最后结果。令Hk表示N×k维观测矩阵,θ
LSE,即
ˆk=HTθkHk
基于Hk的LS误差的最小值为
()
−1
HTkx (8-25)
ˆJmink=x−Hkθk
(ˆ))(x−Hθ (8-26)
T
kk
现将阶次增加为k+1,这将产生一个新的观测矩阵Hk+1,其分块矩阵形式如下
Hk+1=[Hkhk+1]=[N×kN×1] (8-27)
那么,参量的按阶递推LSE为
−1T⊥⎤TT⎡HHHhhPkkkk+1k+1kxˆ⎢θk−⎥T⊥⎡k×1⎤hPh⎢⎥k+1kk+1==⎢T⊥
⎢⎥⎣1×1⎥hk+1Pkx⎦ (8-28) ⎢⎥T⊥
hk+1Pkhk+1⎣⎦
−1
()
ˆk+1θ
T⊥
式中Pk⊥=I−Hk(HT这个kHk)Hk。Pk是x在子空间上的投影矩阵,
子空间与Hk的各列张成的子空间正交。 为了避开计算HTkHk的逆矩阵,令
Dk=HHk
并使用下面的递推公式
(
T
k
)
−1
(8-29)
T
⎡⎤DkHTDkHTkhk+1hk+1HkDkkhk+1
−T⊥⎢Dk+⎥⎡k×kk×1⎤T⊥
hPhhPhk+1kk+1k+1kk+1⎥=⎢Dk+1=⎢T (8-30) hk+1HkDk1⎢⎥⎣1×k1×1⎥⎦−T⊥
⎢⎥T⊥
hPhhPhk+1kk+1k+1kk+1⎦⎣
式中Pk⊥=I−HkDkHTk。
LS误差最小值的阶次递推式为
Jmink+1=Jmink
(h
−
PxT
hk+1Ph
Tk+1
2⊥k⊥kk+1
)
(8-31)
整套算法不需要求逆阵。递推的初值θˆ1、Jmin1和D1
分别由式(8-25)、(8-26)和(8-29)确定,式(8-28)、(8-30)和(8-31)称为按阶递推最小二乘。
下面将这个方法应用于前面的直线拟合举例中。 例8-5 直线拟合
由于s1[n]=A1和s2[n]=A2+B2n(n=0,1,",N−1),所以有
⎡1⎤
⎢1⎥
H1=⎢⎥=1
⎢#⎥⎢⎥⎣1⎦
H2=[H1h2]
⎡0⎤⎢1⎥
⎥。利用式(8-25)和(8-26)计算递推初值 式中h2=⎢
⎢#⎥⎢⎥N−1⎣⎦
−1TTˆˆA1=θ1=(H1H1)H1x=
ˆJ=x−Hθ11 min1
()(x−Hθˆ)=∑(x[n]−)
T
11
n=0
N−1
2
ˆBˆˆ2=A下一步,利用(8-28)计算θ22
[
],即
T
−1TTT⊥⎤⎡HHHhhP111221xˆ⎥⎢θ1−T⊥h2P1h2⎥ˆ2=⎢θ⊥ ⎥⎢hTPx21
⎥⎢T⊥
hPh212⎦⎣
()
其中所需要的项
(H
T
1
H1
)
−1
=
1N
T
P1⊥=I−H1H1H1
()
−1
T
=I−H1
1T
11 N
P1⊥x=x−
111N
T
x=x−
⊥TT
hTPx=hx−21221=
∑nx[n]−∑n
n=0
n=0
N−1
N−1
2
N−1N−1
1T21⎛⎞2hP
h2=hh2−h21=∑n−⎜∑n⎟
NN⎝n=0⎠n=0
T
2
⊥1
T2
()
N−1
⎡1N−1⎡N−1⎤⎤
nnxnn[]−∑⎢⎢∑⎥⎥Nn=0⎣n=0n=0⎦⎥⎢−
2N−1N−1⎥⎢1⎛⎞2N−1
1⎤⎡nn−⎥⎢⎜⎟∑∑ˆN⎝n=0⎠n=0⎥=⎢−N∑nB2⎥ˆ2=⎢θn=0N−1N−1代入得 ⎥ ⎥⎢⎢
ˆnx[n]−nB2⎥⎥⎢⎢∑⎦⎣
n=0n=0⎥⎢
2N−1⎥⎢1⎛N−1⎞2
n−⎜∑n⎟⎥⎢∑N⎝n=0⎠⎥⎢n=0⎦⎣
其中
ˆ=B2
∑∑
N−1n=0
N−1n=0
nx[n]−∑
N−1n=0
n
2
1⎛2n−
N⎜⎝
∑
N−1n=0
⎞
n⎟⎠
化简得
ˆ=B2
∑
N−1n=0
nx[n]−N
N
(N−1)
N
2
2
−12
N−1
6
x[n]+=−∑NN+1n=0N
N−1
12
nx[n] ∑2
N−1n=0
和
N−11N−1ˆˆˆA2=−∑nB2=−B2Nn=02
N−12(2N−1)N−16
x[n]−nx[n] =∑∑NN+1n=0NN+1n=0
这个结果与(
8-24)式相同。
根据(8-31)式LS误差的最小值为
2
Jmin2=Jmin1
h(−
N−1
T
2
Px)
⊥1
⊥
hTP21h2
N−1
2
=Jmin1
⎛⎞
[]nxn−n∑⎜∑⎟
n=0−n=0
2N−1N−1
1⎛⎞2
n−⎜∑n⎟∑N⎝n=0⎠n=0
N−1
N−1
=Jmin1
⎛⎞[]−nxnn∑⎜∑⎟
n=0⎠−⎝n=0
NN2−12
通常,递推过程是连续求解LS问题,直到获得理想阶次为止。
这样不仅获得了理想模型的LSE,同时,也获得了所有低阶次模型的LSE。正如前面所述,这对于模型阶次未知的情况很有用。
8.5 序贯(递推)最小二乘估计
在多数信号处理应用中,都是对连续时间波形进行采样而获得数据,显然,这些数据是按时间顺序获得的,当需要很多的数据时,我们有两种选择,一种是在获得所有有用的数据后一次处理,另一种就是这一节将介绍的,及时处理每一个数据。 特别是,假设已经确定了基于{x[0],x[1],",x[N−1]}的LSE
ˆ,θ现在又观测到了新的数据x[N],能否在不计算
LSE式—(8-10)
ˆ获得基于{x[0],x[1],",x[N]}的LSE? 式,通过修改θ
回答是肯定的。为了区别起见,这种方法称为序贯(递推)最小二乘(sequential least squares)法,前面的最小二乘法称为批量(batch)法。
8.5.1 标量递推最小二乘
仍考虑例8-1直流量信号的估计。已知基于当前观测到的数据的LSE为
ˆ[N−1]=1A
N
ˆ[N]=A
∑x[n]
n=0
N−1
现在又观测到了一个新的数据x[N],则这时的LSE应为
N
1
x[n] ∑N+1n=0
ˆ[N−1]和x[N]求出Aˆ[N]。 现在我们不重新计算这个求和,而通过A
ˆ[NA
1⎛N−1⎞+xnxN[]]=[]∑⎟N+1⎜⎝n=0⎠
N1ˆ[N−1]+=Ax[N] (8-35)N+1N+1ˆ[N−1]+1x[N]−Aˆ[N−1]=A
N+1 (8-36)
()
由此看出,新的估计等于原来的估计加上一个修正项。
ˆ[N−1]是基于很多的 修正项随着N的增大而减小,这反映了A
ˆ[N−1]可以看作数据,因此,应该给以更多的权重。同时,x[N]−A
是预测误差,如果这个误差是零,则不进行修正。否则,修正新的估计。
LS误差的最小值也可得出递推形式。基于截止N−1时刻数据的LS误差的最小值为
ˆ[N−1]2 Jmin[N−1]=∑x[n]−A
n=0N−1
((
)
增加新的数据后的这个误差
Jmin[N]=∑
n=0N
2
ˆx[n]−A[N]
)
整理后得
J
min[N]=Jmin[N−1]+
Nˆ[N−1]2 x[N]−AN+1
()
由此看出,增加新的数据使LS误差的最小值增大,这很容
易解释,这是因为对于每一个新的数据,都使平方误差项增加。
8.5.2 加权递推最小二乘
如果对线性最小二乘问题进行加权,假设加权矩阵W是
[W]ii
=2的对角阵,由(8-15)式得加权LSE为
ˆ[N−1]=A
∑
N−1n=0
N−1
x[n]
σσ
2n2n
∑
1
n=0
不难推出,其递推形式为
ˆ[N]=A
∑σ
n=0
Nn=0
N
x[n]
2n2n
∑σ
1
=
∑σ
n=0
N−1
x[n]
2nN
+1
x[N]
2σN
∑σ
n=0
⎛N−11⎞ˆx[N]
−AN1[]⎜∑2⎟2
σNn=0σn=+N
N
11
2
n
∑σ
n=0
2n
∑σ
n=0
2n
ˆ[N−1]−=A
σN
1ˆ
A[N−1]2
x[N]+
2
σN
∑σ
n=0
N
1
2n
∑σ
n=0
N
1
2n
1
ˆ[N−1]+=A
2
σN
∑σ
n=0
N
1
n
ˆ[N−1])(x[N]−A
(8-37)
正如所料,如果σn2=
σ2,则与前面的结果(8-36)式相同。
式(8-37)中修正项的增益因子取决于对新的采样值的可信
2
度,如果新的采样值是噪声或σN
→∞,则不改变前面的LSE,即
2ˆ[N]=Aˆ[N−1]。相反,如果新的采样值是无噪声的或σN那么,A→0,
ˆ[N]→x[N],即忽略所有前面的采样值。即采样因子表示对新的A
数据采样相对于前面采样值的可信度。
现在我们进一步解释这个结果。如果x[n]=A+w[n],w[n]是具有方差σn2的零均值不相关噪声。那么,根据例6-2,这个LSE就是BLUE,所以有
ˆ[N−1]=varA
()
1
∑
N−1n=0
1
σ
2
n
根据(8-37)式,第N次修正项的增益因子为
1
K[N]=
σN
1
1
2n
σ
n=0
N
=
σN
σ
2N
+
[N−1]varA
=
ˆ[N−1]varA
)
ˆ[N−1]+σ2 (8-38) varAN
(
)
ˆ[N−1]很大,则修正项大,否则,0≤K[N]≤1,由此看出,如果varA
ˆ[N−1]),可推导得出K[N]的递推形式:相反。由于K[N]取决于var(A
(
ˆ[NvarA
(])=
1
∑
N
1
=
1
n=0
σ
2
n
∑
N−1n=0
1
σ
2n
+
2N2N
1
=
1
[N−1]varA
σ
2
N
+
σ
2N
=
ˆ[N−1]σvarA
[N−1]+σvarA
(
)
ˆ[N−1]⎛varA
=⎜1−
[N−1]+
σ⎜varA
⎝
(
)
2N
⎞⎟var⎟⎠
(
ˆ[N−1]A
)
或
ˆ[N]=(1−K[N])varAˆ[N−1]varA (8-39)
()()
利用(8-38)和(8-39)式可递推K[N]。这样以来,在递推增益的过程中,同时也递推了LSE的方差。现将这些结果总结如下: 估计量的递推式:
ˆ[N]=Aˆ[N−1]+K[N]x[N]−Aˆ[N−1] (8-40) A
其中
ˆ[N−1]varA
K[N]=
N−1+σ2varA
()
(
)
(8-41)
N
方差的递推式:
ˆ[N]=(1−K[N])varAˆ[N−1]varA (8-42)
()()
具体递推过程如下: 1. 在开始递推时,令
ˆ[0]=x[0]A
ˆ[0]=σ2 varA0
()
ˆ[1]。ˆ[1]2. 由(8-41)式确定K[1],(8-40)式确定A下次迭代的varA
由(8-42)式确定。 3. 重复第2步。 图8-9示出了A=10和σ
n
2
仿真结果。
=1时的递推LS法的Monte-Carlo计算机
图8-9 WGN中的直流量的递推LS估计
图中看到,随着N的增大,方差和增益趋于零,这是由于
11ˆ[N])= var(A=N
∑
1N+1
n=0
σ
2
n
和
K[N]=
1=1
N+1
+1N
同时也看到,估计近似收敛于真值A=10。最后,给出LS误差的最小值的递推式
2
ˆx[N]−A[N−1]Jmin[N]=Jmin[N−1]+
N−1+σ2 (8-43)
varA
(
)
N
8.5.3 矢量加权递推最小二乘
现将获得的标量参量的递推LSE扩展到矢量参量(附录8C)。考虑W=C−1的加权LS误差准则,其中C表示零均值噪声的协方差矩阵,或
J=(x−Hθ)C−1(x−Hθ)
T
这与BLUE有同样的假设条件。根据(8-16)式这个加权LSE
T−1−1T−1ˆ为 θ=HCHHCx
()
ˆ的协方差,根由于C是噪声的协方差矩阵,所以我们用Cθˆ表示θ
据(6-17)式有
Cθˆ=HCH
(
T−1
)
−1
ˆ的递推 如果C是对角阵或噪声是不相关的,那么可以得到θ
LSE,否则,则不能。现假设满足条件,令
222
C[n]=diagσ0,σ1,",σn
()
⎡H[n−1]⎤⎡n×p⎤H[n]=⎢T=⎢⎥⎥ []nph1×⎣⎦⎣⎦
x[n]=[x[0]x[1]"n]T
ˆ的加权ˆ[n]表示基于x[n]或(n+1)个采样数据的θθ
LSE,那么,批量
LSE估计为
ˆ[n]=HT[n]C−1[n]H[n]θ
ˆ[n]的协方差阵为 θ
()
−1
HT[n]C−1[n]x[n] (8-44)
Cθˆ=Σ[n]=HT[n]C−1[n]H[n]
()
−1
(8-45)
再次强调,C[n]是噪声的协方差,∑[n]是LSE的协方差。
递推LSE为: 估计量的递推式:
ˆ[n]=θˆ[n−1]+K[n]x[n]−hT[n]θˆ[n−1] (8-46) θ
()
式中
K[n]=
σn2
Σ[n−1]h[n]
(8-47)
+hT[n]Σ[n−1]h[n]
方差的递推式:
Σ[n]=I−K[n]hT[n]Σ[n−1] (8-48)
()
这里的增益因子K[n]是p×1矢量,协方差阵Σ[n]是p×p维的。从上述等式中发现,矢量参量的递推LSE不需要求逆阵。
估计量的递推过程总结于图8-10,图中粗箭头表示矢量。递
ˆ[n−1]和推过程:首先使用批量估计(8-44)和(8-45)式确定θ
Σ[n−1],然后,对于n≥p利用递推等式(8-46)-(8-48)式进行
递推。需要指出的是,在计算 ˆ[n−1]θ和Σ[n−1]时,要求
HT[n−1]C−1[n−1]H[n−1]可逆。为保证这一点,H[n−1]的秩必须大于
或等于p,又由于H[n−1]是n×p维矩阵,所以必须n≥p(假设它的所有列是线性独立的)。
图8-10 序贯LSE估计
最后,给出LS误差的最小值的递推式
2
ˆx[n]−h[n]θ[n−1]Jmin[n]=Jmin[n−1]+2 (8-49)
σn+hT[n]Σ[n−1]h[n]
T
()
下面讨论关于LSE的应用问题。LSE通常在下述情况下使用: 1. 噪声的统计特性未知,或噪声的统计特性已知,但不能找到最优的MVU估计。
2. 对于已知噪声统计特性,但渐近最优的MLE实现起来太复杂。
8.6 约束最小二乘估计
有时,我们会遇到未知参数受到限制的LS问题。例如:我们希望估计几个信号的幅度,但预先知道这些信号的幅度中有部分是相等的。于是,利用先验知识,总的被估计参数的数目应当有所减少。对于这个例子,参数是线性相关的,从而引出线性约束的线性最小二乘估计问题。
假设参数θ受到r
Aθ=b (8-50)
其中A是一个r×p矩阵,而且已知b是一个r×1矢量。
举例来说,如果p=2而且已知一个参数是另一个参数的负数,于是约束可以表示为θ1+θ2=0。于是A=[1 1]和b=0。
矩阵A总是假设为满秩的(等于r),这对于约束条件相互独立是必要的。应当认识到在约束LS问题中,实际上只有p−r个独立参数。
为了求出线性约束的LSE,我们使用拉格朗日乘因子。即通过使如下的拉格朗日乘因子最小来确定θˆc(c表示约束LSE),
Jc=(x−Hθ)
T
T
x−Hθ+λ(Aθ−b) ()
其中λ是拉格朗日乘因子的r×1矢量。
展开这个表达式,得
Jc=xTx−2θTHTx+θTHTHx+λTAθ−λTb
取关于θ的梯度并利用(4-3)式得
∂Jc
Jc=−2HTx+2HTHx+ATλ
∂θ
令上式等于零,得
1TT−1ˆθc=(HH)Hx−(HTH)−1ATλ
2
ˆ−(HTH)−1ATλ/2=θ (8-51)
其中θ是无约束的LSE。
ˆ
为了求λ,我们利用(8-50)的约束条件,于是
ˆ=Aθˆ−A(HTH)−1ATλ/2=b Aθc
因此,
ˆ−b) λ/2=[A(HTH)−1AT]-1(Aθ
代入(8-51)式,得解为
ˆ=θˆ−(HTH)−1AT[A(HTH)−1AT]-1(Aθˆ−b)(8-52) θc
ˆ=(HTH)−1HTx。 其中,θ
约束LSE是无约束lSE的修正形式。如θˆ恰好满足约束条件,即Aθ=b,则约束LSE和无约束lSE具有相同的估计量。
例8.8 约束信号 如果信号模型为
⎧θ1⎪
s[n]=⎨θ2
⎪⎩01
n=0n=1n=2
我们观测到的数据为{x[0],x[1],x[2]},那么观测矩阵为
⎡1H=⎢⎢0
⎢⎣0
10⎤⎥⎥⎥⎦
观测到的信号矢量
⎡θ1⎤⎥θs=Hθ=⎢⎢2⎥⎢⎣0⎥⎦
则无约束lSE为
⎡x[0]⎤T−1Tˆθ=(HH)Hx=⎢⎥⎣x[1]⎦
信号估计为
⎡x[0]⎤ˆ=⎢x[1]⎥ˆ=Hθs⎢⎥
⎢⎣0⎥⎦
现假设我们预先知道θ1
=θ2
.根据(8-50)式,有
−1]θ=0
[1
则
A=[1−1] b=0
根据(8-52)式,得到约束lSE为
ˆ=θˆ−(HTH)−1AT[A(HTH)−1Aθc
ˆ−AT(AAT)-1Aθˆ=θ
=[I−A⎡⎢=⎢⎢⎢⎣
T
T
ˆ−b)]-1(Aθ
(AA
T
ˆ)-1A]θ
1⎤
(x[0]+x[1])⎥2
⎥
1
(x[0]+x[1])⎥
⎥⎦2
约束信号估计为
⎡1⎤([0][1])xx+⎢2⎥⎢⎥1ˆ=⎢(x[0]+x[1])⎥ˆc=Hθsc
⎢2⎥⎢⎥0⎢⎥⎢⎥⎣⎦
因为θ1=θ2,正好取两个观测的平均。
8.7 非线性最小二乘估计
LS方法是通过使下式最小来估计模型参数θ
J=(x−s(θ))
T
(x−s(θ))
其中s(θ)是x的信号模型,显然与θ有关。在线性LSE问题中,信号表示为特殊的形式s(θ)=Hθ,从而引出简单的线性LSE。 然而,一般情况下,s(θ)不能表示为这种形式,而是θ的一个N维非线性函数。这种非线性LS问题的形式在统计学上称为非线性回归问题。非线性LS问题的求解必须基于迭代方法。
在讨论非线性LS问题求解方法之前,我们先讲述两个能降低这种问题复杂度的方法,即参数变换和参数分离方法。 8.7.1 参数变换方法
该方法中,首先寻找一个θ的一对一变换,从而得到新空间中的一个线性信号模型。为此,我们令
α=g(θ)
其中g是θ的一个p维函数,它的反函数存在。如果能找到这样一个g满足
s(θ(α))=s(g(α))=Hα
则信号模型与α呈线性关系。
-1
因此,很容易求得α的线性LSE,的非线性LSE可以通过下式
ˆ=g-1(αˆ) θ
ˆ=(HH)Hx。 求得。其中,α
例8.9 正弦信号参数估计 设一个正弦信号模型为
s[n]=Acos(2πf0n+φ)n=0,1,...,N-1
式中f0是已知的,幅度A和相位φ是待估计参量。则LSE在A和φ上通过使下式最小找到。
J(A)=∑(x[n]−Acos(2πf0n+φ))
n=0N−1
2
T−1T
这是一个非线性最小二乘问题。
然而,因为
Acos(2πf0n+φ)=Acosφcos2πf0n−Asinφsin2πf0n
α1=Acosφ
如果令 α=Asinφ
2
信号模型变为
s[n]=α1cos2πf0n−α2sin2πf0n
利用矩阵形式表示,上式为 s=Hα 其中
10⎡⎤
⎢cos2πf⎥sin2fπ00⎥H=⎢⎢⎥##
⎢⎥ cos2f(N1)sin2f(N1)ππ−−00⎣⎦
现在,它与新参数呈线性关系,α的LSE为
ˆ=(HTH)−1HTx α
ˆ,必须求出逆变换g-1(α
ˆ) 为求出θ
φ=arctan(
−α2
α1
)
则非线性的LSE为
⎡⎤
ˆ⎡A⎤⎢⎥ˆ=⎢⎥=θ⎢arctan(−α2⎥ˆφ⎢⎦⎥⎢⎣α1⎥⎣⎦
然而,只有一小部分非线性LSE问题能用本方法解决。
8.7.2 参数分离方法
假设正弦信号模型为
s[n]=Acos2πf0n n=0,1,...,N-1
如果A 和f0均为待估计参量,
J(A,f0)=∑(x[n]−Acos2πf0n)
n=0N−1
2
则是可分离的最小二乘问题。一般而言,可分离的信号模型具有如下的形式 s=H(α)β
⎡α⎤⎡(p-q)×1⎤θ=⎢⎥=⎢⎥×βq1其中, ⎣⎦⎣⎦
这里,H(α)是一个与α有关的N×q矩阵,模型与β呈线性关系但是与α呈非线性关系。因此,LS误差可以求关于β的最小值,于是化简成仅为α的函数。
J(α,β)=(x-H(α)β)T(x-H(α)β)
对于给定的α,使J达到最小的β为
ˆ=(HT(α)H(α))−1HT(α)xβ (8-53)
根据式(8-22),得到LS误差为
ˆ)=xT[I-H(α)(HT(α)H(α))−1HT(α)]xJ(α,β
则问题化为在上求下式的最大值
xTH(α)(HT(α)H(α))−1HT(α)x
例8.10 阻尼指数信号
假设有如下的信号模型
s[n]=A1rn+A2r2n+A3r3n n=0,1,...,N-1
其中未知参数为{A1,A2,A3,r}。已知0
β=[A1,A2,A3]T呈线性关系但是与α=r呈非线性关系。利用式
(8-54),通过在0
xTH(r)(HT(r)H(r))−1HT(r)x
其中
1⎡1
⎢r2
r
H(r)=⎢
⎢##⎢N−12(N−1)
r⎣r
ˆ
一旦求得r,就可以求出幅度的LSE
1⎤
r3⎥⎥#⎥3(N−1)⎥r⎦
ˆ=(HT(rˆ)H(rˆ))−1HT(rˆ)x β
8.8 小结
1. 最小二乘法是选择θ使(8-1)式最小的参量估计,式中的信号取决于θ。
2. 在8.2节提出了线性LS问题和非线性LS问题。
3. 一般的线性LS问题是使(8-9)式最小,从而导出(8-10)式的LSE和(8-11)-(8-13)式的LS误差的最小值。
4. 式(8-14)是加权最小二乘误差准则,由它产生了(8-16)式的加权LSE和(8-17)加权LS误差的最小值。
5. 在8.4节描述了LS的几何解释,从而导出了重要的正交原理。
6. 当矢量参量的维数未知时,或信号模型未知时,阶次递推最小二乘法是非常有用的。它是随着未知参量数目的增加,递推计算最小二乘估计。递推过程总结于(8-28)-(8-31)式。
7. 如果希望随着新的数据的增加及时修正LSE,那么可利用递推法,它是基于前次的估计和新的数据确定LSE。等式(8-46)-(8-48)给出了所需要的计算等式。