最小二乘法溯源追踪

假设我们有四个观测点\((1,6),(2,5),(3,7),(4,10)\)，我们希望找到一个线性 函数\(y = \beta_{1} + \beta_{2}x\) 拟合这四个点。也就是说，我们希望找 到\(\beta_{1},\beta_{2}\) 适配这个系统。另外，我们在这里也假定这个系统 是线性系统，则有：

\begin{eqnarray} \label{eq:2} 6&=& \beta_{1} + \beta_{2} \\ 5&=& \beta_{1} + 2\beta_{2} \\ 7&=& \beta_{1} + 3\beta_{3} \\ 10&=& \beta_{1} + 4\beta_{3} \end{eqnarray}

以上有四个方程，两个未知数。如果系数矩阵的秩大于2，那么这个方程组就是 过定方程组。

最小二乘的目的在于找到\(\beta_{1},\beta_{2}\)使得式 (\ref{eq:3})最小：

\begin{eqnarray} \label{eq:3} S(\beta_{1},\beta_{2}) &=& [6 - (\beta_{1} + \beta_{2})]^{2} \\ &+& [5 - (\beta_{1} + 2\beta_{2})]^{2} \\ &+& [7 - ( \beta_{1} + 3\beta_{3})]^{2} \\ &+& [10 - ( \beta_{1} + 4\beta_{3} )]^{2} \end{eqnarray}

\(S(\beta_{1},\beta_{2})\)对\(\beta_{1},\beta_{2}\)求导，并令其等于零。得到：

\begin{equation} \label{eq:4} \beta_{1} = 3.5 \end{equation} \begin{equation} \label{eq:5} \beta_{2} = 1.4 \end{equation}

最终的拟合结果为：\( y = 3.5 + 1.4x \)。

考虑系统：

\begin{equation} \label{eq:6} \sum_{j = 1}^{n}A_{ij}\beta_{j} = y_{i}, i=1,2,\ldots,m \end{equation}

这个系统有\(m\)个观测值\(y_{i},i=1,\ldots,m\)，\(n\)个未知系数 \(\beta_{1},\ldots,\beta_{n}\)。比如一个房子的属性有面积，单价，卧室数 量，车位数量四个属性，根据这四个属性我们会对房子做一个估价\(y\)。我们 想要拟合一个房价和这四个属性之间的线性关系。对应到式~(\ref{eq:6})

\begin{equation} \label{eq:7} y = \beta_{1}H_{1} + \beta_{2}H_{2} + \beta_{3}H_{3} + \beta_{4}H_{4} \end{equation}

其中\(H_{1}\)代表面积，\(H_{2}\)代表单价，\(H_{3}\)代表卧室数量， \(H_{4}\)代表车位数量，\(\beta_{1},\ldots,\beta_{4}\)代表这四个属性的 权重，\(y\)代表房价。这里我们假定这四个因素和总价\(y\)之间只存在线性关 系。为了得到\(\beta_{1},\ldots,\beta_{4}\)我们需要一些房子的数据表，根 据这个表，我们来拟合这个关系。

对应到式 (\ref{eq:6}) ，我们有\(m\)个观测，每个观测涉及\(n\)个属性。把 式~(\ref{eq:6})写成矩阵的形式有：

\begin{equation} \label{eq:8} \mathbf{Y} = \mathbf{A} \mathbf{\beta} \end{equation}

其中：

\begin{equation*} \mathbf{Y} = [y_{1}, y_{2}, \ldots, y_{m} ] \end{equation*}

\(\mathbf{Y}\)是\(m\times 1\)的列向量。

\begin{equation*} \mathbf{A} = \begin{bmatrix} A_{11} & A_{12} & \ldots & A_{1n} \\ A_{21} & A_{22} & \ldots & A_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ A_{m1} & A_{m2} & \ldots & A_{mn} \end{bmatrix} \end{equation*}

\(\mathbf{A}\)是\(m\times n\)的矩阵，\(m\)代表一共有\(m\)个观测数据， \(n\) 代表一共有\(n\)个属性。

\begin{equation*} \mathbf{\beta} = [\beta_{1}, \beta_{2},\ldots, \beta_{n}] \end{equation*}

\(\mathbf{\beta}\)是\(n\times 1\)的列向量，\(n\)代表\(n\)个属性。 \(\mathbf{\beta}\) 向量代表了\(n\)个属性的权重。

通常，方程~(\ref{eq:8}) 没有解。我们的目的是找到一个最合适的 \(\mathbf{\beta}\) 使得\(\mathbf{Y}\)与\(\mathbf{A}\mathbf{\beta}\)的 差值的平方最小。

\begin{equation} \label{eq:10} \arg\min_{\mathbf{\beta}} S(\mathbf{Y},\mathbf{A}\mathbf{\beta}) \end{equation}

