感谢B站Up 凩子白的讲解视频, 大多数的RLS算法介绍都是从各种专业领域角度讲解的(比如滤波器等角度), 对于缺乏专业背景的同学入门较难, 本文主要是以上提到的视频的文字化, 加入了自己的一些理解, 也许有一些地方不是那么严谨, 不过希望能帮助其他同学快速了解一下RLS算法的思想。
PRELIMINARIES
对于样本数据对儿( x , y ) (\mathbf{x},y)(x ,y ), 其中输入数据向量x = [ x 11 , x 12 , . . . , x 1 m ] T ∈ R m \mathbf{x}=[x_{11},x_{12},…,x_{1m}]^T \in \mathbb{R}^m x =[x 11 ,x 12 ,…,x 1 m ]T ∈R m, 输出样本为y ∈ R y\in \mathbb{R}y ∈R; 使用参数为w \mathbf{w}w 的模型来拟合数据( x , y ) (\mathbf{x},y)(x ,y )之间的真实映射关系; 认为模型w \mathbf{w}w的输出为y y y的估计值y ^ ∈ R \hat{y}\in \mathbb{R}y ^∈R, 满足y ^ ∼ f ( w ; x ) \hat{y} \sim f({\mathbf{w}};\mathbf{x})y ^∼f (w ;x ), 拟合模型满足如下形式
y 1 ^ = w 1 x 11 + w 2 x 12 + . . . w m x 1 m = x 1 T w y 2 ^ = x 2 T w ⋮ y n ^ = x n T w (1) \hat{y_1}=w_1x_{11}+w_2x_{12}+…w_mx_{1m}=\mathbf{x_1^T}\mathbf{{w}}\ \hat{y_2}=\mathbf{x_2^T }\mathbf{{w}}\ \vdots\ \hat{y_n}=\mathbf{x_n^T }\mathbf{{w}} \tag{1}y 1 ^=w 1 x 11 +w 2 x 12 +…w m x 1 m =x 1 T w y 2 ^=x 2 T w ⋮y n ^=x n T w (1 )
最小二乘法的思路, 就是希望近似模型参数w \mathbf{{w}}w在这n n n个样本输入数据X n × m = [ x 1 , x 2 , . . . , x n ] T X_{n\times m}=[\mathbf{x_1},\mathbf{x_2},…,\mathbf{x_n}]^T X n ×m =[x 1 ,x 2 ,…,x n ]T (以后简记为X X X)上得出的估计值y ^ = [ y 1 ^ , y 2 ^ , . . . , y n ^ ] T \hat{\mathbf y}=[\hat{y_1},\hat{y_2},…,\hat{y_n}]^T y ^=[y 1 ^,y 2 ^,…,y n ^]T 与ground truth 输出样本数据y = [ y 1 , y 2 , . . . , y n ] T \mathbf y=[y_1,y_2,…,y_n]^T y =[y 1 ,y 2 ,…,y n ]T 之间的差值平方和最小,即
w = arg m i n ∑ i = 1 n ( y i − y i ^ ) 2 = arg m i n ∑ i = 1 n ( y i − x i T w ) 2 = arg m i n ∥ y − X w ∥ 2 2 = arg m i n E ( w ) (2) {\mathbf{w}} =\arg min \sum\limits_{i=1}^n (y_i-\hat{y_i})^2\ = \arg min \sum\limits_{i=1}^n (y_i-\mathbf{x_i^T }\mathbf{{w}})^2\ =\arg min \begin{Vmatrix} \mathbf{y}-X\mathbf{w} \end{Vmatrix}2^2\ =\arg min\ E(\mathbf{w}) \tag{2}w =ar g min i =1 ∑n (y i −y i ^)2 =ar g min i =1 ∑n (y i −x i T w )2 =ar g min y −X w 2 2 =ar g min E (w )(2 )
误差E ( w ) E(\mathbf{w})E (w )对参数w \mathbf{w}w 求梯度,
∇ w E = ∇ w ∥ y − X w ∥ 2 2 = ∇ w [ ( y − X w ) T ( y − X w ) ] = 2 X T ( y − X w ) (3) \nabla\mathbf{w} E =\nabla_\mathbf{w}\begin{Vmatrix} \mathbf{y}-X\mathbf{w} \end{Vmatrix}2^2\ =\nabla\mathbf{w} \big[(\mathbf{y}-X\mathbf{w})^T(\mathbf{y}-X\mathbf{w})\big] \ =2X^T(\mathbf{y}-X\mathbf{w}) \tag{3}∇w E =∇w y −X w 2 2 =∇w [(y −X w )T (y −X w )]=2 X T (y −X w )(3 )
令∇ w E = 0 \nabla_\mathbf{w} E=0 ∇w E =0, 即可求出
w = X − 1 y = ( X T X ) − 1 X T y (4) \mathbf{w}=X^{-1}\mathbf{y}=(X^TX)^{-1}X^T\mathbf{y} \tag{4}w =X −1 y =(X T X )−1 X T y (4 )
注意, 公式3 {3}3中的y \mathbf y y是样本数据, 将X − 1 X^{-1}X −1表述为( X T X ) − 1 X T (X^TX)^{-1}X^T (X T X )−1 X T的原因是矩阵X X X不一定是n × n n\times n n ×n形状,因此不一定有逆矩阵, 而X T X X^TX X T X的逆是存在的?
当一次性给出所有样本集合( X , y ) (X,\mathbf{y})(X ,y )时, 可以通过公式4 {4}4来直接计算出最优的拟合模型参数w \mathbf{w}w , 然而, 在实际应用中, 这种直接计算法并不常见, 主要是因为公式中求逆部分( X T X ) − 1 (X^TX)^{-1}(X T X )−1 的计算量大, 在样本数据量大时计算量更是明显增大; 另外, 现实生活中,往往出现样本数据可能也并不是一次性给出, 而是不断给出新的样本数据, 以一种数据流的形式给出样本数据, 例如传感器随时间不断读取信号等, 这种情况下利用公式4 {4}4直接计算最优模型参数w {\mathbf{w}}w 就需要每次进行直接计算, 也是不现实的.
因此, 为了利用最小误差平方和原则, 求解在大样本量, 或者数据流情况下的最优模型参数w {\mathbf{w}}w, 一种方法可以将大样本分成多批次(batch), 计算旧模型在新批次样本上的梯度, 不断进行梯度下降来进行迭代求解(也可以将数据流当做一个个batch来梯度更新); 另一种则是解析的方法, 就是这里提到的递归最小二乘解法(RLS).
本质上, 递归最小二乘法RLS和梯度下降、直接计算法一样, 都是为了求解满足最小误差平方和原则的最优模型参数w {\mathbf{w}}w, 只是在实现方式上有所不同.
如之前提到的, RLS的主要应用场景, 是假设输入样本数据X X X在不断添加新数据, 例如X = [ x 1 , x 2 , . . . , x n ] T → X ′ = [ x 1 , x 2 , . . . , x n , x n + 1 ] T X=[\mathbf{x_1},\mathbf{x_2},…,\mathbf{x_n}]^T\rightarrow X’=[\mathbf{x_1},\mathbf{x_2},…,\mathbf{x_n} , \mathbf{x_{n+1}}]^T X =[x 1 ,x 2 ,…,x n ]T →X ′=[x 1 ,x 2 ,…,x n ,x n +1 ]T, y = [ y 1 , y 2 , . . . , y n ] T → y ′ = [ y 1 , y 2 , . . . , y n , y n + 1 ] T \mathbf y=[y_1,y_2,…,y_n]^T\rightarrow\mathbf y’=[y_1,y_2,…,y_n,y_{n+1}]^T y =[y 1 ,y 2 ,…,y n ]T →y ′=[y 1 ,y 2 ,…,y n ,y n +1 ]T, 即, 以一种数据流的形式给定样本; 这种情况下最优模型参数也将发生变化w → w ′ \mathbf{w}\rightarrow \mathbf{w}’w →w ′, 那么如果使用公式4 {4}4 就必须不断一次次计算逆矩阵( X T X ) − 1 (X^TX)^{-1}(X T X )−1, 由于计算逆矩阵非常耗时, 上述的计算方法显然是不实用的, 因此希望找到一种以公式4 {4}4为基础的递归求解新参数w ′ \mathbf{w’}w ′的方法, 使得求解出的新模型w ′ \mathbf{w’}w ′在当前最新的样本集( X ′ , y ′ ) (X’,\mathbf{y’})(X ′,y ′)上仍然满足误差平方和最小原则.
递归最小二乘具体解法
假设, 我们手头已经有了一个在已有样本( X , y ) (X,\mathbf{y})(X ,y ) 上满足最小误差平方和的模型参数w \mathbf{w}w （至于最初的模型参数的获取见下文）, 我们希望找到一种递推公式, 能够得到更新数据前后的参数w → w ′ \mathbf{w}\rightarrow \mathbf{w’}w →w ′之间的关系, 避免一次次重新计算逆矩阵( X T X ) − 1 (X^TX)^{-1}(X T X )−1, 就是RLS算法的主要动机.
对公式4 {4}4 进行分析, 定义R = d e f X T X , z = d e f X T y R\overset{\underset{def}{}}{=} X^TX, \mathbf{z} \overset{\underset{def}{}}{=} X^T\mathbf{{y}}R =d e f X T X ,z =d e f X T y ,则公式4 {4}4可改写为
w = R − 1 ⋅ z (5) \mathbf{w}=R^{-1}\cdot \mathbf{z} \tag{5}w =R −1 ⋅z (5 )
在发生数据更新后, 新的权重矩阵记为w ′ \mathbf{w}’w ′ , 新数据矩阵为X ′ X’X ′ 新矩阵R ′ R’R ′ 公式4 {4}4可更新为
w ′ = R ′ − 1 ⋅ z ′ (6) \mathbf{w’}=R’^{-1}\cdot \mathbf{z’} \tag{6}w ′=R ′−1 ⋅z ′(6 )
在更新数据之后, 公式4 {4}4求解新权重矩阵w ′ \mathbf{w}’w ′的主要计算量在于求逆部分R − 1 R^{-1}R −1, 因此先对矩阵R R R进行计算处理, 根据分块矩阵计算,可以得到更新后矩阵R ′ R’R ′ 与更新前矩阵R R R 之间的递推公式
R ′ = X ′ T X ′ = [ X T ∣ x n + 1 ] [ X x n + 1 T ] = X T X + x n + 1 x n + 1 T = R + x n + 1 x n + 1 T (7) R’ = {X’}^TX’ = [X^T|\mathbf{x_{n+1}}] \begin{bmatrix} X\ \hline {\mathbf{x_{n+1}}}^T \end{bmatrix} =X^TX+\mathbf{x_{n+1}}\mathbf{x_{n+1}}^T = R + \mathbf{x_{n+1}}\mathbf{x_{n+1}}^T \tag{7}R ′=X ′T X ′=[X T ∣x n +1 ][X x n +1 T ]=X T X +x n +1 x n +1 T =R +x n +1 x n +1 T (7 )
在现实中, 新的数据往往比旧数据更有价值, 因此一般为公式7 {7}7添加遗忘因子λ ≤ 1 \lambda \leq 1 λ≤1, 这样越旧的数据在迭代过程中比重就越小, 即
R ′ = λ R + x n + 1 x n + 1 T (8) R’= \lambda R+\mathbf{x_{n+1}}\mathbf{x_{n+1}}^T \tag{8}R ′=λR +x n +1 x n +1 T (8 )
公式8 {8}8表明了矩阵R R R与R ′ R’R ′的迭代关系, 但是并不包含对求逆过程的处理, 我们更希望, 能够获得矩阵R − 1 R^{-1}R −1与R ′ − 1 {R’}^{-1}R ′−1 之间的递推关系. 在计算地推关系前, 需要引入如下引理
Theorem 1 : 如果矩阵A A A可以表示为如下形式
A = B − 1 + C D − 1 C T (9) A= B^{-1} + CD^{-1}C^T \tag{9}A =B −1 +C D −1 C T (9 )
则逆矩阵A − 1 A^{-1}A −1可以表示
A − 1 = B − B C ( D + C T B C ) − 1 C T B (10) A^{-1}=B-BC(D+C^TBC)^{-1}C^TB \tag{10}A −1 =B −BC (D +C T BC )−1 C T B (10 )
将公式9 {9}9,10 {10}10相乘即可证明该引理
对比公式8 {8}8,9 {9}9令A = d e f R ′ , B = d e f ( λ R ) − 1 , C = d e f x n + 1 , D = d e f 1 A\overset{\underset{def}{}}{=} R’, B\overset{\underset{def}{}}{=}(\lambda R)^{-1}, C\overset{\underset{def}{}}{=} \mathbf{x_{n+1}},D \overset{\underset{def}{}}{=} 1 A =d e f R ′,B =d e f (λR )−1 ,C =d e f x n +1 ,D =d e f 1 则根据公式10 {10}10计算得到R ′ − 1 {R’}^{-1}R ′−1 为
R ′ − 1 = ( λ R ) − 1 − ( λ R ) − 1 x n + 1 ( 1 + x n + 1 T ( λ R ) − 1 x n + 1 ) − 1 x n + 1 T ( λ R ) − 1 = 1 λ R − 1 − 1 λ R − 1 x n + 1 1 1 + 1 λ x n + 1 T R − 1 x n + 1 x n + 1 T 1 λ R − 1 = 1 λ R − 1 − 1 λ 2 R − 1 x n + 1 x n + 1 T R − 1 1 + 1 λ x n + 1 T R − 1 x n + 1 = 1 λ R − 1 − 1 λ R − 1 x n + 1 x n + 1 T R − 1 λ + x n + 1 T R − 1 x n + 1 = 1 λ R − 1 − 1 λ R − 1 x n + 1 λ + x n + 1 T R − 1 x n + 1 x n + 1 T R − 1 (11) {R’}^{-1} =(\lambda R)^{-1} – (\lambda R)^{-1}\mathbf{x_{n+1}}\big(1+\mathbf{x_{n+1}}^T(\lambda R)^{-1}\mathbf{x_{n+1}}\big)^{-1}\mathbf{x_{n+1}}^T(\lambda R)^{-1}\ = \frac{1}{\lambda}R^{-1} – \frac{1}{\lambda}R^{-1}\mathbf{x_{n+1}}\frac{1}{1+\frac{1}{\lambda}\mathbf{x_{n+1}}^TR^{-1}\mathbf{x_{n+1}}}\mathbf{x_{n+1}}^T\frac{1}{\lambda}R^{-1}\ =\frac{1}{\lambda}R^{-1} – \frac{\frac{1}{\lambda^2}R^{-1}\mathbf{x_{n+1}}\mathbf{x_{n+1}}^TR^{-1}}{1+\frac{1}{\lambda}\mathbf{x_{n+1}}^TR^{-1}\mathbf{x_{n+1}}}\ =\frac{1}{\lambda}R^{-1} – \frac{\frac{1}{\lambda}R^{-1}\mathbf{x_{n+1}}\mathbf{x_{n+1}}^TR^{-1}}{\lambda+\mathbf{x_{n+1}}^TR^{-1}\mathbf{x_{n+1}}}\ =\frac{1}{\lambda}R^{-1} – \frac{1}{\lambda}\frac{R^{-1}\mathbf{x_{n+1}}}{\lambda+\mathbf{x_{n+1}}^TR^{-1}\mathbf{x_{n+1}}}\mathbf{x_{n+1}}^TR^{-1} \tag{11}R ′−1 =(λR )−1 −(λR )−1 x n +1 (1 +x n +1 T (λR )−1 x n +1 )−1 x n +1 T (λR )−1 =λ1 R −1 −λ1 R −1 x n +1 1 +λ1 x n +1 T R −1 x n +1 1 x n +1 T λ1 R −1 =λ1 R −1 −1 +λ1 x n +1 T R −1 x n +1 λ2 1 R −1 x n +1 x n +1 T R −1 =λ1 R −1 −λ+x n +1 T R −1 x n +1 λ1 R −1 x n +1 x n +1 T R −1 =λ1 R −1 −λ1 λ+x n +1 T R −1 x n +1 R −1 x n +1 x n +1 T R −1 (11 )
公式11 {11}11计算新的逆矩阵R ′ − 1 R’^{-1}R ′−1的过程仅仅需要之前的旧的逆矩阵R − 1 R^{-1}R −1 以及新添加的数据向量x n + 1 \mathbf{x_{n+1}}x n +1 即可, 避免了直接求逆, 因此计算复杂度比直接求逆要小很多.
对公式11 {11}11作进一步简化, 令P ′ = d e f R ′ − 1 , P = d e f R − 1 P’\overset{\underset{def}{}}{=} R’^{-1},P\overset{\underset{def}{}}{=} R^{-1}P ′=d e f R ′−1 ,P =d e f R −1, 定义增益向量k = d e f R − 1 x n + 1 λ + x n + 1 T R − 1 x n + 1 k\overset{\underset{def}{}}{=} \frac{R^{-1}\mathbf{x_{n+1}}}{\lambda+\mathbf{x_{n+1}}^TR^{-1}\mathbf{x_{n+1}}}k =d e f λ+x n +1 T R −1 x n +1 R −1 x n +1 可转变为
P ′ = 1 λ P − 1 λ k ⋅ x n + 1 T P (12) P’ = \frac{1}{\lambda}P – \frac{1}{\lambda}k\cdot \mathbf{x_{n+1}}^TP \tag{12}P ′=λ1 P −λ1 k ⋅x n +1 T P (12 )
需要指出的是, 对公式12 {12}12两侧都右乘向量x n + 1 \mathbf{x}{n+1}x n +1 恰好满足如下关系
P ′ x n + 1 = 1 λ P x n + 1 − 1 λ P x n + 1 x n + 1 T P x n + 1 λ + x n + 1 T P x n + 1 = P x n + 1 λ + x n + 1 T P x n + 1 = k (13) P’\mathbf{x}{n+1} = \frac{1}{\lambda}P\mathbf{x}{n+1} -\frac{\frac{1}{\lambda} P\mathbf{x{n+1}}\mathbf{x_{n+1}}^TP\mathbf{x_{n+1}}}{\lambda+\mathbf{x_{n+1}}^T P\mathbf{x_{n+1}}}\ =\frac{P\mathbf{x_{n+1}}}{\lambda+\mathbf{x_{n+1}}^T P\mathbf{x_{n+1}}}\ =k \tag{13}P ′x n +1 =λ1 P x n +1 −λ+x n +1 T P x n +1 λ1 P x n +1 x n +1 T P x n +1 =λ+x n +1 T P x n +1 P x n +1 =k (13 )
这样, 根据公式12 {12}12就得到了旧逆矩阵P P P与更新后逆矩阵P ′ P’P ′ 之间的递推关系; 重新表示公式6 {6}6为
w ′ = P ′ ⋅ z ′ (14) \mathbf{w’}=P’\cdot \mathbf{z’} \tag{14}w ′=P ′⋅z ′(14 )
对公式6 {6}6中的向量z ′ \mathbf{z’}z ′同样利用分块矩阵计算
z ′ = X ′ T y ′ = [ X T ∣ x n + 1 ] [ y y n + 1 ] = X T y + x n + 1 y n + 1 = z + x n + 1 y n + 1 (15) \mathbf{z’} = {X’}^T\mathbf{{y’}} = [X^T|\mathbf{x_{n+1}}] \begin{bmatrix} \mathbf{{y}} \ \hline {y}{n+1} \end{bmatrix} =X^T\mathbf{{y}}+\mathbf{x{n+1}}{y}{n+1} = \mathbf{z}+\mathbf{x{n+1}}{y}{n+1} \tag{15}z ′=X ′T y ′=[X T ∣x n +1 ][y y n +1 ]=X T y +x n +1 y n +1 =z +x n +1 y n +1 (15 )
添加遗忘因子λ ≤ 1 \lambda\leq 1 λ≤1 ,得到递推公式
z ′ = λ z + x n + 1 y n + 1 (16) \mathbf{z’} =\lambda\mathbf{z}+\mathbf{x{n+1}}{y}_{n+1} \tag{16}z ′=λz +x n +1 y n +1 (16 )
结合公式 12 {12}12,13 {13}13,14 {14}14,16 {16}16,进行多步推导可以得到
w ′ = P ′ ⋅ z ′ = P ′ [ λ z + x n + 1 y n + 1 ] = λ P ′ z + P ′ x n + 1 y n + 1 = λ [ 1 λ P − 1 λ k ⋅ x n + 1 T P ] z + P ′ x n + 1 y n + 1 = P z − k ⋅ x n + 1 T P z + P ′ x n + 1 y n + 1 = P z − k ⋅ x n + 1 T P z + k ⋅ y n + 1 = w − k ( x n + 1 T w − y n + 1 ) (17) \mathbf{w’} =P’\cdot \mathbf{z’}\ =P'[\lambda\mathbf{z}+\mathbf{x_{n+1}}{y}{n+1}] \ =\lambda P’\mathbf{z}+P’\mathbf{x{n+1}}{y}{n+1}\ =\lambda \bigg[\frac{1}{\lambda}P – \frac{1}{\lambda}k\cdot\mathbf{x{n+1}}^T P \bigg]\mathbf{z}+ P’\mathbf{x_{n+1}}{y}{n+1} \ =P\mathbf{z} – k\cdot \mathbf{x{n+1}}^T P\mathbf{z} + P’\mathbf{x_{n+1}}{y}{n+1} \ =P\mathbf{z} – k\cdot\mathbf{x{n+1}}^T P\mathbf{z} + k\cdot {y}{n+1} \ =\mathbf{w}-k(\mathbf{x{n+1}}^T\mathbf{w}-{y}{n+1}) \tag{17}w ′=P ′⋅z ′=P ′[λz +x n +1 y n +1 ]=λP ′z +P ′x n +1 y n +1 =λ[λ1 P −λ1 k ⋅x n +1 T P ]z +P ′x n +1 y n +1 =P z −k ⋅x n +1 T P z +P ′x n +1 y n +1 =P z −k ⋅x n +1 T P z +k ⋅y n +1 =w −k (x n +1 T w −y n +1 )(17 )
注意其中, x n + 1 T w − y n + 1 \mathbf{x{n+1}}^T\mathbf{w}-{y}{n+1}x n +1 T w −y n +1 项中, 模型参数w \mathbf{w}w是旧模型参数,如果定义e = d e f x n + 1 T w − y n + 1 e\overset{\underset{def}{}}{=}\mathbf{x{n+1}}^T\mathbf{w}-{y}_{n+1}e =d e f x n +1 T w −y n +1 则公式17 {17}17可变形为
w ′ = w − k ⋅ e (18) \mathbf{w’}=\mathbf{w}-k\cdot e \tag{18}w ′=w −k ⋅e (18 )
这就是RLS的最终计算目标.
RLS主要描述的是一种推理关系, 不断地在原来的旧最优模型参数上进行迭代得到最新模型参数; 那么最初进行迭代时, 需要一个初始的模型参数, 这个模型参数最好是满足最小平方和误差原则; 公式(5) 通过以上介绍, 可以改写为
w = P ⋅ z (19) \mathbf{w} = P\cdot \mathbf{z} \tag{19}w =P ⋅z (19 )
其中, z = d e f X T y \mathbf{z} \overset{\underset{def}{}}{=} X^T\mathbf{{y}}z =d e f X T y 可通过已有样本计算得出, 初始的P P P 一般取
P = k ⋅ I (20) P=k\cdot I \tag{20}P =k ⋅I (20 )
同时, 初始k k k取一个较大的数（保证P P P不会在递归过程中减小为负).
RLS主要是在误差平方和最小的原则基础上, 提出一种解析的拟合模型参数w \mathbf{w}w的迭代递推公式; 可以实现在新的样本数据到来时, 利用新的样本数据以及旧的最优模型参数来便捷地计算新的满足最小二乘最优模型参数, 从而避免直接计算方法中的逆矩阵运算.
[3] 还有一个忘记了
Original: https://blog.csdn.net/u012057432/article/details/125836404
Author: hymwgk
Title: RLS递归最小二乘法(Recursive Least Squares)
原创文章受到原创版权保护。转载请注明出处：https://www.johngo689.com/614304/
转载文章受原作者版权保护。转载请注明原作者出处！