奇异值分解法计算广义逆

奇异值分解法计算广义逆 线性最小二乘问题的广义逆求解
(丁梁波 整理)
对于任意的n m ⨯方程组：b Ax =
其中⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=mn m n a a a a A 1111 ⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=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 阶正交阵
⎥
⎥⎥
⎥⎥⎥⎥⎥⎦⎤
⎢⎢⎢⎢⎢
⎢⎢⎢⎣
⎡=∑001
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
⎥⎥
⎥
⎥⎥
⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡-=11
0 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. 按下式构造列旋转矩阵
⎥⎥
⎥
⎥⎥
⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡-=11
1 c s s c Q 并计算Q 1 BH 0 3. 构造列旋转矩阵
⎥⎥
⎥⎥
⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡-=1111 c s s c H
并计算Q 1 BH 0H 1以及H 0H 1
4. 构造列旋转矩阵
⎥⎥
⎥⎥
⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡-=1112 c s s c Q
并计算Q 2 Q 1 BH 0H 1以及Q 2 Q 1
5. 按类似（3），（4）的方法构造列旋转矩阵，并计算相应的新矩阵Q i …Q 2 Q 1 BH 0H 1
…
H i-1,
直
到
i=n
，
并
记
1
21
1Q Q Q Q n =，
1101
1-=n H H H H ，11012111-==n n H H BH Q Q Q Q B 即1
1111BH 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 +=
程序说明：
程序共有一个主程序MAIN.FOR和三个主要功能子程序：
MAIN.FOR—主要功能有：方程组的初始化，输出系数矩阵及其广义逆、广义逆法的最小二乘解以及逆的逆对方法进行验证。

BMUA V.FOR—程序的核心部分，奇异值分解程序，输入系数矩阵，输出分解后的U，V，∑
AGMIV.FOR—计算广义逆+A以及方程组的最小二乘解
BGINV.FOR—仅计算广义逆+A
另外还有一个奇异值分解的辅助小程序SSS.FOR。