递归最小二乘滤波器

递归最小二乘滤波器（RLS）属于自适应滤波器，会针对和输入信号有关的加权最小二乘（英语：Weighted least squares）损失函数，递归查找可以使其最小化的系数。此方法和想要减少均方误差的最小均方滤波器（LMS）不同。在推导递归最小二乘滤波器时，会假设输入信号是确定性的，而最小均方滤波及其他算法会假设信号随机。和其他的方法比较起来，递归最小二乘滤波器的收敛速度特别快。但其代价是非常高的运算复杂度。

演进

递归最小二乘滤波器是由卡尔·弗里德里希·高斯发现，但被遗忘无人使用，直到Plackett在1950年发现高斯1821年的著作之后，才再获使用。一般来说，递归最小二乘滤波器可以求得任何可以用自适应滤波器求解的问题。例如，假设信号${\displaystyle d(n)}$从有回声的有噪信道传输，接收到的是

${\displaystyle x(n)=\sum _{k=0}^{q}b_{n}(k)d(n-k)+v(n)}$

其中${\displaystyle v(n)}$代表加性高斯白噪声。RLS滤波器的目的是用${\displaystyle p+1}$阶有限冲激响应（FIR）滤波器${\displaystyle \mathbf {w} }$，还原想要的信号${\displaystyle d(n)}$：

${\displaystyle d(n)\approx \sum _{k=0}^{p}w(k)x(n-k)=\mathbf {w} ^{\mathit {T}}\mathbf {x} _{n}}$

其中${\displaystyle \mathbf {x} _{n}=[x(n)\quad x(n-1)\quad \ldots \quad x(n-p)]^{T}}$是包括${\displaystyle x(n)}$最近${\displaystyle p+1}$次取样的列向量。接收到想要信号的估测值为

${\displaystyle {\hat {d}}(n)=\sum _{k=0}^{p}w_{n}(k)x(n-k)=\mathbf {w} _{n}^{\mathit {T}}\mathbf {x} _{n}}$

其目的是估测滤波器${\displaystyle \mathbf {w} }$的参数，在每个时间${\displaystyle n}$，会将目前的估测值用${\displaystyle \mathbf {w} _{n}}$表示，自适应的最小方差估测值为${\displaystyle \mathbf {w} _{n+1}}$。${\displaystyle \mathbf {w} _{n}}$也是如下的列向量，其转置${\displaystyle \mathbf {w} _{n}^{\mathit {T}}}$则是行向量。矩阵乘法${\displaystyle \mathbf {w} _{n}^{\mathit {T}}\mathbf {x} _{n}}$（也是${\displaystyle \mathbf {w} _{n}}$和${\displaystyle \mathbf {x} _{n}}$的点积）是${\displaystyle {\hat {d}}(n)}$为标量。若${\displaystyle {\hat {d}}(n)-d(n)}$在最小二乘法的概念下，其值是小的，其估测就算是好的。

随着时间演进，会希望避免完全重做最小方差算法来找新估测值${\displaystyle \mathbf {w} _{n+1}}$，会希望可以将新估测值${\displaystyle \mathbf {w} _{n+1}}$用${\displaystyle \mathbf {w} _{n}}$配合其他变量来表示。

递归最小二乘滤波器的好处是不用进行反矩阵运算，因此可以节省计算时间。另一个好处是在其结果之后，提供了类似卡尔曼滤波的直觉信息。

概念

递归最小二乘滤波器背后的概念是在收到新资料时，适当选择滤波器系数${\displaystyle \mathbf {w} _{n}}$、更新滤波器，以及让损失函数 ${\displaystyle C}$最小化。误差信号${\displaystyle e(n)}$和期望信号${\displaystyle d(n)}$之间的关系，可以用以下的负反馈框图来说明：

误差会透过估测值${\displaystyle {\hat {d}}(n)}$，受到滤波器系数影响：

${\displaystyle e(n)=d(n)-{\hat {d}}(n)}$

加权最小方差函数${\displaystyle C}$—希望可以最小化的费用函数—是${\displaystyle e(n)}$的函数，因此也会受到滤波器系数的影响：

${\displaystyle C(\mathbf {w} _{n})=\sum _{i=0}^{n}\lambda ^{n-i}e^{2}(i)}$

其中${\displaystyle 0<\lambda \leq 1}$，是“遗忘因子”（forgetting factor），以指数的方式让较旧的资料有较小的加权。

