带遗忘因子的递推最小二乘法推导

摘要： 首先，已知最小二乘法的解后，推导出解的递推形式；然后通过一个简单的公式y=x+w的形式引入遗忘因子，并改写成递推形式，最后将遗忘因子引入带递推的最小二乘法中。

递推最小二乘法

对多组数据 ${\vec{x}}_i$ 和 $y_i$，满足
$$y_i={\vec{x}}_i^{\mathrm{T}}\vec{\theta }$$
其中 ${\vec{x}}_i$ 是输入数据向量，$y_i$ 是输出数据标量。写成矩阵形式
$$\vec{y}=X\vec{\theta }$$
其中（以3组数据、2个未知参数为例）
$$\vec{y}=\left[\begin{array}{c}{y}_1\\ y_2\\ y_3\end{array}\right],X=\left[\begin{array}{c}{\vec{x}}_1^{\mathrm{T}}\\ {\vec{x}}_2^{\mathrm{T}}\\ {\vec{x}}_3^{\mathrm{T}}\end{array}\right],\theta =\left[\begin{array}{c}{\theta }_1\\ {\theta }_2\end{array}\right],y_i=x_{i1}{\theta }_1+x_{i2}{\theta }_2$$
完整地写作（以3组数据、2个未知参数为例）
$$\left[\begin{array}{c}{y}_1\\ y_2\\ y_3\end{array}\right]=\left[\begin{array}{cc}{x}_{11}& x_{12}\\ x_{21}& x_{22}\\ x_{31}& x_{32}\end{array}\right]\left[\begin{array}{c}{\theta }_1\\ {\theta }_2\end{array}\right]$$
有$k$组数据时记作
$$\begin{gathered} X_k^{\mathrm{T}}=\left[\begin{array}{cccc}{\vec{x}}_1& {\vec{x}}_2& \cdots & {\vec{x}}_k\end{array}\right] \\ {\vec{y}}_k=\left[\begin{array}{cccc}{y}_1& y_2& \cdots & y_k\end{array}\right] \end{gathered}$$
(此时注意区分一些符号，例如 ${\vec{y}}_k$ 与 $y_k$，$X$ 与 $\vec{x}$ 等)
已知最小二乘法的解为
$$\vec{\theta }=(X^{\mathrm{T}}X{)}^{-1}{X}^{\mathrm{T}}\vec{y}\text{\hspace{1em}(1)}$$
令 $P=(X^{\mathrm{T}}X{)}^{-1}$，
则加上递推后
$$\begin{array}{rl}{P}_k^{-1}=& X_k^{\mathrm{T}}{X}_k=\sum _{i=1}^k{\vec{x}}_i{\vec{x}}_i^{\mathrm{T}}\\ =& \sum _{i=1}^{k-1}{\vec{x}}_i{\vec{x}}_i^{\mathrm{T}}+{\vec{x}}_k{\vec{x}}_k^{\mathrm{T}}\\ =& P_{k-1}^{-1}+{\vec{x}}_k{\vec{x}}_k^{\mathrm{T}}\end{array}\text{\hspace{1em}(2)}$$
同理可得
$$X_k^{\mathrm{T}}{\vec{y}}_k=X_{k-1}^{\mathrm{T}}{\vec{y}}_{k-1}+{\vec{x}}_k{y}_k\text{\hspace{1em}(3)}$$
于是代入式(1)得到
$$\begin{array}{rl}{\vec{\theta }}_k=& P_k{X}_k^{\mathrm{T}}{\vec{y}}_k\\ =& P_k(X_{k-1}^{\mathrm{T}}{\vec{y}}_{k-1}+{\vec{x}}_k{y}_k)\\ =& P_k(P_{k-1}^{-1}{\vec{\theta }}_{k-1}+{\vec{x}}_k{y}_k)\\ =& P_k(P_k^{-1}{\vec{\theta }}_{k-1}-{\vec{x}}_k{\vec{x}}_k^{\mathrm{T}}{\theta }_{k-1}+{\vec{x}}_k{y}_k)\\ =& {\vec{\theta }}_{k-1}+P_k{\vec{x}}_k(y_k-{\vec{x}}_k^{\mathrm{T}}{\theta }_{k-1})\\ =& {\vec{\theta }}_{k-1}+(P_{k-1}^{-1}+{\vec{x}}_k{\vec{x}}_k^{\mathrm{T}}{)}^{-1}{\vec{x}}_k(y_k-{\vec{x}}_k^{\mathrm{T}}{\theta }_{k-1})\\ =& {\vec{\theta }}_{k-1}+(P_{k-1}-\frac{P_{k-1}{\vec{x}}_k{\vec{x}}_k^{\mathrm{T}}{P}_{k-1}}{1+{\vec{x}}_k^{\mathrm{T}}{P}_{k-1}{\vec{x}}_k}){\vec{x}}_k(y_k-{\vec{x}}_k^{\mathrm{T}}{\theta }_{k-1})\end{array}$$
其中倒数第3行到倒数第一行的过程
$$P_k=P_{k-1}-\frac{P_{k-1}{\vec{x}}_k{\vec{x}}_k^{\mathrm{T}}{P}_{k-1}}{1+{\vec{x}}_k^{\mathrm{T}}{P}_{k-1}{\vec{x}}_k}$$
用到了下面的矩阵求逆公式及其引理
$$\begin{array}{rl}(A+BCD{)}^{-1}=& A^{-1}-A^{-1}B(D{A}^{-1}B+C^{-1}{)}^{-1}D{A}^{-1}\\ (A+\vec{u}{\vec{u}}^{\text{T}}{)}^{-1}=& A^{-1}-\frac{A^{-1}\vec{u}{\vec{u}}^{\text{T}}{A}^{-1}}{1+{\vec{u}}^{\text{T}}{A}^{-1}\vec{u}}\end{array}$$
令
$$K_k=\frac{P_{k-1}{\vec{x}}_k}{1+{\vec{x}}_k^{\mathrm{T}}{P}_{k-1}{\vec{x}}_k}$$
则
$$\begin{array}{rl}{P}_k=& (I-K_k{\vec{x}}_k^{\mathrm{T}})P_{k-1}\\ P_k{\vec{x}}_k=& P_{k-1}{\vec{x}}_k-K_k{\vec{x}}_k^{\mathrm{T}}{P}_{k-1}{\vec{x}}_k\\ =& \frac{P_{k-1}{\vec{x}}_k(1+{\vec{x}}_k^{\mathrm{T}}{P}_{k-1}{\vec{x}}_k)-P_{k-1}{\vec{x}}_k{\vec{x}}_k^{\mathrm{T}}{P}_{k-1}{\vec{x}}_k}{1+{\vec{x}}_k^{\mathrm{T}}{P}_{k-1}{\vec{x}}_k}\\ =& \frac{P_{k-1}{\vec{x}}_k}{1+{\vec{x}}_k^{\mathrm{T}}{P}_{k-1}{\vec{x}}_k}\\ =& K_k\end{array}$$
总结成递推公式得到
$$\begin{array}{rl}{K}_k=& \frac{P_{k-1}{\vec{x}}_k}{1+{\vec{x}}_k^{\mathrm{T}}{P}_{k-1}{\vec{x}}_k}\\ P_k=& (I-K_k{\vec{x}}_k^{\mathrm{T}})P_{k-1}\\ {\vec{\theta }}_k=& {\vec{\theta }}_{k-1}+K_k(y_k-{\vec{x}}_k^{\mathrm{T}}{\theta }_{k-1})\end{array}$$