其中：

\begin{equation} \label{eq:11} S(\mathbf{Y},\mathbf{A}\mathbf{\beta}) = \sum_{i=1}^{m} | y_{i} - \sum_{j=1}^{n}A_{ij} \beta_{j} |^{2} = \| \mathbf{Y} - \mathbf{A} \mathbf{\beta} \|^{2} \end{equation}

假定\(\mathbf{A}\)的各列是相互独立的，也就是说\(n\)个属性互不相关，通 过正则化方程~(\ref{eq:8})：

\begin{equation} \label{eq:12} \mathbf{A}^{T}\mathbf{Y} = \mathbf{A}^{T}\mathbf{A}\mathbf{\beta} \end{equation}

\(\mathbf{A}^{T}\mathbf{A}\)是一个半正定的矩阵，其逆一定存在（假设 \(\mathbf{A}\) 的各列是相互独立的）。最终，有：

\begin{equation} \label{eq:13} \mathbf{\beta} = (\mathbf{A}^{T}\mathbf{A})^{-1} \mathbf{A}^{T} \mathbf{Y} \end{equation}

式~(\ref{eq:13})就是最小二成法的解。

在上一节，我们不加说明的给出：式~(\ref{eq:10})的解就是~(\ref{eq:12})的 解。现在，我们给出式 (\ref{eq:12})的由来。

首先定义拟合值和观测值的差：

\begin{equation} \label{eq:14} r_{i} = y_{i} - \sum_{j=1}^{n}A_{ij}\beta_{j} \end{equation}

那么：

\begin{equation} \label{eq:15} S = \sum_{i=1}^{m} r_{i}^{2} \end{equation}

为了求得\(S\)最小时的\(\mathbf{\beta}\)，我们需要依\(S\)对 \(\mathbf{\beta}\)求偏导。

\begin{equation} \label{eq:16} \frac{\partial S}{\partial \beta_{j}} = 2\sum_{i=1}^{m}r_{i} \frac{\partial r_{i}}{\partial \beta_{j}}, j = 1,2,\ldots,n \end{equation}

另外：

\begin{equation} \label{eq:17} \frac{\partial r_{i}}{\beta_{j}} = -H_{ij} \end{equation}

综上有：

\begin{equation} \label{eq:18} \frac{\partial S}{\partial \beta_{j}} = 2 \sum_{i=1}^{m}( y_{i} - \sum_{k=1}^{n} A_{ik}\beta_{k})(-A_{ij}) \end{equation}

令：

\begin{equation} \label{eq:19} 2 \sum_{i=1}^{m}( y_{i} - \sum_{k=1}^{n} A_{ik}\beta_{k})(-A_{ij}) =0 \end{equation}

重新组织式 (\ref{eq:19})，有：

\begin{equation} \label{eq:20} \sum_{i=1}^{m}\sum_{k=1}^{n}A_{ij}A_{ik}\beta_{k} = \sum_{i=1}^{m}A_{ij}y_{i} , j = 1,2,\ldots,m \end{equation}

写成矩阵的形式有：

\begin{equation} \label{eq:21} (\mathbf{A}^{T}\mathbf{A}) \mathbf{\beta} = \mathbf{A}^{T} \mathbf{Y} \end{equation}

啊哈，现在我们总算把最小二成准则和对应的方程联系起来了。锦上添花，我们 再给出一种获取式 (\ref{eq:21})的方式。这次我们以矢量的方式给出来。

定义：

\begin{equation} \label{eq:22} S(\mathbf{Y},\mathbf{\beta}) = \|\mathbf{Y} - \mathbf{A}\mathbf{\beta}\|^{2} \end{equation}

展开有：

\begin{equation} \label{eq:23} (\mathbf{Y} - \mathbf{A}\mathbf{\beta})^{T} (\mathbf{Y} - \mathbf{A}\mathbf{\beta}) = \mathbf{Y}^{T}\mathbf{Y} - \mathbf{\beta}^{T}\mathbf{A}^{T}\mathbf{Y}- \mathbf{Y}^{T}\mathbf{A}\mathbf{\beta} + \beta^{T}\mathbf{A}^{T}\mathbf{A}\mathbf{\beta} \end{equation}

对式 (\ref{eq:23}) 就\(\beta\)求导，有：

\begin{equation} \label{eq:24} \mathbf{A}^{T} \mathbf{Y} = (\mathbf{A}^{T}\mathbf{A}) \mathbf{A} \beta \end{equation}

假设有一个电阻，但是阻值没有标明，我们需要通过万用表测量以得到正确的阻 值。但是，万用表也是有误差的，每次测量都会存在误差。假设我们测量了 \(k\)次，则有：