费用函数最小化的方式，是将系数向量${\displaystyle \mathbf {w} _{n}}$中的所有项次微分，并让微分为零

${\displaystyle {\frac {\partial C(\mathbf {w} _{n})}{\partial w_{n}(k)}}=\sum _{i=0}^{n}2\lambda ^{n-i}e(i)\cdot {\frac {\partial e(i)}{\partial w_{n}(k)}}=-\sum _{i=0}^{n}2\lambda ^{n-i}e(i)\,x(i-k)=0\qquad k=0,1,\ldots ,p}$

接着，用以下的误差信号来代替${\displaystyle e(n)}$

${\displaystyle \sum _{i=0}^{n}\lambda ^{n-i}\left[d(i)-\sum _{\ell =0}^{p}w_{n}(\ell )x(i-\ell )\right]x(i-k)=0\qquad k=0,1,\ldots ,p}$

将等式重组如下

${\displaystyle \sum _{\ell =0}^{p}w_{n}(\ell )\left[\sum _{i=0}^{n}\lambda ^{n-i}\,x(i-\ell )x(i-k)\right]=\sum _{i=0}^{n}\lambda ^{n-i}d(i)x(i-k)\qquad k=0,1,\ldots ,p}$

可以表示为以下的矩阵

${\displaystyle \mathbf {R} _{x}(n)\,\mathbf {w} _{n}=\mathbf {r} _{dx}(n)}$

其中${\displaystyle \mathbf {R} _{x}(n)}$是${\displaystyle x(n)}$的加权样本协方差（英语：Sample mean and sample covariance）矩阵，而${\displaystyle \mathbf {r} _{dx}(n)}$是${\displaystyle d(n)}$和${\displaystyle x(n)}$互协方差的等效估计。依照上式可以找到最小化费用函数的系数

${\displaystyle \mathbf {w} _{n}=\mathbf {R} _{x}^{-1}(n)\,\mathbf {r} _{dx}(n)}$

如何选择λ

${\displaystyle \lambda }$越小，旧数据对协方差矩阵的影响越小，让滤波器对最近的数据较敏感，这也会让滤波器的co-efficients比较容易振荡。${\displaystyle \lambda =1}$称为growing window递归最小二乘算法。在实务上，会让${\displaystyle \lambda }$0.98介于1之间^[1]。利用第二型的最大可能区间估测，可以用资料中估测到最佳的${\displaystyle \lambda }$^[2]。

递归算法

以上叙述的结论是可以决定费用函数的参数，使费用函数最小化的方程式。以下则说明如何找到此形式的递归解

${\displaystyle \mathbf {w} _{n}=\mathbf {w} _{n-1}+\Delta \mathbf {w} _{n-1}}$

其中${\displaystyle \Delta \mathbf {w} _{n-1}}$是时间${\displaystyle {n-1}}$的修正因子。首先将互协方差${\displaystyle \mathbf {r} _{dx}(n)}$用${\displaystyle \mathbf {r} _{dx}(n-1)}$来表示

${\displaystyle \mathbf {r} _{dx}(n)}$	${\displaystyle =\sum _{i=0}^{n}\lambda ^{n-i}d(i)\mathbf {x} (i)}$
	${\displaystyle =\sum _{i=0}^{n-1}\lambda ^{n-i}d(i)\mathbf {x} (i)+\lambda ^{0}d(n)\mathbf {x} (n)}$
	${\displaystyle =\lambda \mathbf {r} _{dx}(n-1)+d(n)\mathbf {x} (n)}$

其中${\displaystyle \mathbf {x} (i)}$是${\displaystyle {p+1}}$维的资料向量

${\displaystyle \mathbf {x} (i)=[x(i),x(i-1),\dots ,x(i-p)]^{T}}$

接下来以相似的方式，用${\displaystyle \mathbf {R} _{x}(n-1)}$表示${\displaystyle \mathbf {R} _{x}(n)}$

${\displaystyle \mathbf {R} _{x}(n)}$	${\displaystyle =\sum _{i=0}^{n}\lambda ^{n-i}\mathbf {x} (i)\mathbf {x} ^{T}(i)}$
	${\displaystyle =\lambda \mathbf {R} _{x}(n-1)+\mathbf {x} (n)\mathbf {x} ^{T}(n)}$

为了要找到其系数向量，接下来要关注的是决定性自协方差矩阵的反矩阵。这问题可以使用伍德伯里矩阵恒等式（英语：Woodbury matrix identity）。若

