简介递推最小二乘法RLS的MATLAB实现源码包面向信号处理、控制理论与机器学习等领域的在线参数估计任务适合需要实时更新模型参数的研究者、工程师及高年级学生。与普通最小二乘相比RLS通过迭代更新参数并维护逆协方差矩阵具有收敛速度快、计算效率高的特点特别适合系统辨识、自适应滤波、在线预测等动态场景。压缩包内包含1个.m源文件整体仅2KB代码精简可直接在MATLAB中运行或嵌入现有工程便于逐行剖析算法核心。该实现覆盖初始化、预测值计算、误差求解、逆协方差矩阵更新、Kalman增益计算及参数估计递推等关键环节运行后输出参数估计序列可直观观察参数随数据更新的收敛过程也能对比不同初值或数据条件下的估计效果。已有2529人学习下载对于入门RLS算法、完成课程设计或作为在线辨识算法教学演示均有切实参考价值。1. 先说清楚RLS解决了什么问题做系统辨识、自适应滤波、在线参数估计的人十有八九都会在某个晚上对着递推最小二乘法RLS的公式发愁。我在研究生时第一次在MATLAB里写RLS当时只是照着论文抄代码跑通了也就没再管。后来工作里真正做传感器数据校核和在线辨识才发现这个算法的分量。很多人看教材觉得RLS只是一堆矩阵递推背下来公式就能用。但真要你从零写一份MATLAB源程序或者把算法移植到嵌入式设备上时问题全出来了初始化P矩阵给多大遗忘因子怎么取数据激励不够时为什么会发散等等。先归纳一下RLS到底解决什么问题。传统最小二乘法是批处理模式把所有历史数据组成矩阵求正规方程的解。假设你有N组数据模型有m个参数正规方程里涉及一个m×m的矩阵求逆。数据量少还好办一旦数据不断到来每次都要重新构建矩阵、重新求逆计算量越来越大而且必须保存所有历史数据。这在实时系统里完全不现实。RLS的核心思路是每次来一个新样本用上一时刻的参数估计和协方差矩阵通过一个递推公式修正参数不需要重新计算全局解。它本质上是把全局最小二乘解转换成递归更新解保证在每一步当前估计都等价于基于到当前时刻为止所有数据的加权最小二乘解。这也是为什么RLS的收敛速度远快于LMS因为它每一步都在用完整的二阶统计信息“校正”更新方向而不是像LMS那样只沿着瞬时梯度走一小步。实际工程里RLS最常见的三个出口是系统辨识、自适应滤波、在线数据校核。系统辨识需要在线估计被控对象参数自适应滤波需要跟踪非平稳信号的最优滤波器系数数据校核则是用模型预测值与实际测量值比较判断传感器是否漂移、系统是否出现故障。后面我会给一个数据校核的扩展场景那是RLS特别能出彩的地方。2. 核心公式推导看懂才敢改代码2.1 带遗忘因子的代价函数RLS的出发点是最小化这样一个带权重的代价函数J(θ) ∑_{i1}^{n} λ^{n-i} (d(i) - φ^T(i) θ)^2其中 n 是当前时刻λ 是遗忘因子取值范围一般是0.95到1之间。d(i) 是期望输出φ(i) 是回归向量θ 是要估计的参数向量。λ^{n-i} 的作用是给历史数据加权越旧的数据权重越小越新的数据权重越大。λ1 时就是标准最小二乘所有数据同等对待λ1 时算法具备跟踪时变参数的能力。这个形式和普通最小二乘的唯一区别就是多了一个指数衰减权重。为什么要这么做因为实际系统很少是定常的。比如一个电机负载发生变化它的时间常数和增益会慢慢变。如果我们还用所有历史数据做同等权重的估计旧数据会拖住参数让估计结果反应迟钝。加上遗忘因子相当于算法自动忘掉太久远的信息。2.2 矩阵求逆引理带来的递推形式对代价函数求导并令其为零可以得到正规方程R(n) θ(n) r(n)其中 R(n) ∑_{i1}^n λ^{n-i} φ(i) φ^T(i)r(n) ∑_{i1}^n λ^{n-i} φ(i) d(i)。如果每次都直接求逆那就失去了递推的意义。这里的关键技巧是矩阵求逆引理(A BCD)^{-1} A^{-1} - A^{-1}B(C^{-1} DA^{-1}B)^{-1}DA^{-1}令 P(n) R^{-1}(n)利用 R(n) λR(n-1) φ(n)φ^T(n)可以得到K(n) P(n-1)φ(n) / (λ φ^T(n)P(n-1)φ(n)) θ(n) θ(n-1) K(n)e(n) P(n) (P(n-1) - K(n)φ^T(n)P(n-1)) / λ其中 K(n) 称为增益向量e(n) d(n) - φ^T(n)θ(n-1) 是先验误差。我见过很多人卡在第一步就不愿意往下看了觉得引理推导太诡异。但你不需要记住引理的完整证明只需要知道它把一个高维矩阵求逆问题化成了一个标量除法问题。注意分母 λ φ^T P φ 是一个数不是矩阵所以代码里那行除法根本不涉及矩阵求逆。整个递推过程只涉及矩阵乘法和一个标量除法计算复杂度大约是 O(m²)m 是参数个数远低于每次重新求逆的 O(m³)。2.3 三个递推量的物理含义要真正掌握RLS不能只背公式。θ(n) 很好理解就是当前时刻的参数估计。P(n) 和 K(n) 的含义值得多说几句。P(n) 近似于参数估计误差协方差矩阵。P 大说明当前参数估计很不确定算法对新信息会更敏感P 收敛到很小说明估计已经很可靠算法不太愿意被新样本牵着走。这也是为什么初始化时 P(0) δI 要设一个较大的 δ比如10到1000表示我对初始参数完全没信心要靠数据尽快纠正。K(n) 是增益向量它决定了新样本误差 e(n) 以多大比例修正到参数上。K(n) 和 P(n-1) 成正比和 φ(n) 方向有关。如果本次的回归向量方向恰好是历史数据中信息量比较小的方向K会相对大一些如果这个方向已经被充分激励过K就小。这就是RLS自适应步长的由来它每个方向的修正步长不是均匀的而是根据历史的激励情况自动调节。对比LMSLMS用一个固定的标量 μ 乘梯度方向所以收敛慢而且各方向严重耦合。3. MATLAB源程序实现3.1 最小可运行的RLS函数先给一份干净的MATLAB源程序。不是我为了凑篇幅写的伪代码是真正能用、能跑、能改的直接版本。function [w, e, y_est] rls_identify(x, d, nOrder, lambda, delta) % RLS递推最小二乘辨识实现 % 输入: % x - 输入信号列向量长度N % d - 期望输出信号列向量长度N % nOrder - 模型阶数即参数个数m % lambda - 遗忘因子典型值 0.95~1 % delta - 初始协方差对角元素典型值 10~1000 % 输出: % w - 最终参数估计nOrder×1列向量 % e - 每一步的预测误差N×1列向量 % y_est - 每一步的模型预测输出N×1列向量 N length(x); w zeros(nOrder, 1); P delta * eye(nOrder); e zeros(N, 1); y_est zeros(N, 1); for n nOrder:N phi x(n:-1:n-nOrder1); % 构建回归向量 y_est(n) phi * w; % 先用当前参数预测 e(n) d(n) - y_est(n); % 计算先验误差 K P * phi / (lambda phi * P * phi); % 增益向量 w w K * e(n); % 更新参数 P (P - K * phi * P) / lambda; % 更新协方差矩阵 end end注意这段代码里我直接用了输入信号自身的延迟结构作为回归向量适合 FIR 模型和部分ARX模型辨识。如果你的问题不是时间序列而是普通的多元线性回归只需要把phi的构建方式换成对应的特征向量即可后面的递推核心完全不用动。3.2 逐段说明代码细节第一段初始化里w设为零向量这是最常用的做法表示对初始参数没有先验偏好。P delta * eye(nOrder)是RLS里最关键的初始化。δ不能太小否则后面来的新数据几乎改变不了参数算法会僵住但也不能太大极端情况下会让前几步的增益过大产生很大的参数跳变。循环部分需要注意从nOrder开始。因为构建回归向量需要充足的过去样本前面几个点凑不齐完整的历史数据。实际工程中如果信号初始段很重要你可以采用零填充的方式从第一点开始循环但一般的辨识场景直接从nOrder起步即可。phi x(n:-1:n-nOrder1)这行是在倒序取窗口。为什么倒序因为对于因果系统当前输出受过去输入影响最近的历史对应最大的权重。虽然从最小二乘解的角度不关心中间顺序但保持倒序可以让phi的第一个元素是当前最新输入方便和后续扩展的差分方程模型对应。增益计算的写法P * phi / (lambda phi * P * phi)我特意没用括号把分母括成一个标量再用inv或者\因为MATLAB里标量除法就是最直接最稳定的而且速度最快。这里phi * P * phi的结果是一个标量所以整个表达式没有矩阵求逆。协方差更新P (P - K * phi * P) / lambda注意必须先使用当前的P更新完毕再给下一轮用顺序不能乱。我见过有人把这一步和参数更新写反结果算出来的P和K不自洽算法行为完全不对劲。3.3 测试脚本识别一个未知FIR系统光有函数没有验证脚本代码再漂亮也白搭。下面是一份可以直接运行的测试脚本% test_rls.m clear; clc; close all; % 真实系统参数一个三阶FIR滤波器 b_true [0.8, -0.5, 0.2]; nOrder length(b_true); % 生成输入输出数据加一点噪声 N 500; x randn(N, 1); d filter(b_true, 1, x) 0.01 * randn(N, 1); % RLS参数 lambda 0.99; delta 10; % 调用RLS [w, e, y_est] rls_identify(x, d, nOrder, lambda, delta); % 画参数收敛轨迹 figure; for k 1:nOrder hold on; end脚本后半段画参数收敛轨迹需要保存每一轮的历史参数。理论上改造函数很简单把w在循环里往矩阵里塞即可。我更推荐单独写一个带历史记录版本的调试脚本避免把功能函数和可视化逻辑混在一起。4. 仿真实验与结果分析4.1 系统辨识例子用上面这段代码跑一轮你会发现RLS的收敛速度快得惊人。真实系数是[0.8, -0.5, 0.2]输入是白噪声输出加了方差0.01的高斯噪声在 λ0.99、δ10 的情况下大概前50到100个样本三个参数就从0迅速逼近真实值。大概到第200个样本时估计曲线已经基本稳住了最后的误差不再明显变化。你可以在同一个脚本里打印最终参数和真实参数对比disp([真实系数: , num2str(b_true)]); disp([估计系数: , num2str(w)]);配合好的噪声水平误差量级通常能在10^{-2}以下。注意这里说的误差取决于信噪比。如果噪声加大到0.1最终估计的抖动也会变大但收敛速度基本不受影响。这个性质非常实用因为很多工业场景中不需要等到系统完全平稳就能很快得到一个可用的模型。4.2 和LMS对比一下才理解RLS值在哪里LMS是很多自适应滤波教材的入门算法更新公式是w w μ * phi * e。它的计算量只有O(m)比RLS的O(m²)省不少但收敛性在不同特征值分布下差异很大。我做过一个对比实验同一个三阶系统输入不是白噪声而是有色噪声例如通过低通滤波后的信号LMS如果步长设小收敛非常慢步长设大又可能直接发散。而RLS不管输入信号的相关性强弱都能在200步以内收敛到接近真实值。原因是白噪声激励下系统所有模式都能被激活LMS的梯度方向比较正收敛自然快。一旦输入信号变成窄带或高度相关的信号LMS等效于在特征值散布很大的误差曲面上跑梯度下降步长只能迁就最小特征值方向导致整体收敛慢。RLS通过P矩阵对每一维做归一化相当于把误差曲面掰圆再走所以对输入信号的相关性不敏感。当然代价也有。RLS每一步要更新一个m×m的矩阵内存占用和计算量都大。对大多数嵌入式应用m在10以下时区别不明显m到50以上时你就得认真考虑是不是需要用快速RLS或QR-RLS这类变体了。5. 遗忘因子 λ 与参数初始化的工程心得5.1 如何选λ先给结论再解释遗忘因子λ决定了算法跟踪时变参数的能力和抗噪能力之间的平衡。工程上我经常用以下经验值场景λ建议范围等效记忆长度系统参数固定不变1全部历史数据缓慢漂移传感器漂移0.99 ~ 0.999100~1000个样本中等速度时变负载突变0.96 ~ 0.9925~100个样本快速时变强非平稳0.90 ~ 0.9510~20个样本等效记忆长度是1/(1-λ)意思是历史数据的有效权重衰减到初始的约36.8%时所需的样本数。λ0.98时等效记忆长度约50个样本λ0.99时约100个样本。这个指标能帮你快速判断当前采样率下参数变化时间尺度是100个样本还是1000个样本然后反过来选λ。如果你拿不准先设λ0.99跑一版看收敛后的参数波动。如果参数波动太大说明等效记忆长度太短噪声影响偏大就把λ往上调整到0.995甚至0.999。如果参数跟踪明显滞后比如真实参数已经变了估计值却还慢吞吞跟着旧值走说明记忆太长了把λ往下调。5.2 δ到底怎么设δ是P(0)的对角线元素它代表初始协方差的大小。δ的取值直接影响前期的收敛速度。理论分析表明P(0)δI 等价于在代价函数上加了正则项初始误差构成一个分布。如果δ取0增益K始终是0算法完全不会更新这是新手最容易踩的坑之一。常规做法是δ100左右。对于参数数量小于5的小模型100到1000都能接受前几步的收敛轨迹略有不同最终稳态基本一致。但δ也不能无脑大。我试过把δ设成10^6系统辨识在初始阶段会出现很大的参数跳变如果模型输出同时参与控制这种跳变可能会引入危险的瞬时控制量。一个更稳妥的做法是把δ和实际数据尺度挂钩。先大致估计一下回归向量的数量级。如果φ的最大值在10数量级δ设在10到100之间就够。如果φ来自归一化后的数据δ直接用10。5.3 工程调参清单调参的完整流程我总结成下面几步先用批处理最小二乘算一版参数作为参考判断模型结构是否合理是否欠阶。RLS参数先设λ1δ100跑模拟数据确认递推结果和批处理结果一致。加入遗忘因子从λ0.99开始用历史数据回放观察追踪效果。调整δ观察前50步收敛速度。δ太小则前期收敛慢δ太大会出现明显跳变。最后用一段独立测试数据验证不要只盯着训练段的拟合误差。6. 常见坑与排查记录6.1 算法发散P矩阵爆炸最典型的意外是跑着跑着参数突然冲上天P矩阵变成接近非正定或者数值上很大的矩阵。原因通常是 λ 太小同时又碰上一段数据激励不足。λ太小意味着算法对近期数据过度信任遗忘掉过去积累的信息P矩阵的数值会随递推不断膨胀。这种情况下如果输入信号有一段时间持续为0或者恒定φ^T P φ可能非常小K变得很大一个很小的噪声误差就能把参数推飞。解决办法不是简单减小λ而是检查信号激励质量。如果输入本身持续激励不足RLS再快也是无米之炊。工程上常用的做法是给输入叠加一个小幅随机扰动保证所有模态都被激发。如果需要保持λ很小以跟踪快变可以考虑带死区的更新策略当 |e(n)| 小于某个阈值时不更新P矩阵或者对P矩阵做上下限约束。6.2 参数收敛到错误值有一种情况是RLS收敛得很快但收敛到的值和真实参数不一致。最可能是模型结构不对。比如真正的系统是三阶你却用了二阶模型RLS会把二阶参数调到能最好拟合的位置看起来收敛了但预测效果有问题。另一种可能是噪声模型不匹配如果输出噪声是有色噪声并且和输入相关普通RLS的估计是有偏的。这时候需要扩展的递推最小二乘或者辅助变量法。检查方法很简单对比系统输出的仿真信号和实际输出。如果预测误差平稳且小参数可能还可以如果误差有系统性的低频分量多半是模型结构有问题而不是调参问题。6.3 数值稳定性问题m不超过10时直接用上面代码一般没问题。m大了之后经典的RLS递推会在长距离运行中出现数值退化。原因很简单P矩阵在数学上应该保持对称正定但有限精度计算下舍入误差会逐渐累积导致P不再对称甚至负定。这时候有两种做法。一种是每N步重新对称化P (P P) / 2。这是最简单的修补但不改变根本问题。更彻底的做法是使用平方根滤波或QR分解类RLS它们从设计上保证P的非负定性。如果只是在MATLAB里做仿真先用对称化顶一顶完全够如果要做嵌入式长跑最好从一开始就上平方根形式。6.4 RLS常见问题速查表现象可能原因解决方向参数前几步跳变极大δ太大调小δ参数一直不更新δ0或φ为零检查初始化检查输入信号稳态误差大λ太小增大λ提高输出信噪比参数跟踪滞后λ接近1减小λ缩短记忆长度运行几千步后发散P数值退化对称化或改用平方根RLS收敛但参数不对模型阶次不足提高阶次或改用辅助变量法7. 一个典型的工程延伸传感器数据校核提到RLS数据校核应用场景大多数人第一反应是“这不就是拿模型预测值和实测值对比吗”。实际操作中RLS可以做更细致的工作在线跟踪传感器的增益和零偏及时发现增益突变或漂移。假设某个测量通道的理论关系是 y a·x b其中a是增益b是零偏。传感器老化后a可能缓慢变化b也会漂移。如果直接用固定系数做校核正常变化会被误判为故障。用RLS在线估计a和b就能把“正常漂移”和“异常突变”区分开。代码实现只需要把回归向量改成两维phi [x(n); 1];其他递推完全不变。这里常数项1就是给零偏留出的系数。x(n)是实际激励y(n)是传感器输出。RLS估出的两个参数分别对应增益和零偏。你可以在线观察这两个参数的变化。我做这个场景时有个深刻体会遗忘因子的选择直接决定你能检测多快的故障。如果传感器漂移是缓慢的λ0.995等效记忆500个样本可以平滑跟踪漂移趋势。但如果传感器可能发生突发性跳变λ0.995就太迟钝了跳变发生后几十个样本才能跟上。要检测突变可以把λ设成0.95但付出的代价是估计噪声变大正常工况下参数抖动明显。折中的方案是做一个简单的变化检测连续若干个样本的误差都超过3倍标准差时把解释为“突变发生”然后将λ临时调到0.9快速追踪等估计稳定后再把λ调回去。RLS做数据校核还有一个容易被忽略的好处它天然给出模型预测残差 e(n)。这个残差比原始测量值更适合做报警判据因为它减掉了系统本身的动态。如果系统是一个惯性过程温度设定改变后输出会自然延迟上升直接用输出阈值很容易误报但用RLS模型预测残差动态部分被抵消只有真正偏离模型的行为才会让残差变大。8. 写在最后按我个人的经验学RLS最快的路径就是亲手在MATLAB里写一遍然后用一个带噪声的仿真系统验证再去尝试跟踪随时间变化的参数。代码跑通只是开始真正理解P矩阵、K向量和遗忘因子的相互作用才算把算法内化成自己的工具。之后遇到在线辨识、自适应滤波、数据校核这类问题你第一反应不会是去找工具箱而是直接想“RLS能不能解决、遗忘因子怎么配”。很多项目里一个几十行的RLS循环比复杂的机器学习模型实用得多因为它在嵌入式环境里跑得动在数据流上实时工作而且结果可解释。祝你把这份源程序跑出自己的心得。本文还有配套的精品资源点击获取