奇异值分解

地球物理系反演报告

实验一奇异值分解计算广义逆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

一、用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 ，从而可以根据下列公式计算方程组的最小二乘解：