RLS递归最小二乘法(Recursive Least Squares)

感谢B站Up 凩子白的讲解视频, 大多数的RLS算法介绍都是从各种专业领域角度讲解的(比如滤波器等角度), 对于缺乏专业背景的同学入门较难, 本文希望单纯从数学角度出发，便于缺乏信号处理等领域背景知识的同学理解。本文主要是对以上提到的视频的文字化, 同时也加入了自己的一些理解, 也许有一些地方不是那么严谨, 但希望能帮助大家快速了解一下RLS算法的基本思想。

PRELIMINARIES

最小二乘法

对于样本数据对儿$(\mathbf{x},y)$, 其中输入数据向量$\mathbf{x}=[x_{11},x_{12},...,x_{1m}{]}^T\in {\mathbb{R}}^m$, 输出样本为$y\in \mathbb{R}$; 使用参数为$\mathbf{w}$ 的模型来拟合数据$(\mathbf{x},y)$之间的真实映射关系; 认为模型$\mathbf{w}$的输出为$y$的估计值$\hat{y}\in \mathbb{R}$, 满足$\hat{y}\sim f(\mathbf{w};\mathbf{x})$, 拟合模型满足如下形式
$$\begin{gathered} \widehat{y_1}=w_1{x}_{11}+w_2{x}_{12}+...w_m{x}_{1m}={\mathbf{x}}_1^{\mathbf{T}}\mathbf{w} \\ \widehat{y_2}={\mathbf{x}}_2^{\mathbf{T}}\mathbf{w} \\ ⋮ \\ \widehat{y_n}={\mathbf{x}}_{\mathbf{n}}^{\mathbf{T}}\mathbf{w}\text{\hspace{1em}(1)} \end{gathered}$$
最小二乘法的思路, 就是希望近似模型参数$\mathbf{w}$在这$n$个样本输入数据$X_{n\times m}=[{\mathbf{x}}_1,{\mathbf{x}}_2,...,{\mathbf{x}}_{\mathbf{n}}{]}^T$ (以后简记为$X$)上得出的估计值$\hat{\mathbf{y}}=[\widehat{y_1},\widehat{y_2},...,\widehat{y_n}{]}^T$ 与ground truth 输出样本数据$\mathbf{y}=[y_1,y_2,...,y_n{]}^T$ 之间的差值平方和最小,即
$$\begin{gathered} \mathbf{w}=\arg min\sum _{i=1}^n(y_i-\widehat{y_i}{)}^2 \\ =\arg min\sum _{i=1}^n(y_i-{\mathbf{x}}_{\mathbf{i}}^{\mathbf{T}}\mathbf{w}{)}^2 \\ =\arg min{\Vert \begin{array}{c}\mathbf{y}-X\mathbf{w}\end{array}\Vert }_2^2 \\ =\arg minE(\mathbf{w})\text{\hspace{1em}(2)} \end{gathered}$$
误差$E(\mathbf{w})$对参数$\mathbf{w}$ 求梯度,
$$\begin{gathered} {\nabla }_{\mathbf{w}}E={\nabla }_{\mathbf{w}}{\Vert \begin{array}{c}\mathbf{y}-X\mathbf{w}\end{array}\Vert }_2^2 \\ ={\nabla }_{\mathbf{w}}[(\mathbf{y}-X\mathbf{w}{)}^T(\mathbf{y}-X\mathbf{w})] \\ =2X^T(\mathbf{y}-X\mathbf{w})\text{\hspace{1em}(3)} \end{gathered}$$
令${\nabla }_{\mathbf{w}}E=0$, 即可求出
$$\mathbf{w}=X^{-1}\mathbf{y}=(X^{T}X{)}^{-1}{X}^T\mathbf{y}\text{\hspace{1em}(4)}$$
注意, 公式$3$中的$\mathbf{y}$是样本数据, 将$X^{-1}$表述为$(X^{T}X{)}^{-1}{X}^T$的原因是矩阵$X$不一定是$n\times n$形状,因此不一定有逆矩阵, 而$X^{T}X$的逆是存在的?

