核自适应滤波器的递归最小二乘（RLS）实现与权重更新过程
题目描述
本题目讲解核自适应滤波器的一种经典实现方式：递归最小二乘算法。我们将聚焦于如何将传统的线性递归最小二乘（RLS）算法，通过核技巧映射到高维特征空间，从而处理非线性自适应滤波问题。核心内容包括：问题定义、核函数引入、核空间下的递归最小二乘权重更新推导，以及预测过程。
解题过程
1. 问题定义：非线性自适应滤波
假设我们有一个在线到达的输入-输出数据流：在时间点 \(n\)，我们观测到输入向量 \(\mathbf{u}(n) \in \mathbb{R}^d\) 和期望输出 \(d(n) \in \mathbb{R}\)。
目标：学习一个函数 \(f\)，使得预测输出 \(y(n) = f(\mathbf{u}(n))\) 尽可能接近 \(d(n)\)。
在线性自适应滤波中，我们假设 \(f\) 是线性的：\(y(n) = \mathbf{w}^T \mathbf{u}(n)\)，权重 \(\mathbf{w}\) 通过递归最小二乘（RLS）在线更新。
挑战：实际系统往往是非线性的，我们需要一个非线性函数 \(f\)。
核方法思路：将输入 \(\mathbf{u}\) 通过一个非线性映射 \(\phi(\mathbf{u})\) 映射到高维（可能是无限维）的再生核希尔伯特空间，然后在那个空间中进行线性处理：
\[y(n) = \langle \mathbf{w}, \phi(\mathbf{u}(n)) \rangle \]
其中 \(\mathbf{w}\) 是特征空间中的权重向量。由于 \(\phi(\cdot)\) 可能很复杂（甚至无限维），我们不会显式地计算它，而是通过核技巧来隐式处理。
2. 从线性RLS到核RLS
2.1 线性RLS回顾
线性RLS的目标是：最小化从初始时刻到当前时刻 \(n\) 的加权平方误差和：
\[J(\mathbf{w}) = \sum_{i=1}^{n} \lambda^{n-i} |d(i) - \mathbf{w}^T \mathbf{u}(i)|^2 + \delta \lambda^n \|\mathbf{w}\|^2 \]
其中：
- \(\lambda \in (0, 1]\) 是遗忘因子，让近期误差的权重更大。
- \(\delta > 0\) 是正则化参数（初始时保证逆矩阵可逆）。
其解为：
\[\mathbf{w}(n) = \mathbf{P}(n) \mathbf{r}(n) \]
其中：
- \(\mathbf{P}(n) = \left( \sum_{i=1}^{n} \lambda^{n-i} \mathbf{u}(i)\mathbf{u}(i)^T + \delta \lambda^n \mathbf{I} \right)^{-1}\) 是相关矩阵的逆。
- \(\mathbf{r}(n) = \sum_{i=1}^{n} \lambda^{n-i} d(i) \mathbf{u}(i)\) 是互相关向量。
RLS通过递归方式更新 \(\mathbf{P}(n)\) 和 \(\mathbf{w}(n)\)，避免每次直接求逆，计算复杂度为 \(O(d^2)\)。
2.2 引入核函数
在特征空间中，输入变为 \(\phi(\mathbf{u}(i))\)，权重 \(\mathbf{w}\) 可表示为历史样本的线性组合（由表示定理保证）：
\[\mathbf{w} = \sum_{i=1}^{n} a(i) \phi(\mathbf{u}(i)) \]
代入预测公式：
\[y(n) = \langle \mathbf{w}, \phi(\mathbf{u}(n)) \rangle = \sum_{i=1}^{n} a(i) \langle \phi(\mathbf{u}(i)), \phi(\mathbf{u}(n)) \rangle \]
定义核函数 \(k(\mathbf{u}(i), \mathbf{u}(j)) = \langle \phi(\mathbf{u}(i)), \phi(\mathbf{u}(j)) \rangle\)，则：
\[y(n) = \sum_{i=1}^{n} a(i) k(\mathbf{u}(i), \mathbf{u}(n)) \]
常用核函数：高斯核 \(k(\mathbf{u}, \mathbf{v}) = \exp(-\|\mathbf{u}-\mathbf{v}\|^2 / (2\sigma^2))\)。
目标转化为：求系数向量 \(\mathbf{a}(n) = [a(1), a(2), \dots, a(n)]^T\)，使得预测误差最小。
3. 核RLS的推导
在特征空间中，定义：
- 核矩阵 \(\mathbf{K}(n) \in \mathbb{R}^{n \times n}\)，其中 \([\mathbf{K}(n)]_{ij} = k(\mathbf{u}(i), \mathbf{u}(j))\)，对 \(1 \le i,j \le n\)。
- 期望输出向量 \(\mathbf{d}(n) = [d(1), d(2), \dots, d(n)]^T\)。
则特征空间中的加权正则化最小二乘问题为：
\[\min_{\mathbf{a}(n)} J = \|\mathbf{d}(n) - \mathbf{K}(n)\mathbf{a}(n)\|^2_{\mathbf{\Lambda}(n)} + \delta \lambda^n \mathbf{a}(n)^T \mathbf{K}(n) \mathbf{a}(n) \]
其中 \(\mathbf{\Lambda}(n) = \text{diag}([\lambda^{n-1}, \lambda^{n-2}, \dots, 1])\) 是加权矩阵。
令导数为零，得到解析解：
\[\mathbf{a}(n) = \left( \mathbf{K}(n)^T \mathbf{\Lambda}(n) \mathbf{K}(n) + \delta \lambda^n \mathbf{K}(n) \right)^{-1} \mathbf{K}(n)^T \mathbf{\Lambda}(n) \mathbf{d}(n) \]
如果 \(\mathbf{K}(n)\) 可逆，可简化为：
\[\mathbf{a}(n) = \left( \mathbf{K}(n) \mathbf{\Lambda}(n) \mathbf{K}(n) + \delta \lambda^n \mathbf{K}(n) \right)^{-1} \mathbf{K}(n) \mathbf{\Lambda}(n) \mathbf{d}(n) = \left( \mathbf{K}(n) \mathbf{\Lambda}(n) + \delta \lambda^n \mathbf{I} \right)^{-1} \mathbf{\Lambda}(n) \mathbf{d}(n) \]
但直接求逆计算量 \(O(n^3)\) 随 \(n\) 增长而剧增，不可行。因此需要递归更新。
4. 递归更新公式推导
核心思想：当新数据 \((\mathbf{u}(n), d(n))\) 到达时，我们递归地更新系数向量 \(\mathbf{a}(n)\) 和辅助矩阵 \(\mathbf{P}(n)\)（这里 \(\mathbf{P}(n)\) 是在特征空间中定义，与线性RLS中的类似但不同）。
步骤1：定义并更新 \(\mathbf{P}(n)\)
在特征空间中，定义相关矩阵：
\[\mathbf{\Phi}(n) = [\sqrt{\lambda^{n-1}}\phi(\mathbf{u}(1)), \dots, \phi(\mathbf{u}(n))]^T \]
则加权的自相关矩阵为 \(\mathbf{\Phi}(n)^T \mathbf{\Phi}(n) + \delta \lambda^n \mathbf{I}\)。但注意，我们不会显式计算 \(\phi\)，而是用核函数。实际推导中，我们定义：
\[\mathbf{P}(n) = \left( \mathbf{K}(n) \mathbf{\Lambda}(n) + \delta \lambda^n \mathbf{I} \right)^{-1} \]
其中 \(\mathbf{K}(n)\) 是 \(n \times n\) 的核矩阵，\(\mathbf{\Lambda}(n) = \text{diag}([\lambda^{n-1}, \dots, 1])\).
步骤2：递归更新 \(\mathbf{P}(n)\)
当新数据到来时，核矩阵扩展为：
\[\mathbf{K}(n) = \begin{bmatrix} \lambda \mathbf{K}(n-1) & \lambda^{1/2} \mathbf{k}(n-1) \\ \lambda^{1/2} \mathbf{k}(n-1)^T & k(\mathbf{u}(n), \mathbf{u}(n)) \end{bmatrix} \]
其中 \(\mathbf{k}(n-1) = [k(\mathbf{u}(1), \mathbf{u}(n)), \dots, k(\mathbf{u}(n-1), \mathbf{u}(n))]^T\)。
类似地，加权矩阵为 \(\mathbf{\Lambda}(n) = \begin{bmatrix} \lambda \mathbf{\Lambda}(n-1) & 0 \\ 0 & 1 \end{bmatrix}\).
通过分块矩阵求逆引理，可得到 \(\mathbf{P}(n)\) 的递归公式：
\[\mathbf{P}(n) = \frac{1}{\lambda} \begin{bmatrix} \mathbf{P}(n-1) & 0 \\ 0^T & 0 \end{bmatrix} + \mathbf{q}(n) \mathbf{q}(n)^T / r(n) \]
其中：
- \(r(n) = k(\mathbf{u}(n), \mathbf{u}(n)) + \delta \lambda^n - \mathbf{k}(n-1)^T \mathbf{P}(n-1) \mathbf{k}(n-1) / \lambda\) 是标量。
- \(\mathbf{q}(n) = \begin{bmatrix} -\mathbf{P}(n-1)\mathbf{k}(n-1)/\lambda \\ 1 \end{bmatrix}\).
步骤3：递归更新系数向量 \(\mathbf{a}(n)\)
系数向量 \(\mathbf{a}(n)\) 满足：
\[\mathbf{a}(n) = \mathbf{P}(n) \mathbf{\Lambda}(n) \mathbf{d}(n) \]
通过代入 \(\mathbf{P}(n)\) 的递归形式，可得 \(\mathbf{a}(n)\) 的更新：
\[\mathbf{a}(n) = \begin{bmatrix} \mathbf{a}(n-1) - \mathbf{q}(n) e(n) / r(n) \\ e(n) / r(n) \end{bmatrix} \]
其中 \(e(n) = d(n) - \mathbf{k}(n-1)^T \mathbf{a}(n-1)\) 是先验误差（用旧系数对新输入的预测误差）。
步骤4：预测输出
在时间 \(n+1\)，对新的输入 \(\mathbf{u}(n+1)\)，预测输出为：
\[y(n+1) = \mathbf{k}_n(\mathbf{u}(n+1))^T \mathbf{a}(n) \]
其中 \(\mathbf{k}_n(\mathbf{u}(n+1)) = [k(\mathbf{u}(1), \mathbf{u}(n+1)), \dots, k(\mathbf{u}(n), \mathbf{u}(n+1))]^T\).
5. 算法总结
- 初始化：\(\mathbf{a}(0) = []\)，\(\mathbf{P}(0) = 1/\delta\)（标量），\(n=1\)。
- 循环（对每个新样本 \((\mathbf{u}(n), d(n))\)）：
a. 计算核向量：\(\mathbf{k}(n-1) = [k(\mathbf{u}(i), \mathbf{u}(n))]_{i=1}^{n-1}\)。
b. 计算先验误差：\(e(n) = d(n) - \mathbf{k}(n-1)^T \mathbf{a}(n-1)\)。
c. 计算中间量：
\[ r(n) = k(\mathbf{u}(n), \mathbf{u}(n)) + \delta \lambda^n - \mathbf{k}(n-1)^T \mathbf{P}(n-1) \mathbf{k}(n-1) / \lambda \]
\[ \mathbf{q}(n) = \begin{bmatrix} -\mathbf{P}(n-1)\mathbf{k}(n-1)/\lambda \\ 1 \end{bmatrix} \]
d. 更新 \(\mathbf{P}(n)\)：
\[ \mathbf{P}(n) = \frac{1}{\lambda} \begin{bmatrix} \mathbf{P}(n-1) & 0 \\ 0^T & 0 \end{bmatrix} + \mathbf{q}(n) \mathbf{q}(n)^T / r(n) \]
e. 更新系数向量：
\[ \mathbf{a}(n) = \begin{bmatrix} \mathbf{a}(n-1) - \mathbf{q}(n) e(n) / r(n) \\ e(n) / r(n) \end{bmatrix} \]
- 预测：对新输入 \(\mathbf{u}_{\text{new}}\)，计算 \(y = \mathbf{k}_{\text{all}}(\mathbf{u}_{\text{new}})^T \mathbf{a}(\text{最新})\)。
核心要点
- 核RLS 将线性RLS扩展到非线性，通过核函数隐式地在高维特征空间中操作。
- 通过递归更新 \(\mathbf{P}(n)\) 和 \(\mathbf{a}(n)\)，避免了存储所有历史数据和求大规模矩阵逆，但计算复杂度仍为 \(O(n^2)\) 每步（因需计算核向量与历史所有样本的内积），因此常结合稀疏化技巧（如近似线性依赖）来控制增长。
- 该算法适用于在线、非线性、自适应滤波场景，如时间序列预测、系统辨识、信道均衡等。