最小二乘法，简明公式整理，数学证明，Matlab实现（自写代码、lsqcurvefit函数、fminsearch函数）

批处理最小二乘方法

  考虑线性输入输出系统，其输出值$y$可由输入状态$\mathbf{x}$线性表示，数学模型描述为：
$$y={\mathbf{x}}^T\theta +\xi$$其中，$\xi$为系统白噪声；$\mathbf{x}$和$\theta$为维度相同的列向量，分别表示为$\mathbf{x}={\left[\begin{array}{cccc}{x}_1& x_2& \cdots & x_n\end{array}\right]}^T$，$\theta ={\left[\begin{array}{cccc}{c}_1& c_2& \cdots & c_n\end{array}\right]}^T$。

  当获取到$k$组输入输出信息，所有这些信息可被用于参数辨识，记作输入矩阵$X\left(k\right)={\left[\begin{array}{cccc}\mathbf{x}\left(1\right)& \mathbf{x}\left(2\right)& \cdots & \mathbf{x}\left(k\right)\end{array}\right]}^T$，输出矩阵$Y\left(k\right)={\left[\begin{array}{cccc}y\left(1\right)& y\left(2\right)& \cdots & y\left(k\right)\end{array}\right]}^T$。参数辨识的本质是对$\theta$进行最优估计$\hat{\theta }$，从而实现：
$$\min J={\left(X\hat{\theta }-Y\right)}^T\left(X\hat{\theta }-Y\right)$$

  对上式求取 $\hat{\theta }$ 的偏导数并设置为零向量，$\frac{\partial J}{\partial \hat{\theta }}=2X^{T}X\hat{\theta }-2X^{T}Y=0$。求解得到全局极值点，如下所示。由于$J$的二阶偏导数$\frac{{\partial }^{2}J}{\partial {\hat{\theta }}^2}=2X^{T}X$为正定矩阵，该解必然为全局最小值。
$$\hat{\theta }\left(k\right)={\left(X^{T}X\right)}^{-1}{X}^{T}Y$$

递推最小二乘方法

   随着$k$的增大，上述批处理最小二乘方法需要占用大量的存储资源，且存在求解效率低、实时性差等缺点。递推形式的最小二乘参数辨识可以有效解决这些问题，$\hat{\theta }\left(k\right)$的递归形式可以被表示为：
