解线性最小二乘问题的正交化方法及应用

第7卷　第1期

太原师范学院学报(自然科学版)Vol.7No.1

　　2008年3月　JOURNALOFTAIYUANNORMALUNIVERSITY(NaturalScienceEdition)　Mar.2008

解线性最小二乘问题的正交化方法及应用

曹　蓉

(汕头职业技术学院,广东汕头515041)

　　〔摘要〕　文章讲述在求解线性最小二乘问题时系数矩阵出现病态时,应用Householder变换

把系数矩阵正交三角化,从而求出最小二乘解,以及在实际中的应用.

〔关键词〕　最小二乘问题;正交化方法;奇异矩阵〔文章编号〕　1672-2027(2008)01-0035-04　〔中图分类号〕　O241.2　〔文献标识码〕　A

0　引言

最小二乘问题在数据拟合、参数估计和函数逼近等方面有着广泛的应用.从最优化方法的观点解释最小二乘问题,就是要求矛盾方程组Ab=Y中的参数b在最小二乘意义下的最佳估计值.记矛盾方程组[1](1)Ab=Y

其中A=(aij)为m×n(m>n)矩阵,b=(b1,b2,…,bn)T,Y=(y1,y2,…,yn)T.矛盾方程组(1)一般不存在通常意义下的解,即任何n维向量b,一般Y-Ab≠0.现在应用加权最小二乘法估计参数向量b,即矩阵为:(AWA)b=AWY(2)

其中W=diag(w1,w2,…,wm),wi>0(i=1,2,…,m),若矩阵A=(aij)是列满秩矩阵,即rank(A)=n,得出(2)的系数矩阵AWA是非奇异的,则可以用古典最小二乘求解得b=(AWA)精度,有时ATWA并不一定非奇异.

11

设矛盾方程组(1)中的系数矩阵A=

w1+w2

ATWA=

w1

2

T

^

T

-1

T

T

AWY,但要考虑计算机的

T

00

,于是法方程组(2)中的系数矩阵:

0 000

w1w1

w1

2

w1+w3

2

w1w1w1+w42

假定 =10-3且计算ATWA是在字长为6的10进制浮点数的计算机上进行的,则 i+ i = i+ i10-6

被舍入为 1,这意味着A的最后四行的信息全部丧失,则ATWA的浮点表示fl(ATWA)贮存在计算机内为:

w1w1w1

w1w1w1

w1w1.w这是一个奇异矩阵,即(2)是病态的.也就是说最小二乘估计要求系数矩阵必须是非奇异的,而当系数矩阵出现奇异时,无法用古老的最小二乘估计,为此我们将系数矩阵进行变换,将系数矩阵正交三角化,然后求最小二乘解.

1　Householder变换

Gauss消去法是用左乘一系列初等下三角阵来约化一个矩阵为上三角阵.讨论另一种三角化方法

23

,

36

太原师范学院学报(自然科学版)　　　　　　　　　　　　　　　第7卷　

—Householder变换,它的应用非常广泛,并将应用它解最小二乘问题.

定义[2～3]　Householder矩阵是指形式为I-2uuT的矩阵,其中uTu=1,通常记作H.

根据定义知H有性质:对称性(HT=H),正交性(HTH=I)和对合性(H2=I).

-1T

定理1[2～3]　设u≠0,令 =‖u‖2,则H=I- 是一个Householder矩阵.uu2

Householder矩阵的一个关键是能把一个给定向量的若干个指定的分量变为零,具体地说就是下列定理:定理2　设0≠X=(x1,x2,…,xn)∈Rn,!=±‖X‖,且假定X≠-!e1,则可找到一个Householder矩阵H使HX=-!e1,其中e1=(1,0,…,0)T.

证明　作u=X+!e1,求出 =‖u‖2,则H=1- -1uuT则即为所求的Householder矩阵.

2

T222

因为 =(X+!e1)(X+!e1)=(‖X‖+2!x1+!)=!+!x1

22

11T-1T

所以HX=(I- uu)X=XX=-!e1.

!+!x1

定理2的证明是构造性的,即它叙述了计算u与 的过程,一旦!的正负确定后,就可以求出u1=x1+!,ui=xi(i=2,…,n)以及 =!(!+x1).如果!与x1异号,在计算时就会发生抵消,因此,我们取!=sign(x1)‖X‖.

这样我们将运用Householde变换把系数矩阵化为上梯形矩阵.

