基本最小二乘到递推最小二乘

基本最小二乘（LS）

先导知识：

从函数出发

  假设一个函数$$y=\left[\begin{array}{cccc}{\theta }_1& {\theta }_2& \cdots & {\theta }_n\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_n\end{array}\right]=\mathbf{\theta }X=\sum _{i=1}^n{\theta }_i{x}_i$$

我们约定粗体和大写字母均表示矩阵或者向量。

残差

  现在假设我们有一个这样的函数函数$y=\mathbf{\theta }X$但是我们并不知道$\mathbf{\theta }$的各个量的具体数值。我们只能用一系列（一共$k$组）的
$$X_j=\left[\begin{array}{c}{x}_1^j\\ x_2^j\\ ⋮\\ x_n^j\end{array}\right]\in {\mathbb{R}}^{n\times 1}(j=1,2,\cdots ,k)$$作为输入往这个函数的入口放（这里的$j$并不是次方的意思，而是表示这是第$j$组数据中的某个量）。然后我们得到了$k$组输出
$$Y=\left[\begin{array}{cccc}{y}_1& y_2& \cdots & y_k\end{array}\right]\in {\mathbb{R}}^{1\times k}$$
也即：
$$Y=\mathbf{\theta }\left[\begin{array}{cccc}{X}_1& X_2& \cdots & X_k\end{array}\right]\in ({\mathbb{R}}^{1\times n}\ast {\mathbb{R}}^{n\times k})={\mathbb{R}}^{1\times k}$$
  我们可以估计出一个$\mathbf{\theta }$的取值$\hat{\mathbf{\theta }}$。我们再把$X_j(j=1,2,\cdots ,k)$代入我们估计出的方程模型$y=\hat{\mathbf{\theta }}X$。由于估计的模型不可能完全拟合原始数据，所以我们会得到一组拟合值$\hat{Y}=\hat{\mathbf{\theta }}\left[\begin{array}{cccc}{X}_1& X_2& \cdots & X_j\end{array}\right]\in {\mathbb{R}}^{1\times k}$。
  我们把$\mathbf{e}=Y-\hat{Y}\in {\mathbb{R}}^{1\times k}$称为残差

梳理

  我们在这里将一些常数和符号做一个约定：

• n——函数的输入变量$x_i$的个数，同时也是${\theta }_i$的总数
• k——测试用的数据总组数
• $\mathbf{\Phi }$——$\left[\begin{array}{cccc}{X}_1& X_2& \cdots & X_k\end{array}\right]\in {\mathbb{R}}^{n\times k}$
• $i$——一般作为输入数据及其参数的下标（${\theta }_i{x}_i$）
• $j$——一般作为输入组数数据的下标（$X_j$）

推导

基本思想：

  我们再估计参数$\hat{\mathbf{\theta }}$时，我们自然是希望与客观存在的值$\mathbf{\theta }$尽可能的接近。因此，我们需要引入一个评价我们估计的好坏的一个“标准”。前面提到的残差可以作为一个很好的标准。
  但是我们发现$e\in {\mathbb{R}}^{1\times k}$作为一个向量，不好作为一个直观的标准。同时它的每一项中既有正也有负，自然是不适合直接作为评价标准来使用。所以我们引入一个指标函数$J$$$J=\sum _{k=n+1}^{n+N}{e}^2(k)\stackrel{def}{=}\mathbf{e}\cdot {\mathbf{e}}^T=(\mathbf{Y}-\hat{\mathbf{\theta }}\Phi )(\mathbf{Y}-\hat{\mathbf{\theta }}\Phi {)}^T$$

开始推导

  所以，最小二乘法就是让指标函数$J$最小的参数估计方法。既有：
$$\hat{\mathbf{\theta }}=\underset{\mathbf{\theta }}{\min }J$$

  而$J$取最小值，我们先讨论J为极值时：
$$\frac{\partial J}{\partial \hat{\mathbf{\theta }}}=0$$

$$\Rightarrow \frac{\partial [(\mathbf{Y}-\hat{\mathbf{\theta }}\mathbf{\Phi })(\mathbf{Y}-\hat{\mathbf{\theta }}\mathbf{\Phi }{)}^T]}{\partial \hat{\mathbf{\theta }}}=0$$

$$\Rightarrow -2\mathbf{\Phi }(\mathbf{Y}-\hat{\mathbf{\theta }}\mathbf{\Phi }{)}^T=0$$

$$\mathbf{\Phi }{\mathbf{\Phi }}^{\mathbf{T}}{\hat{\mathbf{\theta }}}^{\mathbf{T}}=\mathbf{\Phi }{\mathbf{Y}}^{\mathbf{T}}$$