${\displaystyle A}$	${\displaystyle =\lambda \mathbf {R} _{x}(n-1)}$ 是 ${\displaystyle (p+1)\times (p+1)}$矩阵
${\displaystyle U}$	${\displaystyle =\mathbf {x} (n)}$ 是 ${\displaystyle (p+1)\times 1}$（列向量）
${\displaystyle V}$	${\displaystyle =\mathbf {x} ^{T}(n)}$ 是 ${\displaystyle 1\times (p+1)}$（行向量）
${\displaystyle C}$	${\displaystyle =\mathbf {I} _{1}}$ 是 ${\displaystyle 1\times 1}$ 单位矩阵

依照伍德伯里矩阵恒等式，可得到下式

${\displaystyle \mathbf {R} _{x}^{-1}(n)}$	${\displaystyle =}$	${\displaystyle \left[\lambda \mathbf {R} _{x}(n-1)+\mathbf {x} (n)\mathbf {x} ^{T}(n)\right]^{-1}}$
	${\displaystyle =}$	${\displaystyle {\dfrac {1}{\lambda }}\left\lbrace \mathbf {R} _{x}^{-1}(n-1)-{\dfrac {\mathbf {R} _{x}^{-1}(n-1)\mathbf {x} (n)\mathbf {x} ^{T}(n)\mathbf {R} _{x}^{-1}(n-1)}{\lambda +\mathbf {x} ^{T}(n)\mathbf {R} _{x}^{-1}(n-1)\mathbf {x} (n)}}\right\rbrace }$
		${\displaystyle }$

为了和标准的文献一致，定义

${\displaystyle \mathbf {P} (n)}$	${\displaystyle =\mathbf {R} _{x}^{-1}(n)}$
	${\displaystyle =\lambda ^{-1}\mathbf {P} (n-1)-\mathbf {g} (n)\mathbf {x} ^{T}(n)\lambda ^{-1}\mathbf {P} (n-1)}$

其中的增益向量${\displaystyle g(n)}$为

${\displaystyle \mathbf {g} (n)}$	${\displaystyle =\lambda ^{-1}\mathbf {P} (n-1)\mathbf {x} (n)\left\{1+\mathbf {x} ^{T}(n)\lambda ^{-1}\mathbf {P} (n-1)\mathbf {x} (n)\right\}^{-1}}$
	${\displaystyle =\mathbf {P} (n-1)\mathbf {x} (n)\left\{\lambda +\mathbf {x} ^{T}(n)\mathbf {P} (n-1)\mathbf {x} (n)\right\}^{-1}}$

在往下推导之前，需要将${\displaystyle \mathbf {g} (n)}$改为以下的形式

${\displaystyle \mathbf {g} (n)\left\{1+\mathbf {x} ^{T}(n)\lambda ^{-1}\mathbf {P} (n-1)\mathbf {x} (n)\right\}}$	${\displaystyle =\lambda ^{-1}\mathbf {P} (n-1)\mathbf {x} (n)}$
${\displaystyle \mathbf {g} (n)+\mathbf {g} (n)\mathbf {x} ^{T}(n)\lambda ^{-1}\mathbf {P} (n-1)\mathbf {x} (n)}$	${\displaystyle =\lambda ^{-1}\mathbf {P} (n-1)\mathbf {x} (n)}$

等式两侧减去左边的第二项，得到

${\displaystyle \mathbf {g} (n)}$	${\displaystyle =\lambda ^{-1}\mathbf {P} (n-1)\mathbf {x} (n)-\mathbf {g} (n)\mathbf {x} ^{T}(n)\lambda ^{-1}\mathbf {P} (n-1)\mathbf {x} (n)}$
	${\displaystyle =\lambda ^{-1}\left[\mathbf {P} (n-1)-\mathbf {g} (n)\mathbf {x} ^{T}(n)\mathbf {P} (n-1)\right]\mathbf {x} (n)}$

配合${\displaystyle \mathbf {P} (n)}$的递归式定义，希望的形式如下

${\displaystyle \mathbf {g} (n)=\mathbf {P} (n)\mathbf {x} (n)}$

此时就可以完成递归，如以上讨论

${\displaystyle \mathbf {w} _{n}}$	${\displaystyle =\mathbf {P} (n)\,\mathbf {r} _{dx}(n)}$
	${\displaystyle =\lambda \mathbf {P} (n)\,\mathbf {r} _{dx}(n-1)+d(n)\mathbf {P} (n)\,\mathbf {x} (n)}$

