在RLS自適應濾波器的實現過程中,難免不涉及矩陣的求逆運算。而求逆操作雙是非常耗時的,一個很自然的想法就是盡可能的避免直接對矩陣進行求逆運算。那么,在RLS自適應濾波器的實現中,有沒有一種方法能避免直接求逆運算呢?答案當然是有的:使用矩陣求逆引理來避免對矩陣進行直接求逆。
這里先對矩陣求逆引理做下介紹,也叫做Woodbury矩陣恆等式(或者稱做Sherman–Morrison formula,這里統一稱矩陣求逆引理)在線性代數中:
\[{\left( {A + UCV} \right)^{ - 1}} = {A^{ - 1}} - {A^{ - 1}}U{\left( {{C^{ - 1}} + V{A^{ - 1}}U} \right)^{ - 1}}V{A^{ - 1}}\]
其中,A、C必須是可逆方陣,U、V可以是方陣也可以不是方陣。先不忙着記公式,這個不是最重要的,首先要弄明白的是,矩陣求逆引理的寫成這個樣子,它要解決的是什么樣的問題,其思想是什么?
矩陣求逆引理要解決的問題是:已知一個矩陣A的逆矩陣,當A矩陣產生了變化時,能不能根據已知的A的逆矩陣,求產生變化后的矩陣的逆。
現在我們知道當然是可以的,這里的變化,指的就是恆等式中的矩陣UCV。
知道了矩陣求逆引理解決的是什么樣的問題,主要思想后,感覺還是有點不得足,最好能證明一下,好讓自己有個“底”,同時滿足下好奇心。知其然同時又知其所以然,有什么不好!這里只對其中一種證明方法做下詳細的介紹,以方便理解,考慮以下方程組:
\[\left\{ {\begin{array}{*{20}{c}}
{AX + UY = I} \\
{VX - {C^{ - 1}}Y = 0} \\
\end{array}} \right.\]
及其矩陣形式
\[\left[ {\begin{array}{*{20}{c}}
A & U \\
V & { - {C^{ - 1}}} \\
\end{array}} \right]\left[ {\begin{array}{*{20}{c}}
X \\
Y \\
\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}
I \\
0 \\
\end{array}} \right]\]
根據第2個方程可知:$Y = CVX$,代入第1個方程得到
\[\left( {A + UCV} \right)X = I \Leftrightarrow X = {\left( {A + UCV} \right)^{ - 1}}\]
根據第一個方程可得:$X = {A^{ - 1}}(I - UY)$,將X的表達式代回第二個方程,得到
\[{C^{ - 1}}Y = V{A^{ - 1}}(I - UY) = V{A^{ - 1}} - V{A^{ - 1}}UY\]
將$V{A^{ - 1}}$與${C^{ - 1}}Y$交換等式兩邊位置
\[V{A^{ - 1}} = {C^{ - 1}}Y + V{A^{ - 1}}UY = \left( {{C^{ - 1}} + V{A^{ - 1}}U} \right)Y\]
可得到Y的表達式$Y = {\left( {{C^{ - 1}} + V{A^{ - 1}}U} \right)^{ - 1}}V{A^{ - 1}}$,這里將Y再次第1個方程$AX + UY = I$,得到
\[\begin{array}{l}
AX + U{\left( {{C^{ - 1}} + V{A^{ - 1}}U} \right)^{ - 1}}V{A^{ - 1}} = I\\
AX = I - U{\left( {{C^{ - 1}} + V{A^{ - 1}}U} \right)^{ - 1}}V{A^{ - 1}}\\
X = {A^{ - 1}}\left[ {I - U{{\left( {{C^{ - 1}} + V{A^{ - 1}}U} \right)}^{ - 1}}V{A^{ - 1}}} \right]\\
X = {A^{ - 1}} - {A^{ - 1}}U{\left( {{C^{ - 1}} + V{A^{ - 1}}U} \right)^{ - 1}}V{A^{ - 1}}\\
{\left( {A + UCV} \right)^{ - 1}} = X = {A^{ - 1}} - {A^{ - 1}}U{\left( {{C^{ - 1}} + V{A^{ - 1}}U} \right)^{ - 1}}V{A^{ - 1}}
\end{array}\]
關於其它證明部分,建議可以參考wiki中對Woodbury矩陣恆等式的介紹,里面有好幾種不同的證明方法(直接證明、代數證明、矩陣分塊推導、LDU分解推導),寫的都非常好,建議喜歡追根究底的朋友認真看下。
另外,矩陣求逆引理更通用的是二項式反轉定理(Binomial inverse theorem),很多的變種都可以視為二項式反轉定理的特殊情況,這些變種icoolmedia還沒有整理完畢,先不一一說明了。
下面就用矩陣求逆引理來解決RLS濾波器中的求逆問題,這里假定讀者有最小二乘法的基礎,首先給出更新濾波器輸入向量相關矩陣的公式和矩陣求逆引理的表達式(下面的推導過程來自於Simon Haykin《自適應濾波器原理》第10章的學習)。
\[{\bf{\Phi }}(n) = \lambda {\bf{\Phi }}(n - 1) + {\bf{u}}(n){{\bf{u}}^{\bf{H}}}(n)\]
\[{{\bf{A}}^{ - 1}} = {\bf{B}} - {\bf{BC}}{\left( {{\bf{D}} + {{\bf{C}}^{\bf{H}}}{\bf{BC}}} \right)^{ - 1}}{{\bf{C}}^{\bf{H}}}{\bf{B}}\]
可以看到,上面的矩陣求逆引理與本文開始里介紹的不同,這也是一個求逆引理的變種。不多做介紹,先記下來就好,后面另在博客中詳述吧。如果設:${\bf{A}} = \Phi (n)$,${{\bf{B}}^{ - 1}} = \lambda \Phi (n - 1)$,${\bf{C}} = {\bf{u}}(n)$,${\bf{D}} = 1$。對波器輸入向量相關矩陣迭代更新公式應用矩陣求逆引理公式,可得
\[{{\bf{\Phi }}^{ - 1}}(n) = {\lambda ^{ - 1}}{{\bf{\Phi }}^{ - 1}}(n - 1) - \frac{{{\lambda ^{ - 2}}{{\bf{\Phi }}^{ - 1}}(n - 1){\bf{u}}(n){{\bf{u}}^{\bf{H}}}(n){{\bf{\Phi }}^{ - 1}}(n - 1)}}{{1 + {\lambda ^{ - 1}}{{\bf{u}}^{\bf{H}}}(n){{\bf{\Phi }}^{ - 1}}(n - 1){\bf{u}}(n)}}\]
如果設
\[{\bf{P}}(n) = {{\bf{\Phi }}^{ - 1}}(n)\]
\[{\bf{k}}(n) = \frac{{{\lambda ^{ - 1}}{\bf{P}}(n){\bf{u}}(n)}}{{1 + {\lambda ^{ - 1}}{{\bf{u}}^{\bf{H}}}(n){{\bf{\Phi }}^{ - 1}}(n - 1){\bf{u}}(n)}}\]
這里對增益k(n)做如下變換,得到其簡化形式:
\[\begin{array}{l}
{\bf{k}}(n) = {\lambda ^{ - 1}}{\bf{P}}(n - 1){\bf{u}}(n) - {\lambda ^{ - 1}}{\bf{k}}(n){{\bf{u}}^{\bf{H}}}(n){\bf{P}}(n - 1){\bf{u}}(n) \\
= \left[ {{\lambda ^{ - 1}}{\bf{P}}(n - 1) - {\lambda ^{ - 1}}{\bf{k}}(n){{\bf{u}}^{\bf{H}}}(n){\bf{P}}(n - 1)} \right]{\bf{u}}(n) \\
= {\bf{P}}(n){\bf{u}}(n) \\
\end{array}\]
或者也可以表示為:${\bf{k}}(n) = {\Phi ^{ - 1}}(n - 1){\bf{u}}(n)$。由線性最小二乘估計可知,濾波器的系數估計表示為:${{\bf{\hat w}}(n) = {\Phi ^{ - 1}}(n){\bf{z}}(n)}$,把增益k(n)代入,得到
\[\begin{array}{l}
{\bf{\hat w}}(n) = {\Phi ^{ - 1}}(n){\bf{z}}(n) = {\bf{P}}(n){\bf{z}}(n) = \lambda {\bf{P}}(n){\bf{z}}(n - 1) + {\bf{P}}(n){\bf{u}}(n){d^*}(n) \\
= \lambda \left[ {{\lambda ^{ - 1}}{\bf{P}}(n - 1) - {\lambda ^{ - 1}}{\bf{k}}(n){{\bf{u}}^{\bf{H}}}(n){\bf{P}}(n - 1)} \right]{\bf{z}}(n - 1) + {\bf{P}}(n){\bf{u}}(n){d^*}(n) \\
= {\bf{P}}(n - 1){\bf{z}}(n - 1) - {\bf{k}}(n){{\bf{u}}^{\bf{H}}}(n){\bf{P}}(n - 1){\bf{z}}(n - 1) + {\bf{P}}(n){\bf{u}}(n){d^*}(n) \\
= {\Phi ^{ - 1}}(n - 1){\bf{z}}(n - 1) - {\bf{k}}(n){{\bf{u}}^{\bf{H}}}(n){\Phi ^{ - 1}}(n - 1){\bf{z}}(n - 1) + {\bf{P}}(n){\bf{u}}(n){d^*}(n) \\
= {\bf{\hat w}}(n - 1) - {\bf{k}}(n){{\bf{u}}^{\bf{H}}}(n){\bf{\hat w}}(n - 1) + {\bf{P}}(n){\bf{u}}(n){d^*}(n) \\
= {\bf{\hat w}}(n - 1) + \left[ {{\bf{k}}(n){d^*}(n) - {\bf{k}}(n){{\bf{u}}^{\bf{H}}}(n){\bf{\hat w}}(n - 1)} \right] \\
= {\bf{\hat w}}(n - 1) + {\bf{k}}(n)\left[ {{d^*}(n) - {{\bf{u}}^{\bf{H}}}(n){\bf{\hat w}}(n - 1)} \right] \\
= {\bf{\hat w}}(n - 1) + {\bf{k}}(n)\xi (n) \\
\end{array}\]
這就是遞歸最小二乘(RLS)濾波器的系數更新公式,可以看到,這個過程,並沒有對輸入向量相關矩陣進行直接求逆,避免了直接求逆運算帶來的復雜性。還有最后一個問題沒有解決。P的初始化問題,這里簡單給出來:${\bf{\hat w}}(n) = 0$,${\bf{P}}(0) = {\delta ^{ - 1}}{\bf{I}}$ 。這里$\delta $是正則化參數,是為了在求解反問題或者最優化問題時解決不適定性,增加抗繞動能力時引入的,這里不做詳細介紹。下面總結下遞歸最小二乘法的實現步驟。
\[\xi (n) = {d^*}(n) - {{\bf{u}}^{\bf{H}}}(n){\bf{\hat w}}(n - 1)\]
\[{\bf{k}}(n) = \frac{{{\lambda ^{ - 1}}{\bf{P}}(n - 1){\bf{u}}(n)}}{{1 + {\lambda ^{ - 1}}{{\bf{u}}^{\bf{H}}}(n){\bf{P}}(n - 1){\bf{u}}(n)}}\]
\[{\bf{\hat w}}(n) = {\bf{\hat w}}(n - 1) + {\bf{k}}(n)\xi (n)\]
\[{\bf{P}}(n) = {\lambda ^{ - 1}}{\bf{P}}(n - 1) - {\lambda ^{ - 1}}{\bf{k}}(n){{\bf{u}}^{\bf{H}}}(n){\bf{P}}(n - 1)\]
終於推導完畢,以上步驟的實現代碼及工程可以在博文下面的QQ交流群的群文件中找到(TestRLS_using the matrix inversion lemma.rar)。本人水平有限,如以上證明和推導過程及代碼如有錯誤之誤,還請大家及時給予批評指正。