查看“︁遞迴最小平方濾波器”︁的源代码

{{NoteTA
|G1 = Math
|G2 = IT
}}
'''遞迴最小平方濾波器'''（RLS）屬於[[自适应滤波器]]，會針對和輸入信號有關的{{link-en|加權最小平方|Weighted least squares}}[[损失函数]]，遞迴尋找可以使其最小化的係數。此方法和想要減少[[均方误差]]的[[最小均方滤波器]]（LMS）不同。在推導遞迴最小平方濾波器時，會假設輸入信號是[[确定性模型|确定性]]的，而最小均方滤波及其他演算法會假設信號{{link-en|隨機|stochastic}}。和其他的方法比較起來，遞迴最小平方濾波器的收斂速度特別快。但其代價是非常高的運算複雜度。
==演進==
遞迴最小平方濾波器是由[[卡爾·弗里德里希·高斯]]發現，但被遺忘無人使用，直到Plackett在1950年發現高斯1821年的著作之後，才再獲使用。一般來說，遞迴最小平方濾波器可以求得任何可以用[[自适应滤波器]]求解的問題。例如，假設信號<math>d(n)</math>從有回聲的[[有噪信道编码定理|有噪信道]]傳輸，接收到的是

:<math>x(n)=\sum_{k=0}^q b_n(k) d(n-k)+v(n)</math>

其中<math>v(n)</math>代表[[加性高斯白噪声]]。RLS濾波器的目的是用<math>p+1</math>階[[有限冲激响应]]（FIR）濾波器<math>\mathbf{w}</math>，還原想要的訊號<math>d(n)</math>：

:<math>d(n) \approx \sum_{k=0}^{p} w(k)x(n-k)=\mathbf{w}^\mathit{T} \mathbf{x}_n</math>

其中<math>\mathbf{x}_n = [x(n)\quad x(n-1)\quad\ldots\quad x(n-p)]^T</math>是包括<math>x(n)</math>最近<math>p+1</math>次取樣的[[行向量與列向量|-{zh-tw:行; zh-cn:列;}-向量]]<!--column-->。接收到想要訊號的估測值為

:<math>\hat{d}(n) = \sum_{k=0}^{p} w_n(k)x(n-k)=\mathbf{w}_n^\mathit{T} \mathbf{x}_n</math>

其目的是估測濾波器<math>\mathbf{w}</math>的參數，在每個時間<math>n</math>，會將目前的估測值用<math>\mathbf{w}_n</math>表示，自適應的最小方差估測值為<math>\mathbf{w}_{n+1}</math>。<math>\mathbf{w}_n</math>也是如下的列向量，其[[转置矩阵|转置]]<math>\mathbf{w}_n^\mathit{T}</math>則是-{zh-tw:列; zh-cn:行;}-<!--row-->向量。[[矩陣乘法]]<math>\mathbf{w}_n^\mathit{T} \mathbf{x}_n</math>（也是<math>\mathbf{w}_n</math>和<math>\mathbf{x}_n</math>的[[點積]]）是<math>\hat{d}(n)</math>為純量。若<math>\hat{d}(n)-d(n)</math>在[[最小二乘法]]的概念下，其值是小的，其估測就算是好的。

隨著時間演進，會希望避免完全重做最小方差演算法來找新估測值<math>\mathbf{w}_{n+1}</math>，會希望可以將新估測值<math>\mathbf{w}_{n+1}</math>用<math>\mathbf{w}_n</math>配合其他變數來表示。

遞迴最小平方濾波器的好處是不用進行反矩陣運算，因此可以節省計算時間。另一個好處是在其結果之後，提供了類似[[卡尔曼滤波]]的直覺資訊。

==概念==
遞迴最小平方濾波器背後的概念是在收到新資料時，適當選擇濾波器係數<math>\mathbf{w}_n</math>、更新濾波器，以及讓[[损失函数]] <math>C</math>最小化。誤差信號<math>e(n)</math>和期望信號<math>d(n)</math>之間的關係，可以用以下的[[负反馈]]方塊圖來說明：

[[File:AdaptiveFilter C.png|500px]]

誤差會透過估測值<math>\hat{d}(n)</math>，受到濾波器係數影響：

:<math>e(n)=d(n)-\hat{d}(n)</math>

