递归最小二乘(RLS)算法详解

本文首发在我的个人博客：宅到没朋友，欢迎大家来玩！

代码下载链接在最下方！

1.最小二乘算法简介

最小二乘算法基于确定性思想，该算法讨论怎样根据有限个观测数据来寻找滤波器的最优解，即求如下图这样具有M个抽头的横向滤波器的最优权向量。

最小二乘算法的目的就是通过最小化误差向量$\vec{e}$的模的平方和，即最小化$J=\Vert \vec{e}{\Vert }^2$，来求解得到最优权向量：

$$\vec{w}=[w_0{w}_1\cdots w_{M-1}]$$

由于：

$$\vec{e}=[e(M)e(M+1)\cdots e(N){]}^H=\vec{b}-\vec{\hat{b}}=\vec{b}-A\vec{w}$$

其中：

• 滤波器阶数：$M$
• 样本点数：$N$
• 误差向量：$\vec{e}=[e(M)e(M+1)\cdots e(N)]$
• 期望向量：$\vec{b}=[d(M)d(M+1)\cdots d(N)]$
• 数据矩阵：$A^H=[\vec{u}(M)\vec{u}(M+1)\cdots \vec{u}(N)]=\left[\begin{array}{cccc}u(M)& u(M+1)& \cdots & u(N)\\ u(M-1)& u(M)& \cdots & u(N-1)\\ ⋮& ⋮& \ddots & ⋮\\ u(1)& u(2)& \cdots & u(N-M+1)\end{array}\right]$

由此可以得到：

$$J=\Vert \vec{e}{\Vert }^2={\vec{e}}^H\vec{e}=(\vec{b}-A\vec{w}{)}^H(\vec{b}-A\vec{w})$$

化简得：
$$J={\vec{b}}^H\vec{b}-{\vec{b}}^{H}A\vec{w}-{\vec{w}}^H{A}^H\vec{b}+{\vec{w}}^H{A}^{H}A\vec{w}$$

要求$J$的最小值，首先要得到$J$关于$\vec{w}$的梯度：

$$▽J=-2A^H\vec{b}+2A^{H}A\vec{w}$$

再令$▽J=0$得：

$$A^{H}A\vec{\hat{w}}=A^H\vec{b}(1)$$

上式被称为确定性正则方程。

当$M<N-M+1$时，若方阵$A^{H}A$是非奇异（可逆）的，则可以得到确定性正则方程的解：

$$\vec{\hat{w}}=(A^{H}A{)}^{-1}{A}^H\vec{b}(2)$$

上式也被称为最小二乘$(LS)$解。

估计向量$\vec{\hat{b}}=A\vec{w}$被称为对期望响应向量$\vec{b}$的最小二乘估计，简称$LS$估计。

2.递归最小二乘$(RLS)$算法

在(2)式中，涉及矩阵$A^{H}A$的求逆运算，运算量较大，递归最小二乘$(RLS)$算法就是用迭代算法代替矩阵求逆达到降低运算量的目的。

将数据矩阵$A^H$做如下扩展：

$$A^H=\left[\begin{array}{cccccc}u(1)& u(2)& \cdots & u(M)& \cdots & u(N)\\ 0& u(1)& \cdots & u(M-1)& \cdots & u(N-1)\\ ⋮& ⋮& \ddots & ⋮& \ddots & ⋮\\ 0& 0& \cdots & u(1)& \cdots & u(N-M+1)\end{array}\right]$$

将期望向量$\vec{b}$做如下扩展：
$$\vec{b}=[d(1)d(2)\cdots d(M)\cdots d(N)]$$

定义输入数据数据的时间相关矩阵：
$$\Phi (N)=A^{H}A=\sum _{i=1}^N\vec{u}(i){\vec{u}}^H(i)$$

定义时间互相关向量为
$$\vec{z}(N)=A^H\vec{b}=\sum _{i=1}^N\vec{u}(i)d^{\ast }(i)$$

于是，如式(1)所示的确定性正则方程变成下式：
$$\Phi (N)\vec{\hat{w}}=\vec{z}(N)(3)$$

(3)式就是N时刻的确定性正则方程。

将N时刻变为任意时刻n，并引入遗忘因子$\lambda (0<\lambda \le 1)$，则$\Phi$、$z$和(3)式变成如下所示：
$$\Phi (n)=\sum _{i=1}^n{\lambda }^{n-i}\vec{u}(i){\vec{u}}^H(i)+\delta {\lambda }^{n}I$$

$$z(n)=\sum _{i=1}^n{\lambda }^{n-i}\vec{u}(i)d^{\ast }(i)$$

$$\Phi (n)\vec{\hat{w}}=\vec{z}(n)$$

注：$\Phi (n)$中引入$\delta {\lambda }^{n}I$是为了使$\Phi (n)$可逆。

令$P(n)={\Phi }^{-1}(n)$，接下来需要的事情就是求出$\vec{\hat{w}}$的迭代公式，推导过程较为繁琐，在此不再展开，直接给出最后的迭代公式。
$$\vec{\hat{w}}(n)=\vec{\hat{w}}(n-1)+\vec{k}(n){\xi }^{\ast }(n)$$
其中：
$$\vec{k}(n)=\frac{{\lambda }^{-1}P(n-1)\vec{u}(n)}{1+{\lambda }^{-1}{\vec{u}}^H(n)P(n-1)\vec{u}(n)}$$
$$\xi (n)=d(n)-\vec{\hat{w}}(n-1)\vec{u}(n)$$
$$P(n)={\lambda }^{-1}P(n-1)-{\lambda }^{-1}\vec{k}(n){\vec{u}}^H(n)P(n-1)$$
可以看到$\vec{\hat{w}}(n)$迭代式中的变量$P(n)、\vec{k}(n))、\xi (n)$都是可以迭代求解的。

2.1 $(RLS)$算法步骤

1. 初始化：
$$P(0)=\delta I，\delta 是小正数$$

$$\vec{\hat{w}}(0)=\vec{0}$$

$$\lambda 取接近于1$$

2. $n=1,2,\cdots ,N$时，做如下迭代运算：
$$\vec{k}(n)=\frac{{\lambda }^{-1}P(n-1)\vec{u}(n)}{1+{\lambda }^{-1}{\vec{u}}^H(n)P(n-1)\vec{u}(n)}$$

$$\xi (n)=d(n)-\vec{\hat{w}}(n-1)\vec{u}(n)$$

$$\vec{\hat{w}}(n)=\vec{\hat{w}}(n-1)+\vec{k}(n){\xi }^{\ast }(n)$$

$$P(n)={\lambda }^{-1}P(n-1)-{\lambda }^{-1}\vec{k}(n){\vec{u}}^H(n)P(n-1)$$

3. 令$n=n+1$，重复步骤2。

2.2 算例

考虑一阶$AR$模型$u(n)=-0.99u(n-1)+v(n)$的线性预测。假设白噪声$v(n)$的方差为${\sigma }_v^2=0.995$，使用抽头数为$M=2$的$FIR$滤波器，用$RLS$算法实现$u(n)$的线性预测，选择遗忘因子$\lambda =0.98$。

2.3 $Matlab$仿真结果

可以看到$RLS$算法收敛的速度很快，收敛特性相当好。

2.4 $Matlab$源码下载

点此下载源码！

欢迎搜索并关注微信公众号“工科南”。