最小二乘法存在的问题（为什么需要递归最小二乘法）

当一次性给出所有样本集合$(X,\mathbf{y})$时, 可以通过公式$4$来直接计算出最优的拟合模型参数$\mathbf{w}$ , 然而, 在实际应用中, 这种直接计算法并不常见, 主要是因为公式中求逆部分$(X^{T}X{)}^{-1}$ 的计算量大, 在样本数据量大时计算量更是明显增大; 另外, 现实生活中,往往出现样本数据可能也并不是一次性给出, 而是不断给出新的样本数据, 以一种数据流的形式给出样本数据, 例如传感器随时间不断读取信号等, 这种情况下利用公式$4$直接计算最优模型参数$\mathbf{w}$ 就需要每次进行直接计算, 也是不现实的.

因此, 为了利用最小误差平方和原则, 求解在大样本量, 或者数据流情况下的最优模型参数$\mathbf{w}$, 一种方法可以将大样本分成多批次(batch), 计算旧模型在新批次样本上的梯度, 不断进行梯度下降来进行迭代求解(也可以将数据流当做一个个batch来梯度更新); 另一种则是解析的方法, 就是这里提到的递归最小二乘解法(RLS).

本质上, 递归最小二乘法RLS和梯度下降、直接计算法一样, 都是为了求解满足最小误差平方和原则的最优模型参数$\mathbf{w}$, 只是在实现方式上有所不同.

递归最小二乘法

如之前提到的, RLS的主要应用场景, 是假设输入样本数据$X$在不断添加新数据, 例如$X=[{\mathbf{x}}_1,{\mathbf{x}}_2,...,{\mathbf{x}}_{\mathbf{n}}{]}^T\to X^{\prime }=[{\mathbf{x}}_1,{\mathbf{x}}_2,...,{\mathbf{x}}_{\mathbf{n}},{\mathbf{x}}_{\mathbf{n}+1}{]}^T$, $\mathbf{y}=[y_1,y_2,...,y_n{]}^T\to {\mathbf{y}}^{\prime }=[y_1,y_2,...,y_n,y_{n+1}{]}^T$, 即, 以一种数据流的形式给定样本; 这种情况下最优模型参数也将发生变化$\mathbf{w}\to {\mathbf{w}}^{\prime }$, 那么如果使用公式$4$ 就必须不断一次次计算逆矩阵$(X^{T}X{)}^{-1}$, 由于计算逆矩阵非常耗时, 上述的计算方法显然是不实用的, 因此希望找到一种以公式$4$为基础的递归求解新参数${\mathbf{w}}^{\prime }$的方法, 使得求解出的新模型${\mathbf{w}}^{\prime }$在当前最新的样本集$(X^{\prime },{\mathbf{y}}^{\prime })$上仍然满足误差平方和最小原则.

递归最小二乘具体解法

假设, 我们手头已经有了一个在已有样本$(X,\mathbf{y})$ 上满足最小误差平方和的模型参数$\mathbf{w}$ （至于最初的模型参数的获取见下文）, 我们希望找到一种递推公式, 能够得到更新数据前后的参数$\mathbf{w}\to {\mathbf{w}}^{\prime }$之间的关系, 避免一次次重新计算逆矩阵$(X^{T}X{)}^{-1}$, 就是RLS算法的主要动机.

对公式$4$ 进行分析, 定义$R\stackrel{\underset{def}{}}{=}{X}^{T}X,\mathbf{z}\stackrel{\underset{def}{}}{=}{X}^T\mathbf{y}$ ,则公式$4$可改写为
$$\mathbf{w}=R^{-1}\cdot \mathbf{z}\text{\hspace{1em}(5)}$$
在发生数据更新后, 新的权重矩阵记为${\mathbf{w}}^{\prime }$ , 新数据矩阵为$X^{\prime }$ 新矩阵$R^{\prime }$ 公式$4$可更新为
$${\mathbf{w}}^{\prime }=R^{\prime -1}\cdot {\mathbf{z}}^{\prime }\text{\hspace{1em}(6)}$$

