带遗忘因子的递推最小二乘法推导
Posted 找不到服务器zhn
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了带遗忘因子的递推最小二乘法推导相关的知识，希望对你有一定的参考价值。
最小二乘法的递推形式、直流信号的遗忘递推形式、遗忘递推最小二乘。
摘要：最小二乘法的递推形式、直流信号的遗忘递推形式、遗忘递推最小二乘。
递推最小二乘法
对多组数据 \\(\\vecx_i\\) 和 \\(y_i\\)，满足
\\[y_i = \\vecx^\\mathrmT_i\\vec\\theta
\\]
其中 \\(\\vecx_i\\) 是输入数据向量，\\(y_i\\) 是输出数据标量。写成矩阵形式
\\[\\vecy = X\\vec\\theta
\\]
其中（以3组数据、2个未知参数为例）
\\[\\vecy = \\beginbmatrix
y_1 \\\\ y_2 \\\\ y_3
\\endbmatrix,
X = \\beginbmatrix
\\vecx^\\mathrmT_1 \\\\
\\vecx^\\mathrmT_2 \\\\
\\vecx^\\mathrmT_3
\\endbmatrix,
\\theta = \\beginbmatrix
\\theta_1 \\\\ \\theta_2
\\endbmatrix,
y_i=x_i1\\theta_1+x_i2\\theta_2
\\]
完整地写作（以3组数据、2个未知参数为例）
\\[\\beginbmatrix
y_1 \\\\ y_2 \\\\ y_3
\\endbmatrix
=\\beginbmatrix
x_11 & x_12 \\\\
x_21 & x_22 \\\\
x_31 & x_32 \\\\
\\endbmatrix
\\beginbmatrix
\\theta_1 \\\\ \\theta_2
\\endbmatrix
\\]
有\\(k\\)组数据时记作
\\[ X^\\mathrmT_k = \\beginbmatrix
\\vecx_1 & \\vecx_2 & \\cdots & \\vecx_k
\\endbmatrix \\\\
\\vecy_k=\\beginbmatrix y_1 & y_2 & \\cdots & y_k\\endbmatrix
\\]
(此时注意区分一些符号，例如 \\(\\vecy_k\\) 与 \\(y_k\\)，\\(X\\) 与 \\(\\vecx\\) 等)
已知最小二乘法的解为
\\[ \\vec\\theta
=(X^\\mathrmTX)^-1
X^\\mathrmT\\vecy \\tag1
\\]
令 \\(P = (X^\\mathrmTX)^-1\\)，
则加上递推后
\\[\\beginaligned
P_k^-1 =& X_k^\\mathrmTX_k
= \\sum_i=1^k\\vecx_i\\vecx^\\mathrmT_i\\\\
=& \\sum_i=1^k-1\\vecx_i\\vecx^\\mathrmT_i
+\\vecx_k\\vecx^\\mathrmT_k \\\\
=& P_k-1^-1 + \\vecx_k\\vecx_k^\\mathrmT
\\endaligned \\tag2\\]
同理可得
\\[ X_k^\\mathrmT\\vecy_k
=X_k-1^\\mathrmT\\vecy_k-1
+\\vecx_ky_k \\tag3
\\]
于是代入式(1)得到
\\[\\beginaligned
\\vec\\theta_k =& P_kX_k^\\mathrmT\\vecy_k \\\\
=& P_k(X_k-1^\\mathrmT\\vecy_k-1
+\\vecx_ky_k) \\\\
=& P_k(P_k-1^-1\\vec\\theta_k-1
+\\vecx_ky_k) \\\\
=& P_k(P_k^-1\\vec\\theta_k-1
-\\vecx_k\\vecx_k^\\mathrmT\\theta_k-1+\\vecx_ky_k) \\\\
=& \\vec\\theta_k-1 + P_k\\vecx_k
(y_k-\\vecx_k^\\mathrmT\\theta_k-1) \\\\
=& \\vec\\theta_k-1
+(P_k-1^-1 + \\vecx_k\\vecx_k^\\mathrmT)^-1
\\vecx_k(y_k-\\vecx_k^\\mathrmT\\theta_k-1) \\\\
=& \\vec\\theta_k-1
+(P_k-1-\\fracP_k-1\\vecx_k\\vecx_k^\\mathrmT
P_k-11+\\vecx_k^\\mathrmTP_k-1\\vecx_k)
\\vecx_k(y_k-\\vecx_k^\\mathrmT\\theta_k-1) \\\\
\\endaligned\\]
其中倒数第3行到倒数第一行的过程
\\[ P_k = P_k-1-\\fracP_k-1\\vecx_k\\vecx_k^\\mathrmT
P_k-11+\\vecx_k^\\mathrmTP_k-1\\vecx_k
\\]
用到了下面的矩阵求逆公式及其引理
\\[\\beginaligned
(A+BCD)^-1 =& A^-1-A^-1 B
(DA^-1B+C^-1)^-1DA^-1 \\\\
(A+\\vecu\\vecu^\\textT)^-1 =& A^-1
-\\fracA^-1\\vecu\\vecu^\\textT A^-1
1+\\vecu^\\textTA^-1\\vecu
\\endaligned\\]
令
\\[ K_k = \\fracP_k-1\\vecx_k
1+\\vecx_k^\\mathrmTP_k-1\\vecx_k
\\]
则
\\[\\beginaligned
P_k =& (I-K_k\\vecx_k^\\mathrmT)P_k-1\\\\
P_k\\vecx_k =& P_k-1\\vecx_k
-K_k\\vecx_k^\\mathrmTP_k-1\\vecx_k \\\\
=& \\fracP_k-1\\vecx_k
(1+\\vecx_k^\\mathrmTP_k-1\\vecx_k)
-P_k-1\\vecx_k\\vecx_k^\\mathrmT
P_k-1\\vecx_k
1+\\vecx_k^\\mathrmTP_k-1\\vecx_k \\\\
=& \\fracP_k-1\\vecx_k
1+\\vecx_k^\\mathrmTP_k-1\\vecx_k \\\\
=& K_k \\\\
\\endaligned\\]
总结成递推公式得到
\\[\\beginaligned
K_k =& \\fracP_k-1\\vecx_k1+\\vecx_k^\\mathrmT
P_k-1\\vecx_k\\\\
P_k =& (I-K_k\\vecx_k^\\mathrmT)P_k-1 \\\\
\\vec\\theta_k =& \\vec\\theta_k-1
+K_k(y_k-\\vecx_k^\\mathrmT\\theta_k-1)
\\endaligned\\]
因为 \\(P=(X^\\textTX)^-1\\)，\\(\\vecx_1\\vecx_1^\\textT\\) 一般很小，所以实际使用时一般把 \\(P_1\\) 初始化为一个接近无穷大的单位阵。
引入遗忘因子
先考虑一个比较简单的情况
\\[x_n=\\theta+w_n
\\]
其中 \\(\\theta\\) 是要估计的参数，\\(w_n\\) 是高斯白噪声，\\(x_n\\) 是观测到的数据，总共有 \\(N\\) 组数据，将 \\(N\\) 组数据取平均得到 \\(\\theta\\) 的一个估计值为
\\[\\hat\\theta=\\frac1N\\sum_n=1^Nx_n
\\]
写成递推形式为
\\[\\hat\\theta_n=\\hat\\theta_n-1+\\frac1n(x_n-\\hat\\theta_n-1)
\\]
例如，
\\[\\beginaligned
\\hat\\theta_1 =& x_1\\\\
\\hat\\theta_2 =& \\hat\\theta_1+\\frac12(x_2-\\hat\\theta_1)
=\\frac12(x_1+x_2) \\\\
\\hat\\theta_3 =& \\hat\\theta_2+\\frac13(x_3-\\hat\\theta_2)
=\\frac13(x_1+x_2+x_3) \\\\
\\endaligned\\]
现在对之前的值乘一个衰减系数 \\(\\lambda\\in(0,1)\\)，也就是
\\[\\beginaligned
\\hat\\theta_2 =& \\frac1\\lambda+1(\\lambda x_1+x_2) \\\\
\\hat\\theta_3 =& \\frac1\\lambda^2+\\lambda+1(\\lambda^2x_1+\\lambda x_2+x_3) \\\\
\\vdots
\\endaligned\\]
写成递推形式为
\\[\\beginaligned
\\hat\\theta_2=&\\frac1\\lambda+1(\\lambda x_1+x_2) \\\\
\\hat\\theta_3=&\\frac1\\lambda^2+\\lambda+1(\\lambda(\\lambda+1)\\hat\\theta_2+x_3) \\\\
=&\\hat\\theta_2+\\frac1\\lambda^2+\\lambda+1(x_3-\\hat\\theta_2)
\\endaligned\\]
当 \\(n\\rightarrow\\infty\\) 时，递推形式可以近似为
\\[\\beginaligned
\\hat\\theta_n =& \\hat\\theta_n-1+\\frac1\\sum_n=1^N\\lambda^n(x_n-\\hat\\theta_n-1) \\\\
=& \\hat\\theta_n-1+(1-\\lambda)(x_n-\\hat\\theta_n-1)
\\endaligned\\]
递推最小二乘中的遗忘因子
类似地，在最小二乘的误差 \\(S_i=||y_i-\\vecx_i^\\textT\\vec\\theta||^2\\) 中加入遗忘因子，得到总误差（以3组数据为例）
\\[J_3=\\lambda^2S_1+\\lambda S_2+S_3
\\]
完整的方程可以写作
\\[\\vecy_3=X_3\\theta \\\\
\\beginbmatrix
\\lambda y_1 \\\\ \\sqrt\\lambday_2 \\\\ y_3
\\endbmatrix
=\\beginbmatrix
\\lambda x_11 & \\lambda x_12 \\\\
\\sqrt\\lambdax_21 & \\sqrt\\lambdax_22 \\\\
x_31 & x_32 \\\\
\\endbmatrix
\\beginbmatrix
\\theta_1 \\\\ \\theta_2
\\endbmatrix
\\]
类似地代入最小二乘的解式(1)，并求式(2)(3)的递推形式
\\[\\beginaligned
X_2^\\textTX_2 =& \\beginbmatrix
\\sqrt\\lambdax_11 & x_21 \\\\
\\sqrt\\lambdax_12 & x_22 \\\\
\\endbmatrix
\\beginbmatrix
\\sqrt\\lambdax_11 & \\sqrt\\lambda x_12 \\\\
x_21 & x_22
\\endbmatrix \\\\
X_3^\\textTX_3 =& \\beginbmatrix
\\lambda x_11 & \\sqrt\\lambdax_21 & x_31 \\\\
\\lambda x_12 & \\sqrt\\lambdax_22 & x_32 \\\\
\\endbmatrix
\\beginbmatrix
\\lambda x_11 & \\lambda x_12 \\\\
\\sqrt\\lambdax_21 & \\sqrt\\lambda x_22 \\\\
x_31 & x_32
\\endbmatrix \\\\
=& \\lambda X_2^\\textTX_2 + \\vecx_3\\vecx_3^\\textT \\\\
X_2^\\textT\\vecy_2 =& \\beginbmatrix
\\sqrt\\lambdax_11 & x_21 \\\\
\\sqrt\\lambdax_12 & x_22 \\\\
\\endbmatrix
\\beginbmatrix
\\sqrt\\lambday_1 \\\\ y_2
\\endbmatrix \\\\
X_3^\\textT\\vecy_3 =& \\beginbmatrix
\\lambda x_11 & \\sqrt\\lambdax_21 & x_31 \\\\
\\lambda x_12 & \\sqrt\\lambdax_22 & x_32 \\\\
\\endbmatrix
\\beginbmatrix
\\lambda y_1 \\\\ \\sqrt\\lambday_2 \\\\ y_3
\\endbmatrix \\\\
=& \\lambda X_2^\\textT\\vecy_2 + \\vecx_3y_3
\\endaligned\\]
即
\\[\\beginaligned
P_k^-1 =& \\lambda P_k-1^-1 + \\vecx_k\\vecx_k^\\textT \\\\
X_k^\\textT\\vecy_k =& \\lambda X_k-1^\\textT\\vecy_k-1 + \\vecx_ky_k
\\endaligned\\]
代入式(1)得到
\\[\\beginaligned
\\vec\\theta_k =& P_kX_k^\\textT\\vecy_k \\\\
=& P_k(\\lambda X_k-1^\\textT\\vecy_k-1 + \\vecx_ky_k) \\\\
=& P_k(\\lambda P_k-1^-1\\vec\\theta_k-1 + \\vecx_ky_k) \\\\
=& P_k(P_k^-1\\vec\\theta_k-1
-\\vecx_k\\vecx_k^\\textT\\theta_k-1+\\vecx_ky_k) \\\\
=& \\vec\\theta_k-1 + P_k\\vecx_k(y_k-\\vecx_k^\\textT\\theta_k-1) \\\\
=& \\vec\\theta_k-1 + (\\lambda P_k-1^-1 + \\vecx_k\\vecx_k^\\textT)^-1
\\vecx_k(y_k-\\vecx_k^\\textT\\theta_k-1) \\\\
=& \\vec\\theta_k-1+(\\fracP_k-1\\lambda
-\\frac\\fracP_k-1\\lambda\\vecx_k\\vecx_k^\\textT
\\fracP_k-1\\lambda1+\\vecx_k^\\textT\\fracP_k-1\\lambda\\vecx_k)
\\vecx_k(y_k-\\vecx_k^\\textT\\theta_k-1) \\\\
=& \\vec\\theta_k-1+(\\fracP_k-1\\lambda
-\\frac1\\lambda\\fracP_k-1\\vecx_k\\vecx_k^\\textTP_k-1
\\lambda + \\vecx_k^\\textTP_k-1\\vecx_k)
\\vecx_k(y_k-\\vecx_k^\\textT\\theta_k-1) \\\\
\\endaligned\\]
令
\\[K_k = \\fracP_k-1\\vecx_k
\\lambda+\\vecx_k^\\mathrmTP_k-1\\vecx_k
\\]
则
\\[\\beginaligned
P_k =& \\frac1\\lambda(I-K_k\\vecx_k^\\mathrmT)P_k-1\\\\
P_k\\vecx_k =& \\frac1\\lambda(P_k-1\\vecx_k
-K_k\\vecx_k^\\mathrmTP_k-1\\vecx_k) \\\\
=& \\frac1\\lambda\\fracP_k-1\\vecx_k
(\\lambda+\\vecx_k^\\mathrmTP_k-1\\vecx_k)
-P_k-1\\vecx_k\\vecx_k^\\mathrmT
P_k-1\\vecx_k
\\lambda+\\vecx_k^\\mathrmTP_k-1\\vecx_k \\\\
=& \\frac1\\lambda\\fracP_k-1\\vecx_k\\lambda
\\lambda+\\vecx_k^\\mathrmTP_k-1\\vecx_k \\\\
=& K_k \\\\
\\endaligned\\]
总结成递推公式得到
\\[\\beginaligned
K_k =& \\fracP_k-1\\vecx_k\\lambda+\\vecx_k^\\mathrmT
P_k-1\\vecx_k\\\\
P_k =& \\frac1\\lambda(I-K_k\\vecx_k^\\mathrmT)P_k-1 \\\\
\\vec\\theta_k =& \\vec\\theta_k-1
+K_k(y_k-\\vecx_k^\\mathrmT\\theta_k-1)
\\endaligned\\]
与不带遗忘因子的公式相比，第3个式子 \\(\\vec\\theta_k\\) 不变，前两个稍微有变化。
一个例子与代码
对方程 \\(y=ax^2+bx+c\\)，输入多组数据，\\(x\\) 取一些高斯噪声，以4组数据为例，
\\[\\beginbmatrix
y_1 \\\\ y_2 \\\\ y_3 \\\\ y_4
\\endbmatrix
=\\beginbmatrix
x_1^2 & x_1 & 1 \\\\
x_2^2 & x_2 & 1 \\\\
x_3^2 & x_3 & 1 \\\\
x_4^2 & x_4 & 1 \\\\
\\endbmatrix
\\beginbmatrix
a \\\\ b \\\\ c
\\endbmatrix
\\]
#include <iostream>
#include <cmath>
#include "zhnmat.hpp"
using namespace zhnmat;
using namespace std;
int main()
constexpr double a = 0.5;
constexpr double b = 1.1;
constexpr double c = 2.1;
double U, V, randx;
double y, lambda=0.5;
Mat K;
Mat vx(3, 1);
Mat P = eye(3)*1e6;
Mat theta(3, 1);
for (int i = 0; i < 10; i++)
U = (rand()+1.0) / (RAND_MAX+1.0);
V = (rand()+1.0) / (RAND_MAX+1.0);
randx = sqrt(-2.0 * log(U))* cos(6.283185307179586477 * V); // `randx`服从标准正态分布
randx *= 2;
vx.set(0, 0, randx*randx);
vx.set(1, 0, randx);
vx.set(2, 0, 1);
y = a*randx*randx + b*randx + c;
K = P * vx / (lambda + (vx.T() * P * vx).at(0, 0));
P = (eye(3) - K * vx.T()) * P / lambda;
theta = theta + K * (y - (vx.T() * theta).at(0, 0));
cout << "theta" << theta << endl;
return 0;
最优化算法-递推最小二乘法
Recursive Least Square(RLS)
最小二乘算法(Least Square)解决的问题是一个多元线性拟合问题:
({a_1,a_2,a_3,...,a_n,b}), 其中(a_i)为自变量, (b)为响应值. 在线系统会不断获得新的观测值({a_1^i,a_2^i,a_3^i,...,a_n^i,b^i}). 在线观测保存所有的观测值并使用经典最小二乘法进行求解, 需要耗费大量的计算资源与内存, 所以需要有一种递推的类似动态规划的方式来在线更新线性模型的参数.
对于线性拟合问题(Ax=b):
[minlimits_x||Ax-b||^2
]
对应的经典最小二乘解可以写成:
[egin{array}{rcl}
A_0 x_0 = b_0 & & x_0=(A_0^TA_0)^{-1}A_0^Tb_0egin{bmatrix}A_0 \ A_1end{bmatrix}x_1 = egin{bmatrix} b_0\b_1 end{bmatrix}& &
x_1=egin{pmatrix}egin{bmatrix}A_0 \ A_1end{bmatrix}^Tegin{bmatrix}A_0 \ A_1end{bmatrix}end{pmatrix}^{-1}egin{bmatrix}A_0 \ A_1end{bmatrix}^Tegin{bmatrix}b_0 \ b_1end{bmatrix}\end{array}
]
接下来需要消去(A_0,b_0)两个变量
[egin{array}{rcl} G_0=A_0^TA_0 &&& x_0=G_0^{-1}A_0^Tb_0end{array} \begin{split}
G_1&=egin{bmatrix}A_0 \ A_1end{bmatrix}^T
egin{bmatrix}A_0 \ A_1end{bmatrix}\ &=G_0+A_1^TA_1 \end{split}
\begin{split}
egin{bmatrix}A_0 \ A_1end{bmatrix}^Tegin{bmatrix}b_0 \ b_1end{bmatrix}
&=A_0^T(A_0x_0)+A_1^Tb_1 \
&=G_0x_0+A_1^Tb_1 \
&=(G_1-A_1^TA_1)x_0 +A_1^Tb_1 &= G_1x_0 + A_1^T(b_1-A_1x_0)
end{split} egin{split}
x_1 &= G_1^{-1}G_1x_0 + G_1^{-1}A_1^T(b_1-A_1x_0) &= x_0 + G_1^{-1}A_1^T(b_1-A_1x_0)
end{split} x_{k+1} = x_k + G_{k+1}^{-1}A_{k+1}^T(b_{k+1}-A_{k+1}x_k)
]
得到线性模型的参数的迭代更新关系(x_{k+1} = x_k + G_{k+1}^{-1}A_{k+1}^T(b_{k+1}-A_{k+1}x_k))
但此时需要对(G_{k+1}=egin{bmatrix}A_k \ A_{k+1}end{bmatrix}^T egin{bmatrix}A_k \ A_{k+1}end{bmatrix}=G_k+A_{k+1}^TA_{k+1}\)求逆, 可以根据Sherman-Morrison-Woodbury formula进一步优化:
[(A+UV)^{-1}= A^{-1}-(A^{-1}U)(I+VA^{-1}U)^{-1}(VA^{-1}) \begin{split}
G_{k+1}^{-1} &= (G_k+A_{k+1}^TA_{k+1})^{-1} &= G_k^{-1}-G_k^{-1}A_{k+1}(I+A_{k+1}^TG_k^{-1}A_{k+1})^{-1}A_{k+1}^TG_k^{-1}G_{k+1}^{-1} equiv P_{k+1}&= P_k - P_kA_{k+1}(I+A_{k+1}^TP_kA_{k+1})^{-1}A_{k+1}^TP_k
end{split}
]
进一步计算((I+A_{k+1}^TG_k^{-1}A_{k+1})^{-1}), 考虑特殊情况(A_{k+1}equiv a_{k+1}^T),即新观测量为单行的向量, (I+A_{k+1}^TG_k^{-1}A_{k+1})退化成标量:
[egin{split}
P_{k+1} &= P_k - P_kA_{k+1}(I+A_{k+1}^TP_kA_{k+1})^{-1}A_{k+1}^TP_k &= P_k - frac{P_ka_{k+1}a_{k+1}^TP_k }{I+a_{k+1}^TP_ka_{k+1}}
end{split}egin{split}
x_{k+1} &= x_k + G_{k+1}^{-1}A_{k+1}^T(b_{k+1}-A_{k+1}x_k) &= x_k + P_{k+1}a_{k+1}(b_{k+1}-a_{k+1}^Tx_k)
end{split}
]
回顾以上推导, 可得如下的初始解和迭代过程
[egin{split}
P_0 &=(A_0^TA_0)^{-1}x_0 &= P_0A_0^Tb_0 P_{k+1} &= P_k - frac{P_ka_{k+1}a_{k+1}^TP_k }{I+a_{k+1}^TP_ka_{k+1}} x_{k+1} &= x_k + P_{k+1}a_{k+1}(b_{k+1}-a_{k+1}^Tx_k)
end{split}
]
Implementation
def RLS(A0: np.ndarray, b0: np.ndarray):
P0 = np.linalg.inv( np.dot(A0.T, A0) )
return P0, np.dot(P0, np.dot(A0.T, b0))
def RLS(Pk: np.ndarray, xk: np.ndarray, a:np.ndarray, b:np.ndarray):
P = Pk - np.dot(np.dot(Pk, a), np.dot(a.T, Pk))/(1+np.dot(np.dot(a.T, P_k), a))
x = xk + np.dot(np.dot(P, a), b - np.dot(a.T, x))
return P, x
Reference
以上是关于带遗忘因子的递推最小二乘法推导的主要内容，如果未能解决你的问题，请参考以下文章