其中，$\mathbf{\Phi }{\mathbf{\Phi }}^{\mathbf{T}}\in {\mathbb{R}}^{n\times n}$。
若其逆矩阵存在，则：
$${\hat{\mathbf{\theta }}}^{\mathbf{T}}=(\mathbf{\Phi }{\mathbf{\Phi }}^{\mathbf{T}}{)}^{-1}\mathbf{\Phi }{\mathbf{Y}}^{\mathbf{T}}$$

$$\hat{\mathbf{\theta }}=\mathbf{Y}{\mathbf{\Phi }}^{\mathbf{T}}(\mathbf{\Phi }{\mathbf{\Phi }}^{\mathbf{T}}{)}^{-1}$$

上述结果只是$J$为极值时的结论，$J$可能是极大值也可能是极小值。我们进一步讨论，要使$J$为极小值的条件为
$$\frac{{\partial }^{2}J}{\partial {\hat{\mathbf{\theta }}}^2}>0$$

$$\frac{\partial J}{\partial \hat{\mathbf{\theta }}}=-2\mathbf{\Phi }(\mathbf{Y}-\mathbf{\Phi }\hat{\mathbf{\theta }}{)}^T$$

$$\Rightarrow \frac{{\partial }^{2}J}{\partial {\hat{\mathbf{\theta }}}^2}=2\mathbf{\Phi }{\mathbf{\Phi }}^{\mathbf{T}}>0$$

也就是 $\mathbf{\Phi }{\mathbf{\Phi }}^{\mathbf{T}}$为正定矩阵。

递推最小二乘法

• 暂留2020/11/4，于2020/11/13完成

背景

  基本最小二乘法（LS）有诸多缺点，例如对于一组动态的数据，每次接收到新数据，就要全部重算一遍。这种重复的计算的成本很大，导致实用性不好。所以对于一组离线数据，基本最小二乘法是适用的。但是如果是实时统计分析一系列数据，那么基本最小二乘法就会遇到计算困难。
  我们希望在获取一个新数据时，可以直接使用该数据和已经计算过的结果进行某种运算，达到“修正”旧结果的目的。

前N个输入输出数据

  为此，我们假设在某个时间点已经获取了$N$组数据：

1. ${\mathbf{\Phi }}_{\mathbf{N}}=\left[\begin{array}{cccc}{X}_1& X_2& \cdots & X_N\end{array}\right]\in {\mathbb{R}}^{n\times N}$（N组输入，每一组有n个分量）
2. ${\mathbf{Y}}_{\mathbf{N}}=\left[\begin{array}{c}{y}_1,y_2,\cdots ,y_N\end{array}\right]=\mathbf{\theta }{\mathbf{\Phi }}_{\mathbf{N}}=\mathbf{\theta }\left[\begin{array}{cccc}{X}_1& X_2& \cdots & X_N\end{array}\right]\in ({\mathbb{R}}^{1\times n}\ast {\mathbb{R}}^{n\times N})={\mathbb{R}}^{1\times N}$（就是前N个输入对应的N个输出）
3. ${\hat{\mathbf{\theta }}}_{\mathbf{N}}={\mathbf{Y}}_{\mathbf{N}}{{\mathbf{\Phi }}_{\mathbf{N}}}^T({\mathbf{\Phi }}_{\mathbf{N}}{\mathbf{\Phi }}_{\mathbf{N}}^{\mathbf{T}}{)}^{-1}$
4. ${\tilde{\mathbf{\theta }}}_{\mathbf{N}}=\mathbf{\theta }-{\hat{\mathbf{\theta }}}_{\mathbf{N}}$
5. $Var\hat{\mathbf{\theta }}={\sigma }^2({\mathbf{\Phi }}_{\mathbf{N}}{\mathbf{\Phi }}_{\mathbf{N}}^{\mathbf{T}}{)}^{-1}$

  我们记${\mathbf{P}}_{\mathbf{N}}=({\mathbf{\Phi }}_{\mathbf{N}}{\mathbf{\Phi }}_{\mathbf{N}}^{\mathbf{T}}{)}^{-1}\in {\mathbb{R}}^{n\times n}$，则
$${\hat{\mathbf{\theta }}}_{\mathbf{N}}={\mathbf{Y}}_{\mathbf{N}}{{\mathbf{\Phi }}_{\mathbf{N}}}^T{\mathbf{P}}_{\mathbf{N}}$$

  现在我们已经拥有了$N$组$I/O$数据，现在我们需要结合第$N+1$个新数据来修正我们的估计${\hat{\mathbf{\theta }}}_{\mathbf{N}}$得到${\hat{\mathbf{\theta }}}_{\mathbf{N}+1}$。这个过程就是一种递推，我们需要得到这种递推的通法。我们记之为：