递推求解矩阵$R^{\prime }$

在更新数据之后, 公式$4$求解新权重矩阵${\mathbf{w}}^{\prime }$的主要计算量在于求逆部分$R^{-1}$, 因此先对矩阵$R$进行计算处理, 根据分块矩阵计算,可以得到更新后矩阵$R^{\prime }$ 与更新前矩阵$R$ 之间的递推公式
$$R^{\prime }={X^{\prime }}^T{X}^{\prime }=[X^T|{\mathbf{x}}_{\mathbf{n}+1}]\left[\begin{array}{c}X\\ {{\mathbf{x}}_{\mathbf{n}+1}}^T\end{array}\right]=X^{T}X+{\mathbf{x}}_{\mathbf{n}+1}{{\mathbf{x}}_{\mathbf{n}+1}}^T=R+{\mathbf{x}}_{\mathbf{n}+1}{{\mathbf{x}}_{\mathbf{n}+1}}^T\text{\hspace{1em}(7)}$$
在现实中, 新的数据往往比旧数据更有价值, 因此一般为公式$7$添加遗忘因子$\lambda \le 1$, 这样越旧的数据在迭代过程中比重就越小, 即
$$R^{\prime }=\lambda R+{\mathbf{x}}_{\mathbf{n}+1}{{\mathbf{x}}_{\mathbf{n}+1}}^T\text{\hspace{1em}(8)}$$

递推求解逆矩阵$R^{\prime -1}$

公式$8$表明了矩阵$R$与$R^{\prime }$的迭代关系, 但是并不包含对求逆过程的处理, 我们更希望, 能够获得矩阵$R^{-1}$与${R^{\prime }}^{-1}$ 之间的递推关系. 在计算地推关系前, 需要引入如下引理

Theorem 1 : 如果矩阵$A$可以表示为如下形式
$$A=B^{-1}+C{D}^{-1}{C}^T\text{\hspace{1em}(9)}$$
​ 则逆矩阵$A^{-1}$可以表示
$$A^{-1}=B-BC(D+C^{T}BC{)}^{-1}{C}^{T}B\text{\hspace{1em}(10)}$$
将公式$9$,$10$相乘即可证明该引理

对比公式$8$,$9$令$A\stackrel{\underset{def}{}}{=}{R}^{\prime },B\stackrel{\underset{def}{}}{=}(\lambda R{)}^{-1},C\stackrel{\underset{def}{}}{=}{\mathbf{x}}_{\mathbf{n}+1},D\stackrel{\underset{def}{}}{=}1$ 则根据公式$10$计算得到${R^{\prime }}^{-1}$ 为
$$\begin{gathered} {R^{\prime }}^{-1}=(\lambda R{)}^{-1}-(\lambda R{)}^{-1}{\mathbf{x}}_{\mathbf{n}+1}(1+{{\mathbf{x}}_{\mathbf{n}+1}}^T(\lambda R{)}^{-1}{\mathbf{x}}_{\mathbf{n}+1}{)}^{-1}{{\mathbf{x}}_{\mathbf{n}+1}}^T(\lambda R{)}^{-1} \\ =\frac{1}{\lambda }{R}^{-1}-\frac{1}{\lambda }{R}^{-1}{\mathbf{x}}_{\mathbf{n}+1}\frac{1}{1+\frac{1}{\lambda }{{\mathbf{x}}_{\mathbf{n}+1}}^T{R}^{-1}{\mathbf{x}}_{\mathbf{n}+1}}{{\mathbf{x}}_{\mathbf{n}+1}}^T\frac{1}{\lambda }{R}^{-1} \\ =\frac{1}{\lambda }{R}^{-1}-\frac{\frac{1}{{\lambda }^2}{R}^{-1}{\mathbf{x}}_{\mathbf{n}+1}{{\mathbf{x}}_{\mathbf{n}+1}}^T{R}^{-1}}{1+\frac{1}{\lambda }{{\mathbf{x}}_{\mathbf{n}+1}}^T{R}^{-1}{\mathbf{x}}_{\mathbf{n}+1}} \\ =\frac{1}{\lambda }{R}^{-1}-\frac{\frac{1}{\lambda }{R}^{-1}{\mathbf{x}}_{\mathbf{n}+1}{{\mathbf{x}}_{\mathbf{n}+1}}^T{R}^{-1}}{\lambda +{{\mathbf{x}}_{\mathbf{n}+1}}^T{R}^{-1}{\mathbf{x}}_{\mathbf{n}+1}} \\ =\frac{1}{\lambda }{R}^{-1}-\frac{1}{\lambda }\frac{R^{-1}{\mathbf{x}}_{\mathbf{n}+1}}{\lambda +{{\mathbf{x}}_{\mathbf{n}+1}}^T{R}^{-1}{\mathbf{x}}_{\mathbf{n}+1}}{{\mathbf{x}}_{\mathbf{n}+1}}^T{R}^{-1}\text{\hspace{1em}(11)} \end{gathered}$$
公式$11$计算新的逆矩阵$R^{\prime -1}$的过程仅仅需要之前的旧的逆矩阵$R^{-1}$ 以及新添加的数据向量${\mathbf{x}}_{\mathbf{n}+1}$即可, 避免了直接求逆, 因此计算复杂度比直接求逆要小很多.

