经典最小二乘与全最小二乘法及其参数估计

理论新探

经典最小二乘与全最小二乘法及其参数估计

王福昌1，曹慧荣2，朱红霞2

（1.中国地震局防灾科技学院基础部，河北三河065201；2.廊坊师范学院数信学院，河北廊坊065000）

摘要：文章对经典的最小二乘和全最小二乘方法的应用背景、原理与算法进行了介绍，给出

了它们在线性模型参数估计中的MATLAB实现；通过计算机仿真说明了在模型中所有变量均具有不可忽略的误差时，全最小二乘法得到的参数估计更接近于真实参数。

关键词：最小二乘；全最小二乘；线性模型中图分类号：O212

文献标识码：A

文章编号：1002-6487（2009）01-0016-02

在模型的参数估计中，最常见的一种拟合准则是经典的最小二乘法，该准则下参数估计特别简单，其假设检验也容易进行。值得指出的是，经典的最小二乘法要求解释变量均为精确无误差的，或其测量误差与模型的因变量的测量误差相比可以忽略不计，所有误差均来自于因变量，这在很多情形下是不能满足的，因而人们提出了考虑解释变量误差模型的全最小二乘（TotalLeastSquares）准则，又叫做正交最小二乘(OrthogonalLeastSquares,OLS)、变量含误差模型(Errors–In-Variables,EV，EIV)、度量误差模型（ErrorMeasurement

n

记X=

ΣΣΣΣΣΣΣ…ΣΣΣΣΣ

11

x11x21xn1

……

x1px

………

…

1x

ΣΣΣΣ2pΣ

ΣΣΣΣΣΣnpΣ

，β=

ΣΣ0ΣΣΣ1ΣΣ…ΣΣΣΣp

bbbΣΣΣΣΣΣΣΣΣΣΣ

，Y=

ΣΣ1ΣΣΣ2ΣΣ…ΣΣΣΣΣn

yy

y

ΣΣΣΣΣΣΣΣΣΣΣΣ

，则问题（1）变为

Q2=Σεi2=||Y-Xβ||2=(Y-Xβ)T(T-Xβ)

i=1

（2）

由向量微分理论可知，对Q2关于向量β求导，令鄣Q=-

2XT(y-Xβ)，由实际意义可知问题的最优解为

赞LS=（XTX）-1XTyβ

（3）

Model）和随机回归模型(RandomRegressor)等等。近年来，全

最小二乘法在统计分析、线性和非线性回归、系统辨识和参数估计及信号处理中有广泛的应用，成为数理统计和数值代数方向的热门问题之一。

本文拟针对线性模型，给出估计参数的算法和实现

赞由于MATLAB是一种向量化编程语言，因此计算系数β

特别简便。在MATLAB语言中使用命令Beta=X\y或者

Beta=pinv(X’*X)*X’*y皆可得到参数估计值。2

线性模型参数估计的全最小二乘法

在经典的最小二乘法中，误差被定义为数据点沿着因变量y轴方向到拟合函数的偏差。而当模型中难以区分因变量和解释变量时，即所有变量都有不可忽略的误差时，使用经典的最小二乘法已经不再合适。这就需要使用考虑了所有变量误差的全最小二乘法。全最小二乘法按照与经典最小二乘法不同的距离标准来估计参数，这时是寻找参数b0,b1,b2,…,

MATLAB源程序，并对两种方法进行计算机仿真，希望研究

结果对使用最小二乘法解决实际问题的科技工作者具有一定的参考价值。

1线性模型参数估计的最小二乘法

已知多元线性模型中有p+1个参数b0,b1,b2,…,bp待定，解

释变量x1,x2,…,xp和因变量y有n次测量值(n>p+1)，假设满足

bp使得数据点(xi1,xi2,…,xip,yi)到超平面距离残差平方和

yi=b0+b1xi1+b2xi2+…+bpxip+εi(i=1,2,…,n)

这是一个超定线性系统，方程可能没有解。在经典的最

QT=Σε=(1+Σb)

2

Ti

i=1

j=1

np

2-1j

i=1

Σ[y-(b+bx+bx+…+bx)]

i

1i1

2i2

pip

n

2

（4）

小二乘法中，假设x1,x2,…,xp的观测值没有误差，寻找参数b0,

取得最小值。如再记Xi=[1,x11,x12,…,x1p]，β0=[b1,b2,…,bp]T，则问题为

n

n

b1,b2,…,bp使得沿着y轴方向的残差平方和

Q2=Σε=Σ[yi-(b0+b1xi1+b2xi2+…+bpxip)]

2i

i=1

i=1

n

n

2

（1）

QT=Σε=(1+ββ0)

2Ti

T0

i=1

-1

i=1

Σ(y-xβ)

i

i

2

（5）

取得最小值。

问题（4）可以看成一个多元函数的最小值问题，可以用无约束优化的有效算法如BFGS算法求解，也可以用下面的

Gauss-Newton迭代法和SVD分解方法求解。2.1

统计与决策２００9年第1期（总第２77期）

用Gauss-Newton迭代法估计参数β

理论新探

軃軃1軃軃軃軃軃2軃軃軃軃軃軃軃軃

y-(b0+b1x11+b2x12+…+bpx1p)

姨1+b+…+b

11

p

向量值函数F(β)=

y-(b0+b2x21+b2x22+…+bpx2p)

姨1+b+…+b姨1+b+…+b

1p

p

yn-(b0+b1xn1+bnxn2+…+bpxnp)

軃軃軃軃軃軃軃軃軃軃軃軃軃軃軃軃軃軃

=

軃軃軃軃軃軃軃軃軃軃軃軃軃軃軃軃

