尧图网站建设 
很多人在学习系统辨识、自适应滤波或者在线参数估计时都会卡在递推最小二乘法Recursive Least Squares, RLS的公式推导上。教材里通常几行带过但实际自己要推一遍或者写代码实现时才发现从“批量求逆”到“递推更新”之间其实隔着不少细节。我最早接触RLS是在做在线辨识的时候当时对着讲义看了好几遍总觉得增益矩阵像个魔法——为什么这个式子加个新息就能不断修正参数为什么协方差矩阵要按那个方式更新后来自己动手推了一遍又把初值、遗忘因子和各种工程变形都玩了一遍才真正觉得这公式“通了”。这篇就把我的推导过程和实践心得完整拆一拆希望能帮你一次看透RLS。1. 从批量最小二乘到逐点递归先弄清楚RLS到底在解什么问题1.1 经典最小二乘闭式解到底长什么样绝大多数人第一次接触最小二乘都是从线性回归开始的。假设系统模型为$$y \varphi^T \theta$$其中 $\varphi$ 是回归向量输入特征$\theta$ 是待估计的参数向量。如果我们拿到的是一批数据 ${(\varphi_i, y_i)}_{i1}^N$那么可以写成矩阵形式$$Y \Phi \theta e$$$Y$ 是输出向量$\Phi$ 的每一行是一个回归向量 $\varphi_i^T$$e$ 是残差。批量最小二乘的优化目标是让残差平方和最小$$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)$。假设你每一秒得到一个样本每秒都要重新算一次求逆硬件功耗和实时性都会吃不消。就算用聪明的矩阵求逆算法反复做全量计算仍然很浪费。第二个尴尬更隐秘批量最小二乘对历史数据一视同仁。如果系统参数本身是缓变的比如飞机飞行时的气动参数随高度变化电机绕组电阻随温度变化那么旧数据对当前时刻的估计其实已经没有太大参考价值了。批量最小二乘只会把旧数据和新数据混在一起最终得到一个“历史平均”的参数根本追不上系统变化。1.3 RLS的核心思路一句话讲完RLS的核心思路其实一句话就能说清楚把上一时刻的参数估计当作基础当新样本到来时用这个样本带来的“新息”innovation对参数进行修正同时通过一个递推式更新协方差矩阵避免显式地重算矩阵求逆。换句话说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的协方差更新中帮我们把一个 $n \times n$ 的矩阵求逆转化成括号内一个标量或小矩阵的求逆。2.2 怎么把引理套进RLS的协方差更新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]$$你会发现原来需要对整个 $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$ 时老数据按指数衰减系统参数变化时估计值能更快跟上。令 $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$ 求导等于零得到最优解$$\hat{\theta}(t) R(t)^{-1} Q(t)$$这就是RLS推导的出发点。接下来要做的事是把 $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节已经推导出来了。把 $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公式中最容易看头晕的地方。我之前卡了很久后来发现只需要盯着 $P(t)$ 的递推式把所有 $P(t)$ 都替换成在第2节得到的结果大部分中间项会自动消掉。记住一个关键点分母是标量所以“求逆”根本不是真正的矩阵求逆而是一次普通除法。3.3 协方差更新公式的两种等价写法有了 $K(t)$$P(t)$ 的更新式可以写成更紧凑的形式。第2节的结果其实等价于$$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$ 是一个秩一矩阵意味着一维新样本对高维协方差矩阵的修正是“秩一更新”。这个结构的几何直觉是只有平行于 $\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最小均方算法收敛更快的原因——它把历史数据的二阶统计信息全部压缩在 $P$ 矩阵中。4. 初值设定、遗忘因子和完整的RLS算法流程4.1 协方差矩阵初值怎么给到底该用大数还是小数代码实现RLS时第一个问题就是 $P(0)$ 怎么设。大多数教材推荐 $P(0) \delta I$$\delta$ 取一个较大的数比如 $100$ 或 $10^3$。这个做法的理由是$P$ 矩阵本身可以理解为参数估计协方差矩阵的近似初始状态下我们对参数几乎一无所知所以把它的“不确定性”设置得很大让前几个样本能快速修正参数。另一种场景是如果你已经有一个比较靠谱的先验参数估计比如上一批次辨识好的参数那 $\delta$ 就可以取小一些比如 $0.1$ 甚至 $0.01$这样初期就不会因为第一个样本就把参数拉飞。还有一种常见做法是用一批小数据先做一次批量最小二乘得到初始值再把对应的协方差矩阵拷贝给 $P(0)$。这个做法最稳但对于在线系统来说不一定有这么多先验数据。4.2 遗忘因子 $\lambda$ 的选择内存长度和跟踪速度的权衡$\lambda$ 的取值直接影响算法对时变系统的跟踪能力。$\lambda 1$ 时算法拥有无限记忆适合参数恒定的系统$\lambda 1$ 时相当于对不同时刻的数据赋了一个指数衰减权重。工程上常把“有效记忆长度”近似为$$N_{\text{eff}} \approx \frac{1}{1 - \lambda}$$比如 $\lambda 0.99$ 大约相当于只有最近100个样本在起作用$\lambda 0.95$ 则只有大约20个样本。$\lambda$ 越小跟踪越快但受噪声影响也越大$\lambda$ 太小时估计方差会明显增加。所以实际使用时需要根据系统的变化速度和噪声水平折中我的经验是先从 $\lambda 0.98$ 左右试起然后观察估计曲线的抖动幅度再逐步调整。宁可让跟踪慢一点也别让估计值抖成心电图。4.3 标准RLS算法流程伪代码形式标准RLS算法每来一个新样本只需要五步初始化 theta zeros(n, 1) P delta * eye(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. 更新协方差 P (P - K * phi^T * P) / lambda 5. 检查 P 是否为对称正定工程上常用 P (P P^T) / 2 强制对称注意第4步的除法 $\lambda$ 是针对标量的可以直接除在矩阵上。如果想避免第5步的额外检查也可以用平方根RLS之类的变形这一点我会在下一章详细说。5. 数值稳定性问题和工程改进措施5.1 为什么P矩阵会变成非正定运行RLS时间长了你可能会发现本来应该正定的协方差矩阵 $P(t)$ 会逐渐失去对称正定性甚至出现负特征值。原因主要有三个第一计算机有限字长误差。RLS递推中反复做 $K(t)\varphi_t^T P(t-1)$ 这种矩阵乘法每一步都有舍入误差长时间累积后可能导致对称性破坏。第二遗忘因子导致“旧信息被指数衰减”当信号激励不足时$P(t)$ 的某些方向会不断被放大或缩小最终变得病态甚至非正定。第三输入 $\varphi_t$ 持续相关或者某一段激励太弱也会让信息矩阵长时间不增长$P$ 就可能在数值上发散。一旦 $P$ 失去正定性参数估计可能会出现剧烈跳变增益向量 $K(t)$ 也可能出现符号异常。建议是在算法里添加监控比如检查 $P$ 的对称性和特征值如果出现非正定及时重置或者强制对称化。5.2 工程上常用的修正策略最简单的修正是每次更新完以后执行$$P_{\text{sym}} \frac{P P^T}{2}$$这个操作不会带来太大成本却能有效避免因非对称导致的累积性误差。其次可以引入正则化项在信息矩阵上叠加一个小量 $\epsilon I$这样 $P^{-1}$ 始终有界。也可以设置一个“更新判定条件”当新息的绝对值特别大时暂时不更新协方差矩阵只更新参数避免异常样本对 $P$ 造成污染。更稳定的做法是使用平方根RLSSquare-Root RLS它把 $P$ 分解为 $P S S^T$递推更新 $S$ 而不是更新 $P$ 本身。因为 $S$ 的特征值都是正的只要 $S$ 不溢出$P$ 就能一直保持正定。平方根RLS的代价是额外多一些三角运算但对长期运行的在线系统来说非常值得。5.3 平方根RLS的基本思想简述平方根RLS的核心是用QR分解的思想更新信息矩阵的平方根。具体来说将 $P(t)^{1/2}$ 作为递推量利用旋转矩阵Givens旋转或Householder变换使得更新后的平方根矩阵保持正定。标准RLS中的分母 $\lambda \varphi_t^T P(t-1) \varphi_t$ 可以看作一个标量在平方根版本中它将合并进一个增广矩阵的QR更新过程。实现细节比较繁琐但如果你的系统要求长时间不间断运行强烈建议不要直接用基础RLS而是上平方根版本。6. 一段可运行的Python实验验证RLS公式有没有推错6.1 实验设计静态参数和时变参数两个场景纸上推了半天还是得用代码验证。这里我做一个最简单的单输入单输出SISO系统辨识实验模型是$$y_t a x_t b n_t$$其中 $a1.2$$b0.8$$n_t$ 是均值为0、方差为0.01的高斯噪声。回归向量取 $\varphi_t [x_t, 1]^T$参数向量 $\theta [a, b]^T$。先测试静态参数下RLS是否收敛到真实值再设计一个时变参数场景比如 $a$ 在第500个样本时从1.2跳到0.5检验遗忘因子能不能帮算法跟上这个突变。6.2 核心代码和结果解读下面是一份非常精简的Python实现直接用基础RLS没有花哨优化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 踩坑记录遗忘因子初值不当引发的发散实验过程中最容易遇到的现象就是前几步参数直接飞上天然后彻底发散。我早先试过把 $\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$ 以后就稳定多了。这类问题在公式推导时完全看不出来只有跑过代码才深有体会。7. 复盘我对RLS公式推导和落地的一些经验推完这一整套RLS公式我最大的感受是数学公式的每个变形都不是孤立的。矩阵求逆引理、遗忘因子、增益矩阵它们其实是同一个目标的不同侧面。你不需要死记公式只需要记住两条主线一条是信息矩阵 $R(t)$ 的递推另一条是参数估计 $\hat{\theta}(t)$ 的递推。所有的RLS变体都是围绕这两条主线做数值稳定性和计算效率上的改进。实际调试中我习惯先在散点图上画出参数的收敛轨迹一旦发现轨迹异常就优先检查输入信号的激励程度——如果输入一直恒定不变任何最小二乘类算法都会在某几个方向上失去可辨识性这不是算法能救的。最后再分享一个小技巧如果系统是慢时变的但又怕 $\lambda$ 太大会跟不上可以用双遗忘因子对输入功率大的样本用较大的 $\lambda$对输入功率小的样本用稍小的 $\lambda$这样既能跟踪突变又不会让噪声把参数抖得厉害。这个思路我也是从现场调参里慢慢摸出来的公式层面看不出来但工程上非常实用。 递推最小二乘RLS公式推导与工程落地详解
很多人在学习系统辨识、自适应滤波或者在线参数估计时都会卡在递推最小二乘法Recursive Least Squares, RLS的公式推导上。教材里通常几行带过但实际自己要推一遍或者写代码实现时才发现从“批量求逆”到“递推更新”之间其实隔着不少细节。我最早接触RLS是在做在线辨识的时候当时对着讲义看了好几遍总觉得增益矩阵像个魔法——为什么这个式子加个新息就能不断修正参数为什么协方差矩阵要按那个方式更新后来自己动手推了一遍又把初值、遗忘因子和各种工程变形都玩了一遍才真正觉得这公式“通了”。这篇就把我的推导过程和实践心得完整拆一拆希望能帮你一次看透RLS。1. 从批量最小二乘到逐点递归先弄清楚RLS到底在解什么问题1.1 经典最小二乘闭式解到底长什么样绝大多数人第一次接触最小二乘都是从线性回归开始的。假设系统模型为$$y \varphi^T \theta$$其中 $\varphi$ 是回归向量输入特征$\theta$ 是待估计的参数向量。如果我们拿到的是一批数据 ${(\varphi_i, y_i)}_{i1}^N$那么可以写成矩阵形式$$Y \Phi \theta e$$$Y$ 是输出向量$\Phi$ 的每一行是一个回归向量 $\varphi_i^T$$e$ 是残差。批量最小二乘的优化目标是让残差平方和最小$$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)$。假设你每一秒得到一个样本每秒都要重新算一次求逆硬件功耗和实时性都会吃不消。就算用聪明的矩阵求逆算法反复做全量计算仍然很浪费。第二个尴尬更隐秘批量最小二乘对历史数据一视同仁。如果系统参数本身是缓变的比如飞机飞行时的气动参数随高度变化电机绕组电阻随温度变化那么旧数据对当前时刻的估计其实已经没有太大参考价值了。批量最小二乘只会把旧数据和新数据混在一起最终得到一个“历史平均”的参数根本追不上系统变化。1.3 RLS的核心思路一句话讲完RLS的核心思路其实一句话就能说清楚把上一时刻的参数估计当作基础当新样本到来时用这个样本带来的“新息”innovation对参数进行修正同时通过一个递推式更新协方差矩阵避免显式地重算矩阵求逆。换句话说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的协方差更新中帮我们把一个 $n \times n$ 的矩阵求逆转化成括号内一个标量或小矩阵的求逆。2.2 怎么把引理套进RLS的协方差更新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]$$你会发现原来需要对整个 $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$ 时老数据按指数衰减系统参数变化时估计值能更快跟上。令 $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$ 求导等于零得到最优解$$\hat{\theta}(t) R(t)^{-1} Q(t)$$这就是RLS推导的出发点。接下来要做的事是把 $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节已经推导出来了。把 $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公式中最容易看头晕的地方。我之前卡了很久后来发现只需要盯着 $P(t)$ 的递推式把所有 $P(t)$ 都替换成在第2节得到的结果大部分中间项会自动消掉。记住一个关键点分母是标量所以“求逆”根本不是真正的矩阵求逆而是一次普通除法。3.3 协方差更新公式的两种等价写法有了 $K(t)$$P(t)$ 的更新式可以写成更紧凑的形式。第2节的结果其实等价于$$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$ 是一个秩一矩阵意味着一维新样本对高维协方差矩阵的修正是“秩一更新”。这个结构的几何直觉是只有平行于 $\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最小均方算法收敛更快的原因——它把历史数据的二阶统计信息全部压缩在 $P$ 矩阵中。4. 初值设定、遗忘因子和完整的RLS算法流程4.1 协方差矩阵初值怎么给到底该用大数还是小数代码实现RLS时第一个问题就是 $P(0)$ 怎么设。大多数教材推荐 $P(0) \delta I$$\delta$ 取一个较大的数比如 $100$ 或 $10^3$。这个做法的理由是$P$ 矩阵本身可以理解为参数估计协方差矩阵的近似初始状态下我们对参数几乎一无所知所以把它的“不确定性”设置得很大让前几个样本能快速修正参数。另一种场景是如果你已经有一个比较靠谱的先验参数估计比如上一批次辨识好的参数那 $\delta$ 就可以取小一些比如 $0.1$ 甚至 $0.01$这样初期就不会因为第一个样本就把参数拉飞。还有一种常见做法是用一批小数据先做一次批量最小二乘得到初始值再把对应的协方差矩阵拷贝给 $P(0)$。这个做法最稳但对于在线系统来说不一定有这么多先验数据。4.2 遗忘因子 $\lambda$ 的选择内存长度和跟踪速度的权衡$\lambda$ 的取值直接影响算法对时变系统的跟踪能力。$\lambda 1$ 时算法拥有无限记忆适合参数恒定的系统$\lambda 1$ 时相当于对不同时刻的数据赋了一个指数衰减权重。工程上常把“有效记忆长度”近似为$$N_{\text{eff}} \approx \frac{1}{1 - \lambda}$$比如 $\lambda 0.99$ 大约相当于只有最近100个样本在起作用$\lambda 0.95$ 则只有大约20个样本。$\lambda$ 越小跟踪越快但受噪声影响也越大$\lambda$ 太小时估计方差会明显增加。所以实际使用时需要根据系统的变化速度和噪声水平折中我的经验是先从 $\lambda 0.98$ 左右试起然后观察估计曲线的抖动幅度再逐步调整。宁可让跟踪慢一点也别让估计值抖成心电图。4.3 标准RLS算法流程伪代码形式标准RLS算法每来一个新样本只需要五步初始化 theta zeros(n, 1) P delta * eye(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. 更新协方差 P (P - K * phi^T * P) / lambda 5. 检查 P 是否为对称正定工程上常用 P (P P^T) / 2 强制对称注意第4步的除法 $\lambda$ 是针对标量的可以直接除在矩阵上。如果想避免第5步的额外检查也可以用平方根RLS之类的变形这一点我会在下一章详细说。5. 数值稳定性问题和工程改进措施5.1 为什么P矩阵会变成非正定运行RLS时间长了你可能会发现本来应该正定的协方差矩阵 $P(t)$ 会逐渐失去对称正定性甚至出现负特征值。原因主要有三个第一计算机有限字长误差。RLS递推中反复做 $K(t)\varphi_t^T P(t-1)$ 这种矩阵乘法每一步都有舍入误差长时间累积后可能导致对称性破坏。第二遗忘因子导致“旧信息被指数衰减”当信号激励不足时$P(t)$ 的某些方向会不断被放大或缩小最终变得病态甚至非正定。第三输入 $\varphi_t$ 持续相关或者某一段激励太弱也会让信息矩阵长时间不增长$P$ 就可能在数值上发散。一旦 $P$ 失去正定性参数估计可能会出现剧烈跳变增益向量 $K(t)$ 也可能出现符号异常。建议是在算法里添加监控比如检查 $P$ 的对称性和特征值如果出现非正定及时重置或者强制对称化。5.2 工程上常用的修正策略最简单的修正是每次更新完以后执行$$P_{\text{sym}} \frac{P P^T}{2}$$这个操作不会带来太大成本却能有效避免因非对称导致的累积性误差。其次可以引入正则化项在信息矩阵上叠加一个小量 $\epsilon I$这样 $P^{-1}$ 始终有界。也可以设置一个“更新判定条件”当新息的绝对值特别大时暂时不更新协方差矩阵只更新参数避免异常样本对 $P$ 造成污染。更稳定的做法是使用平方根RLSSquare-Root RLS它把 $P$ 分解为 $P S S^T$递推更新 $S$ 而不是更新 $P$ 本身。因为 $S$ 的特征值都是正的只要 $S$ 不溢出$P$ 就能一直保持正定。平方根RLS的代价是额外多一些三角运算但对长期运行的在线系统来说非常值得。5.3 平方根RLS的基本思想简述平方根RLS的核心是用QR分解的思想更新信息矩阵的平方根。具体来说将 $P(t)^{1/2}$ 作为递推量利用旋转矩阵Givens旋转或Householder变换使得更新后的平方根矩阵保持正定。标准RLS中的分母 $\lambda \varphi_t^T P(t-1) \varphi_t$ 可以看作一个标量在平方根版本中它将合并进一个增广矩阵的QR更新过程。实现细节比较繁琐但如果你的系统要求长时间不间断运行强烈建议不要直接用基础RLS而是上平方根版本。6. 一段可运行的Python实验验证RLS公式有没有推错6.1 实验设计静态参数和时变参数两个场景纸上推了半天还是得用代码验证。这里我做一个最简单的单输入单输出SISO系统辨识实验模型是$$y_t a x_t b n_t$$其中 $a1.2$$b0.8$$n_t$ 是均值为0、方差为0.01的高斯噪声。回归向量取 $\varphi_t [x_t, 1]^T$参数向量 $\theta [a, b]^T$。先测试静态参数下RLS是否收敛到真实值再设计一个时变参数场景比如 $a$ 在第500个样本时从1.2跳到0.5检验遗忘因子能不能帮算法跟上这个突变。6.2 核心代码和结果解读下面是一份非常精简的Python实现直接用基础RLS没有花哨优化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 踩坑记录遗忘因子初值不当引发的发散实验过程中最容易遇到的现象就是前几步参数直接飞上天然后彻底发散。我早先试过把 $\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$ 以后就稳定多了。这类问题在公式推导时完全看不出来只有跑过代码才深有体会。7. 复盘我对RLS公式推导和落地的一些经验推完这一整套RLS公式我最大的感受是数学公式的每个变形都不是孤立的。矩阵求逆引理、遗忘因子、增益矩阵它们其实是同一个目标的不同侧面。你不需要死记公式只需要记住两条主线一条是信息矩阵 $R(t)$ 的递推另一条是参数估计 $\hat{\theta}(t)$ 的递推。所有的RLS变体都是围绕这两条主线做数值稳定性和计算效率上的改进。实际调试中我习惯先在散点图上画出参数的收敛轨迹一旦发现轨迹异常就优先检查输入信号的激励程度——如果输入一直恒定不变任何最小二乘类算法都会在某几个方向上失去可辨识性这不是算法能救的。最后再分享一个小技巧如果系统是慢时变的但又怕 $\lambda$ 太大会跟不上可以用双遗忘因子对输入功率大的样本用较大的 $\lambda$对输入功率小的样本用稍小的 $\lambda$这样既能跟踪突变又不会让噪声把参数抖得厉害。这个思路我也是从现场调参里慢慢摸出来的公式层面看不出来但工程上非常实用。