第二步是从${\displaystyle \mathbf {r} _{dx}(n)}$的递归式定义开始，接着使用${\displaystyle \mathbf {P} (n)}$的递归式定义，配合调整后的${\displaystyle \mathbf {g} (n)}$，可以得到

${\displaystyle \mathbf {w} _{n}}$	${\displaystyle =\lambda \left[\lambda ^{-1}\mathbf {P} (n-1)-\mathbf {g} (n)\mathbf {x} ^{T}(n)\lambda ^{-1}\mathbf {P} (n-1)\right]\mathbf {r} _{dx}(n-1)+d(n)\mathbf {g} (n)}$
	${\displaystyle =\mathbf {P} (n-1)\mathbf {r} _{dx}(n-1)-\mathbf {g} (n)\mathbf {x} ^{T}(n)\mathbf {P} (n-1)\mathbf {r} _{dx}(n-1)+d(n)\mathbf {g} (n)}$
	${\displaystyle =\mathbf {P} (n-1)\mathbf {r} _{dx}(n-1)+\mathbf {g} (n)\left[d(n)-\mathbf {x} ^{T}(n)\mathbf {P} (n-1)\mathbf {r} _{dx}(n-1)\right]}$

配合${\displaystyle \mathbf {w} _{n-1}=\mathbf {P} (n-1)\mathbf {r} _{dx}(n-1)}$，可以得到以下的更新方程式

${\displaystyle \mathbf {w} _{n}}$	${\displaystyle =\mathbf {w} _{n-1}+\mathbf {g} (n)\left[d(n)-\mathbf {x} ^{T}(n)\mathbf {w} _{n-1}\right]}$
	${\displaystyle =\mathbf {w} _{n-1}+\mathbf {g} (n)\alpha (n)}$

其中${\displaystyle \alpha (n)=d(n)-\mathbf {x} ^{T}(n)\mathbf {w} _{n-1}}$ 是先验误差。将此和后验误差（在滤波器更新后计算的误差）比较

${\displaystyle e(n)=d(n)-\mathbf {x} ^{T}(n)\mathbf {w} _{n}}$

这就找到了修正因子

${\displaystyle \Delta \mathbf {w} _{n-1}=\mathbf {g} (n)\alpha (n)}$

这个结论指出了修正系数直接和误差和增益向量成正比，增益向量会透过加权因子${\displaystyle \lambda }$影响想要的灵敏度，这个结论很符合直觉。

RLS算法摘要

p阶RLS滤波器的算法可以摘要如下

参数：	${\displaystyle p=}$阶数
	${\displaystyle \lambda =}$遗忘因子
	${\displaystyle \delta =\mathbf {P} (0)}$的初始值
开始：	${\displaystyle \mathbf {w} (0)=0}$,
	${\displaystyle x(k)=0,k=-p,\dots ,-1}$,
	${\displaystyle d(k)=0,k=-p,\dots ,-1}$
	${\displaystyle \mathbf {P} (0)=\delta I}$其中${\displaystyle I}$是${\displaystyle p+1}$阶的单位矩阵
计算：	针对${\displaystyle n=1,2,\dots }$
	${\displaystyle \mathbf {x} (n)=\left[{\begin{matrix}x(n)\\x(n-1)\\\vdots \\x(n-p)\end{matrix}}\right]}$
	${\displaystyle \alpha (n)=d(n)-\mathbf {x} ^{T}(n)\mathbf {w} (n-1)}$
	${\displaystyle \mathbf {g} (n)=\mathbf {P} (n-1)\mathbf {x} (n)\left\{\lambda +\mathbf {x} ^{T}(n)\mathbf {P} (n-1)\mathbf {x} (n)\right\}^{-1}}$
	${\displaystyle \mathbf {P} (n)=\lambda ^{-1}\mathbf {P} (n-1)-\mathbf {g} (n)\mathbf {x} ^{T}(n)\lambda ^{-1}\mathbf {P} (n-1)}$
	${\displaystyle \mathbf {w} (n)=\mathbf {w} (n-1)+\,\alpha (n)\mathbf {g} (n)}$.

${\displaystyle P}$的递归依照代数Riccati方程，也类似卡尔曼滤波的结果^[3]。

相关条目

• 自适应滤波器
• 核自适应滤波器（英语：Kernel adaptive filter）
• 最小均方滤波器
• 迫零均衡器（英语：Zero-forcing equalizer）