奇异值分解法计算广义逆

1、奇异值分解法计算广义逆线性最小二乘问题的广义逆求解（丁梁波整理）对于任意的mxn方程组：ax = balla1n其中a=：a mla mn一 xjx = :xnbl b =:m j如果m=n,只要n方阵a非奇异，就有逆阵a，从而得到解x = a/b。然而, 对于m#n的一般情况，a是长方阵，就没有通常的逆阵。不过它仍然可以有相 应于特定方程类型的几种形式的广义逆矩阵，其中适于任何情况的广义逆叫做 penrose广义逆，记为a+。于是，方程的解可以为：由奇异值分解（svd）可以将a分解为:a = u 三vt 其中u , v分别为m ,n阶正交阵1* =仃r0 +。一这样a的广义逆a 4可表示为:

2、a+二vut巴其中41/01一r| 00这样我们可以看出，完成a的奇异值分解后，求解a的广义逆就变得很简单，从而可以方便地求出方程组的最小二乘解。 下面我们说明对矩阵进行奇异值分解的 方法和步骤。通常情况下我们考虑mn时矩阵a的奇异值分解，因为当 mn时，可以将n-m 行补零使其成为方阵后再进行分解。这样我们就将矩阵a的奇异值分解分为两大 步，若干小步如下：一.用householder变换将a约化为双对角矩阵。具体步骤如下：1 .以a的第1列作为v,取i=1 ,按下列式子构造 householder矩p$ q式中hi为q,为了方便以后的说明我们还用 hi表示tt ,vm )hi =i -岭其中

3、,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）式构造householder矩阵h2,右乘qia得到qia h2。3 .取qia h2的第2列为v, i=2,按（1）式构造householder矩阵q2,左乘qia h2,得到q2 qia h2,并将计算q2 qi将其存入qi 04 .取q2 q1a h2的第2行为v, i=3,按（1）式构造householder矩阵h3,右乘 q2 qia h2,得至u q2 q1a h2

4、 h3,并将 h2 h3存入 h2。5 .依次类推，计算出qnqn-1- q1ah2 h3 hn-1为双对角矩阵，并将 qnqn-1 q1存入到q1中,h2 h33hn-1存入到h2中。qnqn-1- q1ah2 h3 - - hn-1为双对角矩阵记为：但匕1ip2y3+ +-00000需要注意的是：当m = n时，只计算到qn-qiah2包山？二.用原点位移qr算法进行迭代，计算所有的奇异值，并最终结合（一）计算 出出u和vc -ss ch 0 =式中c = 11r(2)1.按下式列旋转矩阵ho1*1s = 2 / rt _ b y2 -1 2.:2 . 2 _ 1112 . 22 _ 2

5、2 .4 2 2 1/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 qibh0h1 hi-i , 直 至iji=n , 并 记 q11 = qn q2q1,h 11 =h0h1 hn，b1 =q； =qnq2q1bh0h1hn即

6、b1 =q1bh；6 .判断b1的次对角线元素是否在误差范围内可以认为是0,若是则分解完毕，若否，则将b1作为上面的b重复步骤1, 2, 3, 4, 5, 6。直到bk可以近似看 作是对角阵。即：bk =q：bkhk；记 q = q：；q；q； , h=h11h:h：j则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一程序的核心部分，