一些注意事项

  因为 $P=(X^{\text{T}}X{)}^{-1}$，${\vec{x}}_1{\vec{x}}_1^{\text{T}}$ 一般很小，所以实际使用时一般把 $P_1$ 初始化为一个接近无穷大的单位阵。
  可以看出 $P$ 是个对称矩阵，如果因为数值计算的原因导致 $P$ 不严格对称时，代码里需要加一步：$P\leftarrow (P+P^{\text{T}})/2$。

引入遗忘因子

先考虑一个比较简单的情况
$$x_n=\theta +w_n$$
其中 $\theta$ 是要估计的参数，$w_n$ 是高斯白噪声，$x_n$ 是观测到的数据，总共有 $N$ 组数据，将 $N$ 组数据取平均得到 $\theta$ 的一个估计值为
$$\hat{\theta }=\frac{1}{N}\sum _{n=1}^N{x}_n$$
写成递推形式为
$${\hat{\theta }}_n={\hat{\theta }}_{n-1}+\frac{1}{n}(x_n-{\hat{\theta }}_{n-1})$$
例如，
$$\begin{array}{rl}{\hat{\theta }}_1=& x_1\\ {\hat{\theta }}_2=& {\hat{\theta }}_1+\frac{1}{2}(x_2-{\hat{\theta }}_1)=\frac{1}{2}(x_1+x_2)\\ {\hat{\theta }}_3=& {\hat{\theta }}_2+\frac{1}{3}(x_3-{\hat{\theta }}_2)=\frac{1}{3}(x_1+x_2+x_3)\end{array}$$
现在对之前的值乘一个衰减系数 $\lambda \in (0,1)$，也就是
$$\begin{array}{rl}{\hat{\theta }}_2=& \frac{1}{\lambda +1}(\lambda x_1+x_2)\\ {\hat{\theta }}_3=& \frac{1}{{\lambda }^2+\lambda +1}({\lambda }^2{x}_1+\lambda x_2+x_3)\\ ⋮& \end{array}$$
写成递推形式为
$$\begin{array}{rl}{\hat{\theta }}_2=& \frac{1}{\lambda +1}(\lambda x_1+x_2)\\ {\hat{\theta }}_3=& \frac{1}{{\lambda }^2+\lambda +1}(\lambda (\lambda +1){\hat{\theta }}_2+x_3)\\ =& {\hat{\theta }}_2+\frac{1}{{\lambda }^2+\lambda +1}(x_3-{\hat{\theta }}_2)\end{array}$$
当 $n\to \infty$ 时，递推形式可以近似为
$$\begin{array}{rl}{\hat{\theta }}_n=& {\hat{\theta }}_{n-1}+\frac{1}{{\sum }_{i=0}^{n-1}{\lambda }^i}(x_n-{\hat{\theta }}_{n-1})\\ =& {\hat{\theta }}_{n-1}+(1-\lambda )(x_n-{\hat{\theta }}_{n-1})\end{array}$$

