遞迴最小平方濾波器

遞迴最小平方濾波器（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）