y1-x1β

1(1+βT0β0)y1-x2β

(1+βT0β0)y1-xnβ

1(1+βT0β0)

軃軃軃軃軃軃軃軃軃軃軃軃軃軃軃軃軃軃軃

赞OT=(ATA-σp+12I)-(3)计算过原点的超平面的全最小二乘解β

1

ATb；

赞TLS=[b-Aβ赞(4)平移得到原始问题的全最小二乘回归系数β

，令

OT

赞OT]。;β

MATLAB程序略。

计算机仿真

为直观起见，讨论二维平面上确定拟合直线的情形。假

F(β)=0（6）

2.3

其中O为维列向量，则求问题（6）最小二乘解与求问题（4）的最小值点是等价的，可以通过Gauss-Newton迭代方法求解。

相应的向量值函数F(β)关于向量β的Jacobi矩阵为

軃

軃軃軃軃軃軃軃軃軃軃軃軃軃軃軃

设原始的直线方程为y=b+kx，其中b=1,k=0.5，围绕这条直线产生一组正态分布的随机数据，然后分别用经典最小二乘和全最小二乘法确定回归系数系数，并进行比较。

-b0

T0

2

-x11(1+βT0β0）-b1（y1-x1β)…-x1p(1+βT0β0）-bp（y1-x1β)

(1+ββ0)

T

T0

2

(1+ββ0)-b0

T0

(1+ββ0)

T0

T0

2

DF(β)=鄣F

(1+ββ0)-b0

-x21(1+ββ0）-b1（y1-x2β)…-x2p(1+ββ0）-bp（y1-x2β)

(1+ββ0)

T0

(1+ββ0)T0

(1+βT0β0)

-xn1(1+βT0β0）-b1（y1-xnβ)…-xnp(1+βT0β0）-bp（y1-xnβ)

(1+βT0β0)

(1+βT0β0)

軃

軃軃軃軃軃軃軃軃軃軃軃軃軃軃軃軃軃軃

利用MATLAB编写程序，独立运行50次，得到的仿真结果见表1。

表1

独立运行50次的得到的回归系数比较

截距均值和标准差

经典最小二乘法全最小二乘法

斜率均值和标准差

1.5511±0.69641.0689±0.74750.4436±0.06100.4915±0.0673

估计β的Gauss-Newton迭代算法描述如下：

(1)取初值β(0)，指定精度ε；

(2)计算向量值函数F(β)和函数的Jacobi矩阵Ak=DF

(k)

由于原始的系数为1和0.5，可见全最小二乘法得到回归系数更接近原始系数，只是由它得到的回归系数的标准差稍大一些。

为了进一步比较，在图1中虚线表示原始的直线方程，实点表示围绕这条直线产生一组正态分布的随机数据，实线表示按照经典最小二乘法和全最小二乘法确定的直线。

(β)进行迭代：β=β-(AAk)F(β)，k=0,1,2,…；

(k)

(k+1)

(k)

Tk

-1

(k)

(3)若||β(k+1)-β(k)||

如算法迭代停止，则迭代序列{β}收敛于线性模型参数

(k)

赞TLS。的全最小二乘估计β

尽管这种方法理论上存在着初始解选取和迭代序列发

赞LS为初始值β(0)可散的问题，但实际经验表明多数情形下以β赞TLS。以很快收敛到β的全最小二乘估计β2.2

用SVD分解法估计参数

为便于理解，从几何的观点考察线性拟合问题。对于2维数据来说，要寻找1维的一条直线近似刻画它们的关系；对于3维数据来说，要找2维的一张平面来近似刻画它们的关系。假设有p+1维空间中一组数据，我们要用p维子空间

图1随机仿真31个数据点情形下回归结果比较

从图1中可以看出，全最小二乘法得到的回归直线与原始直线几乎重合，而经典最小二乘法与原始直线差别较大。

(直线、平面或超平面)来拟合p+1维空间中的数据，因此这是

数值分析上的一个降秩逼近问题，可以使用奇异值分解

3结果与讨论

本文讨论了线性模型系数估计的经典最小二乘法和全

(SVD)方法解决。算法如下：

(1)考虑到计算简便和数据点到超平面的距离平方和与

超平面是否通过原点无关，首先对数据平移，记平移后的数

軃

軃軃11軃軃軃軃21軃軃軃…軃軃軃軃n1

最小二乘法，给出了估计参数的算法和MATLAB源程序，对两种方法进行了计算机仿真，发现在模型中所有变量均含有误差时，由全最小二乘法得到回归参数与真实模型参数更为接近。这些算法和程序对于使用最小二乘法解决工程问题的科技工作者具有一定的参考价值。

参考文献：

軃x-x軃…x-x軃x-x12121p軃x-x軃…x-x軃x-x12222p

…

…

…

据为A=

軃x-x軃…x-x軃x-x12n2np

軃

軃p軃軃軃軃p軃軃軃軃軃軃軃p軃

，b=

軃軃軃1軃軃軃軃2軃軃軃…軃軃軃軃n

軃y-y軃y-y

軃y-y

軃

軃軃軃軃軃軃軃軃軃軃軃軃軃

，其中j=1

Σx，j=

ij

i=1

n

1,2,…,p；1

Σy。

i

i=1

n

(2)对M=[Ab]进行奇异值分解，使得M=USVT，其中U、V

为正交阵，S=diag(σ12,σ22,…,σp+12)，σ12≥σ22≥σp+12≥0为对角阵。

[1]程龙生.基于errors-in-variables的预测模型及其应用[J].数理统计与管理,2005,24(2).

(责任编辑/亦民)

统计与决策２００9年第1期（总第２77期）

17