奇異值分解法計(jì)算廣義逆

1、奇異值分解法計(jì)算廣義逆線性最小二乘問題的廣義逆求解（丁梁波整理）對于任意的mxn方程組：Ax = balla1n其中A=：a mla mn一 xjx = :Xnbl b =:"m J如果m=n,只要n方陣A非奇異，就有逆陣A，從而得到解x = A/b。然而, 對于m#n的一般情況，A是長方陣，就沒有通常的逆陣。不過它仍然可以有相 應(yīng)于特定方程類型的幾種形式的廣義逆矩陣，其中適于任何情況的廣義逆叫做 Penrose廣義逆，記為A+。于是，方程的解可以為：由奇異值分解（SVD）可以將A分解為:A = U 三VT 其中U , V分別為m ,n階正交陣1*£ =仃r0 +。一這樣A

2、的廣義逆A 4可表示為:a+二v"ut巴其中41/01一r| 00這樣我們可以看出，完成A的奇異值分解后，求解A的廣義逆就變得很簡單，從而可以方便地求出方程組的最小二乘解。 下面我們說明對矩陣進(jìn)行奇異值分解的 方法和步驟。通常情況下我們考慮m>n時矩陣A的奇異值分解，因?yàn)楫?dāng) m<n時，可以將n-m 行補(bǔ)零使其成為方陣后再進(jìn)行分解。這樣我們就將矩陣A的奇異值分解分為兩大 步，若干小步如下：一.用Householder變換將A約化為雙對角矩陣。具體步驟如下：1 .以A的第1列作為v,取i=1 ,按下列式子構(gòu)造 Householder矩P$ Q式中Hi為Q,為了方便以后的說明我

3、們還用 Hi表示TT ,vm )Hi =I -嶺其中,Ui =(0,0, ,Visign(Vi) v 2，Vi 1,v=(0。,Vi,Vi“Vm)Tm|v(i)|2=(£ v2)1/2 k ±2 .將Qi左乘A得到矩陣Qi A ,并以Qi A的第1行作為v,取i=2,按（1）式構(gòu)造Householder矩陣H2,右乘QiA得到QiA H2。3 .取QiA H2的第2列為v, i=2,按（1）式構(gòu)造Householder矩陣Q2,左乘QiA H2,得到Q2 QiA H2,并將計(jì)算Q2 Qi將其存入Qi 04 .取Q2 Q1A H2的第2行為v, i=3,按（1）式構(gòu)造Hous

4、eholder矩陣H3,右乘 Q2 QiA H2,得至U Q2 Q1A H2 H3,并將 H2 H3存入 H2。5 .依次類推，計(jì)算出QnQn-1- Q1AH2 H3 Hn-1為雙對角矩陣，并將 QnQn-1 Q1存入到Q1中,H2 H33Hn-1存入到H2中。QnQn-1- Q1AH2 H3 - - Hn-1為雙對角矩陣記為：但匕1IP2y3+ +-00000需要注意的是：當(dāng)m = n時，只計(jì)算到Qn-QiAH2包山？二.用原點(diǎn)位移QR算法進(jìn)行迭代，計(jì)算所有的奇異值，并最終結(jié)合（一）計(jì)算 出出U和Vc -ss cH 0 =式中c = 11r(2)1.按下式列旋轉(zhuǎn)矩陣Ho1*1s = 2 /

5、rt _ B y2 -1 2.':2 . 2 _ 1112 . 22 _ 2 2 .4 2 2 1/2n J n1- n n nd nJ 4 n nJ2并將計(jì)算BHo2 .按下式構(gòu)造列旋轉(zhuǎn)矩陣一 c s-s cQ1 =1*J并計(jì)算Q1 BHo3 .構(gòu)造列旋轉(zhuǎn)矩陣11c -s s cH1=1+:.1 一并計(jì)算Q1 BH0H1以及HoH14 .構(gòu)造列旋轉(zhuǎn)矩陣一11c s-s cQ2 二,1+1 一并計(jì)算Q2 Qi BH0H1以及Q2 Qi5 .按類似(3), (4)的方法構(gòu)造列旋轉(zhuǎn)矩陣，并計(jì)算相應(yīng)的新矩陣QlQ2 QiBH0H1 Hi-i , 直 至iji=n , 并 記 Q11 = Qn

6、 Q2Q1,H 11 =H0H1 Hn，B1 =Q； =QnQ2Q1BH0H1Hn即 B1 =Q1BH；6 .判斷B1的次對角線元素是否在誤差范圍內(nèi)可以認(rèn)為是0,若是則分解完畢，若否，則將B1作為上面的B重復(fù)步驟1, 2, 3, 4, 5, 6。直到Bk可以近似看 作是對角陣。即：Bk =Q：BkHk；記 Q = Q：；Q；Q； , H=H11H:H：j則Bk的對角線元素就是矩陣 A的奇異值，即A=U工VT中的工已經(jīng)求得，從上面的過程中我們可以將A按下面的式子進(jìn)行分解：A=Q1QBkHH2對比 A=UEVT , U =Q1QV = H：HT 工=Bk ,這樣我們就完成了矩陣A的奇異值分解，由于U和V都是正交陣，我們能夠得到A的廣義逆A +,從而可以根據(jù)下列公式計(jì)算方程組的最小二乘解：x = A b程序說明：程序共有一個主程序MAIN.FOR和三個主要功能子程序：MAIN.FOR 一主要功能有：方程組的初始化，輸出系數(shù)矩陣及其廣義逆、廣義逆 法的最小二乘解以及逆的逆對方法進(jìn)行驗(yàn)證。BMUAV.FOR一程序的核心部