$$\{\begin{array}{l}\hat{\theta }\left(k\right)=\hat{\theta }\left(k-1\right)+K\left(k\right)\left[y\left(k\right)-{\mathbf{x}}^T\left(k\right)\hat{\theta }\left(k-1\right)\right]\\ K\left(k\right)=\frac{P\left(k-1\right)\mathbf{x}\left(k\right)}{1+{\mathbf{x}}^T\left(k\right)P\left(k-1\right)\mathbf{x}\left(k\right)}\\ P(k)=\left[I-K\left(k\right){\mathbf{x}}^T\left(k\right)\right]P\left(k-1\right)\end{array}$$其中，$K\left(k\right)$为增益向量；$P\left(k\right)={\left(X^{T}X\right)}^{-1}$表示误差的协方差矩阵，其初始值可由${\left(X^{T}X\right)}^{-1}$给出或设置为$\alpha I$，$I$为单位矩阵，$\alpha =10^6\sim 10^{10}$。

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

  递推最小二乘参数辨识方法可以有效缓解存储与计算压力，但依然存在着数据饱和问题。这是最小二乘算法本身属性所决定的，算法中的所有数据被同等对待，这也就意味着新数据所带来的更新效果会被逐渐削弱。这种内在属性使得一般最小二乘算法难以收敛到参数真实值，且无法跟踪运行过程中的参数变化。为解决该问题，可以引入具有遗忘因子的最小二乘参数辨识方法，对数据的时序进行考虑，并将原始优化问题修改为：
$$\min J={\left(X\hat{\theta }-Y\right)}^{T}W\left(X\hat{\theta }-Y\right)$$其中，$W=diag\left({\lambda }^{k-1},{\lambda }^{k-2},...,\lambda ,1\right)$为加权对角矩阵，$\lambda$为遗忘因子且$0.98<\lambda <1$。

  很容易理解，历史观测数据距离当前时刻越远，其所的占权重也会越低，参数辨识的过程会更依赖于近期加入的数据。对上式进行求解，可以得到具有遗忘因子的递推最小二乘参数辨识方法：
$$\{\begin{array}{l}\hat{\theta }\left(k\right)=\hat{\theta }\left(k-1\right)+K\left(k\right)\left[y\left(k\right)-{\mathbf{x}}^T\left(k\right)\hat{\theta }\left(k-1\right)\right]\\ K\left(k\right)=\frac{P\left(k-1\right)\mathbf{x}\left(k\right)}{\lambda +{\mathbf{x}}^T\left(k\right)P\left(k-1\right)\mathbf{x}\left(k\right)}\\ P\left(k\right)=\frac{1}{\lambda }\left[I-K\left(k\right){\mathbf{x}}^T\left(k\right)\right]P\left(k-1\right)\end{array}$$

其中，协方差矩阵$P\left(k\right)={\left(X^{T}WX\right)}^{-1}$，其初始值的选取原则与一般递归最小而成方法保持一致。

Matlab案例分析

自写代码

  假设系统输出$y$与时间$t$满足如下关系表达式，并已知$t=1:100$期间的$y$值，求取系数$c_1,c_2,c_3$的最优表达式：
$$y=c_{1}t+c_2{\left(t-60\right)}^2+c_3\cos \left(\pi t/10\right)+\xi$$
  通过定义$x_1=t$，$x_2={\left(t-60\right)}^2$，$x_3=\cos \left(\pi t/10\right)$，我们完全可以把上述表达式转换为一个标准的线性表达式：
$$y=c_1{x}_1+c_2{x}_2+c_3{x}_3+\xi$$
  于是，我们可以通过递推最小二乘方法对系数进行求解，MATLAB代码如下：

clc; clear

% 构建三个输入项，一共101组
t = 0 : 1 : 100;
x1 = t;
x2 = (t - 60).^2;
x3 = cos(pi.*t / 10);

% 假设系统真实系数
c1 = 2.22;
c2 = 0.05;
c3 = 23.1;

% 得到系统输出，得到对应的101个输出（人为加入了一些噪声）
Y = c1 * x1 + c2 * x2 + c3 * x3  + randn(1,101);

% 绘制一下采样数据
plot(t, Y,'o')
hold on

% 设置初始值，执行递推过程
theta = [0;0;0];
Pk_ = 1e6 * eye(3);
lambda = 0.995;
for i = 1 : 1 : 100
    x = [x1(i); x2(i); x3(i)];
    y = Y(i);
    Kk = Pk_ * x / (lambda + x' * Pk_ * x);
    theta = theta + Kk * (y - x'*theta);
    Pk_ = (1/lambda)*(eye(3) - Kk * x') * Pk_;
end

% 基于最小二乘的结果绘制曲线，并于原始数据做对比
Y_Est = theta(1) * x1 + theta(2) * x2 + theta(3) * x3;

% 绘制出来拟合的曲线，并于采样点进行比较
plot(t, Y,'r')

  最终估算的参数向量theta = [2.2198; 0.0499; 23.0333]，与系统真实参数$c1=2.22$，$c2=0.05$，$c3=23.1$是近似相等的。注意：由于随机噪声的存在。每次跑的结果可能会存在微小的差异。基于估算系数进行曲线拟合，并对比采样数据，结果如下：

matlab之lsqcurvefit函数

  如果只是简单应用，是可以直接调用matlab库函数的。lsqcurvefit函数基于最小二乘方法，进行曲线拟合，该函数内部进行了很多的优化以满足对非线性拟合的需要。

clc; clear;

% 生成一组模拟数据
x = linspace(-10, 10, 100);
y = 2*x.^2 - 3*sin(x) + 0.5*randn(size(x));

% 定义要拟合的函数形式，包括二次项和三角函数项
f = @(p,x) p(1)*x.^2 + p(2)*sin(p(3)*x) + p(4);

% 定义初始参数矩阵，其中包含二次项系数、三角函数的振幅、频率和常数项
p0 = [1, 1, 1, 1]; % 对应真实值为[2, -3, 1, 0]

% 使用最小二乘法拟合曲线，并返回最优参数和拟合误差
[p, resnorm] = lsqcurvefit(f, p0, x, y);

% 生成拟合曲线并绘制
y_fit = f(p, x);
plot(x, y, 'o', x, y_fit, '-')
legend('原始数据', '拟合曲线')

  lsqcurvefit的输入有四项，分别为待拟合函数、待拟合函数参数向量初始值、待拟合函数输入向量、采样结果，且待拟合函数定义的时候，必须将$x$也作为输入的一部分，仿真结果如下所示。感觉用起来还挺费劲的，相比较而言，更推荐下面的fminsearch函数，比较直观一些。

matlab之fminsearch函数

  相比较而言，fminsearch函数更为直观一些，因为其输入只需要包含两项：目标函数、参数向量初始值。而目标函数直接对所有误差进行了包含，只保留了参数向量作为输入项，或者说待优化项。代码如下所示：

clc; clear;

% Generate simulated data
x = linspace(0, 2*pi, 50)';
y = 2*sin(2*x) + 0.5*x.^2 + randn(size(x));

% Define error function
fun = @(c) sum((y - (c(1)*sin(c(2)*x) + c(3)*x.^2 + c(4))).^2);

% Find least-squares solution
c0 = [1 2 1 1]; % 对应真实值为[2 2 0.5 0]
c = fminsearch(fun, c0);

% Plot results
xfit = linspace(0, 2*pi, 100)';
yfit = c(1)*sin(c(2)*xfit) + c(3)*xfit.^2 + c(4);
plot(x, y, 'ko', xfit, yfit, 'b-');
legend('Data', 'Fit');

  仿真结果如下所示：

附录1：递推最小二乘的数学证明

$$\{\begin{array}{rl}\hat{\theta }\left(k+1\right)& ={\left[X_{k+1}^{\mathrm{T}}{X}_{k+1}\right]}^{-1}{X}_{k+1}^{\mathrm{T}}{Y}_{k+1}\\ & ={\left[X_k^{\mathrm{T}}{X}_k+\mathbf{x}\left(k+1\right){\mathbf{x}}^{\mathrm{T}}\left(k+1\right)\right]}^{-1}\left[\mathbf{x}\left(k+1\right)y\left(k+1\right)+X_k^{\mathrm{T}}{Y}_k\right]，见附录1.1\\ & =\left[I-K(k+1){\mathbf{x}}^{\mathrm{T}}(k+1)\right]P(k)\left[\mathbf{x}\left(k+1\right)y\left(k+1\right)+X_k^{\mathrm{T}}{Y}_k\right]，见附录1.2\\ & =\left[I-K(k+1){\mathbf{x}}^{\mathrm{T}}(k+1)\right]\left[P\left(k\right)\mathbf{x}\left(k+1\right)y\left(k+1\right)+\hat{\theta }\left(k\right)\right]\\ & =\left[I-K\left(k+1\right){\mathbf{x}}^{\mathrm{T}}(k+1)\right]\hat{\theta }\left(k\right)+\left[I-K(k+1){\mathbf{x}}^{\mathrm{T}}(k+1)\right]P\left(k\right)\mathbf{x}\left(k+1\right)y\left(k+1\right)\\ & =\left[I-K(k+1){\mathbf{x}}^{\mathrm{T}}(k+1)\right]\hat{\theta }\left(k\right)+\left[I-\frac{P(k)\mathbf{x}(k+1)}{1+{\mathbf{x}}^{\mathrm{T}}(k+1)P(k)\mathbf{x}(k+1)}{\mathbf{x}}^{\mathrm{T}}(k+1)\right]P\left(k\right)\mathbf{x}\left(k+1\right)y\left(k+1\right)\\ & =\left[I-K(k+1){\mathbf{x}}^{\mathrm{T}}(k+1)\right]\hat{\theta }\left(k\right)+\frac{P\left(k\right)\mathbf{x}\left(k+1\right)}{1+{\mathbf{x}}^{\mathrm{T}}(k+1)P(k)\mathbf{x}(k+1)}y\left(k+1\right)\left[1+{\mathbf{x}}^{\mathrm{T}}(k+1)P(k)\mathbf{x}(k+1)\right]-\frac{P(k)\mathbf{x}(k+1)}{1+{\mathbf{x}}^{\mathrm{T}}(k+1)P(k)\mathbf{x}(k+1)}{\mathbf{x}}^{\mathrm{T}}\left(k+1\right)P\left(k\right)\mathbf{x}\left(k+1\right)y\left(k+1\right)\\ & =\left[I-K(k+1){\mathbf{x}}^{\mathrm{T}}(k+1)\right]\hat{\theta }\left(k\right)+K\left(k+1\right)y\left(k+1\right)\left[1+{\mathbf{x}}^{\mathrm{T}}(k+1)P(k)\mathbf{x}(k+1)\right]-K\left(k+1\right){\mathbf{x}}^{\mathrm{T}}\left(k+1\right)P\left(k\right)\mathbf{x}\left(k+1\right)y\left(k+1\right)\\ & =\hat{\theta }\left(k\right)-K(k+1){\mathbf{x}}^{\mathrm{T}}(k+1)\hat{\theta }\left(k\right)+K\left(k+1\right)y\left(k+1\right)+K\left(k+1\right){\left[{\mathbf{x}}^{\mathrm{T}}(k+1)P(k)\mathbf{x}(k+1)\right]}_{1\times 1}y\left(k+1\right)-K\left(k+1\right){\mathbf{x}}^{\mathrm{T}}\left(k+1\right)P\left(k\right)\mathbf{x}\left(k+1\right)y\left(k+1\right)\\ & =\hat{\theta }\left(k\right)-K(k+1){\mathbf{x}}^{\mathrm{T}}(k+1)\hat{\theta }\left(k\right)+K\left(k+1\right)y\left(k+1\right)\\ & =\hat{\theta }\left(k\right)+K(k+1)\left[y\left(k+1\right)-{\mathbf{x}}^{\mathrm{T}}(k+1)\hat{\theta }\left(k\right)\right]\end{array}$$

附录1.1：数学证明步骤一

  根据批处理最小二乘方法：
$$\hat{\theta }\left(k\right)={\left(X_k^T{X}_k\right)}^{-1}{X}_k^T{Y}_k$$  进一步，在下一刻的估计可以表示为：
$$\hat{\theta }\left(k+1\right)={\left[X_{k+1}^{\mathrm{T}}{X}_{k+1}\right]}^{-1}{X}_{k+1}^{\mathrm{T}}{Y}_{k+1}$$其中，
$$X_{k+1}=\left[\begin{array}{c}{X}_k\\ {\mathbf{x}}^{\mathrm{T}}\left(k+1\right)\end{array}\right],Y_{k+1}=\left[\begin{array}{c}{Y}_k\\ y\left(k+1\right)\end{array}\right]$$  于是有
$$\{\begin{array}{rl}\hat{\theta }\left(k+1\right)& ={\left[X_{k+1}^{\mathrm{T}}{X}_{k+1}\right]}^{-1}{X}_{k+1}^{\mathrm{T}}{Y}_{k+1}\\ & ={\left[X_k^{\mathrm{T}}{X}_k+\mathbf{x}\left(k+1\right){\mathbf{x}}^{\mathrm{T}}\left(k+1\right)\right]}^{-1}\left[\mathbf{x}\left(k+1\right)y\left(k+1\right)+X_k^{\mathrm{T}}{Y}_k\right]\end{array}$$

附录1.2：数学证明步骤二

  为了证明这一步，首先介绍一个数学引理，这个稍后会用到。
【引理】设$A$，$C$和$A+BCD$均为非奇异方阵，则有
$$\{\begin{array}{l}(A+BCD{)}^{-1}=A^{-1}-A^{-1}B{\left(C^{-1}+D{A}^{-1}B\right)}^{-1}D{A}^{-1}\\ 特例(C=I)，(A+BD{)}^{-1}=A^{-1}-A^{-1}B{\left(I+D{A}^{-1}B\right)}^{-1}D{A}^{-1}\end{array}$$  该引理的证明过程如下：
$$\{\begin{array}{l}(A+BCD)\left[A^{-1}-A^{-1}B{\left(C^{-1}+D{A}^{-1}B\right)}^{-1}D{A}^{-1}\right]\\ =I+BCD{A}^{-1}-B\left(C^{-1}+D{A}^{-1}B\right)D{A}^{-1}-BCD{A}^{-1}B{\left(C^{-1}+D{A}^{-1}B\right)}^{-1}D{A}^{-1}\\ =I+BCD{A}^{-1}-B\left(I+CD{A}^{-1}B\right){\left(C^{-1}+D{A}^{-1}B\right)}^{-1}D{A}^{-1}\\ =I+BCD{A}^{-1}-BC\left(C^{-1}+D{A}^{-1}B\right){\left(C^{-1}+D{A}^{-1}B\right)}^{-1}D{A}^{-1}\\ =I\end{array}$$
  上面已经证明到：
$$\{\begin{array}{rl}\hat{\theta }\left(k+1\right)& ={\left[X_{k+1}^{\mathrm{T}}{X}_{k+1}\right]}^{-1}{X}_{k+1}^{\mathrm{T}}{Y}_{k+1}\\ & ={\left[X_k^{\mathrm{T}}{X}_k+\mathbf{x}\left(k+1\right){\mathbf{x}}^{\mathrm{T}}\left(k+1\right)\right]}^{-1}\left[\mathbf{x}\left(k+1\right)y\left(k+1\right)+X_k^{\mathrm{T}}{Y}_k\right]\end{array}$$对其中的求逆部分进行变换：
$$\left[\begin{array}{c}(A+BD{)}^{-1}=A^{-1}-A^{-1}B{\left(I+D{A}^{-1}B\right)}^{-1}D{A}^{-1}\\ \Downarrow \\ {\left[X^{\mathrm{T}}X+\mathbf{x}(k+1){\mathbf{x}}^{\mathrm{T}}(k+1)\right]}^{-1}={\left(X^{\mathrm{T}}X\right)}^{-1}-{\left(X^{\mathrm{T}}X\right)}^{-1}\mathbf{x}\left(k+1\right){\left[I+{\mathbf{x}}^{\mathrm{T}}\left(k+1\right){\left(X^{\mathrm{T}}X\right)}^{-1}\mathbf{x}\left(k+1\right)\right]}^{-1}{\mathbf{x}}^{\mathrm{T}}\left(k+1\right){\left(X^{\mathrm{T}}X\right)}^{-1}\end{array}\right]$$

  定义$P\left(k\right)={\left(X_k^{\mathrm{T}}{X}_k\right)}^{-1}$，由于${\mathbf{x}}^{\mathrm{T}}\left(k+1\right){\left(X^{\mathrm{T}}X\right)}^{-1}\mathbf{x}\left(k+1\right)$为1×1的标量，基于上式可以得到：
$$\{\begin{array}{rl}P\left(k+1\right)& ={\left[X_{k+1}^{\mathrm{T}}{X}_{k+1}\right]}^{-1}\\ P\left(k+1\right)& ={\left[X^{\mathrm{T}}X+\mathbf{x}(k+1){\mathbf{x}}^{\mathrm{T}}(k+1)\right]}^{-1}\\ & ={\left(X^{\mathrm{T}}X\right)}^{-1}-{\left(X^{\mathrm{T}}X\right)}^{-1}\mathbf{x}\left(k+1\right){\left[1+{\mathbf{x}}^{\mathrm{T}}\left(k+1\right){\left(X^{\mathrm{T}}X\right)}^{-1}\mathbf{x}\left(k+1\right)\right]}^{-1}{\mathbf{x}}^{\mathrm{T}}\left(k+1\right){\left(X^{\mathrm{T}}X\right)}^{-1}\\ & =P\left(k\right)-P\left(k\right)\frac{\mathbf{x}(k+1){\mathbf{x}}^{\mathrm{T}}(k+1)}{1+{\mathbf{x}}^{\mathrm{T}}\left(k+1\right)P\left(k\right)\mathbf{x}\left(k+1\right)}P\left(k\right)\\ & =\left[I-K(k+1){\mathbf{x}}^{\mathrm{T}}(k+1)\right]P(k)，递归的第三个方程\end{array}$$其中，$K(k+1)$定义表达式为：
$$K(k+1)=\frac{P(k)\mathbf{x}(k+1)}{1+{\mathbf{x}}^{\mathrm{T}}(k+1)P(k)\mathbf{x}(k+1)}，递归的第二个方程$$  于是有：
$$\{\begin{array}{rl}\hat{\theta }\left(k+1\right)& ={\left[X_{k+1}^{\mathrm{T}}{X}_{k+1}\right]}^{-1}{X}_{k+1}^{\mathrm{T}}{Y}_{k+1}\\ & ={\left[X_k^{\mathrm{T}}{X}_k+\mathbf{x}\left(k+1\right){\mathbf{x}}^{\mathrm{T}}\left(k+1\right)\right]}^{-1}\left[\mathbf{x}\left(k+1\right)y\left(k+1\right)+X_k^{\mathrm{T}}{Y}_k\right]\\ & =\left[I-K(k+1){\mathbf{x}}^{\mathrm{T}}(k+1)\right]P(k)\left[\mathbf{x}\left(k+1\right)y\left(k+1\right)+X_k^{\mathrm{T}}{Y}_k\right]\end{array}$$