对公式$11$作进一步简化, 令$P^{\prime }\stackrel{\underset{def}{}}{=}{R}^{\prime -1},P\stackrel{\underset{def}{}}{=}{R}^{-1}$, 定义增益向量$k\stackrel{\underset{def}{}}{=}\frac{R^{-1}{\mathbf{x}}_{\mathbf{n}+1}}{\lambda +{{\mathbf{x}}_{\mathbf{n}+1}}^T{R}^{-1}{\mathbf{x}}_{\mathbf{n}+1}}$ 可转变为
$$P^{\prime }=\frac{1}{\lambda }P-\frac{1}{\lambda }k\cdot {{\mathbf{x}}_{\mathbf{n}+1}}^{T}P\text{\hspace{1em}(12)}$$
需要指出的是, 对公式$12$两侧都右乘向量${\mathbf{x}}_{n+1}$恰好满足如下关系
$$\begin{gathered} P^{\prime }{\mathbf{x}}_{n+1}=\frac{1}{\lambda }P{\mathbf{x}}_{n+1}-\frac{\frac{1}{\lambda }P{\mathbf{x}}_{\mathbf{n}+1}{{\mathbf{x}}_{\mathbf{n}+1}}^{T}P{\mathbf{x}}_{\mathbf{n}+1}}{\lambda +{{\mathbf{x}}_{\mathbf{n}+1}}^{T}P{\mathbf{x}}_{\mathbf{n}+1}} \\ =\frac{P{\mathbf{x}}_{\mathbf{n}+1}}{\lambda +{{\mathbf{x}}_{\mathbf{n}+1}}^{T}P{\mathbf{x}}_{\mathbf{n}+1}} \\ =k\text{\hspace{1em}(13)} \end{gathered}$$
这样, 根据公式$12$就得到了旧逆矩阵$P$与更新后逆矩阵$P^{\prime }$ 之间的递推关系; 重新表示公式$6$为
$${\mathbf{w}}^{\prime }=P^{\prime }\cdot {\mathbf{z}}^{\prime }\text{\hspace{1em}(14)}$$

递推求${\mathbf{z}}^{\prime }$