$${\hat{\mathbf{\theta }}}_{\mathbf{N}+1}=f({\hat{\mathbf{\theta }}}_{\mathbf{N}},{\mathbf{X}}_{N+1},y_{N+1})$$

开始递推

注意到：
$${\hat{\mathbf{\theta }}}_{N+1}={\mathbf{Y}}_{N+1}{\mathbf{\Phi }}_{N+1}^T{\mathbf{P}}_{N+1}$$

先分析${\mathbf{P}}_{N+1}\in {\mathbb{R}}^{n\times n}$
$${\mathbf{P}}_{N+1}=({\mathbf{\Phi }}_{\mathbf{N}+1}{\mathbf{\Phi }}_{\mathbf{N}+1}^{\mathbf{T}}{)}^{-1}$$

其中
$${\mathbf{\Phi }}_{\mathbf{N}+1}{\mathbf{\Phi }}_{\mathbf{N}+1}^{\mathbf{T}}=\left[\begin{array}{cccc}{X}_1& X_2& \cdots & X_{N+1}\end{array}\right]\left[\begin{array}{c}{X}_1^T\\ X_2^T\\ ⋮\\ X_{N+1}^T\end{array}\right]$$
$$=\sum _{i=1}^{N+1}{X}_i{X}_i^T=\sum _{i=1}^N{X}_i{X}_i^T+X_{N+1}{X}_{N+1}^T={\mathbf{\Phi }}_{\mathbf{N}}{\mathbf{\Phi }}_{\mathbf{N}}^{\mathbf{T}}+X_{N+1}{X}_{N+1}^T$$

故而
$${\mathbf{P}}_{N+1}^{-1}={\mathbf{P}}_N^{-1}+X_{N+1}{X}_{N+1}^T$$

