最优估计准则与方法（4）最小二乘估计(LS)_学习笔记

前言

最小二乘估计（Least Squres Estimation，LSE） 是非常重要的参数(静态)估计方法，是典型的批量处理（Batch Processing） 方法。自1795年诞生于伟大数学家高斯（Gauss）以来，至今200多年来广泛应用于回归分析、曲线拟合和数据建模等领域，一直作为各行业工程界的主要技术手段。本文将详细介绍最优估计领域中的最小二乘估计（LS），即线性最小二乘估计（LLS）。

背景介绍

图1 Gauss与Ceres
最小二乘法（Least Squares Method） 的首次清晰简洁的阐述是由阿德里安·马里·勒让德（Adrien Marie Legendre）于 1805 年发表的，描述为一种用于将线性方程与数据拟合的代数方法1809 年，卡尔·弗里德里希·高斯(Carl Friedrich Gauss)发表了计算天体轨道的方法，他声称自己早在1795年就掌握了最小二乘法。这自然引发与勒让德的优先权之争，但高斯成功超越了前者。通过将最小二乘法与概率原理及正态分布联系起来，提出最小二乘估计，在此过程中发明了正态分布[1]。高斯方法的强大之处在早期的运用中就已有所体现，当时它被用于预测新发现的小行星（矮行星）“凯瑞斯（Ceres）”的未来位置。1801 年 1 月 1 日，意大利天文学家朱塞佩·皮亚齐（Giuseppe Piazzi）发现了凯瑞斯，并能够追踪其轨迹 40 天，之后它便在太阳的强光照射下消失不见了。基于这些数据，天文学家们希望在凯瑞斯从太阳后面出来之前确定其位置，而无需解决开普勒关于行星运动的复杂非线性方程。只有由 24 岁的匈牙利天文学家弗朗茨·亚瑟·冯·扎克（Franz Xaver von Zach）运用高斯提供的最小二乘估计法分析得出的预测，才成功地找到了凯瑞斯的位置[1]。
根据模型参数与观测数据的关系，最小二乘法分为 线性最小二乘法（Linear Least Squares） 和 非线性最小二乘法（Nonlinear Least Squares）。最优估计理论中研究的 最小二乘估计（LS）为线性最小二乘估计（LLS），包括 古典最小二乘估计（LS）、 加权最小二乘估计（WLS）和递推最小二乘估计（RLS）。关于 非线性最小二乘法（NLLS） 是数值优化理论的研究内容，是工程中解决非线性系统问题的主要方法。

线性参数估计问题描述

设$X$为$n$维未知参数向量，$Z$为$k$维观测向量，表示经过$k$组实验观测得到的观测值向量，其中元素$z_i$表示第i次观测实验得到的观测值，显然其是1维观测标量，$V$为$k$维观测噪声向量，其中元素$v_i$表示第i次观测实验的观测噪声，显然其是1维噪声标量。一般情况下$k>n$且希望$k$比$n$大得多。单次观测值为多维的情况将在后续其他Blog中讨论。观测实验依据的自变量为$\theta$，则将观测量$z_i$表示为关于$\theta$的未知函数$f(\theta ,X)$：
$$\begin{array}{r}{z}_i=f(\theta ,X)=\sum _{j=1}^n\left[x_j{h}_{i,j}(\theta )\right]+v_i=x_1{h}_{i,1}(\theta )+x_2{h}_{i,2}(\theta )+\cdots +x_n{h}_{i,n}(\theta )+v_i\end{array}\text{\hspace{1em}(1)}$$
其中
$$\begin{array}{r}X=\left[\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_n\end{array}\right]Z=\left[\begin{array}{c}{z}_1\\ z_2\\ ⋮\\ z_k\end{array}\right]V=\left[\begin{array}{c}{v}_1\\ v_2\\ ⋮\\ v_k\end{array}\right]\end{array}$$
式(1)中$h_{i,j}(\theta )$表示第$i$次观测第$j$个基函数，常用为多项式、三角函数或自然指数函数形式：
$$\begin{array}{rl}{h}_{i,j}(\theta )& ={\theta }^{j-1}\\ h_{i,j}(\theta )& =sin(j\theta )\\ h_{i,j}(\theta )& =exp({\lambda }_j\theta )\end{array}$$
其中，${\lambda }_j$为自然数指数参数。
当观测实验进行，上述基函数均可根据$\theta$求得。令$h_i=\left[\begin{array}{cccc}{h}_{i,1}(\theta )& h_{i,2}(\theta )& \cdots & h_{i,n}(\theta )\end{array}\right]$且为已知，其为$n$维常向量，将式(1)改写为：
$$\begin{array}{r}Z=HX+V\end{array}\text{\hspace{1em}(2)}$$
其中，$H$为参数向量$X$到观测向量$Z$的$k\times n$维转移矩阵：
$$\begin{array}{r}H=\left[\begin{array}{c}{h}_1\\ h_2\\ ⋮\\ h_k\end{array}\right]=\left[\begin{array}{cccc}{h}_{1,1}(\theta )& h_{1,2}(\theta )& \cdots & h_{1,n}(\theta )\\ h_{2,1}(\theta )& h_{2,2}(\theta )& \cdots & h_{2,n}(\theta )\\ ⋮& ⋮& \ddots & ⋮\\ h_{k,1}(\theta )& h_{k,2}(\theta )& \cdots & h_{k,n}(\theta )\end{array}\right]\end{array}$$
显然，观测向量$Z$与被估参数向量$X$存在线性关系，依据最优准则求对$X$的估计值$\hat{X}$是一个线性参数估计问题，自然对应线性最小二乘估计（LLS）。