加權最小方差函數<math>C</math>—希望可以最小化的費用函數—是<math>e(n)</math>的函數，因此也會受到濾波器係數的影響：
:<math>C(\mathbf{w}_n)=\sum_{i=0}^n\lambda^{n-i}e^2(i)</math>
其中<math>0<\lambda\le 1</math>，是「遺忘因子」（forgetting factor），以指數的方式讓較舊的資料有較小的加權。

費用函數最小化的方式，是將係數向量<math>\mathbf{w}_{n}</math>中的所有項次微分，並讓微分為零
:<math>\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</math>
接著，用以下的誤差信號來代替<math>e(n)</math>
:<math>\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</math>
將等式重組如下
:<math>\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</math>
可以表示為以下的矩陣
:<math>\mathbf{R}_{x}(n)\,\mathbf{w}_{n}=\mathbf{r}_{dx}(n)</math>
其中<math>\mathbf{R}_{x}(n)</math>是<math>x(n)</math>的加權{{link-en|樣本平均數和樣本協方差|Sample mean and sample covariance|樣本協方差}}矩陣，而<math>\mathbf{r}_{dx}(n)</math>是<math>d(n)</math>和<math>x(n)</math>[[互协方差]]的等效估計。依照上式可以找到最小化費用函數的係數
:<math>\mathbf{w}_{n}=\mathbf{R}_{x}^{-1}(n)\,\mathbf{r}_{dx}(n)</math>