$$\Rightarrow \{\begin{array}{l}{\mathbf{P}}_{N+1}=({\mathbf{P}}_N^{-1}+X_{N+1}{X}_{N+1}^T{)}^{-1}\\ {\mathbf{P}}_N=({\mathbf{P}}_{N+1}^{-1}-X_{N+1}{X}_{N+1}^T{)}^{-1}\end{array}$$

再分析${\mathbf{Y}}_N{\mathbf{\Phi }}_N^T\in {\mathbb{R}}^{1\times n}$
$${\mathbf{Y}}_{N+1}{\mathbf{\Phi }}_{N+1}^T=\left[\begin{array}{cccc}{y}_1& y_2& \cdots & y_{N+1}\end{array}\right]\left[\begin{array}{c}{X}_1^T\\ X_2^T\\ ⋮\\ X_{N+1}^T\end{array}\right]$$

$$=\sum _{i=1}^{N+1}{y}_i{X}_i^T=\sum _{i=1}^N{y}_i{X}_i^T+y_{N+1}\cdot X_{N+1}^T={\mathbf{Y}}_N{\mathbf{\Phi }}_N^T+y_{N+1}\cdot X_{N+1}^T$$

  我们得到了关于${\hat{\mathbf{\theta }}}_{N+1}$的两个部分的递推式，我们将其代入到${\hat{\mathbf{\theta }}}_{N+1}$中：
$${\hat{\mathbf{\theta }}}_{N+1}={\mathbf{Y}}_{N+1}{\mathbf{\Phi }}_{N+1}^T{\mathbf{P}}_{N+1}=({\mathbf{Y}}_N{\mathbf{\Phi }}_N^T+y_{N+1}\cdot X_{N+1}^T)\cdot ({\mathbf{P}}_N^{-1}+X_{N+1}{X}_{N+1}^T{)}^{-1}$$

$$={\mathbf{Y}}_N{\mathbf{\Phi }}_N^T\cdot {\mathbf{P}}_{N+1}+y_{N+1}\cdot X_{N+1}^T\cdot {\mathbf{P}}_{N+1}$$

又因为：
$${\hat{\mathbf{\theta }}}_N={\mathbf{Y}}_N{\mathbf{\Phi }}_N^T{\mathbf{P}}_N\Rightarrow {\hat{\mathbf{\theta }}}_N{\mathbf{P}}_N^{-1}={\mathbf{Y}}_N{\mathbf{\Phi }}_N^T$$

所以
$${\hat{\mathbf{\theta }}}_{N+1}={\hat{\mathbf{\theta }}}_N{\mathbf{P}}_N^{-1}\cdot {\mathbf{P}}_{N+1}+y_{N+1}\cdot X_{N+1}^T\cdot {\mathbf{P}}_{N+1}$$

$$={\hat{\mathbf{\theta }}}_N\cdot ({\mathbf{P}}_{N+1}^{-1}-X_{N+1}{X}_{N+1}^T){\mathbf{P}}_{N+1}+y_{N+1}\cdot X_{N+1}^T\cdot {\mathbf{P}}_{N+1}$$

$$={\hat{\mathbf{\theta }}}_N-{\hat{\mathbf{\theta }}}_N{X}_{N+1}{X}_{N+1}^T{\mathbf{P}}_{N+1}+y_{N+1}\cdot X_{N+1}^T\cdot {\mathbf{P}}_{N+1}$$

$$={\hat{\mathbf{\theta }}}_N+(y_{N+1}-{\hat{\mathbf{\theta }}}_N{X}_{N+1})\cdot X_{N+1}^T{\mathbf{P}}_{N+1}$$

我们令${\mathbf{K}}_{N+1}=X_{N+1}^T{\mathbf{P}}_{N+1}$，${\epsilon }_{N+1}=y_{N+1}-{\hat{\mathbf{\theta }}}_N{X}_{N+1}$。

  综上，我们得到了${\hat{\mathbf{\theta }}}_{N+1}={\hat{\mathbf{\theta }}}_N+$修正量 的形式：

• ${\mathbf{P}}_{N+1}=({\mathbf{P}}_N^{-1}+X_{N+1}{X}_{N+1}^T{)}^{-1}$
• ${\mathbf{K}}_{N+1}=X_{N+1}^T{\mathbf{P}}_{N+1}$
• ${\epsilon }_{N+1}=y_{N+1}-{\hat{\mathbf{\theta }}}_N{X}_{N+1}$
• ${\hat{\mathbf{\theta }}}_{N+1}={\hat{\mathbf{\theta }}}_N+{\epsilon }_{N+1}{\mathbf{K}}_{N+1}$

递推优化

  我们发现，${\mathbf{P}}_{N+1}=({\mathbf{P}}_N^{-1}+X_{N+1}{X}_{N+1}^T{)}^{-1}$需要求逆，很麻烦。我们引入公式：
$$[A+BCD{]}^{-1}=A^{-1}-A^{-1}B[C^{-1}+D{A}^{-1}B{]}^{-1}D{A}^{-1}$$

$${\mathbf{P}}_{N+1}=({\mathbf{P}}_N^{-1}+X_{N+1}{X}_{N+1}^T{)}^{-1}$$

$$Note:A={\mathbf{P}}_N^{-1},B={\mathbf{X}}_{N+1},C=1,D={\mathbf{X}}_{N+1}^T$$

得
$${\mathbf{P}}_{N+1}={\mathbf{P}}_N-\frac{{\mathbf{P}}_N{\mathbf{X}}_{N+1}{\mathbf{X}}_{N+1}^T{\mathbf{P}}_N}{1+{\mathbf{X}}_{N+1}^T{\mathbf{P}}_N{\mathbf{X}}_{N+1}}$$

结论

对于函数：
$$y=\left[\begin{array}{cccc}{\theta }_1& {\theta }_2& \cdots & {\theta }_n\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_n\end{array}\right]=\mathbf{\theta }X=\sum _{i=1}^n{\theta }_i{x}_i$$
的递推最小二乘估计：

• ¹${\mathbf{P}}_{N+1}={\mathbf{P}}_N-\frac{{\mathbf{P}}_N{\mathbf{X}}_{N+1}{\mathbf{X}}_{N+1}^T{\mathbf{P}}_N}{1+{\mathbf{X}}_{N+1}^T{\mathbf{P}}_N{\mathbf{X}}_{N+1}}$
• ${\mathbf{K}}_{N+1}=X_{N+1}^T{\mathbf{P}}_{N+1}$
• ${\epsilon }_{N+1}=y_{N+1}-{\hat{\mathbf{\theta }}}_N{X}_{N+1}$
• ${\hat{\mathbf{\theta }}}_{N+1}={\hat{\mathbf{\theta }}}_N+{\epsilon }_{N+1}{\mathbf{K}}_{N+1}$

我们需要的是一个估计初值${\hat{\mathbf{\theta }}}_0$和${\mathbf{P}}_0$。下面给出常用的初值取值的方法

Matlab 示例

代码部分

为了进一步加深理解，这里附上一段Matlab实现的递推最小二乘代码。

%递推最小二乘示例代码
%完成时间：2020/11/14
%作者：lamphungry
%原创代码，供大家参考
clc;clear;
times=2;                                    %重复4个周期
x1=0;x0=0;u1=0;u0=0;N=50+2;                 %定义状态参数的4个初值(取0)，并定义状态参数的总个数
a1=1.5;a2=-0.7;b1=1;b2=0.5;                 %初始化系统参数(这是我们需要拟合求的参数))
f=@(x1,x0,u1,u0) a1*x1+a2*x0+b1*u1+b2*u0;   %递推函数(这是我们需要拟合的线性多元方程)
x=zeros(1,times*N);x(1:2)=[x0 x1];          %预分配空间并初始化前两个值