\begin{eqnarray} \label{eq:26} y_{1}&=& x + n_{1} \\ y_{2}&=& x + n_{2} \\ \vdots & & \vdots \\ y_{k}&=& x + n_{k} \\ \end{eqnarray}

矩阵形式有：

\begin{equation} \label{eq:27} \begin{bmatrix} y_{1} \\ y_{2} \\ \vdots \\ y_{k} \end{bmatrix} = \begin{bmatrix} 1 \\ 1 \\ \vdots \\ 1 \end{bmatrix} x + \begin{bmatrix} n_{1} \\ n_{2} \\ \vdots \\ n_{k} \end{bmatrix} \end{equation}

根据\(\hat{x} = (A^{T}A)^{-1} A^{T}y \)，则：

\begin{equation} \label{eq:28} \hat{x} = \frac{1}{k} \sum_{i=1}^{k}y_{i} \end{equation}

可以看出来，在这个例子中，最小二乘法结果与多次测量求均值的直观感觉一 致。更进一步：

\begin{equation} \label{eq:29} \hat{x}= \frac{1}{k}\sum_{i=1}^{k} y_{i} = x +\frac{1}{k}\sum_{i}^{k}n_{i} \end{equation}

在通信系统里\(y = x + n\)是AWGN信道的模型。多测测量等效于多次发送同一 信息，然后经过不同的噪声污染。多次发送求平均之后，信号功率没有变化，噪 声功率却变为原来的 \(\frac{1}{n}\)。噪声是误差的来源。在统计学中，误 差不可能被消除，只能被降低。多次测量求均值的过程就是降低误差的过程，即 降低噪声的过程。在通信系统中，噪声不能被消除，只能被抑制。无论在通信系 统还是统计学领域，其背后的数学原理都是一样的。

注意在这里我们的目标是获取\(x\)的最优估计，\(x\)是我们要估计的量。之前 我们给出来的曲线拟合的例子里，我们要获取的是多项式的系数。虽然场景不同， 但是其背后的数学原理相同。

之所以会出现加权的最小二乘法，是因为并不是每次的观测可靠度都是一样的。 比如之前的测量电阻阻值的例子，有可能有些测量是用昂贵的仪器测量的，有些 测量是用便宜的仪器测量的。这两类仪器的测量结果当然不同，我们不能对其平 等看待。所以就要对那些可靠的结果加上较大的权重，对那些稍微不可靠的结果 加上较小的权重。

我们先把这个电阻的例子放在一边。考虑加权的最小二乘法的数学描述。假设 \(x\)是一个矢量，每个元素都是常数。系统方程为：

\begin{equation} \label{eq:30} \begin{bmatrix} y_{1} \\ y_{2} \\ \ldots \\ y_{k} \end{bmatrix} = \begin{bmatrix} H_{11} & H_{12} & \ldots & H_{1n} \\ H_{21} & H_{22} & \ldots & H_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ H_{k1} & H_{k2} & \ldots & H_{kn} \end{bmatrix} \begin{bmatrix} x_{1} \\ x_{2} \\ \ldots \\ x_{n} \end{bmatrix} + \begin{bmatrix} n_{1} \\ n_{2} \\ \ldots \\ n_{k} \end{bmatrix} \end{equation}

且\(E(n_{i}^{2}) = \sigma_{i}^{2}\)， 即\(n_{i}\)是均值为零，方差为 \(\sigma_{i}^{2}\)的高斯噪声。假设\(n_{i},i=1,\ldots,k\)是相互独立的， 则有：

\begin{equation} \label{eq:32} R = E[\mathbf{n}\mathbf{n}^{T}] = \begin{bmatrix} \sigma_{1}^{2} & \ldots & 0 \\ \vdots & & \vdots \\ 0 &\ldots & \sigma_{k}^{2} \end{bmatrix} \end{equation}

不同的\(\sigma_{i}^{2}\)意味着其对应的测量结果的权重不同。 \(\sigma_{i}^{2}\)越小，则权重越大；反之，越小。

那么现在我们的损失函数（优化目标）是：

\begin{equation} \label{eq:34} J = \epsilon_{y_{1}}^{2}/ \sigma_{1}^{2} + \ldots + \epsilon_{y_{k}}^{2} / \sigma_{k}^{2} \end{equation}

其中\(\epsilon_{y_{i}} = y_{i} - \sum_{j=1}^{n} H_{ij}x_{j}\). 所以 \(J\)又可以被写为：

\begin{eqnarray} \label{eq:36} J&=& \epsilon_{y}^{T}R^{-1}\epsilon_{y} \\ &=& (y-H\hat{x})^{T} R^{-1} (y-H \hat{x}) \\ &=& (y^{T}R^{-1} - \hat{x}^{T}H^{T}R^{-1})(y-H \hat{x}) \\ &=& y^{T}R^{-1}y - y^{T}R^{-1} H \hat{x} - x^{T}H^{T}R^{-1} y + \hat{x}^{T} H^{T}R^{-1} H \hat{x} \end{eqnarray}

