【自适应滤波】LMS算法

现在我们需要处理的是一个变化的环境（非平稳，Non-Stationary）的目标估计：X1,X2,…,Xn→θn(n:Time)X_1,X_2,…,X_n o heta_n quad (n: ext{Time})X1​,X2​,…,Xn​→θn​(n:Time)，我们面临两个问题：

数据能否有效地描绘变化的目标？Metric的定义也随时间发生变化： E∣θ^n−θn∣2=εnmathbb E|hat heta_n- heta_n|^2=varepsilon_nE∣θ^n​−θn​∣2=εn​

采用优化求解，优化通常使用的手段是迭代：x0→x1→x2→…→xn→n→∞xopt=arg⁡min⁡xf(x)x_0 o x_1 o x_2 o … o x_nstackrel{n oinfty}{ o}x_{opt}=argmin_{x}f(x)x0​→x1​→x2​→…→xn​→n→∞xopt​=argminx​f(x)

而迭代通常使用的方式是用泰勒展开来观察函数局部：(One Order): f(x+Δ)=f(x)+∇f(x)T⋅Δ+O(∣Δ∣) ext{(One Order): }f(x+Delta)=f(x)+
abla f(x)^T·Delta+O(|Delta|)(One Order): f(x+Δ)=f(x)+∇f(x)T⋅Δ+O(∣Δ∣)

ΔDeltaΔ决定了迭代搜索的方向，而我们希望朝着使f(x)f(x)f(x)最小的方向：min⁡Δ(∇f(x)T⋅Δ)min_{Delta}(
abla f(x)^T·Delta)minΔ​(∇f(x)T⋅Δ)

柯西不等式告诉我们：Δ=−∇f(x)Delta=-
abla f(x)Δ=−∇f(x)

柯西不等式的核心结论之一是：向量内积的取值范围被两个向量的模长约束，即

于是我们找到了一阶泰勒展开下函数局部下降最快的方向：xk+1=xk−∇f(xk)x_{k+1}=x_k-
abla f(x_k)xk+1​=xk​−∇f(xk​)，该方法称之为最速下降（Steepest Descent）或者梯度法（Gradient）。

有了一阶，自然还可以展开成二阶：(Second Order): f(x+Δ)=f(x)+∇f(x)T⋅Δ+12ΔT∇2f(x)Δ+O(∣∣Δ∣∣2),（∇2f是海森矩阵） ext{(Second Order): }f(x+Delta)=f(x)+
abla f(x)^T·Delta+frac{1}{2}Delta^T
abla^2f(x)Delta+O(||Delta||^2), （
abla^2f是海森矩阵）(Second Order): f(x+Δ)=f(x)+∇f(x)T⋅Δ+21​ΔT∇2f(x)Δ+O(∣∣Δ∣∣2),（∇2f是海森矩阵）

找到下降最快的方向：min⁡Δ(∇f(x)T⋅Δ+12ΔT∇2f(x)Δ)  ⟹  求梯度∇f(x)+∇2f(x)Δ=0  ⟹  Δ=−(∇2f(x))−1∇f(x)min_{Delta}(
abla f(x)^T·Delta+frac{1}{2}Delta^T
abla^2f(x)Delta)stackrel{求梯度}implies
abla f(x)+
abla^2f(x)Delta=0impliesDelta=-(
abla^2 f(x))^{-1}
abla f(x)minΔ​(∇f(x)T⋅Δ+21​ΔT∇2f(x)Δ)⟹求梯度​∇f(x)+∇2f(x)Δ=0⟹Δ=−(∇2f(x))−1∇f(x)

于是乎迭代可以写成：xk+1=xk−(∇2f(xk))−1∇f(xk)x_{k+1}=x_k-(
abla^2 f(x_k))^{-1}
abla f(x_k)xk+1​=xk​−(∇2f(xk​))−1∇f(xk​)，该方法称之为牛顿法（Newton）。

我们准备好来优化：min⁡wεnmin_{old w}varepsilon_nminw​εn​

假定：θ^n=∑k=1nwkXk=wTXhat heta_n=sum_{k=1}^nw_k X_k=old w^TXθ^n​=∑k=1n​wk​Xk​=wTX
其中：w=[w1,w2,…,wn]T，X=[X1,X2,…,Xn]Told w=[w_1,w_2,…,w_n]^T，X=[X_1,X_2,…,X_n]^Tw=[w1​,w2​,…,wn​]T，X=[X1​,X2​,…,Xn​]T