===如何選擇λ===
<math>\lambda</math>越小，舊數據對協方差矩陣的影響越小，讓濾波器對最近的數據較敏感，這也會讓濾波器的co-efficients比較容易振盪。<math>\lambda=1</math>稱為growing window遞迴最小平方演算法。在實務上，會讓<math>\lambda</math>0.98介於1之間<ref>Emannual C. Ifeacor, Barrie W. Jervis. Digital signal processing: a practical approach, second edition. Indianapolis: Pearson Education Limited, 2002, p. 718</ref>。利用第二型的最大可能區間估測，可以用資料中估測到最佳的<math>\lambda</math><ref>Steven Van Vaerenbergh, Ignacio Santamaría, Miguel Lázaro-Gredilla [http://gtas.unican.es/files/pub/stkrls_mlsp2012.pdf "Estimation of the forgetting factor in kernel recursive least squares"], 2012 IEEE International Workshop on Machine Learning for Signal Processing, 2012, accessed June 23, 2016.</ref>。
<!-- someone might like to include a diagram to show said fluctuations -->
==遞迴演算法==
以上敘述的結論是可以決定費用函數的參數，使費用函數最小化的方程式。以下則說明如何找到此形式的遞迴解
:<math>\mathbf{w}_{n}=\mathbf{w}_{n-1}+\Delta\mathbf{w}_{n-1}</math>
其中<math>\Delta\mathbf{w}_{n-1}</math>是時間<math>{n-1}</math>的修正因子。首先將互協方差<math>\mathbf{r}_{dx}(n)</math>用<math>\mathbf{r}_{dx}(n-1)</math>來表示
:{|
|-
|<math>\mathbf{r}_{dx}(n)</math>
|<math>=\sum_{i=0}^{n}\lambda^{n-i}d(i)\mathbf{x}(i)</math>
|-
|
|<math>=\sum_{i=0}^{n-1}\lambda^{n-i}d(i)\mathbf{x}(i)+\lambda^{0}d(n)\mathbf{x}(n)</math>
|-
|
|<math>=\lambda\mathbf{r}_{dx}(n-1)+d(n)\mathbf{x}(n)</math>
|}
其中<math>\mathbf{x}(i)</math>是<math>{p+1}</math>維的資料向量
:<math>\mathbf{x}(i)=[x(i), x(i-1), \dots , x(i-p) ]^{T}</math>
接下來以相似的方式，用<math>\mathbf{R}_{x}(n-1)</math>表示<math>\mathbf{R}_{x}(n)</math>
:{|
|-
|<math>\mathbf{R}_{x}(n)</math>
|<math>=\sum_{i=0}^{n}\lambda^{n-i}\mathbf{x}(i)\mathbf{x}^{T}(i)</math>
|-
|
|<math>=\lambda\mathbf{R}_{x}(n-1)+\mathbf{x}(n)\mathbf{x}^{T}(n)</math>
|}
為了要找到其係數向量，接下來要關注的是決定性自協方差矩陣的反矩陣。這問題可以使用{{link-en|伍德伯里矩陣恆等式|Woodbury matrix identity}}。若
:{|
|-
|<math>A</math>
|<math>=\lambda\mathbf{R}_{x}(n-1)</math> 是 <math>(p+1) \times (p+1)</math>矩陣
|-
|<math>U</math>
|<math>=\mathbf{x}(n)</math> 是 <math>(p+1) \times 1</math>（-{zh-tw:行; zh-cn:列;}-<!--column-->向量）
|-
|<math>V</math>
|<math>=\mathbf{x}^{T}(n)</math> 是 <math>1 \times (p+1)</math>（-{zh-tw:列; zh-cn:行;}-<!--row-->向量）
|-
|<math>C</math>
|<math>=\mathbf{I}_1</math> 是 <math>1 \times 1</math> [[單位矩陣]]
|}
依照伍德伯里矩陣恆等式，可得到下式
:{|
|-
|<math>\mathbf{R}_{x}^{-1}(n)</math>
|<math>=</math>
|<math>\left[\lambda\mathbf{R}_{x}(n-1)+\mathbf{x}(n)\mathbf{x}^{T}(n)\right]^{-1}</math>
|-
|
|<math>=</math>
|<math>\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</math>
|-
|
|
|<math></math>

|}
為了和標準的文獻一致，定義
:{|
|-
|<math>\mathbf{P}(n)</math>
|<math>=\mathbf{R}_{x}^{-1}(n)</math>
|-
|
|<math>=\lambda^{-1}\mathbf{P}(n-1)-\mathbf{g}(n)\mathbf{x}^{T}(n)\lambda^{-1}\mathbf{P}(n-1)</math>
|}
其中的增益向量<math>g(n)</math>為
:{|
|-
|<math>\mathbf{g}(n)</math>
|<math>=\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}</math>
|-
|
|<math>=\mathbf{P}(n-1)\mathbf{x}(n)\left\{\lambda+\mathbf{x}^{T}(n)\mathbf{P}(n-1)\mathbf{x}(n)\right\}^{-1}</math>
|}
在往下推導之前，需要將<math>\mathbf{g}(n)</math>改為以下的形式
:{|
|-
|<math>\mathbf{g}(n)\left\{1+\mathbf{x}^{T}(n)\lambda^{-1}\mathbf{P}(n-1)\mathbf{x}(n)\right\}</math>
|<math>=\lambda^{-1}\mathbf{P}(n-1)\mathbf{x}(n)</math>
|-
|<math>\mathbf{g}(n)+\mathbf{g}(n)\mathbf{x}^{T}(n)\lambda^{-1}\mathbf{P}(n-1)\mathbf{x}(n)</math>
|<math>=\lambda^{-1}\mathbf{P}(n-1)\mathbf{x}(n)</math>
|}
等式兩側減去左邊的第二項，得到
:{|
|-
|<math>\mathbf{g}(n)</math>
|<math>=\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)</math>
|-
|
|<math>=\lambda^{-1}\left[\mathbf{P}(n-1)-\mathbf{g}(n)\mathbf{x}^{T}(n)\mathbf{P}(n-1)\right]\mathbf{x}(n)</math>
|}
配合<math>\mathbf{P}(n)</math>的遞迴式定義，希望的形式如下
:<math>\mathbf{g}(n)=\mathbf{P}(n)\mathbf{x}(n)</math>
此時就可以完成遞迴，如以上討論
:{|
|-
|<math>\mathbf{w}_{n}</math>
|<math>=\mathbf{P}(n)\,\mathbf{r}_{dx}(n)</math>
|-
|
|<math>=\lambda\mathbf{P}(n)\,\mathbf{r}_{dx}(n-1)+d(n)\mathbf{P}(n)\,\mathbf{x}(n)</math>
|}
第二步是從<math>\mathbf{r}_{dx}(n )</math>的遞迴式定義開始，接著使用<math>\mathbf{P}(n)</math>的遞迴式定義，配合調整後的<math>\mathbf{g}(n)</math>，可以得到
:{|
|-
|<math>\mathbf{w}_{n}</math>
|<math>=\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)</math>
|-
|
|<math>=\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)</math>
|-
|
|<math>=\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]</math>
|}
配合<math>\mathbf{w}_{n-1}=\mathbf{P}(n-1)\mathbf{r}_{dx}(n-1)</math>，可以得到以下的更新方程式
:{|
|-
|<math>\mathbf{w}_{n}</math>
|<math>=\mathbf{w}_{n-1}+\mathbf{g}(n)\left[d(n)-\mathbf{x}^{T}(n)\mathbf{w}_{n-1}\right]</math>
|-
|
|<math>=\mathbf{w}_{n-1}+\mathbf{g}(n)\alpha(n)</math>
|}
其中<math>\alpha(n)=d(n)-\mathbf{x}^{T}(n)\mathbf{w}_{n-1}</math>
是[[先驗與後驗|先驗]]誤差。將此和後驗誤差（在濾波器更新後計算的誤差）比較

