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