地球物理系反演报告
实验一奇异值分解计算广义逆G+
专业：地球物理学
姓名：
学号：
指导教师：邵广周
实验一 奇异值分解计算广义逆G +
一、基本原理
对于任意的n m ⨯方程组：b Ax =
其中⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=mn m n a a a a A
1
111
⎥⎥
⎥⎦
⎤
⎢⎢⎢⎣⎡=n x x x 1 ⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=m b b b 1 如果n m =，只要n 方阵A 非奇异，就有逆阵1-A ，从而得到解b A x 1-=。

然而，对于n m ≠的一般情况，A 是长方阵，就没有通常的逆阵。

不过它仍然可以有相应于特定方程类型的几种形式的广义逆矩阵，其中适于任何情况的广义逆叫做Penrose 广义逆，记为+A 。

于是，方程的解可以为：
b A x +=
由奇异值分解（SVD ）可以将A 分解为：
T V U A ∑=
其中U ，V 分别为m ,n 阶正交阵
⎥
⎥⎥⎥⎥⎥⎥⎥⎦⎤
⎢⎢⎢⎢⎢
⎢⎢⎢⎣
⎡=∑00
1
r
σσ 这样A 的广义逆+A 可表示为：
T U V A 1-+∑=
其中
⎥⎦⎤⎢⎣⎡∑=∑--0001
1
r
⎥⎥⎥⎦⎤⎢⎢⎢⎣
⎡=∑---1111r r
σσ
这样我们可以看出，完成A 的奇异值分解后，求解A 的广义逆就变得很简单，从而可以方便地求出方程组的最小二乘解。

下面我们说明对矩阵进行奇异值分解的方法和步骤。

通常情况下我们考虑m>n 时矩阵A 的奇异值分解，因为当m<n 时，可以将n-m 行补零使其成为方阵后再进行分解。

这样我们就将矩阵A 的奇异值分解分为两大步，若干小步如下：
一、用Householder 变换将A 约化为双对角矩阵。

具体步骤如下：
1. 以A 的第1列作为v ，取i=1，按下列式子构造Householder 矩阵Q 式中i H 为Q ，为了方便以后的说明我们还用i H 表示
2
/122
)
(1)(12
)(22
)(),,,,0,0(),,,)(,,0,0()
1(∑=++==+=-
=m
i k k i T m i i i T m i i i i i i
T i i i v v
v v v v v v v v sign v u u u u I H 其中，
2. 将Q 1左乘A 得到矩阵Q 1 A ，并以Q 1 A 的第1行作为v ，取i=2，按（1）式构造Householder 矩阵H 2， 右乘Q 1A 得到Q 1A H 2。

3. 取Q 1A H 2的第2列为v ，i=2，按（1）式构造Householder 矩阵Q 2，左乘Q 1A H 2，得到Q 2 Q 1A H 2，并将计算Q 2 Q 1将其存入Q 1。

4. 取Q 2 Q 1A H 2的第2行为v ，i=3，按（1）式构造Householder 矩阵H 3，右乘Q 2 Q 1A H 2，得到Q 2 Q 1A H 2 H 3，并将H 2 H 3存入H 2。

5. 依次类推，计算出Q n Q n-1…Q 1AH 2 H 3…H n-1为双对角矩阵，并将Q n Q n-1…Q 1存入到Q 1中，H 2 H 3…H n-1存入到H 2 中。

Q n Q n-1…Q 1AH 2 H 3…H n-1为双对角矩阵记为：
⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎢
⎢⎢⎢⎢⎣
⎡=-00
01
3
22
1
n n n B βγβγβγβ 需要注意的是：当n m =时，只计算到Q n-1…Q 1AH 2 H 3…H n-2
二、用原点位移QR 算法进行迭代，计算所有的奇异值，并最终结合（一）计算出出U 和V 。

1. 按下式列旋转矩阵H 0
⎥⎥
⎥
⎥⎥
⎥⎦
⎤
⎢⎢⎢⎢⎢⎢⎣⎡-=110 c s s c H （2）
式中
()
(
)
[]
2
/121
2
2
2
1
21222121222
122112
2212142
121//-----+--+-+++==-=+===n n n n n n n n n n r r
s r c γ
γγβγβγβγβσγβξσβξξξξξ
并将计算BH 0
2. 按下式构造列旋转矩阵并计算Q 1 BH 0
⎥⎥
⎥
⎥⎥
⎥⎦
⎤⎢⎢⎢⎢⎢⎢⎣⎡-=11
1 c s s c Q
3. 构造列旋转矩阵并计算Q 1 BH 0H 1以及H 0H 1
⎥⎥
⎥⎥
⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡-=1111 c s s c H
4. 构造列旋转矩阵并计算Q 2 Q 1 BH 0H 1以及Q 2 Q 1
⎥⎥
⎥⎥
⎥⎥⎥⎥⎦
⎤
⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡-=1112 c s s c Q
5. 按类似（3），（4）的方法构造列旋转矩阵，并计算相应的新矩阵Q i …Q 2 Q 1
BH 0H 1…H i-1,直到i=n ，并记1211Q Q Q Q n =，1101
1-=n H H H H ，110121
11-==n n H H BH Q Q Q Q B ，即1111
1BH Q B = 6. 判断B 1的次对角线元素是否在误差范围内可以认为是0，若是则分解完毕，若否，则将B 1作为上面的B 重复步骤1，2，3，4，5，6。

直到B k 可以近
似看作是对角阵。

即：1
1111-----=k k k k k k H B Q B
记112211Q Q Q Q k k --=，112211--=k k H H H H
则B k 的对角线元素就是矩阵A 的奇异值，即T V U A ∑=中的∑已经求得，从上面的过程中我们可以将A 按下面的式子进行分解：
21HH QB Q A k =
对比T V U A ∑=，k T T B H H V Q Q U =∑==，，21，这样我们就完成了矩阵A 的奇
异值分解，由于U 和V 都是正交阵，我们能够得到A 的广义逆+A ，从而可以根据下列公式计算方程组的最小二乘解：
b A x +=
二、程序设计及结果分析
1、本次实验给的是109⨯的方程组：b Ax =，具体数值如下：
1001001001601001001016001001001
1611100000023000111000
260000
0011
1291.414000 1.4140000.41438.4850000.4140000.414
037.0713⎡⎤⎡⎤⎡⎤
⎢⎥⎢⎥⎢⎥⎢
⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢
⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥
⎢⎥⎢⎥⎢⎥=⎢
⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥
⎢⎥⎢⎥⎢⎢⎥⎢⎥⎢⎢⎥⎢⎥⎢⎣
⎦⎣⎦⎣⎦
⎥⎥⎥ 其中，计算精度为：000001.0=EPS 。

2、计算结果
123456789 1.000350.999651.000002.000352.000001.999652.999293.000353.00035x x x x x x x x x ⎡⎤⎡⎤
⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥
⎣
⎦⎣⎦
3、结果分析
由以上结果分析可知，奇异值分解对求解欠定矩阵是很有用的。

计算结果与
真实模型值相符。

且误差较小。

可见，奇异值分解对于求解奇异矩阵还是很有用的，具有很
好的可行性，其结果也符合方程组。

通过本次实验，我对奇异值分解
有了更进一步的理解，了解了其基本原理及具体求解步骤。