递推最小二乘中的遗忘因子

类似地，在最小二乘的误差 $S_i=||y_i-{\vec{x}}_i^{\text{T}}\vec{\theta }|{|}^2$ 中加入遗忘因子，得到总误差（以3组数据为例）
$$J_3={\lambda }^2{S}_1+\lambda S_2+S_3$$
完整的方程可以写作
$$\begin{gathered} {\vec{y}}_3=X_3\theta \\ \left[\begin{array}{c}\lambda y_1\\ \sqrt{\lambda }{y}_2\\ y_3\end{array}\right]=\left[\begin{array}{cc}\lambda x_{11}& \lambda x_{12}\\ \sqrt{\lambda }{x}_{21}& \sqrt{\lambda }{x}_{22}\\ x_{31}& x_{32}\end{array}\right]\left[\begin{array}{c}{\theta }_1\\ {\theta }_2\end{array}\right] \end{gathered}$$
类似地代入最小二乘的解式(1)，并求式(2)(3)的递推形式
$$\begin{array}{rl}{X}_2^{\text{T}}{X}_2=& \left[\begin{array}{cc}\sqrt{\lambda }{x}_{11}& x_{21}\\ \sqrt{\lambda }{x}_{12}& x_{22}\end{array}\right]\left[\begin{array}{cc}\sqrt{\lambda }{x}_{11}& \sqrt{\lambda }{x}_{12}\\ x_{21}& x_{22}\end{array}\right]\\ X_3^{\text{T}}{X}_3=& \left[\begin{array}{ccc}\lambda x_{11}& \sqrt{\lambda }{x}_{21}& x_{31}\\ \lambda x_{12}& \sqrt{\lambda }{x}_{22}& x_{32}\end{array}\right]\left[\begin{array}{cc}\lambda x_{11}& \lambda x_{12}\\ \sqrt{\lambda }{x}_{21}& \sqrt{\lambda }{x}_{22}\\ x_{31}& x_{32}\end{array}\right]\\ =& \lambda X_2^{\text{T}}{X}_2+{\vec{x}}_3{\vec{x}}_3^{\text{T}}\\ X_2^{\text{T}}{\vec{y}}_2=& \left[\begin{array}{cc}\sqrt{\lambda }{x}_{11}& x_{21}\\ \sqrt{\lambda }{x}_{12}& x_{22}\end{array}\right]\left[\begin{array}{c}\sqrt{\lambda }{y}_1\\ y_2\end{array}\right]\\ X_3^{\text{T}}{\vec{y}}_3=& \left[\begin{array}{ccc}\lambda x_{11}& \sqrt{\lambda }{x}_{21}& x_{31}\\ \lambda x_{12}& \sqrt{\lambda }{x}_{22}& x_{32}\end{array}\right]\left[\begin{array}{c}\lambda y_1\\ \sqrt{\lambda }{y}_2\\ y_3\end{array}\right]\\ =& \lambda X_2^{\text{T}}{\vec{y}}_2+{\vec{x}}_3{y}_3\end{array}$$
即
$$\begin{array}{rl}{P}_k^{-1}=& \lambda P_{k-1}^{-1}+{\vec{x}}_k{\vec{x}}_k^{\text{T}}\\ X_k^{\text{T}}{\vec{y}}_k=& \lambda X_{k-1}^{\text{T}}{\vec{y}}_{k-1}+{\vec{x}}_k{y}_k\end{array}$$
代入式(1)得到
$$\begin{array}{rl}{\vec{\theta }}_k=& P_k{X}_k^{\text{T}}{\vec{y}}_k\\ =& P_k(\lambda X_{k-1}^{\text{T}}{\vec{y}}_{k-1}+{\vec{x}}_k{y}_k)\\ =& P_k(\lambda P_{k-1}^{-1}{\vec{\theta }}_{k-1}+{\vec{x}}_k{y}_k)\\ =& P_k(P_k^{-1}{\vec{\theta }}_{k-1}-{\vec{x}}_k{\vec{x}}_k^{\text{T}}{\theta }_{k-1}+{\vec{x}}_k{y}_k)\\ =& {\vec{\theta }}_{k-1}+P_k{\vec{x}}_k(y_k-{\vec{x}}_k^{\text{T}}{\theta }_{k-1})\\ =& {\vec{\theta }}_{k-1}+(\lambda P_{k-1}^{-1}+{\vec{x}}_k{\vec{x}}_k^{\text{T}}{)}^{-1}{\vec{x}}_k(y_k-{\vec{x}}_k^{\text{T}}{\theta }_{k-1})\\ =& {\vec{\theta }}_{k-1}+(\frac{P_{k-1}}{\lambda }-\frac{\frac{P_{k-1}}{\lambda }{\vec{x}}_k{\vec{x}}_k^{\text{T}}\frac{P_{k-1}}{\lambda }}{1+{\vec{x}}_k^{\text{T}}\frac{P_{k-1}}{\lambda }{\vec{x}}_k}){\vec{x}}_k(y_k-{\vec{x}}_k^{\text{T}}{\theta }_{k-1})\\ =& {\vec{\theta }}_{k-1}+(\frac{P_{k-1}}{\lambda }-\frac{1}{\lambda }\frac{P_{k-1}{\vec{x}}_k{\vec{x}}_k^{\text{T}}{P}_{k-1}}{\lambda +{\vec{x}}_k^{\text{T}}{P}_{k-1}{\vec{x}}_k}){\vec{x}}_k(y_k-{\vec{x}}_k^{\text{T}}{\theta }_{k-1})\end{array}$$
令
$$K_k=\frac{P_{k-1}{\vec{x}}_k}{\lambda +{\vec{x}}_k^{\mathrm{T}}{P}_{k-1}{\vec{x}}_k}$$
则
$$\begin{array}{rl}{P}_k=& \frac{1}{\lambda }(I-K_k{\vec{x}}_k^{\mathrm{T}})P_{k-1}\\ P_k{\vec{x}}_k=& \frac{1}{\lambda }(P_{k-1}{\vec{x}}_k-K_k{\vec{x}}_k^{\mathrm{T}}{P}_{k-1}{\vec{x}}_k)\\ =& \frac{1}{\lambda }\frac{P_{k-1}{\vec{x}}_k(\lambda +{\vec{x}}_k^{\mathrm{T}}{P}_{k-1}{\vec{x}}_k)-P_{k-1}{\vec{x}}_k{\vec{x}}_k^{\mathrm{T}}{P}_{k-1}{\vec{x}}_k}{\lambda +{\vec{x}}_k^{\mathrm{T}}{P}_{k-1}{\vec{x}}_k}\\ =& \frac{1}{\lambda }\frac{P_{k-1}{\vec{x}}_k\lambda }{\lambda +{\vec{x}}_k^{\mathrm{T}}{P}_{k-1}{\vec{x}}_k}\\ =& K_k\end{array}$$
总结成递推公式得到
$$\begin{array}{rl}{K}_k=& \frac{P_{k-1}{\vec{x}}_k}{\lambda +{\vec{x}}_k^{\mathrm{T}}{P}_{k-1}{\vec{x}}_k}\\ P_k=& \frac{1}{\lambda }(I-K_k{\vec{x}}_k^{\mathrm{T}})P_{k-1}\\ {\vec{\theta }}_k=& {\vec{\theta }}_{k-1}+K_k(y_k-{\vec{x}}_k^{\mathrm{T}}{\theta }_{k-1})\end{array}$$
与不带遗忘因子的公式相比，第3个式子 ${\vec{\theta }}_k$ 不变，前两个稍微有变化。

一个例子与代码

  对方程 $y=a{x}^2+bx+c$，输入多组数据，$x$ 取一些高斯噪声，以4组数据为例，
$$\left[\begin{array}{c}{y}_1\\ y_2\\ y_3\\ y_4\end{array}\right]=\left[\begin{array}{ccc}{x}_1^2& x_1& 1\\ x_2^2& x_2& 1\\ x_3^2& x_3& 1\\ x_4^2& x_4& 1\end{array}\right]\left[\begin{array}{c}a\\ b\\ c\end{array}\right]$$
代码使用了自用的矩阵运算库zhnmat，见 自用的矩阵运算库zhnmat使用说明

#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;
}