对公式$6$中的向量${\mathbf{z}}^{\prime }$同样利用分块矩阵计算
$${\mathbf{z}}^{\prime }={X^{\prime }}^T{\mathbf{y}}^{\prime }=[X^T|{\mathbf{x}}_{\mathbf{n}+1}]\left[\begin{array}{c}\mathbf{y}\\ y_{n+1}\end{array}\right]=X^T\mathbf{y}+{\mathbf{x}}_{\mathbf{n}+1}{y}_{n+1}=\mathbf{z}+{\mathbf{x}}_{\mathbf{n}+1}{y}_{n+1}\text{\hspace{1em}(15)}$$
添加遗忘因子$\lambda \le 1$ ,得到递推公式
$${\mathbf{z}}^{\prime }=\lambda \mathbf{z}+{\mathbf{x}}_{\mathbf{n}+1}{y}_{n+1}\text{\hspace{1em}(16)}$$

递推求${\mathbf{w}}^{\prime }$

结合公式 $12$,$13$,$14$,$16$,进行多步推导可以得到
$$\begin{gathered} {\mathbf{w}}^{\prime }=P^{\prime }\cdot {\mathbf{z}}^{\prime } \\ =P^{\prime }[\lambda \mathbf{z}+{\mathbf{x}}_{\mathbf{n}+1}{y}_{n+1}] \\ =\lambda P^{\prime }\mathbf{z}+P^{\prime }{\mathbf{x}}_{\mathbf{n}+1}{y}_{n+1} \\ =\lambda [\frac{1}{\lambda }P-\frac{1}{\lambda }k\cdot {{\mathbf{x}}_{\mathbf{n}+1}}^{T}P]\mathbf{z}+P^{\prime }{\mathbf{x}}_{\mathbf{n}+1}{y}_{n+1} \\ =P\mathbf{z}-k\cdot {{\mathbf{x}}_{\mathbf{n}+1}}^{T}P\mathbf{z}+P^{\prime }{\mathbf{x}}_{\mathbf{n}+1}{y}_{n+1} \\ =P\mathbf{z}-k\cdot {{\mathbf{x}}_{\mathbf{n}+1}}^{T}P\mathbf{z}+k\cdot y_{n+1} \\ =\mathbf{w}-k({{\mathbf{x}}_{\mathbf{n}+1}}^T\mathbf{w}-y_{n+1})\text{\hspace{1em}(17)} \end{gathered}$$
注意其中, ${{\mathbf{x}}_{\mathbf{n}+1}}^T\mathbf{w}-y_{n+1}$ 项中, 模型参数$\mathbf{w}$是旧模型参数,如果定义$e\stackrel{\underset{def}{}}{=}{{\mathbf{x}}_{\mathbf{n}+1}}^T\mathbf{w}-y_{n+1}$ 则公式$17$可变形为

$${\mathbf{w}}^{\prime }=\mathbf{w}-k\cdot e\text{\hspace{1em}(18)}$$
这就是RLS的最终计算目标.

关于初始化

RLS主要描述的是一种推理关系, 不断地在原来的旧最优模型参数上进行迭代得到最新模型参数; 那么最初进行迭代时, 需要一个初始的模型参数, 这个模型参数最好是满足最小平方和误差原则; 公式(5) 通过以上介绍, 可以改写为
$$\mathbf{w}=P\cdot \mathbf{z}\text{\hspace{1em}(19)}$$
其中, $\mathbf{z}\stackrel{\underset{def}{}}{=}{X}^T\mathbf{y}$ 可通过已有样本计算得出, 初始的$P$ 一般取
$$P=k\cdot I\text{\hspace{1em}(20)}$$
同时, 初始$k$取一个较大的数（保证$P$不会在递归过程中减小为负).

总结

RLS主要是在误差平方和最小的原则基础上, 提出一种解析的拟合模型参数$\mathbf{w}$的迭代递推公式; 可以实现在新的样本数据到来时, 利用新的样本数据以及旧的最优模型参数来便捷地计算新的满足最小二乘最优模型参数, 从而避免直接计算方法中的逆矩阵运算.

参考

[1] [知识梳理-04] Recursive Least Squares 递归最小二乘法 RLS_哔哩哔哩_bilibili

[2] 线性回归与递归最小二乘算法 (R.L.S algorithm) - 简书 (jianshu.com)

[3] 还有一个忘记了