%定义输入，使用随机0-1序列，总数为N-1
u=(idinput(N)'+1)/2;u(1:2)=[u0 u1];
u=repmat(u,1,times);
T=1;n=1:times*N;t=n*T;

%定义受噪声的输出和理论输出，引入高斯白噪声,\delta_v^2=0.01^2
delta_v=0.01;
z=zeros(1,times*N);z_ori=zeros(1,times*N);
z(1:2)=x(1:2)+delta_v*randn(1,2);			%掺杂噪声
z_ori(1:2)=x(1:2);							%无噪声

for i=3:times*N
   x(i)=f(x(i-1),x(i-2),u(i-1),u(i-2));
   z(i-2)=x(i)+delta_v*randn(1);
   z_ori(i-2)=x(i);
end

figure('Name','噪声输出和理论输出');
plot(t,z,t,z_ori);
legend('受噪声的输出','理论输出');

n=4;%待估计的数值为a1,a2,b1,b2.共4个,输入为也为4个分量,共times*N组
in1=x(2:end-1);
in2=x(1:end-2);
in3=u(2:end-1);
in4=u(1:end-2);

P_ori=1e6*eye(4,4);     %充分大的数字乘以单位矩阵
theta_ori=[0 0 0 0];   %初始参数估计值设为0
X_ori=[in1(1);in2(1);in3(1);in4(1)];
Phi_ori=[X_ori];
Err=zeros(1,times*N);
Theta=zeros(4,times*N);

for i=2:times*N-2
    %求Phi_N
    %Phi=zeros(4,i);
    X=[in1(i);in2(i);in3(i);in4(i)];
    Phi=[Phi_ori X];
    
    %求P_{N+1}
    P=P_ori-(P_ori*X*X'*P_ori)/(1+X'*P_ori*X);
    %P=(P_ori^-1+X*X')^-1;
    
    %求K_N
    K=X'*P;
    
    %求\varepsilon_N
    varepsilon=z(i)-theta_ori*X;
    Err(i)=varepsilon^2;
    
    %求新值
    theta=theta_ori+varepsilon*K;
    Theta(:,i)=theta;
    
    %进行下一轮
    Phi_ori=Phi;
    P_ori=P;
    theta_ori=theta;
    X_ori=X;
end
% theta

%绘图部分
ff=@(theta,x1,x0,u1,u0) theta*[x1;x0;u1;u0];
xx=zeros(1,times*N);xx(1:2)=[x0 x1];
zz_ori=zeros(1,times*N);
Ori_theta=[1.5 -0.7 1 0.5];
%利用拟合估计的参数再进行递推求值,得到我们的拟合系统
for i=3:times*N
   xx(i)=ff(theta,xx(i-1),xx(i-2),u(i-1),u(i-2));
   zz_ori(i-2)=xx(i);
end

hold on;plot(t,zz_ori);legend('受噪声的输出','理论输出','仿真结果');

figure('Name','输入');
plot(t,u);

figure('Name','误差值');
plot(t,Err);

figure('Name','参数');hold on;
plot(repmat([0;times*N],1,4),repmat(Ori_theta,2,1),'LineWidth',2);
plot(repmat(t',1,4),repmat(Theta',1,1));

Matlab结果：

可以发现，效果拟合非常好。

---

1. ${\mathbf{P}}_{\mathbf{N}}=({\mathbf{\Phi }}_{\mathbf{N}}{\mathbf{\Phi }}_{\mathbf{N}}^{\mathbf{T}}{)}^{-1}\in {\mathbb{R}}^{n\times n}$² ↩︎

2. ${\mathbf{\Phi }}_N$——$\left[\begin{array}{cccc}{X}_1& X_2& \cdots & X_N\end{array}\right]\in {\mathbb{R}}^{n\times N}$ ↩︎