很多人在學習系統(tǒng)辨識、自適應濾波或者在線參數估計時都會卡在遞推最小二乘法Recursive Least Squares, RLS的公式推導上。教材里通常幾行帶過但實際自己要推一遍或者寫代碼實現時才發(fā)現從“批量求逆”到“遞推更新”之間其實隔著不少細節(jié)。我最早接觸RLS是在做在線辨識的時候當時對著講義看了好幾遍總覺得增益矩陣像個魔法——為什么這個式子加個新息就能不斷修正參數為什么協(xié)方差矩陣要按那個方式更新后來自己動手推了一遍又把初值、遺忘因子和各種工程變形都玩了一遍才真正覺得這公式“通了”。這篇就把我的推導過程和實踐心得完整拆一拆希望能幫你一次看透RLS。1. 從批量最小二乘到逐點遞歸先弄清楚RLS到底在解什么問題1.1 經典最小二乘閉式解到底長什么樣絕大多數人第一次接觸最小二乘都是從線性回歸開始的。假設系統(tǒng)模型為$$y \varphi^T \theta$$其中 $\varphi$ 是回歸向量輸入特征$\theta$ 是待估計的參數向量。如果我們拿到的是一批數據 ${(\varphi_i, y_i)}_{i1}^N$那么可以寫成矩陣形式$$Y \Phi \theta e$$$Y$ 是輸出向量$\Phi$ 的每一行是一個回歸向量 $\varphi_i^T$$e$ 是殘差。批量最小二乘的優(yōu)化目標是讓殘差平方和最小$$J(\theta) (Y - \Phi\theta)^T (Y - \Phi\theta)$$對 $\theta$ 求導并令其等于零可以得到著名的閉式解$$\hat{\theta} (\Phi^T \Phi)^{-1} \Phi^T Y$$這個式子本身非常簡潔但它把“過去所有時刻的數據”當作一個整體來看待。每次來一個樣本理論上你都要重新構造 $\Phi$ 和 $Y$然后重新算一次 $(\Phi^T\Phi)^{-1}$。當數據量不大時沒問題但一旦數據源源不斷涌過來這種批量解法就變得不太現實了。1.2 數據源源不斷時批處理有兩個尷尬之處第一個尷尬是計算量。$\Phi^T\Phi$ 的維度是參數個數 $n \times n$求逆的復雜度大約是 $O(n^3)$。假設你每一秒得到一個樣本每秒都要重新算一次求逆硬件功耗和實時性都會吃不消。就算用聰明的矩陣求逆算法反復做全量計算仍然很浪費。第二個尷尬更隱秘批量最小二乘對歷史數據一視同仁。如果系統(tǒng)參數本身是緩變的比如飛機飛行時的氣動參數隨高度變化電機繞組電阻隨溫度變化那么舊數據對當前時刻的估計其實已經沒有太大參考價值了。批量最小二乘只會把舊數據和新數據混在一起最終得到一個“歷史平均”的參數根本追不上系統(tǒng)變化。1.3 RLS的核心思路一句話講完RLS的核心思路其實一句話就能說清楚把上一時刻的參數估計當作基礎當新樣本到來時用這個樣本帶來的“新息”innovation對參數進行修正同時通過一個遞推式更新協(xié)方差矩陣避免顯式地重算矩陣求逆。換句話說RLS只記住兩個量當前參數估計值 $\hat{\theta}(t-1)$以及一個能代表歷史信息累積的矩陣 $P(t-1)$。每來一個新數據只做幾次矩陣乘法和一次標量除法就能得到新的 $\hat{\theta}(t)$ 和 $P(t)$。這既解決了計算量問題又可以通過遺忘因子靈活控制歷史數據的權重。2. 矩陣求逆引理RLS推導中最關鍵的一塊跳板2.1 引理本身和它的證明思路RLS公式能“化簡”成可遞推的形式核心依賴于一個線性代數工具——矩陣求逆引理Matrix Inversion Lemma也叫 Sherman-Morrison-Woodbury 公式。它的標準形式是$$(A B C D)^{-1} A^{-1} - A^{-1} B (C^{-1} D A^{-1} B)^{-1} D A^{-1}$$這里要求 $A$ 和 $C$ 都可逆。初看這個式子可能覺得很抽象但我建議你把它理解成一個“升級版”的分配率如果中間那個 $BCD$ 的乘積項很小你當然可以直接把 $A$ 的逆提取出來展開即使不展開成無窮級數也能用這個等式把大矩陣的求逆轉化成小矩陣的求逆。證明這個引理并不難。只要驗證右邊乘以 $(A BCD)$ 等于單位矩陣 $I$或者把等式左右兩邊同時乘開利用 $A^{-1}A I$ 和 $C^{-1}C I$就能逐步化簡得到恒等式。更直觀的記憶方法是“加了一項就減一項括號里是倒數之和”。這個公式的價值在于它會出現在RLS的協(xié)方差更新中幫我們把一個 $n \times n$ 的矩陣求逆轉化成括號內一個標量或小矩陣的求逆。2.2 怎么把引理套進RLS的協(xié)方差更新RLS會維護一個矩陣 $P(t)$它實際是信息矩陣 $R(t) \sum_{i1}^t \lambda^{t-i} \varphi_i \varphi_i^T$ 的逆。當新樣本到達后新的信息矩陣可以寫成遞推形式$$R(t) \lambda R(t-1) \varphi_t \varphi_t^T$$我們的目標是直接遞推 $P(t) R^{-1}(t)$而不是重新求逆。這時候令$$A \lambda R(t-1), \quad B \varphi_t, \quad C 1, \quad D \varphi_t^T$$帶入矩陣求逆引理得到$$R(t)^{-1} \frac{1}{\lambda} \left[ R(t-1)^{-1} - \frac{R(t-1)^{-1} \varphi_t \varphi_t^T R(t-1)^{-1}}{\lambda \varphi_t^T R(t-1)^{-1} \varphi_t} \right]$$你會發(fā)現原來需要對整個 $R(t)$ 求逆的問題現在只需要計算一個標量分母 $\lambda \varphi_t^T P(t-1) \varphi_t$。這就是RLS計算效率高的數學根源。3. 一步一步推出RLS的三個遞推方程3.1 目標函數和加權最小二乘的構建要推RLS先得有帶遺忘因子 $\lambda$ 的目標函數。為什么加遺忘因子因為希望越靠近當前時刻的數據權重越大。定義$$J_t(\theta) \sum_{i1}^{t} \lambda^{t-i} \left( y_i - \varphi_i^T \theta \right)^2$$其中 $0 \lambda \le 1$。當 $\lambda 1$ 時所有歷史數據權重一樣相當于普通最小二乘當 $\lambda 1$ 時老數據按指數衰減系統(tǒng)參數變化時估計值能更快跟上。令 $R(t) \sum_{i1}^t \lambda^{t-i} \varphi_i \varphi_i^T$$Q(t) \sum_{i1}^t \lambda^{t-i} \varphi_i y_i$顯然有遞推關系$$R(t) \lambda R(t-1) \varphi_t \varphi_t^T$$$$Q(t) \lambda Q(t-1) \varphi_t y_t$$讓 $J_t$ 對 $\theta$ 求導等于零得到最優(yōu)解$$\hat{\theta}(t) R(t)^{-1} Q(t)$$這就是RLS推導的出發(fā)點。接下來要做的事是把 $R^{-1}(t)$ 和 $\hat{\theta}(t)$ 都改寫成前一刻值的遞推。3.2 增益矩陣 $K(t)$ 的推導我們想找到形如 $\hat{\theta}(t) \hat{\theta}(t-1) K(t)\left(y_t - \varphi_t^T \hat{\theta}(t-1)\right)$ 的更新方程。這里的 $K(t)$ 稱為增益矩陣或增益向量。直接代入$$\hat{\theta}(t) P(t) Q(t)$$$$Q(t) \lambda Q(t-1) \varphi_t y_t$$又有 $\hat{\theta}(t-1) P(t-1) Q(t-1)$也就是 $Q(t-1) P^{-1}(t-1)\hat{\theta}(t-1)$。把這幾個式子串起來$$\hat{\theta}(t) P(t)\left[\lambda P^{-1}(t-1)\hat{\theta}(t-1) \varphi_t y_t\right]$$關鍵在于 $P(t) R^{-1}(t)$而 $P(t)$ 和 $P(t-1)$ 之間的遞推關系在第2節(jié)已經推導出來了。把 $P(t)$ 的表達式代入并利用 $\lambda P^{-1}(t-1) P(t)$ 這個組合經過整理后就能得到標準的增益矩陣定義$$K(t) \frac{P(t-1)\varphi_t}{\lambda \varphi_t^T P(t-1)\varphi_t}$$這一步是RLS公式中最容易看頭暈的地方。我之前卡了很久后來發(fā)現只需要盯著 $P(t)$ 的遞推式把所有 $P(t)$ 都替換成在第2節(jié)得到的結果大部分中間項會自動消掉。記住一個關鍵點分母是標量所以“求逆”根本不是真正的矩陣求逆而是一次普通除法。3.3 協(xié)方差更新公式的兩種等價寫法有了 $K(t)$$P(t)$ 的更新式可以寫成更緊湊的形式。第2節(jié)的結果其實等價于$$P(t) \frac{1}{\lambda}\left[P(t-1) - K(t)\varphi_t^T P(t-1)\right]$$也可以寫成$$P(t) \frac{1}{\lambda}\left(I - K(t)\varphi_t^T\right)P(t-1)$$這兩種寫法本質相同只是前者更便于觀察“減掉一項”的含義。$K(t)\varphi_t^T$ 是一個秩一矩陣意味著一維新樣本對高維協(xié)方差矩陣的修正是“秩一更新”。這個結構的幾何直覺是只有平行于 $\varphi_t$ 的方向估計的不確定性才會被明顯壓縮其他方向的變化相對較小。3.4 參數更新公式的另一種理解新息加權參數更新的標準形式是$$\hat{\theta}(t) \hat{\theta}(t-1) K(t)\left(y_t - \varphi_t^T \hat{\theta}(t-1)\right)$$括號里的 $y_t - \varphi_t^T \hat{\theta}(t-1)$ 就是新息表示“用現有模型預測的輸出”和“真實輸出”的誤差。增益 $K(t)$ 則告訴我們應該用多大比例來修正參數。當 $P(t-1)$ 較大時說明歷史信息積累不足$\varphi_t$ 方向的增益也會比較大當前新息對參數修正的幅度就大反之當參數已經收斂得很好$P(t-1)$ 很小增益自然變小新信息對參數的影響也變弱。這也是RLS比LMS最小均方算法收斂更快的原因——它把歷史數據的二階統(tǒng)計信息全部壓縮在 $P$ 矩陣中。4. 初值設定、遺忘因子和完整的RLS算法流程4.1 協(xié)方差矩陣初值怎么給到底該用大數還是小數代碼實現RLS時第一個問題就是 $P(0)$ 怎么設。大多數教材推薦 $P(0) \delta I$$\delta$ 取一個較大的數比如 $100$ 或 $10^3$。這個做法的理由是$P$ 矩陣本身可以理解為參數估計協(xié)方差矩陣的近似初始狀態(tài)下我們對參數幾乎一無所知所以把它的“不確定性”設置得很大讓前幾個樣本能快速修正參數。另一種場景是如果你已經有一個比較靠譜的先驗參數估計比如上一批次辨識好的參數那 $\delta$ 就可以取小一些比如 $0.1$ 甚至 $0.01$這樣初期就不會因為第一個樣本就把參數拉飛。還有一種常見做法是用一批小數據先做一次批量最小二乘得到初始值再把對應的協(xié)方差矩陣拷貝給 $P(0)$。這個做法最穩(wěn)但對于在線系統(tǒng)來說不一定有這么多先驗數據。4.2 遺忘因子 $\lambda$ 的選擇內存長度和跟蹤速度的權衡$\lambda$ 的取值直接影響算法對時變系統(tǒng)的跟蹤能力。$\lambda 1$ 時算法擁有無限記憶適合參數恒定的系統(tǒng)$\lambda 1$ 時相當于對不同時刻的數據賦了一個指數衰減權重。工程上常把“有效記憶長度”近似為$$N_{\text{eff}} \approx \frac{1}{1 - \lambda}$$比如 $\lambda 0.99$ 大約相當于只有最近100個樣本在起作用$\lambda 0.95$ 則只有大約20個樣本。$\lambda$ 越小跟蹤越快但受噪聲影響也越大$\lambda$ 太小時估計方差會明顯增加。所以實際使用時需要根據系統(tǒng)的變化速度和噪聲水平折中我的經驗是先從 $\lambda 0.98$ 左右試起然后觀察估計曲線的抖動幅度再逐步調整。寧可讓跟蹤慢一點也別讓估計值抖成心電圖。4.3 標準RLS算法流程偽代碼形式標準RLS算法每來一個新樣本只需要五步初始化 theta zeros(n, 1) P delta * eye(n) 循環(huán) for each t: 1. 計算增益向量 K P * phi / (lambda phi^T * P * phi) 2. 計算新息 e y - phi^T * theta 3. 更新參數 theta theta K * e 4. 更新協(xié)方差 P (P - K * phi^T * P) / lambda 5. 檢查 P 是否為對稱正定工程上常用 P (P P^T) / 2 強制對稱注意第4步的除法 $\lambda$ 是針對標量的可以直接除在矩陣上。如果想避免第5步的額外檢查也可以用平方根RLS之類的變形這一點我會在下一章詳細說。5. 數值穩(wěn)定性問題和工程改進措施5.1 為什么P矩陣會變成非正定運行RLS時間長了你可能會發(fā)現本來應該正定的協(xié)方差矩陣 $P(t)$ 會逐漸失去對稱正定性甚至出現負特征值。原因主要有三個第一計算機有限字長誤差。RLS遞推中反復做 $K(t)\varphi_t^T P(t-1)$ 這種矩陣乘法每一步都有舍入誤差長時間累積后可能導致對稱性破壞。第二遺忘因子導致“舊信息被指數衰減”當信號激勵不足時$P(t)$ 的某些方向會不斷被放大或縮小最終變得病態(tài)甚至非正定。第三輸入 $\varphi_t$ 持續(xù)相關或者某一段激勵太弱也會讓信息矩陣長時間不增長$P$ 就可能在數值上發(fā)散。一旦 $P$ 失去正定性參數估計可能會出現劇烈跳變增益向量 $K(t)$ 也可能出現符號異常。建議是在算法里添加監(jiān)控比如檢查 $P$ 的對稱性和特征值如果出現非正定及時重置或者強制對稱化。5.2 工程上常用的修正策略最簡單的修正是每次更新完以后執(zhí)行$$P_{\text{sym}} \frac{P P^T}{2}$$這個操作不會帶來太大成本卻能有效避免因非對稱導致的累積性誤差。其次可以引入正則化項在信息矩陣上疊加一個小量 $\epsilon I$這樣 $P^{-1}$ 始終有界。也可以設置一個“更新判定條件”當新息的絕對值特別大時暫時不更新協(xié)方差矩陣只更新參數避免異常樣本對 $P$ 造成污染。更穩(wěn)定的做法是使用平方根RLSSquare-Root RLS它把 $P$ 分解為 $P S S^T$遞推更新 $S$ 而不是更新 $P$ 本身。因為 $S$ 的特征值都是正的只要 $S$ 不溢出$P$ 就能一直保持正定。平方根RLS的代價是額外多一些三角運算但對長期運行的在線系統(tǒng)來說非常值得。5.3 平方根RLS的基本思想簡述平方根RLS的核心是用QR分解的思想更新信息矩陣的平方根。具體來說將 $P(t)^{1/2}$ 作為遞推量利用旋轉矩陣Givens旋轉或Householder變換使得更新后的平方根矩陣保持正定。標準RLS中的分母 $\lambda \varphi_t^T P(t-1) \varphi_t$ 可以看作一個標量在平方根版本中它將合并進一個增廣矩陣的QR更新過程。實現細節(jié)比較繁瑣但如果你的系統(tǒng)要求長時間不間斷運行強烈建議不要直接用基礎RLS而是上平方根版本。6. 一段可運行的Python實驗驗證RLS公式有沒有推錯6.1 實驗設計靜態(tài)參數和時變參數兩個場景紙上推了半天還是得用代碼驗證。這里我做一個最簡單的單輸入單輸出SISO系統(tǒng)辨識實驗模型是$$y_t a x_t b n_t$$其中 $a1.2$$b0.8$$n_t$ 是均值為0、方差為0.01的高斯噪聲?；貧w向量取 $\varphi_t [x_t, 1]^T$參數向量 $\theta [a, b]^T$。先測試靜態(tài)參數下RLS是否收斂到真實值再設計一個時變參數場景比如 $a$ 在第500個樣本時從1.2跳到0.5檢驗遺忘因子能不能幫算法跟上這個突變。6.2 核心代碼和結果解讀下面是一份非常精簡的Python實現直接用基礎RLS沒有花哨優(yōu)化import numpy as np import matplotlib.pyplot as plt np.random.seed(42) N 1000 x np.random.randn(N) a_true np.ones(N) * 1.2 b_true np.ones(N) * 0.8 a_true[500:] 0.5 y a_true * x b_true 0.1 * np.random.randn(N) theta np.zeros(2) P 1000 * np.eye(2) lam 0.98 theta_history [] for t in range(N): phi np.array([x[t], 1.0]) K P phi / (lam phi P phi) e y[t] - phi theta theta theta K * e P (P - np.outer(K, phi P)) / lam theta_history.append(theta.copy()) theta_history np.array(theta_history) plt.plot(theta_history[:,0], labela_hat) plt.plot(theta_history[:,1], labelb_hat) plt.axhline(1.2, colorgray, ls--) plt.axhline(0.5, colorgray, ls--) plt.legend() plt.show()運行結果你會看到前200步左右參數快速收斂到真實值附近當 $a$ 在第500步突變時RLS會在一小段延遲后重新逼近新的真值。這就是遺忘因子的作用。如果把 $\lambda$ 設成1你會看到突變后參數幾乎不動需要很久才能慢慢扭過去。6.3 踩坑記錄遺忘因子初值不當引發(fā)的發(fā)散實驗過程中最容易遇到的現象就是前幾步參數直接飛上天然后徹底發(fā)散。我早先試過把 $\lambda$ 設為0.9初始 $P(0)$ 設為 $1000I$結果第一個樣本就把參數修正得過猛因為 $K(1) P(0)\varphi_1 / (0.9 \varphi_1^T P(0)\varphi_1)$ 的分母雖然很大但分子也很大如果 $\varphi_1$ 的模很小導致增益變得巨大。解決辦法是適當減小 $P(0)$ 的初值或者約束參數更新范圍。這不是公式錯而是初值和遺忘因子搭配不當。換成 $\lambda 0.98$、$\delta 100$ 以后就穩(wěn)定多了。這類問題在公式推導時完全看不出來只有跑過代碼才深有體會。7. 復盤我對RLS公式推導和落地的一些經驗推完這一整套RLS公式我最大的感受是數學公式的每個變形都不是孤立的。矩陣求逆引理、遺忘因子、增益矩陣它們其實是同一個目標的不同側面。你不需要死記公式只需要記住兩條主線一條是信息矩陣 $R(t)$ 的遞推另一條是參數估計 $\hat{\theta}(t)$ 的遞推。所有的RLS變體都是圍繞這兩條主線做數值穩(wěn)定性和計算效率上的改進。實際調試中我習慣先在散點圖上畫出參數的收斂軌跡一旦發(fā)現軌跡異常就優(yōu)先檢查輸入信號的激勵程度——如果輸入一直恒定不變任何最小二乘類算法都會在某幾個方向上失去可辨識性這不是算法能救的。最后再分享一個小技巧如果系統(tǒng)是慢時變的但又怕 $\lambda$ 太大會跟不上可以用雙遺忘因子對輸入功率大的樣本用較大的 $\lambda$對輸入功率小的樣本用稍小的 $\lambda$這樣既能跟蹤突變又不會讓噪聲把參數抖得厲害。這個思路我也是從現場調參里慢慢摸出來的公式層面看不出來但工程上非常實用。 遞推最小二乘RLS公式推導與工程落地詳解
很多人在學習系統(tǒng)辨識、自適應濾波或者在線參數估計時都會卡在遞推最小二乘法Recursive Least Squares, RLS的公式推導上。教材里通常幾行帶過但實際自己要推一遍或者寫代碼實現時才發(fā)現從“批量求逆”到“遞推更新”之間其實隔著不少細節(jié)。我最早接觸RLS是在做在線辨識的時候當時對著講義看了好幾遍總覺得增益矩陣像個魔法——為什么這個式子加個新息就能不斷修正參數為什么協(xié)方差矩陣要按那個方式更新后來自己動手推了一遍又把初值、遺忘因子和各種工程變形都玩了一遍才真正覺得這公式“通了”。這篇就把我的推導過程和實踐心得完整拆一拆希望能幫你一次看透RLS。1. 從批量最小二乘到逐點遞歸先弄清楚RLS到底在解什么問題1.1 經典最小二乘閉式解到底長什么樣絕大多數人第一次接觸最小二乘都是從線性回歸開始的。假設系統(tǒng)模型為$$y \varphi^T \theta$$其中 $\varphi$ 是回歸向量輸入特征$\theta$ 是待估計的參數向量。如果我們拿到的是一批數據 ${(\varphi_i, y_i)}_{i1}^N$那么可以寫成矩陣形式$$Y \Phi \theta e$$$Y$ 是輸出向量$\Phi$ 的每一行是一個回歸向量 $\varphi_i^T$$e$ 是殘差。批量最小二乘的優(yōu)化目標是讓殘差平方和最小$$J(\theta) (Y - \Phi\theta)^T (Y - \Phi\theta)$$對 $\theta$ 求導并令其等于零可以得到著名的閉式解$$\hat{\theta} (\Phi^T \Phi)^{-1} \Phi^T Y$$這個式子本身非常簡潔但它把“過去所有時刻的數據”當作一個整體來看待。每次來一個樣本理論上你都要重新構造 $\Phi$ 和 $Y$然后重新算一次 $(\Phi^T\Phi)^{-1}$。當數據量不大時沒問題但一旦數據源源不斷涌過來這種批量解法就變得不太現實了。1.2 數據源源不斷時批處理有兩個尷尬之處第一個尷尬是計算量。$\Phi^T\Phi$ 的維度是參數個數 $n \times n$求逆的復雜度大約是 $O(n^3)$。假設你每一秒得到一個樣本每秒都要重新算一次求逆硬件功耗和實時性都會吃不消。就算用聰明的矩陣求逆算法反復做全量計算仍然很浪費。第二個尷尬更隱秘批量最小二乘對歷史數據一視同仁。如果系統(tǒng)參數本身是緩變的比如飛機飛行時的氣動參數隨高度變化電機繞組電阻隨溫度變化那么舊數據對當前時刻的估計其實已經沒有太大參考價值了。批量最小二乘只會把舊數據和新數據混在一起最終得到一個“歷史平均”的參數根本追不上系統(tǒng)變化。1.3 RLS的核心思路一句話講完RLS的核心思路其實一句話就能說清楚把上一時刻的參數估計當作基礎當新樣本到來時用這個樣本帶來的“新息”innovation對參數進行修正同時通過一個遞推式更新協(xié)方差矩陣避免顯式地重算矩陣求逆。換句話說RLS只記住兩個量當前參數估計值 $\hat{\theta}(t-1)$以及一個能代表歷史信息累積的矩陣 $P(t-1)$。每來一個新數據只做幾次矩陣乘法和一次標量除法就能得到新的 $\hat{\theta}(t)$ 和 $P(t)$。這既解決了計算量問題又可以通過遺忘因子靈活控制歷史數據的權重。2. 矩陣求逆引理RLS推導中最關鍵的一塊跳板2.1 引理本身和它的證明思路RLS公式能“化簡”成可遞推的形式核心依賴于一個線性代數工具——矩陣求逆引理Matrix Inversion Lemma也叫 Sherman-Morrison-Woodbury 公式。它的標準形式是$$(A B C D)^{-1} A^{-1} - A^{-1} B (C^{-1} D A^{-1} B)^{-1} D A^{-1}$$這里要求 $A$ 和 $C$ 都可逆。初看這個式子可能覺得很抽象但我建議你把它理解成一個“升級版”的分配率如果中間那個 $BCD$ 的乘積項很小你當然可以直接把 $A$ 的逆提取出來展開即使不展開成無窮級數也能用這個等式把大矩陣的求逆轉化成小矩陣的求逆。證明這個引理并不難。只要驗證右邊乘以 $(A BCD)$ 等于單位矩陣 $I$或者把等式左右兩邊同時乘開利用 $A^{-1}A I$ 和 $C^{-1}C I$就能逐步化簡得到恒等式。更直觀的記憶方法是“加了一項就減一項括號里是倒數之和”。這個公式的價值在于它會出現在RLS的協(xié)方差更新中幫我們把一個 $n \times n$ 的矩陣求逆轉化成括號內一個標量或小矩陣的求逆。2.2 怎么把引理套進RLS的協(xié)方差更新RLS會維護一個矩陣 $P(t)$它實際是信息矩陣 $R(t) \sum_{i1}^t \lambda^{t-i} \varphi_i \varphi_i^T$ 的逆。當新樣本到達后新的信息矩陣可以寫成遞推形式$$R(t) \lambda R(t-1) \varphi_t \varphi_t^T$$我們的目標是直接遞推 $P(t) R^{-1}(t)$而不是重新求逆。這時候令$$A \lambda R(t-1), \quad B \varphi_t, \quad C 1, \quad D \varphi_t^T$$帶入矩陣求逆引理得到$$R(t)^{-1} \frac{1}{\lambda} \left[ R(t-1)^{-1} - \frac{R(t-1)^{-1} \varphi_t \varphi_t^T R(t-1)^{-1}}{\lambda \varphi_t^T R(t-1)^{-1} \varphi_t} \right]$$你會發(fā)現原來需要對整個 $R(t)$ 求逆的問題現在只需要計算一個標量分母 $\lambda \varphi_t^T P(t-1) \varphi_t$。這就是RLS計算效率高的數學根源。3. 一步一步推出RLS的三個遞推方程3.1 目標函數和加權最小二乘的構建要推RLS先得有帶遺忘因子 $\lambda$ 的目標函數。為什么加遺忘因子因為希望越靠近當前時刻的數據權重越大。定義$$J_t(\theta) \sum_{i1}^{t} \lambda^{t-i} \left( y_i - \varphi_i^T \theta \right)^2$$其中 $0 \lambda \le 1$。當 $\lambda 1$ 時所有歷史數據權重一樣相當于普通最小二乘當 $\lambda 1$ 時老數據按指數衰減系統(tǒng)參數變化時估計值能更快跟上。令 $R(t) \sum_{i1}^t \lambda^{t-i} \varphi_i \varphi_i^T$$Q(t) \sum_{i1}^t \lambda^{t-i} \varphi_i y_i$顯然有遞推關系$$R(t) \lambda R(t-1) \varphi_t \varphi_t^T$$$$Q(t) \lambda Q(t-1) \varphi_t y_t$$讓 $J_t$ 對 $\theta$ 求導等于零得到最優(yōu)解$$\hat{\theta}(t) R(t)^{-1} Q(t)$$這就是RLS推導的出發(fā)點。接下來要做的事是把 $R^{-1}(t)$ 和 $\hat{\theta}(t)$ 都改寫成前一刻值的遞推。3.2 增益矩陣 $K(t)$ 的推導我們想找到形如 $\hat{\theta}(t) \hat{\theta}(t-1) K(t)\left(y_t - \varphi_t^T \hat{\theta}(t-1)\right)$ 的更新方程。這里的 $K(t)$ 稱為增益矩陣或增益向量。直接代入$$\hat{\theta}(t) P(t) Q(t)$$$$Q(t) \lambda Q(t-1) \varphi_t y_t$$又有 $\hat{\theta}(t-1) P(t-1) Q(t-1)$也就是 $Q(t-1) P^{-1}(t-1)\hat{\theta}(t-1)$。把這幾個式子串起來$$\hat{\theta}(t) P(t)\left[\lambda P^{-1}(t-1)\hat{\theta}(t-1) \varphi_t y_t\right]$$關鍵在于 $P(t) R^{-1}(t)$而 $P(t)$ 和 $P(t-1)$ 之間的遞推關系在第2節(jié)已經推導出來了。把 $P(t)$ 的表達式代入并利用 $\lambda P^{-1}(t-1) P(t)$ 這個組合經過整理后就能得到標準的增益矩陣定義$$K(t) \frac{P(t-1)\varphi_t}{\lambda \varphi_t^T P(t-1)\varphi_t}$$這一步是RLS公式中最容易看頭暈的地方。我之前卡了很久后來發(fā)現只需要盯著 $P(t)$ 的遞推式把所有 $P(t)$ 都替換成在第2節(jié)得到的結果大部分中間項會自動消掉。記住一個關鍵點分母是標量所以“求逆”根本不是真正的矩陣求逆而是一次普通除法。3.3 協(xié)方差更新公式的兩種等價寫法有了 $K(t)$$P(t)$ 的更新式可以寫成更緊湊的形式。第2節(jié)的結果其實等價于$$P(t) \frac{1}{\lambda}\left[P(t-1) - K(t)\varphi_t^T P(t-1)\right]$$也可以寫成$$P(t) \frac{1}{\lambda}\left(I - K(t)\varphi_t^T\right)P(t-1)$$這兩種寫法本質相同只是前者更便于觀察“減掉一項”的含義。$K(t)\varphi_t^T$ 是一個秩一矩陣意味著一維新樣本對高維協(xié)方差矩陣的修正是“秩一更新”。這個結構的幾何直覺是只有平行于 $\varphi_t$ 的方向估計的不確定性才會被明顯壓縮其他方向的變化相對較小。3.4 參數更新公式的另一種理解新息加權參數更新的標準形式是$$\hat{\theta}(t) \hat{\theta}(t-1) K(t)\left(y_t - \varphi_t^T \hat{\theta}(t-1)\right)$$括號里的 $y_t - \varphi_t^T \hat{\theta}(t-1)$ 就是新息表示“用現有模型預測的輸出”和“真實輸出”的誤差。增益 $K(t)$ 則告訴我們應該用多大比例來修正參數。當 $P(t-1)$ 較大時說明歷史信息積累不足$\varphi_t$ 方向的增益也會比較大當前新息對參數修正的幅度就大反之當參數已經收斂得很好$P(t-1)$ 很小增益自然變小新信息對參數的影響也變弱。這也是RLS比LMS最小均方算法收斂更快的原因——它把歷史數據的二階統(tǒng)計信息全部壓縮在 $P$ 矩陣中。4. 初值設定、遺忘因子和完整的RLS算法流程4.1 協(xié)方差矩陣初值怎么給到底該用大數還是小數代碼實現RLS時第一個問題就是 $P(0)$ 怎么設。大多數教材推薦 $P(0) \delta I$$\delta$ 取一個較大的數比如 $100$ 或 $10^3$。這個做法的理由是$P$ 矩陣本身可以理解為參數估計協(xié)方差矩陣的近似初始狀態(tài)下我們對參數幾乎一無所知所以把它的“不確定性”設置得很大讓前幾個樣本能快速修正參數。另一種場景是如果你已經有一個比較靠譜的先驗參數估計比如上一批次辨識好的參數那 $\delta$ 就可以取小一些比如 $0.1$ 甚至 $0.01$這樣初期就不會因為第一個樣本就把參數拉飛。還有一種常見做法是用一批小數據先做一次批量最小二乘得到初始值再把對應的協(xié)方差矩陣拷貝給 $P(0)$。這個做法最穩(wěn)但對于在線系統(tǒng)來說不一定有這么多先驗數據。4.2 遺忘因子 $\lambda$ 的選擇內存長度和跟蹤速度的權衡$\lambda$ 的取值直接影響算法對時變系統(tǒng)的跟蹤能力。$\lambda 1$ 時算法擁有無限記憶適合參數恒定的系統(tǒng)$\lambda 1$ 時相當于對不同時刻的數據賦了一個指數衰減權重。工程上常把“有效記憶長度”近似為$$N_{\text{eff}} \approx \frac{1}{1 - \lambda}$$比如 $\lambda 0.99$ 大約相當于只有最近100個樣本在起作用$\lambda 0.95$ 則只有大約20個樣本。$\lambda$ 越小跟蹤越快但受噪聲影響也越大$\lambda$ 太小時估計方差會明顯增加。所以實際使用時需要根據系統(tǒng)的變化速度和噪聲水平折中我的經驗是先從 $\lambda 0.98$ 左右試起然后觀察估計曲線的抖動幅度再逐步調整。寧可讓跟蹤慢一點也別讓估計值抖成心電圖。4.3 標準RLS算法流程偽代碼形式標準RLS算法每來一個新樣本只需要五步初始化 theta zeros(n, 1) P delta * eye(n) 循環(huán) for each t: 1. 計算增益向量 K P * phi / (lambda phi^T * P * phi) 2. 計算新息 e y - phi^T * theta 3. 更新參數 theta theta K * e 4. 更新協(xié)方差 P (P - K * phi^T * P) / lambda 5. 檢查 P 是否為對稱正定工程上常用 P (P P^T) / 2 強制對稱注意第4步的除法 $\lambda$ 是針對標量的可以直接除在矩陣上。如果想避免第5步的額外檢查也可以用平方根RLS之類的變形這一點我會在下一章詳細說。5. 數值穩(wěn)定性問題和工程改進措施5.1 為什么P矩陣會變成非正定運行RLS時間長了你可能會發(fā)現本來應該正定的協(xié)方差矩陣 $P(t)$ 會逐漸失去對稱正定性甚至出現負特征值。原因主要有三個第一計算機有限字長誤差。RLS遞推中反復做 $K(t)\varphi_t^T P(t-1)$ 這種矩陣乘法每一步都有舍入誤差長時間累積后可能導致對稱性破壞。第二遺忘因子導致“舊信息被指數衰減”當信號激勵不足時$P(t)$ 的某些方向會不斷被放大或縮小最終變得病態(tài)甚至非正定。第三輸入 $\varphi_t$ 持續(xù)相關或者某一段激勵太弱也會讓信息矩陣長時間不增長$P$ 就可能在數值上發(fā)散。一旦 $P$ 失去正定性參數估計可能會出現劇烈跳變增益向量 $K(t)$ 也可能出現符號異常。建議是在算法里添加監(jiān)控比如檢查 $P$ 的對稱性和特征值如果出現非正定及時重置或者強制對稱化。5.2 工程上常用的修正策略最簡單的修正是每次更新完以后執(zhí)行$$P_{\text{sym}} \frac{P P^T}{2}$$這個操作不會帶來太大成本卻能有效避免因非對稱導致的累積性誤差。其次可以引入正則化項在信息矩陣上疊加一個小量 $\epsilon I$這樣 $P^{-1}$ 始終有界。也可以設置一個“更新判定條件”當新息的絕對值特別大時暫時不更新協(xié)方差矩陣只更新參數避免異常樣本對 $P$ 造成污染。更穩(wěn)定的做法是使用平方根RLSSquare-Root RLS它把 $P$ 分解為 $P S S^T$遞推更新 $S$ 而不是更新 $P$ 本身。因為 $S$ 的特征值都是正的只要 $S$ 不溢出$P$ 就能一直保持正定。平方根RLS的代價是額外多一些三角運算但對長期運行的在線系統(tǒng)來說非常值得。5.3 平方根RLS的基本思想簡述平方根RLS的核心是用QR分解的思想更新信息矩陣的平方根。具體來說將 $P(t)^{1/2}$ 作為遞推量利用旋轉矩陣Givens旋轉或Householder變換使得更新后的平方根矩陣保持正定。標準RLS中的分母 $\lambda \varphi_t^T P(t-1) \varphi_t$ 可以看作一個標量在平方根版本中它將合并進一個增廣矩陣的QR更新過程。實現細節(jié)比較繁瑣但如果你的系統(tǒng)要求長時間不間斷運行強烈建議不要直接用基礎RLS而是上平方根版本。6. 一段可運行的Python實驗驗證RLS公式有沒有推錯6.1 實驗設計靜態(tài)參數和時變參數兩個場景紙上推了半天還是得用代碼驗證。這里我做一個最簡單的單輸入單輸出SISO系統(tǒng)辨識實驗模型是$$y_t a x_t b n_t$$其中 $a1.2$$b0.8$$n_t$ 是均值為0、方差為0.01的高斯噪聲?；貧w向量取 $\varphi_t [x_t, 1]^T$參數向量 $\theta [a, b]^T$。先測試靜態(tài)參數下RLS是否收斂到真實值再設計一個時變參數場景比如 $a$ 在第500個樣本時從1.2跳到0.5檢驗遺忘因子能不能幫算法跟上這個突變。6.2 核心代碼和結果解讀下面是一份非常精簡的Python實現直接用基礎RLS沒有花哨優(yōu)化import numpy as np import matplotlib.pyplot as plt np.random.seed(42) N 1000 x np.random.randn(N) a_true np.ones(N) * 1.2 b_true np.ones(N) * 0.8 a_true[500:] 0.5 y a_true * x b_true 0.1 * np.random.randn(N) theta np.zeros(2) P 1000 * np.eye(2) lam 0.98 theta_history [] for t in range(N): phi np.array([x[t], 1.0]) K P phi / (lam phi P phi) e y[t] - phi theta theta theta K * e P (P - np.outer(K, phi P)) / lam theta_history.append(theta.copy()) theta_history np.array(theta_history) plt.plot(theta_history[:,0], labela_hat) plt.plot(theta_history[:,1], labelb_hat) plt.axhline(1.2, colorgray, ls--) plt.axhline(0.5, colorgray, ls--) plt.legend() plt.show()運行結果你會看到前200步左右參數快速收斂到真實值附近當 $a$ 在第500步突變時RLS會在一小段延遲后重新逼近新的真值。這就是遺忘因子的作用。如果把 $\lambda$ 設成1你會看到突變后參數幾乎不動需要很久才能慢慢扭過去。6.3 踩坑記錄遺忘因子初值不當引發(fā)的發(fā)散實驗過程中最容易遇到的現象就是前幾步參數直接飛上天然后徹底發(fā)散。我早先試過把 $\lambda$ 設為0.9初始 $P(0)$ 設為 $1000I$結果第一個樣本就把參數修正得過猛因為 $K(1) P(0)\varphi_1 / (0.9 \varphi_1^T P(0)\varphi_1)$ 的分母雖然很大但分子也很大如果 $\varphi_1$ 的模很小導致增益變得巨大。解決辦法是適當減小 $P(0)$ 的初值或者約束參數更新范圍。這不是公式錯而是初值和遺忘因子搭配不當。換成 $\lambda 0.98$、$\delta 100$ 以后就穩(wěn)定多了。這類問題在公式推導時完全看不出來只有跑過代碼才深有體會。7. 復盤我對RLS公式推導和落地的一些經驗推完這一整套RLS公式我最大的感受是數學公式的每個變形都不是孤立的。矩陣求逆引理、遺忘因子、增益矩陣它們其實是同一個目標的不同側面。你不需要死記公式只需要記住兩條主線一條是信息矩陣 $R(t)$ 的遞推另一條是參數估計 $\hat{\theta}(t)$ 的遞推。所有的RLS變體都是圍繞這兩條主線做數值穩(wěn)定性和計算效率上的改進。實際調試中我習慣先在散點圖上畫出參數的收斂軌跡一旦發(fā)現軌跡異常就優(yōu)先檢查輸入信號的激勵程度——如果輸入一直恒定不變任何最小二乘類算法都會在某幾個方向上失去可辨識性這不是算法能救的。最后再分享一個小技巧如果系統(tǒng)是慢時變的但又怕 $\lambda$ 太大會跟不上可以用雙遺忘因子對輸入功率大的樣本用較大的 $\lambda$對輸入功率小的樣本用稍小的 $\lambda$這樣既能跟蹤突變又不會讓噪聲把參數抖得厲害。這個思路我也是從現場調參里慢慢摸出來的公式層面看不出來但工程上非常實用。