自适应 Adaptive
文章目录
1. 最小二乘的两个难题
1.1 概述
我们从最小二乘开始。最小二乘有这样的基本模型。
假设我们有这样的数据
Z = ( z 1 , z 2 , . . . , z n ) T Z = (z_1,z_2,...,z_n)^T Z=(z1,z2,...,zn)T
我们希望用这组数据去估计我们的目标。我们希望能够找到一组数据,来极小化我们的目标函数。最终我们得到的解就是最小二乘解。
ω
=
(
ω
1
,
ω
2
,
.
.
.
,
ω
n
)
T
m
i
n
E
(
d
−
Z
T
ω
)
2
⇒
ω
L
S
\omega = (\omega_1,\omega_2,...,\omega_n)^T \\ min E(d-Z^T \omega)^2 \Rightarrow \omega_{LS}
ω=(ω1,ω2,...,ωn)TminE(d−ZTω)2⇒ωLS
其实我们得到最小二乘解的时候,没有遇到过多的问题。但是实际上,我们可能会面临两方面的困难。
首先我们要面对的问题是,目标函数可能具有自身的复杂性。因为我们是以均方误差作为目标函数的,并且我们采用的是线性估计。均方和线性使得我们在求解目标函数的时候没有遇到很多问题。我们也可以通过引入对角加载项的方法,让最小二乘进行正则化。但是,我们之前做的种种努力,都没有脱离均方和线性的范围。
但是实际情况是,我们的目标函数可能并不是均方误差。这就使得,我们的目标函数求梯度,可能没有解析解。
J ( ω ) ⇒ ∇ J ( ω ) = 0 ⇒ ω J(\omega) \Rightarrow \nabla J(\omega) = 0 \cancel \Rightarrow \omega J(ω)⇒∇J(ω)=0⇒ ω
其次,我们可能会面临非稳态的环境。这种非稳态来源于两部分。首先,可能我们的目标函数本身是在变的。其次,我们的数据也是时变的,其统计特性会随着时间的改变而发生改变。这样的话,我们就不能一劳永逸的求出估计了。我们就需要在每一个时刻,都进行重新求解。
一旦引入了这两个问题,我们之前做的努力,就有待商榷了。
1.2 Complexity of Objective Function
首先我们来看待第一个问题。如果目标函数得不得解析解怎么办。
这个时候我们可以引入数值优化的方法。
1.2.1 最速下降法
最经典的就是最速下降法。Steepest Descent
最速下降法是一种迭代的方法,我们从某一个初值开始,逐次进行迭代,每一次迭代都进行某种计算。我们希望我们的数据能够收敛到最优值。
ω 0 → ω 1 → ω 2 → . . . → ω n \omega_0 \rightarrow \omega_1 \rightarrow \omega_2 \rightarrow... \rightarrow \omega_n ω0→ω1→ω2→...→ωn
ω n → ω o p t ⇒ ∇ J ( ω o p t ) = 0 \omega_n \rightarrow \omega_{opt} \Rightarrow \nabla J(\omega_{opt}) = 0 ωn→ωopt⇒∇J(ωopt)=0
迭代的话,就会有数值方程,我们假设已经有了ωk了,我们想知道ωk和ωk+1之间有什么联系。
ω k + 1 = g ( ω k ) \omega_{k+1} = g(\omega_k) ωk+1=g(ωk)
我们可以通过泰勒展开来寻找这种关系
J ( ω k + 1 ) = J ( ω k ) + ( ∇ J ( ω k ) ) T ( ω k + 1 − ω k ) + o ( ∣ ∣ ω k + 1 − ω k ∣ ∣ 2 ) J(\omega_{k+1}) = J(\omega_k) + (\nabla J(\omega_k))^T(\omega_{k+1} - \omega_k) +o(||\omega_{k+1}- \omega_k||^2) J(ωk+1)=J(ωk)+(∇J(ωk))T(ωk+1−ωk)+o(∣∣ωk+1−ωk∣∣2)
我们假定高阶无穷小可以忽略,只考虑线性的部分。
J ( ω k + 1 ) = J ( ω k ) + ( ∇ J ( ω k ) ) T ( ω k + 1 − ω k ) J(\omega_{k+1}) = J(\omega_k) + (\nabla J(\omega_k))^T(\omega_{k+1} - \omega_k) J(ωk+1)=J(ωk)+(∇J(ωk))T(ωk+1−ωk)
我们想让这个式子尽可能的小,也就是
m i n ( ∇ J ( ω k ) ) T ( ω k + 1 − ω k ) min(\nabla J(\omega_k))^T(\omega_{k+1} - \omega_k) min(∇J(ωk))T(ωk+1−ωk)
根据柯西不等式等号成立条件
ω
k
+
1
−
ω
k
=
−
μ
∇
(
ω
k
)
\omega_{k+1} - \omega_k = -\mu \nabla(\omega_k)
ωk+1−ωk=−μ∇(ωk)
我们就到了迭代关系
ω
k
+
1
=
ω
k
−
μ
∇
(
ω
k
)
(
1
)
\omega_{k+1} = \omega_k-\mu \nabla(\omega_k) \quad\quad(1)
ωk+1=ωk−μ∇(ωk)(1)
1.2.2 牛顿法
最速下降法是一个一阶的优化方法,如果我们不满足于一阶的优化方法,我们想继续做二阶的。我们可以进行二阶展开
J ( ω k + 1 ) = J ( ω k ) + ( ∇ J ( ω k ) ) T ( ω k + 1 − ω k ) + 1 2 ( ω k + 1 − ω k ) T ∇ 2 J ( ω k ) ( ω k + 1 − ω k ) J(\omega_{k+1}) = J(\omega_k) + (\nabla J(\omega_k))^T(\omega_{k+1} - \omega_k) +\frac{1}{2} (\omega_{k+1} - \omega_k)^T \nabla^2 J(\omega_k) (\omega_{k+1} - \omega_k) J(ωk+1)=J(ωk)+(∇J(ωk))T(ωk+1−ωk)+21(ωk+1−ωk)T∇2J(ωk)(ωk+1−ωk)
其中二阶项是海森矩阵
∇ 2 J ( ω k ) ⇒ Hessian Matrix \nabla^2 J(\omega_k) \Rightarrow \text{ Hessian Matrix} ∇2J(ωk)⇒ Hessian Matrix
一阶展开只有线性的部分。二阶展开以后,我们得到了一个二次型。现在,如果我们想要让目标函数尽可能小,我们只要对后面的部分进行求导即可
m i n [ ( ∇ J ( ω k ) ) T ( ω k + 1 − ω k ) + 1 2 ( ω k + 1 − ω k ) T ∇ 2 J ( ω k ) ( ω k + 1 − ω k ) ] min[(\nabla J(\omega_k))^T(\omega_{k+1} - \omega_k) +\frac{1}{2} (\omega_{k+1} - \omega_k)^T \nabla^2 J(\omega_k) (\omega_{k+1} - \omega_k)] min[(∇J(ωk))T(ωk+1−ωk)+21(ωk+1−ωk)T∇2J(ωk)(ωk+1−ωk)]
∇ ω k + 1 ( ( ∇ J ( ω k ) ) T ( ω k + 1 − ω k ) + 1 2 ( ω k + 1 − ω k ) T ∇ 2 J ( ω k ) ( ω k + 1 − ω k ) ) = 0 ∇ J ( ω k ) + 2 ∗ 1 2 ∇ 2 J ( ω k ) ∗ ( ω k + 1 − ω k ) = 0 ω k + 1 = ω k − μ ∇ ( ∇ 2 J ( ω k ) ) − 1 J ( ω k ) ( 2 ) \nabla \omega_{k+1}((\nabla J(\omega_k))^T(\omega_{k+1} - \omega_k) +\frac{1}{2} (\omega_{k+1} - \omega_k)^T \nabla^2 J(\omega_k) (\omega_{k+1} - \omega_k)) =0\\ \nabla J(\omega_k) + 2*\frac{1}{2}\nabla^2 J(\omega_k)*(\omega_{k+1} - \omega_k) =0 \\ \omega_{k+1} = \omega_k-\mu \nabla (\nabla^2 J(\omega_k))^{-1}J(\omega_k) \quad\quad(2) ∇ωk+1((∇J(ωk))T(ωk+1−ωk)+21(ωk+1−ωk)T∇2J(ωk)(ωk+1−ωk))=0∇J(ωk)+2∗21∇2J(ωk)∗(ωk+1−ωk)=0ωk+1=ωk−μ∇(∇2J(ωk))−1J(ωk)(2)
这里需要注意一下
∇ x ( x T A x ) = ( A + A T ) x \nabla_x(x^T A x) = (A +A^T)x ∇x(xTAx)=(A+AT)x
这里求导出来了两倍是因为hessian矩阵是对称的,转置等于其本身。
这个二阶优化的方法叫做牛顿法。因此,如果我们的目标函数解不出解析解,我们就只能通过数值计算的方法进行逼近了。
1.3 Non-Stationary Environment
我们再来说说第二个问题,我们让下标k代表时间。现在我们的数据和目标都是随时间变化的了。并且求解的系数也是个时变的参数了。
k → T i m e Z ( k ) = ( z 1 ( k ) , . . . , z n ( k ) ) T d → d ( k ) ω → ω ( k ) k \rightarrow Time \\ Z(k) = (z_1(k),...,z_n(k))^T \\ d \rightarrow d(k) \\ \omega \rightarrow \omega(k) k→TimeZ(k)=(z1(k),...,zn(k))Td→d(k)ω→ω(k)
现在研究的对象全都是与时间相关的了,我们重写一下目标函数,然后进行求解就行
m i n E ( d ( k ) − Z T ( k ) ω ( k ) ) 2 ⇒ ω o p t ( k ) min E(d(k) - Z^T(k) \omega(k))^2 \Rightarrow \omega_{opt}(k) minE(d(k)−ZT(k)ω(k))2⇒ωopt(k)
2. 最小均方自适应滤波器LMS
2.1 问题引入
我们重新审视一下我们提出的有关最小二乘的两个问题。
假设我们的研究对象是非平稳的,如果优化函数还是均方的,我们要求目标函数可以简单的用最小二乘解。因为虽然ω是随着时间变化的,但是我们计算最小二乘解只需要一步就能求解完成。
但是,如果我们的目标函数不是均方的,很复杂,而且还是非平稳的,目标函数会与时间相关,这个时候怎么办呢?因为,首先,我们的最优解是一步解不出来的,需要通过不断的迭代。而且,最优解又是与时间相关的,会随时变化。即使,我们可以假设,最优解的时间变化是个慢过程,计算是个快过程,计算机计算性能足够高,能够在最优解状态改变之前得到结果。如果这个假设成立,会对计算机产生极大的压力,而且这样的算法也不具备普适性,没有办法得到广泛的应用。
因此,我们需要用其他办法来解决这个既不平稳,目标函数又复杂的问题。事实上,如果我们能够让迭代的尺度和时间的尺度联系起来,这个问题就能得到解决了。
这就是我们今天要思考的核心问题。
我们一开始做的事情是这样的。
在一个状态下,对采样得到的数据,通过不断的迭代,获得最优的参数估计。只要没找到最优解,我们的迭代就不会结束。然后进入下一个状态,继续采样,对采样数据不断迭代,得到新状态的最优估计参数,然后依次类推。
现在,我们不这么做了。我们每次获得一个采样点就迭代一次。采样一次以后,目标函数就变了,我们在新的目标函数下继续迭代一次,是一种阶梯状的迭代方式。这么做其实更有利于对变化的环境进行追踪,因为虽然我们没有在原系统中得到最优解,但是系统马上就发生改变了,得到最优解也没有什么意义。这就是今天要探讨的自适应问题。
2.2 LMS
2.2.1 模型建立
现在,我们要让我们的算法去适应环境的变化。我们写一下新的目标函数
J ( ω ( k ) ) = E ( d ( k ) − Z T ( k ) ω ( k ) ) 2 J(\omega(k)) = E(d(k) - Z^T(k)\omega (k))^2 J(ω(k))=E(d(k)−ZT(k)ω(k))2
我们来求一下最速下降
∇
ω
(
J
(
ω
(
k
)
)
)
=
−
2
Z
(
k
)
∗
(
d
(
k
)
−
Z
T
(
k
)
ω
(
k
)
)
=
−
2
Z
(
k
)
e
(
k
)
\nabla_{\omega}(J(\omega(k))) = -2Z(k)*(d(k) - Z^T(k)\omega (k)) = -2Z(k)e(k)
∇ω(J(ω(k)))=−2Z(k)∗(d(k)−ZT(k)ω(k))=−2Z(k)e(k)
其中e(k)表示估计误差
e ( k ) = d ( k ) − Z T ( k ) ω ( k ) ⇒ error e(k) = d(k) - Z^T(k)\omega (k) \Rightarrow \text{error} e(k)=d(k)−ZT(k)ω(k)⇒error
目前我们使用的符号k还都是表示时间,现在我们混用符号了。
ω
k
+
1
=
ω
k
−
μ
∇
J
(
ω
k
)
\omega_{k+1} = \omega_{k} - \mu \nabla J(\omega_k)
ωk+1=ωk−μ∇J(ωk)
现在我们要把计算得到的梯度进行代入
ω k + 1 = ω k − μ Z ( k ) e ( k ) \omega_{k+1} = \omega_{k} - \mu Z(k)e(k) \\ ωk+1=ωk−μZ(k)e(k)
注意,原本上,ω的角标k是迭代次数,其他的k表示时间,在经过我们的符号统一之后,我们给k赋予了时间的含义,也就是每次采样一次,就发生一次迭代,然后新的迭代是用来估计新状态的最优参数的。其实,我们每次采样就迭代一次的想法是正确的,因为每次采样之后,系统就发生改变了,就算我们把最优参数算出来,也没有多大的意义。因为系统马上就变了,约等于四九年入国军。
ω k + 1 = ω k − μ Z ( k ) ( d ( k ) − Z T ( k ) ω ( k ) ) ( i ) \omega_{k+1}= \omega_{k} - \mu Z(k)(d(k) - Z^T(k)\omega (k)) \quad\quad(i) ωk+1=ωk−μZ(k)(d(k)−ZT(k)ω(k))(i)
(i)式就是非常有名的LMS自适应滤波器。(LMS = Least Mean Square)
2.2.2 LMS细节剖析
其实自适应的目的就是为了把迭代和自身的发展实现统一。
这里其实需要注意一下,因为看起来有些地方有些问题。就是我们在求梯度的时候,并没有写成期望的形式
J ( ω ( k ) ) = E ( d ( k ) − Z T ( k ) ω ( k ) ) 2 J(\omega(k)) = E(d(k) - Z^T(k)\omega (k))^2 J(ω(k))=E(d(k)−ZT(k)ω(k))2
我们求梯度是这么写的,也就是没有加期望。
∇
ω
(
J
(
ω
(
k
)
)
)
=
−
2
Z
(
k
)
∗
(
d
(
k
)
−
Z
T
(
k
)
ω
(
k
)
)
=
−
2
Z
(
k
)
e
(
k
)
\nabla_{\omega}(J(\omega(k))) = -2Z(k)*(d(k) - Z^T(k)\omega (k)) = -2Z(k)e(k)
∇ω(J(ω(k)))=−2Z(k)∗(d(k)−ZT(k)ω(k))=−2Z(k)e(k)
这样是不是错了呢?实际上,并没有,因为我们在做LMS的时候,每次只采用1次,也就是说Z只代表了一个数据。一个数据是没有办法求期望的,因此就索性不写期望符号了。
这里我们使用了一个小技巧,就是令瞬时值代替了期望值。 instant value -> expectation value
LMS是个革命性的算法,里面有两个非常重要的思想
- LMS首先将迭代和时间的角标完成了统一,把采样和优化完成了统一
- 其次LMS使用瞬时值取代了平均值
后面人们发现,这两个做法是非常聪明的,这个算法也得到了非常广泛的应用。因为实现极其简单,并且效果出人意料的好。
这个算法来自与widrow,这个人做了第一个感知器,而且用模拟电路实现了。对神经网络的发展功不可没。而且这个人谦虚的说,自己不懂数学。
2.3 LMS与学习曲线
2.3.1 学习曲线的建立
所谓学习曲线,就是误差随着时间的变化规律,变化规律可能会有很多形式,比如线性的,指数的等等
我们希望知道LMS的误差曲线是什么情况。我们来计算两件事情。
( 1 ) J ( k ) = E ( d ( k ) − Z T ( k ) ω k ) 2 ( 2 ) ξ ( k ) = E ( ω k − ω 0 ) 2 (1) \quad\quad J(k) = E(d(k) - Z^T(k) \omega_k)^2 \\ (2) \quad\quad\quad\quad\quad\quad \xi(k) = E(\omega_k - \omega_0)^2 (1)J(k)=E(d(k)−ZT(k)ωk)2(2)ξ(k)=E(ωk−ω0)2
其中,我们假设
d
(
k
)
=
Z
T
(
k
)
ω
0
+
e
0
(
k
)
(
i
i
)
d(k) = Z^T(k) \omega_0 +e_0(k) \quad\quad(ii)
d(k)=ZT(k)ω0+e0(k)(ii)
e
0
(
k
)
:
White Noise
E
(
e
0
(
k
)
e
0
(
n
)
)
=
δ
k
n
e_0(k):\text{White Noise} \\ E(e_0(k)e_0(n)) = \delta_{kn}
e0(k):White NoiseE(e0(k)e0(n))=δkn
这是一个简化的模型,我们假设目标值d(k)就是最优估计和白噪声。
通常来说,我们的目标值d(k)应该具有这样的形式
d ( k ) = Z T ( k ) ω ( k ) d(k) = Z^T(k) \omega(k) d(k)=ZT(k)ω(k)
这个模型其实是比较难处理的,我们就认为线性系数有一个最优的值可以去追寻。但是,如果直接把线性系数认为是不变的显然不太合理,于是我们又通过增加白噪声的方式,增加了一个加性扰动。
因此我们的两个计算式子的计算意义如下
- 第一个式子是在计算估计值和目标值相差多少
- 第二个式子考察的就是,随着时间k的变化,我们的估计系数与最优的系数之间的距离是怎么减少的。这是一种跟踪的性能,我们跟踪系数,克服白噪声的影响,不断的逼近我们的最优系数。
2.3.2 估计系数误差的递推
我们假设对估计的系数和最优的系数之间的距离是e,并且k时刻之间的距离是e(k)。这个距离就是误差
e
(
k
)
=
ω
k
−
ω
0
(
i
i
i
)
e(k) = \omega_k - \omega_0 \quad\quad(iii)
e(k)=ωk−ω0(iii)
现在我们想寻求k+1时刻的距离,并建立递推关系
e ( k + 1 ) = ω k + 1 − ω 0 ( i v ) e(k+1) = \omega_{k+1} - \omega_0 \quad\quad(iv) e(k+1)=ωk+1−ω0(iv)
我们引入LMS,也就是把(i)式代入(iv)式
e ( k + 1 ) = ω k − μ Z ( k ) ( d ( k ) − Z T ( k ) ω k ) − ω 0 ( v ) e(k+1)= \omega_{k} - \mu Z(k)(d(k) - Z^T(k)\omega_k) - \omega_0 \quad\quad(v) e(k+1)=ωk−μZ(k)(d(k)−ZT(k)ωk)−ω0(v)
我们再引入d(k)的近似模型,就是把(ii)代入(v)
e ( k + 1 ) = ω k − ω 0 − μ Z ( k ) ( Z T ( k ) ω 0 + e 0 ( k ) − Z T ( k ) ω k ) = ω k − ω 0 + μ Z ( k ) Z ( k ) T ( ω k − ω 0 ) − μ Z ( k ) e 0 ( k ) = ( I + μ Z ( k ) Z ( k ) T ) ( ω k − ω 0 ) − μ Z ( k ) e 0 ( k ) ( v i ) e(k+1)= \omega_{k}- \omega_0 - \mu Z(k)( Z^T(k) \omega_0 +e_0(k) - Z^T(k)\omega_k) \\ =\omega_{k}- \omega_0 + \mu Z(k)Z(k)^T(\omega_k - \omega_0) - \mu Z(k)e_0(k) \\ = (I+\mu Z(k)Z(k)^T)(\omega_k - \omega_0)- \mu Z(k)e_0(k) \quad\quad(vi) e(k+1)=ωk−ω0−μZ(k)(ZT(k)ω0+e0(k)−ZT(k)ωk)=ωk−ω0+μZ(k)Z(k)T(ωk−ω0)−μZ(k)e0(k)=(I+μZ(k)Z(k)T)(ωk−ω0)−μZ(k)e0(k)(vi)
再把(iii)代入(vi)可得
e ( k + 1 ) = ( I + μ Z ( k ) Z ( k ) T ) e ( k ) − μ Z ( k ) e 0 ( k ) ( v i i ) e(k+1) = (I+\mu Z(k)Z(k)^T)e(k) - \mu Z(k)e_0(k) \quad\quad(vii) e(k+1)=(I+μZ(k)Z(k)T)e(k)−μZ(k)e0(k)(vii)
我们就得到了误差的递推式子(vii)
但是通过观测,我们发现这个式子中,每一个变量都是随机变量。如果递推系数本身是个确定值的话,递推比较容易进行,但是现在,递推系数也是个随机变量。我们希望能够把某些变量确定下来
当μ足够小的时候,有这样的式子成立
R ≈ Z ( k ) Z T ( k ) μ < < 1 R \approx Z(k)Z^T(k) \\ \mu << 1 R≈Z(k)ZT(k)μ<<1
这样做是有意义的。因为μ是步长,如果步长足够小,就可以对Z进行近似。
ω k + 1 = ω k − μ ∇ J \omega_{k+1} = \omega_k - \mu \nabla J ωk+1=ωk−μ∇J
加入近似条件以后(vii)就变成了
e ( k + 1 ) = ( I + μ R ) e ( k ) − μ Z ( k ) e 0 ( k ) e(k+1) = (I+ \mu R)e(k) - \mu Z(k) e_0(k) e(k+1)=(I+μR)e(k)−μZ(k)e0(k)
但是这个矩阵结构仍然非常复杂,我们可以对R做特征分解
R = Q T Λ Q R = Q^T \Lambda Q R=QTΛQ
e ( k + 1 ) = ( I + μ Q T Λ Q ) e ( k ) − μ Z ( k ) e 0 ( k ) e(k+1) = (I +\mu Q^T \Lambda Q)e(k) - \mu Z(k) e_0(k) \\ e(k+1)=(I+μQTΛQ)e(k)−μZ(k)e0(k)
因为特征矩阵是正交矩阵,满足
Q T ∗ Q = Q ∗ Q T = I Q^T*Q = Q*Q^T = I QT∗Q=Q∗QT=I
所以,我们又可以把式子写成
e ( k + 1 ) = Q T ( I + μ Λ ) Q e ( k ) − μ Z ( k ) e 0 ( k ) Q e ( k + 1 ) = ( I + μ Λ ) Q e ( k ) − μ Q Z ( k ) e 0 ( k ) [ 1 ] e(k+1) = Q^T (I +\mu\Lambda) Qe(k) - \mu Z(k) e_0(k) \\ Qe(k+1) = (I +\mu\Lambda)Qe(k) - \mu QZ(k) e_0(k) \quad\quad[1] e(k+1)=QT(I+μΛ)Qe(k)−μZ(k)e0(k)Qe(k+1)=(I+μΛ)Qe(k)−μQZ(k)e0(k)[1]
这里做一次变量代换,令
v ( k ) = Q e ( k ) ϕ ( k ) = Q Z ( k ) e 0 ( k ) [ 2 ] v(k) = Q e(k) \\ \phi(k) = QZ(k)e_0(k) \quad\quad[2] v(k)=Qe(k)ϕ(k)=QZ(k)e0(k)[2]
则[1]可以写成
v ( k + 1 ) = ( I + μ Λ ) v ( k ) − μ ϕ ( k ) [ 3 ] v(k+1) = (I +\mu \Lambda)v(k) - \mu \phi(k) \quad\quad[3] v(k+1)=(I+μΛ)v(k)−μϕ(k)[3]
由于v(k)和φ(k)可以写成列向量的形式
v ( k ) = ( v 1 ( k ) , . . . , v n ( k ) ) T ϕ ( k ) = ( ϕ 1 ( k ) , . . . , ϕ n ( k ) ) T v(k) = (v_1(k),...,v_n(k))^T \\ \phi(k) = (\phi_1(k),...,\phi_n(k))^T v(k)=(v1(k),...,vn(k))Tϕ(k)=(ϕ1(k),...,ϕn(k))T
所以有
v
i
(
k
+
1
)
=
(
1
+
μ
λ
i
)
v
i
(
k
)
−
μ
ϕ
i
(
k
)
[
4
]
v_i(k+1) = (1+ \mu \lambda_i)v_i(k) - \mu \phi_i(k) \quad\quad[4]
vi(k+1)=(1+μλi)vi(k)−μϕi(k)[4]
我们可以进行递推了
v
i
(
k
)
=
(
1
+
μ
λ
i
)
k
v
i
(
0
)
−
∑
j
=
0
k
−
1
(
1
+
μ
λ
i
)
k
−
1
−
j
(
μ
ϕ
i
(
j
)
)
[
5
]
v_i(k) = (1+ \mu \lambda_i)^k v_i(0) - \sum_{j=0}^{k-1}(1+ \mu \lambda_i)^{k-1-j}(\mu \phi_i(j)) \quad\quad[5]
vi(k)=(1+μλi)kvi(0)−j=0∑k−1(1+μλi)k−1−j(μϕi(j))[5]
这里的这个v是对我们的误差进行某种正交变换的结果,并不影响误差的模
我们发现误差随时间的规律是个指数变化,如果我们要让误差减小,必须让底数小于1。我们来计算一下这个步长的条件
∣ 1 + μ λ i ∣ < 1 |1+ \mu \lambda_i| <1 ∣1+μλi∣<1
⇒ − 1 < 1 + μ λ i < 1 ⇒ − 2 < μ λ i < 0 \Rightarrow -1<1+\mu\lambda_i <1 \\ \Rightarrow -2 < \mu \lambda_i <0 ⇒−1<1+μλi<1⇒−2<μλi<0
因为正交阵的特征值一定大于0,所以我们得到了
⇒ − 2 λ i < μ < 0 ∣ μ ∣ < 2 λ i , ∀ 1 ≤ i ≤ n ⇒ ∣ μ ∣ < 2 λ m a x \Rightarrow -\frac{2}{\lambda_i} <\mu <0 \\ |\mu| < \frac{2}{\lambda_i},\forall 1\leq i \leq n \\ \Rightarrow |\mu| < \frac{2}{\lambda_{max}} ⇒−λi2<μ<0∣μ∣<λi2,∀1≤i≤n⇒∣μ∣<λmax2
因此,μ一定小于最大的特征值分之2
2.3.3 估计目标误差的递推
我们研究的误差是两个。第一个是估计距离估计目标还有多远,第二个是估计系数距离真正的系数还有多远。我们刚才已经完成了第二个误差的计算,下面我们要开始计算第一个误差了。我们希望把这个两个误差能够联系在一起
d ( k ) − Z T ( k ) ω k = d ( k ) − Z T ( k ) ( ω k − ω 0 + ω 0 ) = d ( k ) − Z T ( k ) ω 0 − Z T e ( k ) = d 0 − Z T e ( k ) d(k) - Z^T(k) \omega_k = d(k) - Z^T(k)(\omega_k - \omega_0 + \omega_0) \\ = d(k)- Z^T(k) \omega_0 - Z^T e(k) \\ =d_0 - Z^T e(k) d(k)−ZT(k)ωk=d(k)−ZT(k)(ωk−ω0+ω0)=d(k)−ZT(k)ω0−ZTe(k)=d0−ZTe(k)
其中,d0表示最优的系数产生的估计误差,也就是白噪声e0(k)
d 0 = d ( k ) − Z T ( k ) ω 0 = e 0 ( k ) d_0 = d(k)- Z^T(k) \omega_0 = e_0(k) d0=d(k)−ZT(k)ω0=e0(k)
因此,目标误差可以写做
E ( d ( k ) − Z T ( k ) ω k ) 2 = E ( e 0 ( k ) − Z T ( k ) e ( k ) ) 2 = E [ e 0 ( k ) 2 ] − 2 ∗ E [ e 0 ( k ) ∗ Z T ( k ) ∗ e ( k ) ] + E [ Z T ( k ) e ( k ) ∗ e ( k ) T Z ( k ) ] < 1 > E(d(k) - Z^T(k) \omega_k)^2 = E(e_0(k) - Z^T(k) e(k))^2 \\ = E[e_0(k)^2] - 2 *E[e_0(k) * Z^T(k) * e(k)] + E[Z^T(k) e(k) * e(k)^T Z(k)] \quad\quad<1> E(d(k)−ZT(k)ωk)2=E(e0(k)−ZT(k)e(k))2=E[e0(k)2]−2∗E[e0(k)∗ZT(k)∗e(k)]+E[ZT(k)e(k)∗e(k)TZ(k)]<1>
我们发现这个式子可以写成三部分。
- 第一部分:是最优系数产生的误差,也就是最小误差,是个常数
- 第二部分:是交叉项,因为白噪声和后面的乘积都不相关,这个交叉项实际是0
- 第三部分:是我们需要认真研究的部分
重写一下式<1>
E ( d ( k ) − Z T ( k ) ω k ) 2 = J m i n + E [ Z T ( k ) e ( k ) ∗ e ( k ) T Z ( k ) ] < 2 > E(d(k) - Z^T(k) \omega_k)^2 = J_{min} + E[Z^T(k) e(k) * e(k)^T Z(k)] \quad\quad<2> E(d(k)−ZT(k)ωk)2=Jmin+E[ZT(k)e(k)∗e(k)TZ(k)]<2>
我们把第三项单独拿出来分析
E [ Z T ( k ) e ( k ) ∗ e ( k ) T Z ( k ) ] < 3 > E[Z^T(k) e(k) * e(k)^T Z(k)] \quad\quad<3> E[ZT(k)e(k)∗e(k)TZ(k)]<3>
我们发现,ZTe和eTZ都是个标量。标量乘以标量还是个标量。由于我们知道,标量的迹就是它本身,所以我们可以把式子变成迹的样子。
E [ Z T ( k ) e ( k ) ∗ e ( k ) T Z ( k ) ] = E [ T r ( Z T ( k ) e ( k ) ∗ e ( k ) T Z ( k ) ) ] < 4 > E[Z^T(k) e(k) * e(k)^T Z(k)] =E[Tr(Z^T(k) e(k) * e(k)^T Z(k))] \quad\quad<4> E[ZT(k)e(k)∗e(k)TZ(k)]=E[Tr(ZT(k)e(k)∗e(k)TZ(k))]<4>
矩阵的迹就是矩阵对角线主元之和,也等于特征值之和
T r ( A ) = ∑ k = 1 n A k k = ∑ k = 1 n λ k Tr(A)= \sum_{k=1}^n A_{kk} = \sum_{k=1}^n \lambda_{k} Tr(A)=k=1∑nAkk=k=1∑nλk
迹具有这样的性质,求迹的矩阵可以交换乘积顺序
T r ( A B ) = T r ( B A ) Tr(AB) = Tr(BA) Tr(AB)=Tr(BA)
因此<4>式可以进行这样的变形
E [ Z T ( k ) e ( k ) ∗ e ( k ) T Z ( k ) ] = E [ T r ( Z T ( k ) e ( k ) e ( k ) T ∗ Z ( k ) ) ] = E [ T r ( Z ( k ) ∗ Z T ( k ) e ( k ) e ( k ) T ) ] = T r [ E ( Z ( k ) Z T ( k ) ∗ e ( k ) e ( k ) T ) ] E[Z^T(k) e(k) * e(k)^T Z(k)] =E[Tr(Z^T(k) e(k) e(k)^T *Z(k))] \\ = E[Tr(Z(k)*Z^T(k) e(k) e(k)^T)] \\ = Tr[E(Z(k)Z^T(k) *e(k) e(k)^T)] E[ZT(k)e(k)∗e(k)TZ(k)]=E[Tr(ZT(k)e(k)e(k)T∗Z(k))]=E[Tr(Z(k)∗ZT(k)e(k)e(k)T)]=Tr[E(Z(k)ZT(k)∗e(k)e(k)T)]
相当于前三项作为一个整体,与第四项交换了乘积顺序。
下面我们要引入玄妙的一步。我们假设Z和e是独立的。也就是前后两部分能分开算。这个假设其实是比较离谱的,但是最终得到的结果却是比较贴近的。
E
[
Z
T
(
k
)
e
(
k
)
∗
e
(
k
)
T
Z
(
k
)
]
=
T
r
[
E
(
Z
(
k
)
Z
T
(
k
)
∗
e
(
k
)
e
(
k
)
T
)
]
=
T
r
[
E
(
(
Z
(
k
)
Z
T
(
k
)
)
∗
E
(
e
(
k
)
e
(
k
)
T
)
]
=
T
r
(
R
E
(
k
)
)
<
5
>
E[Z^T(k) e(k) * e(k)^T Z(k)] = Tr[E(Z(k)Z^T(k) *e(k) e(k)^T)] \\ = Tr[E((Z(k)Z^T(k))*E(e(k) e(k)^T)] \\ = Tr(RE(k)) \quad\quad <5>
E[ZT(k)e(k)∗e(k)TZ(k)]=Tr[E(Z(k)ZT(k)∗e(k)e(k)T)]=Tr[E((Z(k)ZT(k))∗E(e(k)e(k)T)]=Tr(RE(k))<5>
其中
R
=
E
(
(
Z
(
k
)
Z
T
(
k
)
)
E
(
k
)
=
E
(
e
(
k
)
e
(
k
)
T
)
R=E((Z(k)Z^T(k)) \\ E(k) = E(e(k) e(k)^T)
R=E((Z(k)ZT(k))E(k)=E(e(k)e(k)T)
把<5>代入<2>中得
J ( k ) = J m i n + T r ( R E ( k ) ) < 6 > J(k) = J_{min} + Tr(RE(k)) \quad\quad<6> J(k)=Jmin+Tr(RE(k))<6>
我们希望把J和v联系在一起
T r ( R E ( k ) ) = T r [ E ( Z ( k ) Z T ( k ) e ( k ) ∗ e ( k ) T ) ] = T r [ E ( e ( k ) T ∗ Z ( k ) Z T ( k ) e ( k ) ) ] = T r [ E ( e ( k ) T ∗ E ( Z ( k ) Z T ( k ) ) ∗ e ( k ) ) ] = T r [ E ( e ( k ) T ∗ R ∗ e ( k ) ) ] = T r [ E ( e ( k ) T ∗ Q T Λ Q ∗ e ( k ) ) ] = T r [ E ( v T ( k ) ∗ Λ ∗ v ( k ) ) ] = T r [ E ( ∑ i = 1 k ( λ i ∗ ∣ v i ( k ) ∣ 2 ) ) ] Tr(RE(k)) = Tr[E(Z(k)Z^T(k)e(k) *e(k)^T)] \\ = Tr[E(e(k)^T*Z(k)Z^T(k)e(k))] \\ =Tr[E(e(k)^T*E(Z(k)Z^T(k))*e(k))] \\ = Tr[E(e(k)^T *R *e(k))] \\ =Tr[E(e(k)^T *Q^T \Lambda Q *e(k))] \\ = Tr [E(v^T(k)*\Lambda* v(k))] \\ = Tr[E(\sum_{i=1}^k(\lambda_i* |v_i(k)|^2))] Tr(RE(k))=Tr[E(Z(k)ZT(k)e(k)∗e(k)T)]=Tr[E(e(k)T∗Z(k)ZT(k)e(k))]=Tr[E(e(k)T∗E(Z(k)ZT(k))∗e(k))]=Tr[E(e(k)T∗R∗e(k))]=Tr[E(e(k)T∗QTΛQ∗e(k))]=Tr[E(vT(k)∗Λ∗v(k))]=Tr[E(i=1∑k(λi∗∣vi(k)∣2))]
因为已经是标量相加的结果了,迹可以去掉了
T r [ E ( ∑ i = 1 k ( λ i ∗ ∣ v i ( k ) ∣ 2 ) ) ] = ∑ i = 1 k ( λ i ∗ E ( ∣ v i ( k ) ∣ 2 ) ) < 7 > Tr[E(\sum_{i=1}^k(\lambda_i* |v_i(k)|^2))] = \sum_{i=1}^k(\lambda_i* E(|v_i(k)|^2)) \quad\quad<7> Tr[E(i=1∑k(λi∗∣vi(k)∣2))]=i=1∑k(λi∗E(∣vi(k)∣2))<7>
<7>代入<6>中
J ( k ) = J m i n + T r ( R E ( k ) ) = J m i n + ∑ i = 1 k ( λ i ∗ E ( ∣ v i ( k ) ∣ 2 ) ) < 8 > J(k) = J_{min}+ Tr(RE(k)) = J_{min} + \sum_{i=1}^k(\lambda_i* E(|v_i(k)|^2)) \quad\quad<8> J(k)=Jmin+Tr(RE(k))=Jmin+i=1∑k(λi∗E(∣vi(k)∣2))<8>
我们对v继续进行处理,代入[5]式的结果,因为我们在计算系数误差的时候已经得到了v的具体表达式了
E ( ∣ v i ( k ) ∣ 2 ) = E [ ( 1 + μ λ i ) k v i ( 0 ) − ∑ j = 0 k − 1 ( 1 + μ λ i ) k − 1 − j ( μ ϕ i ( j ) ) ] 2 = ( 1 + μ λ i ) 2 k ∣ v i ( 0 ) ∣ 2 − 2 ∗ ( 1 + μ λ i ) k v i ( 0 ) ∗ E [ ∑ j = 0 k − 1 ( 1 + μ λ i ) k − 1 − j ( μ ϕ i ( j ) ) ] + E [ ∑ j = 0 k − 1 ∑ s = 0 k − 1 ( 1 + μ λ i ) k − 1 − j ( 1 + μ λ i ) k − 1 − s ( μ ϕ i ( j ) ) ( μ ϕ i ( s ) ) ] < 9 > E(|v_i(k)|^2) = E[(1+ \mu \lambda_i)^k v_i(0) - \sum_{j=0}^{k-1}(1+ \mu \lambda_i)^{k-1-j}(\mu \phi_i(j)) ]^2 \\ = (1+ \mu \lambda_i)^{2k} |v_i(0)|^2 - 2*(1+ \mu \lambda_i)^k v_i(0)*E[\sum_{j=0}^{k-1}(1+ \mu \lambda_i)^{k-1-j}(\mu \phi_i(j))] \\+ E[\sum_{j=0}^{k-1}\sum_{s=0}^{k-1}(1+ \mu \lambda_i)^{k-1-j}(1+ \mu \lambda_i)^{k-1-s}(\mu \phi_i(j))(\mu \phi_i(s))] \quad\quad<9> E(∣vi(k)∣2)=E[(1+μλi)kvi(0)−j=0∑k−1(1+μλi)k−1−j(μϕi(j))]2=(1+μλi)2k∣vi(0)∣2−2∗(1+μλi)kvi(0)∗E[j=0∑k−1(1+μλi)k−1−j(μϕi(j))]+E[j=0∑k−1s=0∑k−1(1+μλi)k−1−j(1+μλi)k−1−s(μϕi(j))(μϕi(s))]<9>
首先,我们来看一下交叉项,由于e0与其他的东西不相关,求期望是0,所以交叉项是0
E ( μ ϕ i ( j ) ) = E ( Q μ Z ( j ) e 0 ( j ) ) = 0 < 10 > E(\mu \phi_i(j)) = E(Q \mu Z(j)e_0(j)) = 0 \quad\quad<10> E(μϕi(j))=E(QμZ(j)e0(j))=0<10>
后面部分中
E
(
ϕ
i
(
j
)
ϕ
i
(
s
)
)
=
0
,
j
=
s
E(\phi_i(j) \phi_i(s)) = 0, j \cancel= s
E(ϕi(j)ϕi(s))=0,j=
s
因此
E ( μ 2 ϕ i ( j ) ϕ i ( s ) ) = μ 2 Q T ∗ Q E ( e 0 ( j ) ) 2 E ( Z i ( j ) ) 2 = μ 2 E ( e 0 ( j ) ) 2 E ( Z i ( j ) ) 2 < 11 > E(\mu^2\phi_i(j) \phi_i(s)) = \mu^2 Q^T*Q E(e_0(j))^2 E(Z_i(j))^2 \\ = \mu^2 E(e_0(j))^2 E(Z_i(j))^2 \quad\quad<11> E(μ2ϕi(j)ϕi(s))=μ2QT∗QE(e0(j))2E(Zi(j))2=μ2E(e0(j))2E(Zi(j))2<11>
这里面由于Q是正交阵,乘积为1
把<10>和<11>代入 <9>中
E ( ∣ v i ( k ) ∣ 2 ) = ( 1 + μ λ i ) 2 k ∣ v i ( 0 ) ∣ 2 + ∑ j = 0 k − 1 ( 1 + μ λ i ) 2 ( k − 1 − j ) ∗ μ 2 E ( e 0 ( j ) ) 2 E ( Z i ( j ) ) 2 < 12 > E(|v_i(k)|^2) = (1+ \mu \lambda_i)^{2k} |v_i(0)|^2 + \sum_{j=0}^{k-1}(1+ \mu \lambda_i)^{2(k-1-j)}*\mu^2 E(e_0(j))^2 E(Z_i(j))^2 \quad\quad<12> E(∣vi(k)∣2)=(1+μλi)2k∣vi(0)∣2+j=0∑k−1(1+μλi)2(k−1−j)∗μ2E(e0(j))2E(Zi(j))2<12>
其中
J
m
i
n
=
E
[
(
e
0
(
j
)
)
2
]
<
13
>
J_{min} = E[(e_0(j))^2] \quad\quad<13>
Jmin=E[(e0(j))2]<13>
我们继续做个假设,这个假设影响不是非常大,就是结果看起来比较简洁
E [ ( Z i ( j ) ) 2 ] = λ i < 14 > E[(Z_i(j))^2] = \lambda_i \quad\quad<14> E[(Zi(j))2]=λi<14>
把<14>和<13>代入<12>中可得
E ( ∣ v i ( k ) ∣ 2 ) = ( 1 + μ λ i ) 2 k ∣ v i ( 0 ) ∣ 2 + μ 2 J m i n λ i ∑ j = 0 k − 1 ( 1 + μ λ i ) 2 ( k − 1 − j ) < 15 > E(|v_i(k)|^2) = (1+ \mu \lambda_i)^{2k} |v_i(0)|^2 +\mu^2 J_{min} \lambda_i \sum_{j=0}^{k-1}(1+ \mu \lambda_i)^{2(k-1-j)} \quad\quad<15> E(∣vi(k)∣2)=(1+μλi)2k∣vi(0)∣2+μ2Jminλij=0∑k−1(1+μλi)2(k−1−j)<15>
我们发现<15>后面部分是个等比数列求和问题,因此,我们可以把<15>继续写成
E ( ∣ v i ( k ) ∣ 2 ) = ( 1 + μ λ i ) 2 k ∣ v i ( 0 ) ∣ 2 + μ 2 ∗ λ i ∗ J m i n ∗ 1 − ( 1 + μ λ i ) 2 k 1 − ( 1 + μ λ i ) 2 = ( 1 + μ λ i ) 2 k [ ∣ v i ( 0 ) ∣ 2 − μ 2 ∗ λ i ∗ J m i n 1 − ( 1 + μ λ i ) 2 ] + μ 2 ∗ λ i ∗ J m i n 1 − ( 1 + μ λ i ) 2 < 16 > E(|v_i(k)|^2) = (1+ \mu \lambda_i)^{2k} |v_i(0)|^2 +\mu^2 *\lambda_i*J_{min} * \frac{1-(1+\mu \lambda_i)^{2k}}{1-(1+\mu \lambda_i)^{2}} \\ = (1+ \mu \lambda_i)^{2k}[|v_i(0)|^2 - \frac{\mu^2 *\lambda_i*J_{min}}{1-(1+\mu \lambda_i)^{2}}] + \frac{\mu^2 *\lambda_i*J_{min}}{1-(1+\mu \lambda_i)^{2}} \quad\quad<16> E(∣vi(k)∣2)=(1+μλi)2k∣vi(0)∣2+μ2∗λi∗Jmin∗1−(1+μλi)21−(1+μλi)2k=(1+μλi)2k[∣vi(0)∣2−1−(1+μλi)2μ2∗λi∗Jmin]+1−(1+μλi)2μ2∗λi∗Jmin<16>
其中
μ 2 ∗ λ i ∗ J m i n 1 − ( 1 + μ λ i ) 2 = μ 2 ∗ λ i ∗ J m i n − 2 μ λ i − μ 2 λ i 2 = μ J m i n − 2 − μ λ i < 17 > \frac{\mu^2 *\lambda_i*J_{min}}{1-(1+\mu \lambda_i)^{2}} = \frac{\mu^2 *\lambda_i*J_{min}}{-2\mu\lambda_i - \mu^2 \lambda_i^2} = \frac{\mu J_{min}}{-2-\mu \lambda_i } \quad\quad<17> 1−(1+μλi)2μ2∗λi∗Jmin=−2μλi−μ2λi2μ2∗λi∗Jmin=−2−μλiμJmin<17>
<17>代入<16>得
E ( ∣ v i ( k ) ∣ 2 ) = [ ∣ v i ( 0 ) ∣ 2 + μ J m i n 2 + μ λ i ] ( 1 + μ λ i ) 2 k − μ J m i n 2 + μ λ i E(|v_i(k)|^2) = [|v_i(0)|^2 +\frac{\mu J_{min}}{2+\mu \lambda_i } ] (1+ \mu \lambda_i)^{2k} - \frac{\mu J_{min}}{2+\mu \lambda_i} E(∣vi(k)∣2)=[∣vi(0)∣2+2+μλiμJmin](1+μλi)2k−2+μλiμJmin
最终我们得到了对目标估计的学习曲线方程
{ E ( ∣ v i ( k ) ∣ 2 ) = [ ∣ v i ( 0 ) ∣ 2 + μ J m i n 2 + μ λ i ] ( 1 + μ λ i ) 2 k − μ J m i n 2 + μ λ i J ( k ) = J m i n + ∑ i = 1 k ( λ i ∗ E ( ∣ v i ( k ) ∣ 2 ) ) \begin{cases} E(|v_i(k)|^2) = [|v_i(0)|^2 +\frac{\mu J_{min}}{2+\mu \lambda_i } ] (1+ \mu \lambda_i)^{2k} - \frac{\mu J_{min}}{2+\mu \lambda_i} \\ J(k) = J_{min} + \sum_{i=1}^k(\lambda_i* E(|v_i(k)|^2)) \end{cases} {E(∣vi(k)∣2)=[∣vi(0)∣2+2+μλiμJmin](1+μλi)2k−2+μλiμJminJ(k)=Jmin+∑i=1k(λi∗E(∣vi(k)∣2))
3. 小结
我们为了求得LMS自适应滤波器的学习曲线,做了很多很多的假设,看似LMS非常不靠谱,因为做的假设往往很离谱。但是实际上,LMS效果非常好,并且学习曲线与我们做了很多假设的理论计算基本吻合。
自适应这个思想,成功解决了迭代优化和非稳态两件事情,实现了采用时间和迭代时间两个时间尺度的统一。我们要构造一个自适应算法可能并不是很难。但是要注意,我们做自适应算法的时候,要做两件事情,首先要做仿真,然后要做理论上的分析,计算这个自适应算法的学习曲线。在计算学习曲线的时候,要敢于假设,把最后学习曲线的计算公式变得非常简洁漂亮才行。 算法+仿真+数学分析才是一个完整的研究。