但有个问题，随着 nnn 变大，wold ww的长度会变长。需要改一下，取最近的M（与n无关）M（与n无关）M（与n无关）个数据：min⁡wεn=E[∣θn−∑k=0M−1wkXn−k∣2]min_{old w} varepsilon_n=mathbb E[| heta_n-sum_{k=0}^{M-1}w_kX_{n-k}|^2]minw​εn​=E[∣θn​−∑k=0M−1​wk​Xn−k​∣2]

在角标符号上，为防止 优化迭代的角标（当前迭代轮数） 跟 数据角标（时间） 之间意义混淆，需做一点区分：min⁡wε(n)=E[∣θ(n)−∑k=0M−1wkX(n−k)∣2]min_{old w} varepsilon(n)=mathbb E[| heta(n)-sum_{k=0}^{M-1}w_kX(n-k)|^2]minw​ε(n)=E[∣θ(n)−∑k=0M−1​wk​X(n−k)∣2]

而实际迭代的是：w(0)→w(1)→⋅⋅⋅→w(n)→n→∞wopt(n)old w^{(0)} o old w^{(1)} o··· o old w^{(n)}stackrel{n o infty}{ o}old w_{opt}(n)w(0)→w(1)→⋅⋅⋅→w(n)→n→∞wopt​(n)

不难发现，woptold w_{opt}wopt​随着nnn在变，这使得我们在某一时刻迭代太多次数的意义并不大。一种比较好的做法是每一个时间步都只做一次迭代，即：把整体迭代的角标和时间等同起来，如图所示：

这个做法就叫自适应（Adaptive）。

遵循梯度法，我们接着来做这个∇wε(n)
abla_{old w}varepsilon(n)∇w​ε(n)：

得到梯度法结果：w(n+1)=w(n)−∇wε(n)=w(n)−2E[θ(n)X(n)]+2E[X(n)XT(n)]w(n)old w(n+1)=old w(n)-
abla_{old w}varepsilon(n)=w(n)-2mathbb E[ heta(n)X(n)]+2mathbb E[X(n)X^T(n)]old w(n)w(n+1)=w(n)−∇w​ε(n)=w(n)−2E[θ(n)X(n)]+2E[X(n)XT(n)]w(n)

值得注意的是，我们还需要一个 步长（Step Size）μmuμ 来控制优化的实际性能，以避免不同初始条件下可能发生的震荡情景，如图所示：

再来写我们的迭代：w(n+1)=w(n)−μ∇wε(n)=w(n)−μE[(θ(n)−wT(n)X(n))⋅X(n)]=w(n)−μE[e(n)⋅X(n)]old w(n+1)=old w(n)-mu
abla_{old w}varepsilon(n)=old w(n)-mumathbb E[( heta(n)-old w^T(n) X(n))·X(n)]=old w(n)-mumathbb E[e(n)·X(n)]w(n+1)=w(n)−μ∇w​ε(n)=w(n)−μE[(θ(n)−wT(n)X(n))⋅X(n)]=w(n)−μE[e(n)⋅X(n)]
其中：e(n)=θ(n)−wT(n)X(n)e(n)= heta(n)-old w^T(n) X(n)e(n)=θ(n)−wT(n)X(n)

我们发现 E[e(n)⋅X(n)]mathbb E[e(n)·X(n)]E[e(n)⋅X(n)] 的期望并不好做，本质原因还是源于场景的非平稳特性。我们只能够把期望去掉不做，即由瞬时值替代原先的期望计算（Widrow）。

w(n+1)=w(n)−μ[e(n)⋅X(n)]old w(n+1)=old w(n)-mu[e(n)·X(n)]w(n+1)=w(n)−μ[e(n)⋅X(n)]，该方法即为LMS算法。

LMS算法包含了三层思想：

用优化而非解线性方程迭代的index和时间等同用瞬时值取代期望值

LMS算法开创了 随机梯度下降（SGD，Stochastic Gradient Descent） 的时代。

下面我们需要从数学上证明一下为什么步长的合适选取，能够帮助收敛？

说到收敛，直观地，我们希望误差能够趋于零：ε(n)→n→∞0varepsilon(n)stackrel{n o infty}{ o}0ε(n)→n→∞0

