用Python从零实现带遗忘因子的递推最小二乘法（附完整代码与调参指南）
用Python从零实现带遗忘因子的递推最小二乘法（附完整代码与调参指南）
在动态系统中，参数往往会随时间变化。传统的最小二乘法对所有历史数据一视同仁，无法适应这种场景。本文将带你用Python从零实现带遗忘因子的递推最小二乘法(RLS)，这种算法能够"遗忘"旧数据，更关注近期观测，特别适合处理时变系统。
1. 算法核心原理
遗忘因子λ（0 < λ ≤ 1）决定了历史数据的衰减速度。当λ=1时退化为普通递推最小二乘；λ越小，对旧数据的遗忘速度越快。算法通过三个关键公式实现递推更新：
K_k = P_{k-1} @ x_k / (λ + x_k.T @ P_{k-1} @ x_k) P_k = (I - K_k @ x_k.T) @ P_{k-1} / λ θ_k = θ_{k-1} + K_k * (y_k - x_k.T @ θ_{k-1})其中P矩阵的初始化很关键。由于初始时缺乏数据，通常设P₀=αI（α取较大值如1e6），相当于给初始参数一个宽松的置信区间。
2. Python完整实现
我们先构建一个完整的RLS类，包含初始化、更新和预测功能：
import numpy as np class ForgettingRLS: def __init__(self, n_features, lambda_=0.95): self.n_features = n_features self.lambda_ = lambda_ # 遗忘因子 self.theta = np.zeros(n_features) # 参数向量 self.P = 1e6 * np.eye(n_features) # 初始协方差矩阵 def update(self, x, y): x = np.asarray(x).reshape(-1, 1) K = self.P @ x / (self.lambda_ + x.T @ self.P @ x).item() self.P = (np.eye(self.n_features) - K @ x.T) @ self.P / self.lambda_ error = y - (x.T @ self.theta).item() self.theta += K.flatten() * error return self.theta def predict(self, x): x = np.asarray(x).reshape(-1, 1) return (x.T @ self.theta).item()这个实现使用了NumPy进行矩阵运算，注意几个关键点：
- 使用
reshape(-1,1)确保输入向量是列向量 item()方法将1x1矩阵转换为标量- 矩阵乘法优先使用
@运算符
3. 实战演示：时变系统参数追踪
我们用一个逐渐变化的线性系统来测试算法效果。真实参数θ会随时间缓慢变化，观察RLS能否有效跟踪：
import matplotlib.pyplot as plt # 生成测试数据 np.random.seed(42) n_samples = 200 X = np.random.randn(n_samples, 2) # 输入特征 true_theta = np.array([1.0, -0.5]) # 初始真实参数 y = np.zeros(n_samples) for i in range(n_samples): # 参数缓慢变化 true_theta += 0.02 * np.random.randn(2) y[i] = X[i] @ true_theta + 0.1 * np.random.randn() # 使用不同λ值进行测试 lambdas = [0.9, 0.95, 0.99] estimates = {f"λ={lam}": [] for lam in lambdas} for lam in lambdas: rls = ForgettingRLS(n_features=2, lambda_=lam) for i in range(n_samples): rls.update(X[i], y[i]) estimates[f"λ={lam}"].append(rls.theta.copy())可视化结果时，我们发现λ=0.95在跟踪速度和稳定性之间取得了较好平衡：
plt.figure(figsize=(12, 6)) for i, lam in enumerate(lambdas, 1): est = np.array(estimates[f"λ={lam}"]) plt.subplot(2, 2, i) plt.plot(est[:, 0], label="θ₁估计") plt.plot(est[:, 1], label="θ₂估计") plt.title(f"λ={lam}参数估计轨迹") plt.legend() plt.tight_layout()4. 调参指南与常见问题
4.1 遗忘因子选择
λ的选择需要权衡：
- λ接近1：记忆性强，适合缓慢变化系统
- λ较小：适应性强，但对噪声更敏感
推荐采用以下策略确定λ：
| λ值范围 | 适用场景 | 优点 | 缺点 |
|---|---|---|---|
| 0.98-1.0 | 准静态系统 | 高精度 | 跟踪滞后 |
| 0.95-0.98 | 缓慢变化 | 平衡性好 | - |
| 0.9-0.95 | 快速变化 | 响应快 | 噪声敏感 |
| <0.9 | 剧烈变化 | 极快适应 | 不稳定 |
提示：可以从λ=0.95开始，根据实际效果微调
4.2 数值稳定性处理
长期运行可能导致P矩阵失去正定性。可定期进行以下维护：
# 对称化处理 self.P = (self.P + self.P.T) / 2 # 重置策略（当迹过大时） if np.trace(self.P) > 1e8: self.P = 1e6 * np.eye(self.n_features)4.3 常见问题排查
参数发散：
- 检查λ是否过小
- 确认输入数据是否标准化
- 尝试减小P的初始值
跟踪滞后：
- 适当减小λ值
- 检查系统变化速度是否超出算法能力
矩阵奇异：
- 确保输入数据有足够激励
- 添加小的正则化项：
P += 1e-6 * np.eye(n_features)
5. 高级应用：非线性系统处理
对于非线性系统，可通过基函数扩展实现：
def rbf_features(x, centers): """径向基函数特征扩展""" return np.exp(-0.5 * np.sum((x - centers)**2, axis=1)) # 使用示例 centers = np.linspace(-3, 3, 10) x = np.array([1.5]) phi = rbf_features(x, centers) # 转换为高维特征 rls = ForgettingRLS(n_features=len(centers)) rls.update(phi, y_observed)这种处理方式使RLS能够逼近非线性函数，同时保持在线更新的优势。