第3章 最小二乘法（LS）与递推最小二乘法（RLS）

说到系统辨识，最小二乘法是绕不开的基石。我个人习惯把它叫做「参数估计的万金油」——简单、直观、数学上漂亮。你想想看，一堆带噪声的观测数据，怎么从中把系统的真实参数揪出来？最小二乘法就是干这个的。

这一章我们分两块走：先搞定批量最小二乘，再升级到递推版本。后者在机器人实时控制里特别实用，我当年做机械臂在线辨识时，RLS 几乎是标配。

3.1 批量最小二乘原理

先从一个最简单的场景说起。假设你有一个线性系统：

y = θ₁ * x₁ + θ₂ * x₂ + ... + θₙ * xₙ + e

其中 y 是输出，x 是输入，θ 是待估计的参数，e 是噪声。我们做了 m 次观测，写成矩阵形式：

Y = X * Θ + E

这里 Y 是 m×1 的输出向量，X 是 m×n 的输入矩阵（也叫设计矩阵），Θ 是 n×1 的参数向量，E 是噪声向量。

最小二乘法的核心思想很简单：让误差的平方和最小。定义损失函数：

J(Θ) = (Y - XΘ)ᵀ (Y - XΘ)

对 Θ 求导，令导数为零，得到正规方程：

XᵀX Θ = XᵀY

解出来就是：

Θ̂ = (XᵀX)⁻¹ XᵀY

这就是批量最小二乘的闭式解。嗯，这里要注意：XᵀX 必须可逆，也就是说输入数据要足够丰富，不能有线性相关。

3.2 矩阵推导的直观理解

你可能觉得公式推导有点抽象。我换个角度解释：

XΘ 是模型预测的输出，Y 是实际测量值。我们希望预测值尽可能接近真实值。说白了，就是在所有可能的 Θ 中，找一个让预测误差最小的那个。

从几何角度看，XΘ 是 X 的列空间中的一个向量，Y 是实际观测向量。最小二乘解其实就是 Y 在 X 列空间上的正交投影。这个投影向量就是 XΘ̂，而残差 Y - XΘ̂ 垂直于 X 的列空间。

这个几何直觉很重要。我在做机器人动力学参数辨识时，经常用这个视角检查数据质量——如果残差有明显的结构（不是白噪声），说明模型结构有问题，或者数据没覆盖到足够的激励空间。

3.3 Python 实现：批量最小二乘

直接上代码。我们用 numpy 实现一个简单的批量 LS 估计器：

import numpy as np
import matplotlib.pyplot as plt

def batch_ls(X, Y):
    """
    批量最小二乘估计
    X: 输入矩阵 (m x n)
    Y: 输出向量 (m x 1)
    返回: 参数估计向量 (n x 1)
    """
    # 检查 XᵀX 是否可逆
    XtX = X.T @ X
    if np.linalg.matrix_rank(XtX) < X.shape[1]:
        raise ValueError("XᵀX 奇异，输入数据不满足 PE 条件")
    
    # 计算参数估计
    theta_hat = np.linalg.inv(XtX) @ X.T @ Y
    return theta_hat

# 模拟数据
np.random.seed(42)
n_samples = 100
true_theta = np.array([2.5, -1.3, 0.8])

# 生成输入数据（满足 PE 条件）
X = np.random.randn(n_samples, 3)
# 生成输出数据，加入噪声
Y = X @ true_theta + 0.5 * np.random.randn(n_samples)

# 估计参数
theta_est = batch_ls(X, Y)
print("真实参数:", true_theta)
print("估计参数:", theta_est)
print("估计误差:", np.linalg.norm(theta_est - true_theta))

运行结果应该接近真实值。噪声越小，估计越准。我在实际项目中一般会先跑一遍批量 LS，看看数据质量，再决定要不要上更复杂的辨识方法。

3.4 递推最小二乘的动机

批量 LS 有个明显的缺点：每次有新数据，都要重新计算全部数据。如果数据量很大（比如每秒采样 1000 次），计算量会越来越大，内存占用也越来越高。

更关键的是，在机器人实时控制中，你不可能等收集完所有数据再辨识。你需要边运行边更新参数。比如机械臂的负载变化了，控制器要立刻调整参数，不能停下来重新算一遍。

这就是递推最小二乘（RLS）的用武之地。RLS 的核心思想是：用上一时刻的参数估计，结合当前的新数据，在线更新参数。每次更新只需要很小的计算量，O(n²) 级别，n 是参数个数。

我曾经在一个移动机器人项目里用 RLS 在线辨识轮胎-地面摩擦参数。机器人一边跑，一边更新参数，控制器根据最新参数调整转向策略。效果比固定参数好得多。

3.5 RLS 算法推导