:<math>e(n)=d(n)-\mathbf{x}^{T}(n)\mathbf{w}_n</math>

這就找到了修正因子
:<math>\Delta\mathbf{w}_{n-1}=\mathbf{g}(n)\alpha(n)</math>

這個結論指出了修正係數直接和誤差和增益向量成正比，增益向量會透過加權因子<math>\lambda</math>影響想要的靈敏度，這個結論很符合直覺。

==RLS演算法摘要==
''p''階RLS濾波器的演算法可以摘要如下
{|
|-
| 參數：|| <math>p=</math>階數
|-
|  || <math>\lambda=</math>遺忘因子
|-
|  || <math>\delta= \mathbf{P}(0)</math>的初始值
|-
|開始：|| <math>\mathbf{w}(0)=0</math>,
|-
| ||<math>x(k)=0, k=-p,\dots,-1</math>,
|-
| ||<math>d(k)=0, k=-p,\dots,-1</math>
|-
| ||<math>\mathbf{P}(0)=\delta I</math>其中<math>I</math>是<math>p+1</math>階的[[單位矩陣]]
|-
|計算：|| 針對<math>n=1,2,\dots </math>
|-
|||
<math> \mathbf{x}(n) = 
\left[
\begin{matrix}
x(n)\\
x(n-1)\\
\vdots\\
x(n-p)
\end{matrix}
\right]
</math>
|-
|||<math> \alpha(n) = d(n)-\mathbf{x}^T(n)\mathbf{w}(n-1)</math>
|-
|||<math>\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}</math>
|-
|||<math>\mathbf{P}(n)=\lambda^{-1}\mathbf{P}(n-1)-\mathbf{g}(n)\mathbf{x}^{T}(n)\lambda^{-1}\mathbf{P}(n-1)</math>
|-
|||<math> \mathbf{w}(n) = \mathbf{w}(n-1)+\,\alpha(n)\mathbf{g}(n)</math>.
|}

<math>P</math>的遞迴依照[[代數Riccati方程]]，也類似[[卡尔曼滤波]]的結果<ref>Welch, Greg and Bishop, Gary [http://www.cs.unc.edu/~welch/media/pdf/kalman_intro.pdf "An Introduction to the Kalman Filter"], Department of Computer Science, University of North Carolina at Chapel Hill, September 17, 1997, accessed July 19, 2011.</ref>。
==相關條目==
*[[自适应滤波器]]
*{{link-en|核自適應濾波器|Kernel adaptive filter}}
*[[最小均方滤波器]]
*{{link-en|迫零均衡器|Zero-forcing equalizer}}

==書目==
{{Reflist}}
* {{Cite book|author=Hayes, Monson H.|title=Statistical Digital Signal Processing and Modeling|url=https://archive.org/details/statisticaldigit0000haye|chapter=9.4: Recursive Least Squares|page=[https://archive.org/details/statisticaldigit0000haye/page/n560 541]|publisher=Wiley|year=1996|isbn=0-471-59431-8}}
* Simon Haykin, ''Adaptive Filter Theory'', Prentice Hall, 2002, {{ISBN|0-13-048434-2}}
* M.H.A Davis, R.B. Vinter, ''Stochastic Modelling and Control'', Springer, 1985, {{ISBN|0-412-16200-8}}
* Weifeng Liu, Jose Principe and Simon Haykin, ''Kernel Adaptive Filtering: A Comprehensive Introduction'', John Wiley, 2010, {{ISBN|0-470-44753-2}}
* R.L.Plackett, ''Some Theorems in Least Squares'', Biometrika, 1950, 37, 149–157, {{ISSN|0006-3444}}
* C.F.Gauss, ''Theoria combinationis observationum erroribus minimis obnoxiae'', 1821, Werke, 4. Gottinge

[[Category:数字信号处理]]
[[Category:滤波器理论]]
[[Category:統計信號處理]]