设A=[aij]∈Rm×n,m≥n,rankA=r>0,且A的前r列线性无关.保留H矩阵的性质,寻找一系列H矩阵H1,H2,…,Hk使HkHk-1…H2H1A为一个上梯形矩阵.记A=[a1,a2,…,an]=[a(11),a(21),…,a(n1)],其中a(j1)=(a1j,a2j,…,amj),j=1,2,…,n,另ei表示m-k+1维单位坐标向量,它的第i个分量是1,其余分量都是0.

m

(1)

(1)

(1)

(k)

第一步,令u1=a,- e,其中 1=-sign(a11)!1,!1=(∑a),H1=Im-b1u1u1,其中i=1

(1)

1

(1)11

2i1

-1

T

122(1)

1- 1 11b1=2‖u1‖=

1

a12

(2)

…… …

a1n

(2)

则H1A=[H1a,H1a,…,Ha]=[ e,a,…, ]=

(1)1(1)2(1)1n(1)11(2)2(2)n

0

2)a(22a(22)n

2)a(mn

.

假设进行了k-1步,得到Householder变换H1,H2,…,Hk-1,使Ak=Hk-1Hk-2…H1A=

A11

(k)

a(m2)2

,(k)

0A22m-k+1k-1　n-k+1

～

A12k-1

(k)

(k)

k)k)

其中A(11是上三角阵,假定A(22=

～

ak,…,an

～

(k)

,其中a(jk)=

～～

akj,…,amj

m

(k)

～

(k)T

,j=k,…,n.

1

(k)(k)(k)2-1T~第k步,令uk=a(kk)- ke1,其中 k=-sign(akk)!k,!k=(∑(aik),Hk=Im-k+1-bkukuk,其中i=k

～(k)(k)22(k)~bk=‖ak- ke1‖= k- kakk,则Hk=diag(Ik-1,Hk)仍为H矩阵,于是Ak+1=HkAk=

2

最后得到HkHk-1,…H2H1A=

R　R1

=

U

k)A(11

k)A(12

.~kA(22k)

H

,其中R为k×k阶上三角矩阵,R1为k×(n-k)阶矩阵,U

0　00

=R,R为k×n阶上梯形矩阵.记Q=H1H2…Hk-1Hk(3)

U

则(3)式可写成A=Q.特别,若r=rankA=n,则当m>n时,经n步可将矩阵A化为一个上三角阵,得

到A

R为n×n阶上三角阵;当m=n时,经n-1步可将A化为一个上三角阵,得到A=QR.

　第1期　　　　　　　　　　　　　曹　蓉:解线性最小二乘问题的正交化方法及应用

37

for　　j=1∶n

[v,#]=house(A(j∶m,j))

T

A(j∶m,j∶n)=(Im-j+1-#vv)A(j∶m,j∶n)d(j)=#

if　j<m

　A(j+1∶m,j)=v(2∶m-j+1)end

end

function:[v,#]=house(x)

n=length(x)

=‖x‖∞　　x=x/

T

!=x(2∶n)x(2∶n)

v(1)=1∶v(2∶n)=x(2∶n)

if!=0　#=0else

　 =x(1)+!　if　x(1)≤0

　　v(1)=x(1)- else

　v(1)=-!/(x(1)+ )

end

22

#=2v(1)/(!+v(1));v=v/v(1)

end

2　解最小二乘问题的正交化方法

求解线性最小二乘问题的新算法,,而直接从矛盾方程组Ab=Y入手.常应用

R

Householder变换把系数矩阵A正交三角化,使QA=,其中R为n阶上三角阵,0为(m-n)×n的零矩阵,

由(3)式得到的一个m.mQY相应地分块成n维向量c与(m-n)维向量d,即QY=ccRc-Rb.于是Qr=QY-QAb=-(4)b=dd0　d

因Q是正交矩阵,所以‖r‖2=‖Qr‖2=‖c-Rb‖2+‖d‖2,若选择b,c-Rb=0(5)

00

那么‖r‖2将达到极小值,此时由(4)可得‖r‖2=‖d‖2,且从(4)有Qr=,故r=QT.

dd

由(5)可知,n阶上三角形方程组的解b就是最小二乘解,它是非常容易求解的.这样,可描述正交化方法求最小二乘解的计算步骤为:

1)对m×n矩阵A建立QR分解A=QR;2)计算C=QY,形成方程组Rb=c;

3)用回代法解Rb=c得出最小二乘解b.

正交化方法因其数值 …… 此处隐藏：3689字，全部文档内容请下载后查看。喜欢就下载吧 ……