但事实上，我们首先需要判断这样一件事：w(n)→n→∞w0(n)（假设w0(n)是n时刻最优的，即w0(n)=E[X(n)XT(n)]−1E[X(n)θ(n)]）old w(n)stackrel{n o infty}{ o}old w_0(n)quad （假设old w_0(n)是n时刻最优的，即old w_0(n)=mathbb E[X(n)X^T(n)]^{-1}mathbb E[X(n) heta(n)]）w(n)→n→∞w0​(n)（假设w0​(n)是n时刻最优的，即w0​(n)=E[X(n)XT(n)]−1E[X(n)θ(n)]）

第一步，两边都减去w0(n)old w_0(n)w0​(n)：w(n+1)−w0(n)=w(n)−w0(n)−μ[θ(n)−wTX(n)]⋅X(n)old w(n+1)-old w_0(n)=old w(n)-old w_0(n)-mu[ heta(n)-old w^T X(n)]·X(n)w(n+1)−w0​(n)=w(n)−w0​(n)−μ[θ(n)−wTX(n)]⋅X(n)

令r(n+1)=w(n+1)−w0(n)r(n+1)=old w(n+1)-old w_0(n)r(n+1)=w(n+1)−w0​(n)，则能化简成：

差分方程解的渐进特性：an+1=Aan+B，{an→0}  ⟺  {∣A∣<1}a_{n+1}=Aa_n+B，quad{a_n o0} iff {|A|<1}an+1​=Aan​+B，{an​→0}⟺{∣A∣<1}

我们发现[I+μX(n)XT(n)]和r(n)left[I+mu X(n)X^T(n)
ight]和r(n)[I+μX(n)XT(n)]和r(n) 原则上讲是不独立的，所以期望无法拆开。

我们只好做出一个独立假设：[I+μX(n)XT(n)]和r(n)独立left[I+mu X(n)X^T(n)
ight]和r(n)独立[I+μX(n)XT(n)]和r(n)独立

因此

令：E(n)=E[r(n+1)]E(n)=mathbb E[r(n+1)]E(n)=E[r(n+1)]，则方程变为：

则方程写为：

我们希望随着 n→∞，Ek→0n o infty，E_k o 0n→∞，Ek​→0，那么

我们再回过头来看 ε(n)→0varepsilon(n) o0ε(n)→0这件事情：

Tr(AB)=Tr(BA) ext{Tr}(AB)= ext{Tr}(BA)Tr(AB)=Tr(BA)

令：R(n)=E[r(n)rT(n)]，K(n)=E[X(n)XT(n)]R(n)=mathbb E[r(n)r^T(n)]，K(n)=mathbb E[X(n)X^T(n)]R(n)=E[r(n)rT(n)]，K(n)=E[X(n)XT(n)]则

自适应这种思想在信号处理当中占据很特殊的地位，因为我们所处理的信号多数情况下都假设具有平稳性，而平稳性在现实当中是很难存在的。因此在不平稳的情况下，我们希望在性能和时效二者之间取一个平衡。

LMS算法导出过程：

假设一组数据：X1,X2,…,XnX_1,X_2,…,X_nX1​,X2​,…,Xn​（角标 nnn 代表着时刻）目标：dnd_ndn​要估计的参数：θn heta_nθn​评估指标为均方距离：L(Xn,dn,θn)=E(dn−θTXn)2⇒∇θL=−2E[(dn−θTXn)Xn]L(X_n,d_n, heta_n)=mathbb E(d_n- heta^TX_n)^2Rightarrow
abla_{ heta}L=-2mathbb E[(d_n- heta^T X_n)X_n]L(Xn​,dn​,θn​)=E(dn​−θTXn​)2⇒∇θ​L=−2E[(dn​−θTXn​)Xn​]构造优化梯度算法：θ^n+1=θ^n−μ∇θL(Xn,dn,θ^n)hat heta_{n+1}=hat heta_n-mu
abla_{ heta}L(X_n,d_n,hat heta_n)θ^n+1​=θ^n​−μ∇θ​L(Xn​,dn​,θ^n​)（此时角标 nnn 既作为优化迭代指标又作为时刻，故称之为自适应）把求导的期望去掉不做，就得到了所谓的 LMS：θ^n+1=θ^n−μenXn,en=dn−θ^nTXnhat heta_{n+1}=hat heta_n-mu e_nX_n,quad e_n=d_n-hat heta_n^TX_nθ^n+1​=θ^n​−μen​Xn​,en​=dn​−θ^nT​Xn​