奇异值分解法计算广义逆

1、奇异值分解法计算广义逆线性最小二乘问题的广义逆求解丁梁波整理对于任意的mxn方程组：Ax = bal1aln其中A二 ： I ：a mla mn一 XiX = :Xn一 bib =:"m J如果m=n,只要n方阵A非奇异,就有逆阵A,从而得到解x = A/b.然而, 对于m#n的一般情况,A是长方阵,就没有通常的逆阵.不过它仍然可以有相 应于特定方程类型的几种形式的广义逆矩阵,其中适于任何情况的广义逆叫做 Penrose广义逆,记为A+.于是,方程的解可以为：由奇异值分解SVD可以将A分解为:A = U 三VT 其中U , V分别为m ,n阶正交阵*仃,0 +0 一这样A的广义逆A

2、4可表示为：a+二v"ut其中00仃1这样我们可以看出,完成A的奇异值分解后,求解A的广义逆就变得很简单,从而可以方便地求出方程组的最小二乘解. 下面我们说明对矩阵进行奇异值分解的 方法和步骤.通常情况下我们考虑m>n时矩阵A的奇异值分解,由于当 m<n时,可以将n-m 行补零使其成为方阵后再进行分解.这样我们就将矩阵A的奇异值分解分为两大 步,假设干小步如下：一.用Householder变换将A约化为双对角矩阵.具体步骤如下：1.以A的第1列作为v,取i=1 ,按以下式子构造 Householder矩P$ Q式中Hi为Q,为了方便以后的说明我们还用 Hi表示THi =I

3、 -岭T ,Vm )打：其中,Ui =(0,0,Vi sign(Vi) v 2,Vi 1,V=(0.,Vi,Vi1Vm)Tm|v(i)2=(£ v2)1/2 k ±2 .将Qi左乘A得到矩阵Qi A ,并以Qi A的第1行作为v,取i=2,按1式构造Householder矩阵H2,右乘QiA得到QiA H2.3 .取Q1A H2的第2列为v, i=2,按1式构造Householder矩阵Q2,左乘QiA H2,得到Q2 Q1A H2,并将计算Q2 Qi将其存入Qi 04 .取Q2 Q1A H2的第2行为v, i=3,按1式构造Householder矩阵H3,右乘 Q2 Q1

4、A H2,得至U Q2 Q1A H2 H3,并将 H2 H3存入 H2.5 .依次类推,计算出QnQn-1Q1AH2 H3-Hn-1为双对角矩阵,并将 QnQn-1Qi存入到Qi中,H2 H3Hn-1存入到H2中.QnQn-1 QiAH 2 H3 - Hn-1为双对角矩阵记为：但匕1|P2y3+ +-00000需要注意的是：当m=n时,只计算到Qn-1- Q1AH2 Ha - Hn-2二.用原点位移QR算法进行迭代,计算所有的奇异值,并最终结合一计算 出出U和Vc -ss cH 0 =式中c = 1 /r1.按下式列旋转矩阵H0(2)1*1s = 2 / rt _ B y2 -1 2.'

5、;:2.2_ 1112 .22 _ 2 2 .4 2 21/2n J n1- n n nd nJ 4 n nJ2并将计算BHo2 .按下式构造列旋转矩阵一 c s-s cQ1 =1*J并计算Q1 BHo3 .构造列旋转矩阵11c -s s cH1=1+:.1 一并计算Q1 BH0H1以及HoH14 .构造列旋转矩阵一11c s-s cQ2 =,1+1 一并计算Q2 Qi BH0H1以及Q2 Qi5 .按类似(3), (4)的方法构造列旋转矩阵,并计算相应的新矩阵QlQ2 Q1BH0H1 Hi-1 , 直 至 iji=n , 并 记 Q11 = Qn Q2Q1,H 11 =H0H1 Hn/,B1

6、 =Q； =QnQ2Q1BH0H1Hn即 B1 =Q；BH；6 .判断B1的次对角线元素是否在误差范围内可以认为是0,假设是那么分解完毕,假设否,那么将B1作为上面的B重复步骤1, 2, 3, 4, 5, 6.直到Bk可以近似看 作是对角阵.即：Bk =Q：B-Hkk；iEQ =Q QfQi , H =H11H； H：；那么Bk的对角线元素就是矩阵 A的奇异值,即A=U工VT中的工已经求得,从上面的过程中我们可以将A按下面的式子进行分解：A=Q1QBkHH2比照 A=UEVT , U =Q1QV = H：Ht £ = Bk,这样我们就完成了矩阵A的奇异值分解,由于U和V都是正交阵,我们能够得到A的广义逆A +,从而可以根据以下公式计算方程组的最小二乘解：x = A b程序说明：程序共有一个主程序MAIN.FOR和三个主要功能子程序：MAIN.FOR 一主要功能有：方程组的初始化,输出系数矩阵及其广义逆、广义逆 法的最小二乘解以及逆的逆对方法进行验证.BMUAV.FOR一程序的核心局部