上式以\(\hat{x}\)求微分，令其等于零，有：

\begin{equation} \label{eq:37} \frac{\partial J}{\partial \hat{x}} = - y^{T}R^{-1}H + \hat{x}^{T}H^{T}R^{-1} H= 0 \end{equation}

进而，有：

\begin{equation} \label{eq:38} \hat{x} = (H^{T}R^{-1}H)^{-1} H^{T}R^{-1}y \end{equation}

注意这个解要求\(R^{-1}\)存在，也就是说每一次测量都要求\(y_{i}\)受点噪 声干扰。但是，也可以看出来，如果这个干扰很小的话，会出现\(R^{-1}\)对角 线上的值非常大，在实际应用中会出现溢出现象。

现在让我们回到测量电阻的那个例子。系统方程仍然不变：

\begin{equation} \label{eq:39} \begin{bmatrix} y_{1} \\ y_{2} \\ \vdots \\ y_{k} \end{bmatrix} = \begin{bmatrix} 1 \\ 1 \\ \vdots \\ 1 \end{bmatrix} x + \begin{bmatrix} n_{1} \\ n_{2} \\ \vdots \\ n_{k} \end{bmatrix} \end{equation}

考虑噪声的协方差矩阵：

\begin{equation} \label{eq:40} R = \mathrm{diag}(\sigma_{1}^{2},\ldots,\sigma_{k}^{2}) \end{equation}

根据 式 (\ref{eq:38})，有：

\begin{eqnarray} \label{eq:41} \hat{x} &=& (H^{T}R^{-1} H )^{-1} H^{T}R^{-1} y \\ &=& ( \sum_{i=1}^{k} 1/\sigma_{i}^{2} )^{-1} \sum_{j=1}^{k} \frac{y_{j}}{\sigma_{j}^{2}} \end{eqnarray}

之前，我们对最小二乘法的原理和正则方程进行了推导。在之前推导最小二乘 的过程中，我们假定不会有新的观测数据，所以系数矩阵都是固定大小的。这 样的方法有个问题：如果有新的数据到来，那么之前的计算需要重新来过。

在本节，我们探讨如何迭代的更新\(\hat{x}\)，即：当我们已经根据\(k-1\)个 估计数据得到了\(\hat{x}\)的一个估计，现在又来了一个数据\(y_{k}\)，如 何不用重新计算 (\ref{eq:38})更新\(\hat{x}\).

一个线性迭代估计可以表示为：

\begin{eqnarray} \label{eq:42} y_{k}&=&H_{k}x + n_{k} \\ \hat{x}_{k} &=& \hat{x}_{k-1} + K_{k}(y_{k} - H_{k}\hat{x}_{k-1}) \end{eqnarray}

我们根据\(\hat{x}_{k-1}\) 和\(y_{k}\)更新\(\hat{x}_{k}\)。\(K_{k}\)是 估计算子的增益矩阵，将在后文推导其计算过程。\(y_{k} - H_{k}\hat{x}_{k-1}\)是修正项。如果\(K_{k}=0\) 或者修正项为零，则 \(\hat{x}_{k} = \hat{x}_{k-1}\)。在计算\(K_{k}\)之前，我们首先考虑迭代 估计的误差期望：

\begin{eqnarray} \label{eq:43} E[\epsilon_{x,k}] &=& E[x - \hat{x}_{k}] \\ &=&E[x - \hat{x}_{k-1} - K_{k}(y_{k} - H_{k}\hat{x}_{k-1}) ] \\ &=&E[x - \hat{x}_{k-1} - K_{k}(H_{k}x + n_{k} - H_{k}\hat{x}_{k-1}) ] \\ &=&E[\epsilon_{x,k-1} - K_{k}H_{k} (x - \hat{x}_{k-1}) - K_{k}n_{k}] \\ &=&E[(\epsilon_{x,k-1} - K_{k}H_{k}\epsilon_{x,k-1}) - K_{k}n_{k}] \\ &=&(I- K_{k}H_{k})E[\epsilon_{x,k-1}] - K_{k}E[n_{k}] \end{eqnarray}

如果\(E[n_{k}] = 0\)，\(E[\epsilon_{x,k-1}] = 0\)，那么 \(E[\epsilon_{x,k}] = 0\) 。也就是说，如果估计噪声是零均值的，并且初始 的估计\(x_{0}= x\)，那么此后所有的估计都是\(x\)。从这点看来，式 (\ref{eq:42})给出来的估计是无偏的。并且这个特性与\(K_{k}\)无关，这是一 个无偏估计算子该有的样子。接下来我们计算\(K_{k}\)。