最小二乘法溯源追踪

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

6=β1+β25=β1+2β27=β1+3β310=β1+4β3

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

最小二乘的目的在于找到β1,β2使得式 (5)最小：

S(β1,β2)=[6−(β1+β2)]2+[5−(β1+2β2)]2+[7−(β1+3β3)]2+[10−(β1+4β3)]2

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

β1=3.5 β2=1.4

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

考虑系统：

n∑j=1Aijβj=yi,i=1,2,…,m

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

y=β1H1+β2H2+β3H3+β4H4

其中H1代表面积，H2代表单价，H3代表卧室数量， H4代表车位数量，β1,…,β4代表这四个属性的 权重，y代表房价。这里我们假定这四个因素和总价y之间只存在线性关 系。为了得到β1,…,β4我们需要一些房子的数据表，根 据这个表，我们来拟合这个关系。

对应到式 (11) ，我们有m个观测，每个观测涉及n个属性。把 式~(11)写成矩阵的形式有：

Y=Aβ

其中：

Y=[y1,y2,…,ym]

Y是m×1的列向量。

A=[A11A12…A1nA21A22…A2n⋮⋮⋱⋮Am1Am2…Amn]

A是m×n的矩阵，m代表一共有m个观测数据， n 代表一共有n个属性。

β=[β1,β2,…,βn]

β是n×1的列向量，n代表n个属性。 β 向量代表了n个属性的权重。

通常，方程~(13) 没有解。我们的目的是找到一个最合适的 β 使得Y与Aβ的 差值的平方最小。

argminβS(Y,Aβ)

其中：

S(Y,Aβ)=m∑i=1|yi−n∑j=1Aijβj|2=‖Y−Aβ‖2

假定A的各列是相互独立的，也就是说n个属性互不相关，通 过正则化方程~(13)：

ATY=ATAβ

ATA是一个半正定的矩阵，其逆一定存在（假设 A 的各列是相互独立的）。最终，有：

β=(ATA)−1ATY

式~(17)就是最小二成法的解。

在上一节，我们不加说明的给出：式~(14)的解就是~(16)的 解。现在，我们给出式 (16)的由来。

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

ri=yi−n∑j=1Aijβj

那么：

S=m∑i=1r2i

为了求得S最小时的β，我们需要依S对 β求偏导。

∂S∂βj=2m∑i=1ri∂ri∂βj,j=1,2,…,n

另外：

∂riβj=−Hij

综上有：

∂S∂βj=2m∑i=1(yi−n∑k=1Aikβk)(−Aij)

令：

2m∑i=1(yi−n∑k=1Aikβk)(−Aij)=0

重新组织式 (23)，有：

m∑i=1n∑k=1AijAikβk=m∑i=1Aijyi,j=1,2,…,m

写成矩阵的形式有：

(ATA)β=ATY

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

定义：

S(Y,β)=‖Y−Aβ‖2

展开有：

(Y−Aβ)T(Y−Aβ)=YTY−βTATY−YTAβ+βTATAβ

对式 (27) 就β求导，有：

ATY=(ATA)Aβ

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

y1=x+n1y2=x+n2⋮⋮yk=x+nk

矩阵形式有：

[y1y2⋮yk]=[11⋮1]x+[n1n2⋮nk]

根据ˆx=(ATA)−1ATy，则：

ˆx=1kk∑i=1yi

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

ˆx=1kk∑i=1yi=x+1kk∑ini

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

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

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

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

[y1y2…yk]=[H11H12…H1nH21H22…H2n⋮⋮⋱⋮Hk1Hk2…Hkn][x1x2…xn]+[n1n2…nk]

且E(n2i)=σ2i， 即ni是均值为零，方差为 σ2i的高斯噪声。假设ni,i=1,…,k是相互独立的， 则有：

R=E[nnT]=[σ21…0⋮⋮0…σ2k]

不同的σ2i意味着其对应的测量结果的权重不同。 σ2i越小，则权重越大；反之，越小。

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

J=ϵ2y1/σ21+…+ϵ2yk/σ2k

其中ϵyi=yi−∑nj=1Hijxj. 所以 J又可以被写为：

J=ϵTyR−1ϵy=(y−Hˆx)TR−1(y−Hˆx)=(yTR−1−ˆxTHTR−1)(y−Hˆx)=yTR−1y−yTR−1Hˆx−xTHTR−1y+ˆxTHTR−1Hˆx

上式以ˆx求微分，令其等于零，有：

∂J∂ˆx=−yTR−1H+ˆxTHTR−1H=0

进而，有：

ˆx=(HTR−1H)−1HTR−1y

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

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

[y1y2⋮yk]=[11⋮1]x+[n1n2⋮nk]

考虑噪声的协方差矩阵：

R=diag(σ21,…,σ2k)

根据 式 (44)，有：

ˆx=(HTR−1H)−1HTR−1y=(k∑i=11/σ2i)−1k∑j=1yjσ2j

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

在本节，我们探讨如何迭代的更新ˆx，即：当我们已经根据k−1个 估计数据得到了ˆx的一个估计，现在又来了一个数据yk，如 何不用重新计算 (44)更新ˆx.

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

yk=Hkx+nkˆxk=ˆxk−1+Kk(yk−Hkˆxk−1)

我们根据ˆxk−1 和yk更新ˆxk。Kk是 估计算子的增益矩阵，将在后文推导其计算过程。yk−Hkˆxk−1是修正项。如果Kk=0 或者修正项为零，则 ˆxk=ˆxk−1。在计算Kk之前，我们首先考虑迭代 估计的误差期望：

E[ϵx,k]=E[x−ˆxk]=E[x−ˆxk−1−Kk(yk−Hkˆxk−1)]=E[x−ˆxk−1−Kk(Hkx+nk−Hkˆxk−1)]=E[ϵx,k−1−KkHk(x−ˆxk−1)−Kknk]=E[(ϵx,k−1−KkHkϵx,k−1)−Kknk]=(I−KkHk)E[ϵx,k−1]−KkE[nk]

如果E[nk]=0，E[ϵx,k−1]=0，那么 E[ϵx,k]=0 。也就是说，如果估计噪声是零均值的，并且初始 的估计x0=x，那么此后所有的估计都是x。从这点看来，式 (49)给出来的估计是无偏的。并且这个特性与Kk无关，这是一 个无偏估计算子该有的样子。接下来我们计算Kk。