这里讨论下超定方程组的矛盾：当$k=n$时，线性方程组有唯一精确解，但当$k>n$，线性方程数大于未知被估参数向量的维度，线性方程组变成线性超定方程组，其解不唯一。最小二乘法的思想是需求统计意义上的近似解，使线性超定方程组中各方程能得到近似相等。

古典最小二乘估计（Classical Least Squares Estimation, CLSE）

最小二乘估计（LS）准则 为：使残差平方和达到最小。
代价函数如下：
$$\begin{array}{r}J={\hat{E}}^T\hat{E}=(Z-H\hat{X}{)}^T(Z-H\hat{X})=\sum _{i=1}^k{\hat{e}}_i^2=\sum _{i=1}^k(z_i-h_i\hat{X}{)}^2=min\end{array}\text{\hspace{1em}(3)}$$
其中，${\hat{e}}_i$为第$i$次观测的残差（Residual Error），$\hat{E}$为$k$维残差向量有：
$$\begin{array}{rl}{\hat{e}}_i& =z_i-h_i\hat{X}\\ \hat{E}& =Z-H\hat{X}=\left[\begin{array}{c}{\hat{e}}_1\\ {\hat{e}}_2\\ ⋮\\ {\hat{e}}_k\end{array}\right]\end{array}\text{\hspace{1em}(4)(5)}$$
最小二乘估计（LS）方法
根据式(3)进行对如下代价函数进行最小化：
$$\begin{array}{rl}J& =(Z-H\hat{X}{)}^T(Z-H\hat{X})\end{array}\text{\hspace{1em}(6)}$$
令$J$对$\hat{X}$求偏导，并令其为0，有：
$$\begin{array}{rl}\frac{\partial J}{\partial \hat{X}}& =0\\ \frac{\partial J}{\partial (Z-H\hat{X})}\frac{\partial (Z-H\hat{X})}{\partial \hat{X}}& =0\\ -2H^T(Z-H\hat{X})& =0\\ \hat{X}& =(H^{T}H{)}^{-1}{H}^{T}Z\end{array}\text{\hspace{1em}(7)}$$
再由$\frac{{\partial }^{2}J}{\partial {\hat{X}}^2}=2H^{T}H>0$，为$\hat{X}$为被估参数向量$X$的最小二乘估计，显然其是观测向量$Z$的线性估计。
此时，由式(6)代价函数的最小值为：
$$\begin{array}{rl}{J}_{min}& =(Z-H\hat{X}{)}^T(Z-H\hat{X})\\ & =Z^T(I-H(H^{T}H{)}^{-1}{H}^T{)}^T(I-H(H^{T}H{)}^{-1}{H}^T)Z\end{array}\text{\hspace{1em}(8)}$$
最小二乘估计（LS）无偏性
令估计误差为$\tilde{X}$，定义为被估参数向量$X$与估计值向量$\hat{X}$的偏差，有：
$$\begin{array}{rl}\tilde{X}& =X-\hat{X}\\ & =(H^{T}H{)}^{-1}{H}^{T}HX-(H^{T}H{)}^{-1}{H}^{T}Z\\ & =(H^{T}H{)}^{-1}{H}^T(HX-Z)\\ & =-(H^{T}H{)}^{-1}{H}^{T}V\end{array}\text{\hspace{1em}(9)(10)}$$
估计误差$\tilde{X}$的数学期望为：
$$\begin{array}{rl}E[\tilde{X}]& =E[X-\hat{X}]\\ & =E[-(H^{T}H{)}^{-1}{H}^{T}V]\\ & =-(H^{T}H{)}^{-1}{H}^{T}E[V]\end{array}\text{\hspace{1em}(11)(12)}$$
由式(12)可知，如果观测噪声$V$为白噪声，即$E[V]=0$，则最小二乘估计$\hat{X}$为无偏线性估计。在该无偏估计情况下，估计误差$\hat{X}$的方差矩阵与估计量$\hat{X}$的均方误差矩阵相等，即：
$$\begin{array}{rl}Var(\tilde{X})& =E[(\tilde{X}-E[\tilde{X}])(\tilde{X}-E[\tilde{X}]{)}^T]\\ & =E[\tilde{X}{\tilde{X}}^T]-(E[\tilde{X}]{)}^2\\ & =E[\tilde{X}{\tilde{X}}^T]-0\\ & =E[\tilde{X}{\tilde{X}}^T]\\ MSE(\hat{X})& =E[(X-\hat{X})(X-\hat{X}{)}^T]\\ & =E[\tilde{X}{\tilde{X}}^T]\end{array}\text{\hspace{1em}(13)(14)}$$
方差矩阵描述估计误差向量偏离其均值的离散程度，而估计量均方误差矩阵描述估计误差平方和均值偏离0的离散程度，二者在无偏估计的情况下等价，说明估计误差越接近0，估计精度越高。
$$\begin{array}{rl}Var(\tilde{X})& =MSE(\hat{X})\\ & =E[\tilde{X}{\tilde{X}}^T]\\ & =E[(-(H^{T}H{)}^{-1}{H}^{T}V)(-(H^{T}H{)}^{-1}{H}^{T}V{)}^T]\\ & =(H^{T}H{)}^{-1}{H}^{T}E[V{V}^T]H^T(H^{T}H{)}^{-1}\\ & =(H^{T}H{)}^{-1}{H}^{T}R{H}^T(H^{T}H{)}^{-1}\end{array}\text{\hspace{1em}(15)(16)}$$
其中，$R$为观测噪声向量$V$的方差矩阵：
$$\begin{array}{rl}R& =\left[\begin{array}{cccc}{\sigma }_1^2& 0& \cdots & 0\\ 0& {\sigma }_2^2& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & {\sigma }_k^2\end{array}\right]\end{array}$$
由式(8)和(16)可知，即使在无偏估计前提下，二者并不相等。因此，最小二乘无偏估计只能保证残差平方和最小，但不能保证估计误差方差最小。

综上，根据最小二乘估计原理，做如下总结：

1. 求最小二乘估计量$\hat{X}$不需要任何观测噪声向量$V$的任何统计信息；
2. 最小二乘估计的无偏性取决于噪声向量$V$的数学期望，如$V$为白噪声，即为无偏估计；
3. 无论是否具备无偏性，最小二乘估计只能保证残差平方和最小而不是估计误差方差最小。

参考文献

[1] Least squares
https://en.wikipedia.org/wiki/Least_squares
[2] 《最优估计理论》，周凤歧，2009，高等教育出版社。
[3] 《最优估计理论》，刘胜，张红梅著，2011，科学出版社。