RLS 的推导从批量 LS 的递推形式出发。假设我们有 k-1 个数据点的估计结果 Θ̂(k-1)，现在来了第 k 个数据点 (x(k), y(k))。

定义协方差矩阵 P(k) = (X(k)ᵀ X(k))⁻¹，其中 X(k) 是前 k 个数据组成的矩阵。可以推导出递推公式：

增益向量:   K(k) = P(k-1) x(k) / (1 + x(k)ᵀ P(k-1) x(k))
参数更新:   Θ̂(k) = Θ̂(k-1) + K(k) [y(k) - x(k)ᵀ Θ̂(k-1)]
协方差更新: P(k) = [I - K(k) x(k)ᵀ] P(k-1)

这个公式的直观理解是：

• 新息：y(k) - x(k)ᵀ Θ̂(k-1) 是当前预测误差
• 增益：K(k) 决定了新息对参数更新的权重
• 协方差：P(k) 反映了参数估计的不确定性，随着数据增多逐渐减小

3.6 Python 实现：递推最小二乘

下面是一个完整的 RLS 实现，包含遗忘因子（forgetting factor）的支持。遗忘因子 λ 让算法能跟踪时变参数，λ 越接近 1，对旧数据「遗忘」得越慢。

class RecursiveLeastSquares:
    """
    递推最小二乘估计器
    支持遗忘因子，适用于时变系统
    """
    def __init__(self, n_params, lambda_=0.99, P_init=1000.0):
        self.n = n_params
        self.lambda_ = lambda_  # 遗忘因子，0 < λ ≤ 1
        self.theta = np.zeros(n_params)  # 参数估计
        self.P = P_init * np.eye(n_params)  # 协方差矩阵初始化
        
    def update(self, x, y):
        """
        单步更新
        x: 输入向量 (n,)
        y: 输出标量
        """
        x = x.reshape(-1, 1)  # 转为列向量
        
        # 计算增益
        Px = self.P @ x
        denom = self.lambda_ + (x.T @ Px).item()
        K = Px / denom
        
        # 计算新息
        innovation = y - (x.T @ self.theta).item()
        
        # 更新参数
        self.theta = self.theta + K.flatten() * innovation
        
        # 更新协方差（带遗忘因子）
        self.P = (self.P - K @ x.T @ self.P) / self.lambda_
        
        return self.theta.copy()
    
    def estimate(self):
        return self.theta.copy()

# 使用示例：跟踪时变参数
rls = RecursiveLeastSquares(n_params=2, lambda_=0.95)

# 模拟时变系统
true_theta = np.array([1.0, 0.5])
theta_history = []
for t in range(200):
    # 参数在 t=100 时突变
    if t == 100:
        true_theta = np.array([2.0, -0.5])
    
    x = np.random.randn(2)
    y = x @ true_theta + 0.1 * np.random.randn()
    
    theta_est = rls.update(x, y)
    theta_history.append(theta_est.copy())

# 检查跟踪效果
print("最终估计:", rls.estimate())
print("真实参数:", true_theta)

注意遗忘因子的选择：λ=0.99 适合慢时变系统，λ=0.95 对突变响应更快但估计方差更大。我一般从 0.98 开始调，看实际跟踪效果再微调。

3.7 在线参数估计实战要点

把 RLS 用到实际系统时，有几个坑要避开：

问题	现象	解决方法
协方差爆炸	P 矩阵元素持续增长，增益过大	限制 P 的最大特征值，或使用带死区的更新
参数漂移	参数缓慢偏离真实值	加入遗忘因子，或使用带约束的 RLS
数值发散	P 矩阵失去正定性	使用平方根滤波或 UD 分解
激励不足	参数更新停滞	注入持续激励信号，或暂停更新

我曾经在一个项目中遇到协方差爆炸的问题。当时 RLS 跑了几分钟后，参数突然跳变到离谱的值。排查后发现是输入信号有一段长时间保持不变，导致 P 矩阵不断增长。后来加了协方差重置机制——如果 P 的最大特征值超过阈值，就把它重置为初始值。问题解决了。

3.8 知识体系总览

下面这张图总结了本章的核心逻辑：从批量 LS 到 RLS 的演进路径，以及各自的应用场景。

从图中可以看出，批量 LS 和 RLS 本质上是同一个思想的不同实现形式。批量 LS 适合离线分析，RLS 适合在线实时应用。实际项目中，我经常先用批量 LS 做一次离线辨识，得到初始参数，然后切换到 RLS 在线微调。这样既保证了初始精度，又能适应系统变化。

好了，这一章的内容就到这里。最小二乘法虽然基础，但吃透了它，后面学更高级的辨识方法（比如极大似然估计、子空间辨识）会轻松很多。记住：理解原理比背公式重要，动手实现比看代码重要。建议你把上面的代码跑一遍，改改参数，看看效果。